Jean-Jacques Delannoy

Symboles 𝜎

AHP Analytical hierarchic process AMC Analyse multi-critères AP Aquifère potential/ Aquifer potential CA Confined aquifers CaPA Canadian Precipitation Analysis CMC Canadian Meteorological Center CP Composante principale DEM Digital elevation model DFHN Distance from hydrological network DGW Deep groundwater DMBF Distance to a major bedrock fault DND Drainage network density

Analyse spatiale pour la hiérarchisation de bassins versants du milieu boréal québécois selon leur dépendance à l'eau souterraine.

Les deux premiers articles ont été rédigés en anglais tandis que le troisième article et le reste de la thèse sont rédigés en français. Compte tenu de la cotutelle internationale, les travaux de terrain ont pu être complétés grâce aux efforts conjoints du doctorant, de la direction de la recherche, du personnel professionnel de l'UQAT et de la communauté locale. Les méthodes et analyses utilisées dans cette thèse ont été initiées par le doctorant avec la concertation et l'approbation de la direction de la recherche.

En tant que premier auteur, le doctorant a été responsable de l'ensemble des phases permettant la réalisation de l'étude. La contribution d'Éric Rosa, directeur de recherche à l'UQAT, a concerné le financement du projet, la révision des trois articles ainsi que la collecte, l'analyse et l'interprétation des données nécessaires au 2 e article. L'apport de Didier Graillot, directeur de recherche à l'EMSE, a plus spécifiquement porté sur la relecture et révision des trois articles avant et après leur évaluation par les pairs. Frédéric Paran, co-directeur à l'EMSE, a contribué aux aspects liés à la collecte et l'analyse des données au sein des systèmes d'information géographique, à la mise en place de l'analyse multicritère et à la relecture des articles. Vincent Cloutier, codirecteur à l'UQAT, a fourni un support à la l'interprétation des données et à la relecture des trois articles. Finalement, Dylan Mayappo, collaborateur, a participé à la préparation et à la relecture du 1 er article et à la collecte des données du 2 e article. REMERCIEMENTS Pendant les quatre dernières années, plusieurs personnes ont participé, de près ou de loin, à mener à terme ce travail de thèse. Sans leur présence et leur accompagnement, la finalité de ce projet aurait été tout autre. Je tiens à remercier mon directeur de thèse à l'UQAT, Éric Rosa, pour sa confiance, sa disponibilité et son indéfectible appui tout au long du projet. Son soutien scientifique m'a permis de concrétiser ce projet dont les premières brides ont été élaborées ensemble en camping dans le Bas-Saint-Laurent. Je remercie mon directeur de thèse à l'EMSE, Didier Graillot, pour avoir embarqué dans le projet; la taille du territoire à étudier ayant éveillé sa curiosité scientifique. Son accueil chaleureux, son accessibilité et sa confiance m'ont accompagné jusqu'à la toute fin. Je tiens à remercier mes deux co-encadrant, Vincent Cloutier, pour son soutien scientifique et administratif ainsi que Frédéric Paran pour ses commentaires toujours constructifs au développement du projet. Un gros merci doit aussi être rendu à tous ceux qui m'ont accompagné et remplacé pour les travaux de terrain. Je pense notamment à Daniel, Guillaume, Dylan, et sans oublier mon ami et collègue

AP values calculated based on the MCA-AHP applied on a 490,000 km 2 territory (Québec, Canada)………………………….…..

3.1

Monitored watersheds, data acquisition and extent of the aquifer type over study area. b) Location and bedrock geology of the study area………………………………………………..……………….…..

3.2

Conceptual representation of the water cycle and annual pattern of the water budget components ……………………….……………….……

3.3

Conceptual representation of the annual geochemical cycle of river discharge…………………………………………….….……………..

3.4

Monthly cumulative water fluxes estimated from remote sensing data and spatial models……………………………………..………………

3.5

Temporal variations in river discharge………………….…………….

3.6

Box plot of geochemical data (𝛿 2 𝐻 -𝛿 18 𝑂 and σ) of groundwater, water from river, monthly composite precipitation, meltwaters from the snowpack, and meltwaters from the river ice cores.……………… xiv 4.12 Représentation cartographique de la projection des sous-bassins versants sur l'axe de la composante principale 2 comparant leur potentiel dépendance à la ressource en eau souterraine (PDES)……… 141 4.13 Résultats de l'analyse de classification hiérarchiques réalisée sur les résultats de l'ACP formant quatre groupes de bassins versants et distribution des regroupements déterminant le niveau de la vigilance proposée……………………………………………………………… 143 RÉSUMÉ Quantifier la ressource en eau stockée dans les aquifères peu profonds et déterminer leurs interactions avec l'eau de surface est essentiel pour mieux évaluer la sensibilité de l'eau souterraine et des écosystèmes qui en dépendent découlant du changement climatique et des impacts de l'activité humaine. Cependant, l'évaluation de la ressource en eau est complexe sur des territoires vastes éloignés où les données préexistantes sont rares et limitées. Cette étude se concentre sur l'utilisation d'approches SIG à grande échelle pour produire de nouvelles informations sur les ressources en eaux souterraines peu profondes. La zone d'étude couvre 185 000 km 2 englobant de vastes bassins versants boréaux du bassin de la Baie-James (Canada). Une analyse spatiale multicritères (AMC) basée sur des données géologiques et hydrogéologiques a été couplée à des calculs géométriques afin d'évaluer le volume des aquifères granulaires peu profonds et d'identifier les zones où la protection des eaux souterraines doit être prioritaire. Des données satellitaires et des modèles spatiaux combinés à des données in situ, telles que des mesures de débit et la composition géochimique de l'eau (isotopes stables [δ 2 H-δ 18 O] et conductivité électrique), ont ensuite été utilisés pour calculer le bilan hydrologique terrestre et estimer les composantes du cycle de l'eau. Enfin, un nouvel indice calculé à partir de 290 images thermiques du satellite Landsat 8, a été créé pour révéler l'emplacement de potentiels occurrences d'écosystèmes dépendant des eaux souterraines (POEDES).

Les résultats ont suggéré que les aquifères potentiels (AP) les plus importants sont principalement associés à des formations fluvioglaciaires. Ils contiennent en moyenne 8,5 m de sédiments saturés pour une évaluation globale d'environ 40 km 3 ± 10%. La contribution annuelle des eaux souterraines au débit des rivières est généralement inférieure à 12% des précipitations. La contribution de l'eau provenant des pluies et de la fonte des neiges est de 30 à 61% et de 18 à 40%, respectivement, avec des ratios évaporation sur les apports de 2 à 10%. Le débit présente des cycles saisonniers avec une contribution plus élevée d'eau de fonte des neiges pendant la crue printanière et de l'eau de pluie à l'automne, une composition évoluant vers celle des eaux souterraines pendant la période de glace et une contribution relativement stable de l'eau souterraine tout au long de l'année. Les terres proches des cours d'eau dans les vallées principales présentent les plus haut POEDES. Une analyse en composantes principales des caractéristiques des bassins versants révèle des corrélations significatives entre les paramètres des eaux souterraines estimés au cours de cette étude. Une topologie de bassins versants est proposée à partir des deux premières composantes qui expliquent 77% de la variance. Les bassins versants au sud du territoire dépendent plus des écoulements souterrains que ceux situés au nord qui montrent une surface plus imperméable. Quatre regroupements de bassins versants identifient les secteurs où une vigilance accrue des eaux souterraines est nécessaire pour préserver la ressource en eau et la biodiversité qui dépend des eaux souterraines.

ABSTRACT

Quantifying the water resources stored in shallow aquifers and their hydrological exchanges with surface waters is critical to better evaluate the sensitivity of groundwater resources and groundwater dependent ecosystems (GDEs) deriving from climate change and human impacts. The evaluation of water resources is still challenging in vast and remote areas, where pre-existing data are sparse and limited. This study focuses on the use of large-scale GIS-based approaches to produce new insights into shallow groundwater resources. The study area covers 185,000 km 2 encompassing vast boreal watersheds of the James Bay basin (Canada). A spatial multicriteria analysis (MCA) based on geological and hydrogeological data was coupled with geometric calculations in order to evaluate the volume of shallow granular aquifers and to identify areas where groundwater protection should be prioritized. Satellite data and spatial models combined with in situ data, such as discharge and geochemical data (stable isotopes [δ 2 H-δ 18 O] and electrical conductivity of water), were then used to calculate terrestrial water budget and estimate water cycle components. Finally, a new index based on the stability of the Temperature-Vegetation dryness index, calculated from 290 images of the satellite Landsat 8, revealed location of potential occurrence of groundwater dependent ecosystem (POGDE).

Results revealed that the highest Aquifer Potential are sparsely distributed and mostly associated with glaciofluvial formations. They contain an average of 8.5 m of saturated sediment for approximately 40 km 3 ± 10% of water over 7% of the territory. The groundwater contribution to river discharge is generally below 12% of watershed precipitation and reach 2.3 km 3 over study area. The rainfall and snowmelt contributed is in proportion of 30-61% and 18-40%, respectively with evaporation over inflow ratios of 2-10%. The drained water presents seasonal patterns with highest snowmelt contribution during the spring flood, rainwater dominating in fall, a composition evolving towards that of groundwater during ice-on period and groundwater providing a stable baseflow contribution throughout the year. Furthermore, isotopic mass balance suggest that plant transpire mostly water from rainfall. The land near streams in main valleys present highest POGDE. A Principal Component Analysis of watershed characteristics revealed significant correlations between groundwater parameters estimated during this study. A watershed topology is proposed from the two first components which explain 77% of the variance. From that one, the southern watershed mainly depends on groundwater flows compared to the northern one, which had more impermeable surface. Four clusters of watersheds locate where increased groundwater vigilance and preoccupation are required to preserve water resources and the current biodiversity that depends on groundwater.

CHAPITRE I INTRODUCTION GÉNÉRALE

1.1 Enjeux globaux de l'eau souterraine La disponibilité de la ressource en eau souterraine constitue un enjeu de société mondiale. Les estimations actuellement disponibles suggèrent que le volume total d'eau souterraine à l'échelle globale se situerait entre 8000 et 23 400 × 10 3 km 3 [START_REF] Oki | Global Hydrological Cycles and World Water Resources[END_REF][START_REF] Schlesinger | Chapter 10 -The Global Water Cycle[END_REF]Gleeson et al. 2015;[START_REF] Abbott | Human domination of the global water cycle absent from depictions and perceptions[END_REF], ce qui correspond approximativement à une tranche d'eau maximale de 170 m étalée sur la surface des continents de la Terre (138 × 10 6 km 2 ; Figure 1.1a). Ce volume représente 48% de la réserve d'eau douce mondiale [START_REF] Oki | Global Hydrological Cycles and World Water Resources[END_REF]. Considérant l'inaccessibilité de l'eau maintenue dans la glace, les aquifères deviennent les plus grandes réserves d'eau douce potentiellement disponibles pour l'approvisionnement en eau potable, bien qu'une portion significative de cette eau se situe dans les aquifères profonds et n'est pas exploitable. Ce n'est que 1,5%

(350 × 10 3 km 3 ) du volume total qui est considéré renouvelable (<50 ans; Gleeson et al. 2015). Aujourd'hui, l'eau souterraine permet l'approvisionnement en eau potable de 36% de population mondiale [START_REF] Döll | Impact of water withdrawals from groundwater and surface water on continental water storage variations[END_REF][START_REF] Taylor | Ground water and climate change[END_REF]. Elle permet aussi de répondre aux besoins écologiques d'un large éventail de milieux naturels [START_REF] Alley | Flow and Storage in Groundwater Systems[END_REF][START_REF] Hayashi | Effects of Ground Water Exchange on the Hydrology and Ecology of Surface Water[END_REF][START_REF] Bertrand | Review: From multi-scale conceptualization to a classification system for inland groundwaterdependent ecosystems[END_REF]Klove et al. 2014a;[START_REF] Eamus | Groundwater Dependent Ecosystems: Classification, Identification Techniques and Threats[END_REF]. Toutefois, certaines régions subissent des stress hydriques majeurs [START_REF] Konikow | Groundwater depletion: A global problem[END_REF][START_REF] Rodell | Satellite-based estimates of groundwater depletion in India[END_REF][START_REF] Feng | Evaluation of groundwater depletion in North China using the Gravity Recovery and Climate Experiment (GRACE) data and ground-based measurements[END_REF][START_REF] Famiglietti | The global groundwater crisis[END_REF]) et les perturbations liées à l'activité humaine engendrent des impacts à long terme sur la ressource [START_REF] Döll | Impact of water withdrawals from groundwater and surface water on continental water storage variations[END_REF][START_REF] Foster | Groundwater-a global focus on the 'local resource[END_REF] sociales et politiques différentes entre les régions influencent son accessibilité et nuancent la perception de la notion d'épuisement [START_REF] Gleeson | Global Groundwater Sustainability, Resources, and Systems in the Anthropocene[END_REF].

L'anthropisation du territoire influence la dynamique naturelle du cycle de l'eau et les rétroactions entre les bassins versants à toutes les échelles [START_REF] Gordon | Human modification of global water vapor flows from the land surface[END_REF][START_REF] Hoekstra | The water footprint of humanity[END_REF][START_REF] Abbott | Human domination of the global water cycle absent from depictions and perceptions[END_REF]. La régulation des débits et la mise en place de seuils et de barrages fragmentent l'écoulement des eaux de surface [START_REF] Dynesius | Fragmentation and Flow Regulation of River Systems in the Northern Third of the[END_REF][START_REF] Nilsson | Fragmentation and Flow Regulation of the World&#039[END_REF] et diminuent la connectivité du réseau hydrologique [START_REF] Grill | Mapping the world's free-flowing rivers[END_REF]. Il en résulte une modification des interactions entre l'eau souterraine et l'eau de surface au sein des hydrosystèmes [START_REF] Amoros | Les concepts d'hydrosystème et de secteur fonctionnel dans l'analyse des systèmes fluviaux à l'échelle des écocomplexes[END_REF]) rendant les milieux naturels dépendant de la proximité de l'eau souterraine plus vulnérables. La croissance démographique mondiale, le développement d'activités anthropiques et les impacts des changements climatiques planétaires anticipés engendrent une pression sans précédent sur la ressource en eau souterraine, certes renouvelable, mais dont la disponibilité est inégale.

Perspectives canadiennes et québécoises

La quantité d'eau souterraine au Canada est estimée à 70 × 10 3 km 3 [START_REF] Rivera | International Year of Planet Earth 3. Groundwater Sustainable Devlopment in Canada -Emergin Issues[END_REF]) ce qui représente approximativement une épaisseur de 8 m d'eau étalée sur l'ensemble du territoire national. En milieu boréal et au sein des bassins versants du Bouclier canadien, la précarité de la ressource en eau souterraine n'est pas encore préoccupante [START_REF] Oki | Global Hydrological Cycles and World Water Resources[END_REF]Kanae 2006, Döll 2009;[START_REF] Mekonnen | Four billion people facing severe water scarcity[END_REF]. Toutefois, elle n'est pas uniformément répartie et elle est localisée généralement dans des aquifères peu profonds (Figure 1.1b). La perméabilité du socle rocheux du Bouclier canadien est observable jusqu'à 300 m de profondeur, mais plus significative dans la tranche superficielle (< 75 m; [START_REF] Rouleau | Aperçu de l'hydrogéologie en socle précambrien au Québec et des problématiques minières[END_REF]. Ce sont les sédiments meubles mis en place lors du dernier retrait glaciaire qui forment les aquifères les plus productifs et qui comportent l'eau présentant la meilleure qualité [START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF]. L'interaction entre les aquifères rocheux et le réseau hydrographique est toutefois souvent considérée faible (Gleeson et al. 2009). Au Canada, on note près de 30% de la population qui s'approvisionne au sein d'aquifères régionaux [START_REF] Rivera | Overview of the groundwater resources of Canada[END_REF]. Sur le Bouclier canadien, dans certains territoires ruraux et communautés autochtones, c'est pratiquement l'ensemble de la population qui dépend de l'eau souterraine pour son approvisionnement en eau potable [START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF]. Sa qualité peut toutefois considérablement varier entre une eau de très bonne qualité maintenue dans des aquifères granulaires [START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF]) et une eau issue d'aquifères fracturés présentant certaines contaminations naturelles, telles qu'à l'arsenic [START_REF] Bondu | Occurrence of geogenic contaminants in private wells from a crystalline bedrock aquifer in western Quebec, Canada: Geochemical sources and health risks[END_REF]). Dans la perspective de préserver l'eau souterraine et les écosystèmes qui en dépendent, la gestion du territoire et des ressources naturelles (ex. :

granulats, minerai et forêt) doit être adaptée. La délimitation des aquifères granulaires ainsi que la quantification de leur réserve et de leurs interactions avec les eaux de surface sont nécessaires. En Abitibi-Témiscamingue, les aquifères granulaires ont un volume total représentant une épaisseur moyenne de 1,2 m sur un territoire de près de 20 000 km 2 (≈ 24 km 3 ; [START_REF] Nadeau | Stratigraphic sequence map for groundwater assessment and protection of unconsolidated aquifers: A case example in the Abitibi-Témiscamingue region, Québec, Canada[END_REF]. À titre comparatif, dans le contexte géologique des Basses-Terres du Saint-Laurent, plus homogène d'un point de vue de la géométrie des aquifères, une épaisseur moyenne de 2,0 m sur un territoire de près de 1500 km 2 a été estimée (≈ 3 km 3 ; [START_REF] Rivera | Overview of the groundwater resources of Canada[END_REF]. Considérant que le contexte géologique du Bouclier canadien présente une topographie irrégulière, l'acquisition de connaissances sur la dispersion et l'étendue des aquifères est devenue nécessaire afin de proposer des modulations de gestion dans l'exploitation des ressources naturelles [START_REF] Nadeau | A GIS-based approach for supporting groundwater protection in eskers: Application to sand and gravel extraction activities in Abitibi-Témiscamingue, Quebec, Canada[END_REF].

Les projets associés au Programme d'Acquisition de Connaissances sur les Eaux

Souterraines (PACES) du ministère de l'Environnement du Québec ont permis de mettre en évidence la grande diversité des contextes hydrogéologiques dans le sud de la province du Québec [START_REF] Larocque | Results from the Quebec Groundwater Knowledge Acquisition Program[END_REF]. Cette base de connaissance a notamment permis de simuler les effets des changements climatiques sur la ressource en eau souterraine [START_REF] Larocque | A review of simulated climate change impacts on groundwater resources in Eastern Canada[END_REF]. Sur ces territoires, où la majeure partie de la population habite, les informations associées aux PACES s'avèrent essentielles pour assurer la protection et la saine gestion des aquifères, notamment dans la perspective de l'aménagement du territoire. Toutefois, au nord du 49 e parallèle, les connaissances relatives à l'eau souterraine sont encore embryonnaires. Un projet de caractérisation est en cours depuis 2018 dans le secteur de la Rivière Eastmain, dans les Basse-Terre de la Baie-James, mais ce secteur constitue une faible proportion de l'ensemble du milieu boréal du Québec. Il s'agit pourtant d'un territoire susceptible d'héberger une ressource en eau souterraine hautement stratégique et vulnérable face à l'augmentation de l'activité humaine. souterraineeau de surface [START_REF] Messier | Physical and Biological Modifications in Response to La Grande Hydroelectric Complex[END_REF][START_REF] Lepage | Salinity intrusion in the Eastmain River estuary following a major reduction of freshwater input[END_REF][START_REF] Schetagne | Suivi environnemental du complexe La Grande[END_REF][START_REF] Métivier | Monitoring saltwater intrusion in Rupert Bay, Québec, Canada, after the partial diversion of a major tributary[END_REF]. La régulation du débit et l'emmagasinement d'eau au sein de réservoirs réduisent l'amplitude des crues printanières et nivèlent les débits tout au long de l'année [START_REF] Woo | Streamflow hydrology in the boreal region under the influences of climate and human interference[END_REF]. Ces effets se superposent aux instabilités hydroclimatiques, contribuant à la modification accrue de l'hydrologie des bassins versants. Les débits observés au sein des rivières du bassin hydrologique des baies de James, d'Hudson et d'Ungava semblent répondre à des variations hydroclimatiques systématiques [START_REF] Déry | Characteristics and Trends of River Discharge into Hudson, James, and Ungava Bays, 1964-2000[END_REF]2009), mais le manque de suivis à long terme empêche de déceler des tendances significatives [START_REF] Déry | Interannual variability and interdecadal trends in Hudson Bay streamflow[END_REF].

Néanmoins, les résultats obtenus depuis le début du siècle dans cette région y suggèrent une intensification des flux du bilan hydrologique [START_REF] Déry | Interannual variability and interdecadal trends in Hudson Bay streamflow[END_REF]. Cette observation semble cohérente avec les impacts prédits par les scénarios de changement climatique sur le territoire (Hennigs et Bleau 2018). L'électrification du transport et la demande croissante anticipée en énergie (Hydro-Québec 2020) rendent de plus en plus probable le développement de nouveaux projets hydroélectriques dans cette région. L'intérêt pour les matières premières entraîne aussi l'émergence de plusieurs projets miniers (MERN 2020a) au nord-ouest du Québec. L'industrie du bois est aussi mise en valeur et constitue une activité en continuel développement. Il existe un intérêt des instances gouvernementales pour développer l'activité économique sur le territoire (SPN 2020), ce qui risque de fragiliser certains milieux naturels.

Les Basses-Terres de la Baie-James constituent l'Eeyou Istchee, territoire de la Nation Crie. Malgré la faible densité de population sur ce vaste territoire de plus de 412 000 km 2 , les impacts sur l'écoulement de l'eau de surface, notamment par la réduction significative du débit de certaines rivières (≈ 90 % dans le cas de la rivière Eastmain), rend la population de plus en plus dépendante de l'eau souterraine. La modification significative des débits a bouleversé le rôle que détenait certaines rivières dans l'approvisionnement en eau et en nourriture de la population riveraine et dans ses interactions au sein de l'hydrosystème. L'intensification du développement anthropique dans ce secteur justifie la nécessité de mieux intégrer la protection de l'eau souterraine comme un critère supplémentaire à prendre en compte dans le développement durable du territoire. Cette protection de la ressource en eau souterraine cadre notamment avec les préoccupations 1) sociales, par le maintien des habitudes culturelles des populations et 2) environnementales, par la préservation de la qualité des milieux naturels dépendant de l'eau souterraine. À ce jour, le milieu boréal du Québec constitue un des rares écosystèmes planétaires où l'écoulement naturel des cours d'eau au sein de certains bassins versants peut encore être observé [START_REF] Déry | Characteristics and Trends of River Discharge into Hudson, James, and Ungava Bays, 1964-2000[END_REF][START_REF] Grill | Mapping the world's free-flowing rivers[END_REF]. Il existe donc toujours une réelle opportunité d'adapter la planification territoire pour préserver ces milieux naturels et maintenir leurs services écosystémiques.

Les Basses-Terres de la Baie-James sont aussi un territoire marqué par une forte présence de milieux humides [START_REF] Tarnocai | Peatlands of Canada : Natural Resources Canada[END_REF]. Les fonctions hydrologiques de ces tourbières peuvent être diversifiées et dépendantes de leur environnement de mise en place (hydrogéologique, géologique, topographique et hydrologique; [START_REF] Rosa | A graphical approach for documenting peatland hydrodiversity and orienting land management strategies[END_REF]. Bien que ces tourbières soient majoritairement ombrotrophes et alimentées par les apports issus des précipitations, la localisation des milieux minérotrophes alimentés par des apports d'eau souterraine est peu documentée. Sans cette base de connaissance, la préservation de la biodiversité exclusive à ces milieux est compromise. À ce jour, il n'existe pas d'approche intégrée permettant d'identifier spécifiquement les émergences d'eau souterraine des grands bassins versants qui favoriserait l'identification des milieux naturels dépendant de l'eau souterraine.

1.4 Enjeux des études hydrogéologiques régionales L'évaluation du volume et de l'étendue des aquifères à l'échelle régionale (>100 000 km 2 ) est un véritable défi. Considérant l'immensité des bassins versants peu accessibles par voie terrestre, l'acquisition de données sur la ressource en eau souterraine (géochimie de l'eau, niveau d'eau, etc.) est complexe, onéreuse et chronophage. De plus, les outils conventionnels d'investigation de l'eau souterraine, tels que les puits et les piézomètres, fournissent des informations ponctuelles sur la ressource. Selon la complexité du contexte hydrogéologique du secteur étudié, une quantité significative de ces données est nécessaire pour apporter des informations fiables sur la dynamique d'écoulement des aquifères régionaux. Sur de vastes territoires peu instrumentés, il apparaît peu réaliste d'envisager implanter suffisamment de puits et de sites d'observation in situ pour obtenir une couverture exhaustive de la ressource. À cette échelle, peu d'études sont réalisées pour examiner les interactions entre les écoulements d'eau de surface et l'eau souterraine [START_REF] Ward | The evolution and state of interdisciplinary hyporheic research[END_REF]. Des modèles hydrogéologiques régionaux élaborés à l'échelle subcontinentale (1 km × 1 km sur 200 000 km 2 ) ont montré des résultats prometteurs sur la base de données produites à l'échelle globale [START_REF] Sutanudjaja | Large-scale groundwater modeling using global datasets: a test case for the Rhine-Meuse basin[END_REF]. Néanmoins, cette approche demande une haute performance des systèmes informatiques et un haut niveau de compétence technique.

La capacité de superposer des données spatialisées de sources diverses dans les systèmes d'informations géographiques (SIG) fournit des informations pertinentes et complémentaires pour identifier des réserves en eau exploitables, délimiter des secteurs favorables à la recharge artificielle, suivre des panaches de contaminations, établir des plans de protection, définir les directions d'écoulement, délimiter les zones de recharge et analyser et suivre la ressource en eau [START_REF] Jha | Groundwater management and development by integrated remote sensing and geographic information systems: prospects and constraints[END_REF]). Dans le domaine de l'aménagement du territoire, les SIG sont aussi largement utilisés, suggérant que le développement d'une approche basée sur des analyses spatiales pour quantifier et localiser les réserves en eaux souterraines serait adaptée. Elle permettrait aussi de fournir un haut niveau de détails malgré le vaste territoire d'étude [START_REF] Malczewski | GIS and Multicriteria decision analysis[END_REF]).

1.5 Objectifs de la thèse Aujourd'hui, l'information disponible sur l'interaction entre l'eau souterraine et l'eau

Abstract

This study focuses on the development of two GIS-based approaches that are used jointly to evaluate the groundwater resources associated with granular aquifers in shield environments. The first approach is a multi-criteria analysis (MCA) using an analytical hierarchic process (AHP) based on geological and hydrogeological data for ranking the probability of finding readily available groundwater resources in a specific territory.

The second approach relies on GIS-based geometric calculations that were developed for evaluating the extent and volume of aquifers. The approaches are applied on a 100 × 100 m grid in a 185,000-km 2 area corresponding to watersheds of the James Bay area in Quebec, Canada. The MCA-AHP approach revealed that the unconfined granular aquifers that present the highest aquifer potential (AP) are sparsely distributed and mostly associated with glaciofluvial formations such as the Harricana and Sakami moraines. The geometric calculations approach allowed for estimating that the total volume of groundwater stored in the unconfined granular aquifers reaches approximately 40 km 3 . When used jointly, the two approaches reveal that the shallow unconfined aquifers that require increased groundwater protection account for approximately 5% of the territory. In areas of confined granular aquifers, the highest APs are located in river valleys and lowlands. A sensitivity analysis conducted on the MCA-AHP approach revealed that the grid size does not significantly affect the results.

Therefore, the approach was expanded northward, to a 490,000-km 2 territory reaching the Ungava Bay area. The proposed method could be adapted and applied in other shield areas.

Introduction

Groundwater accounts for approximately 48% of the freshwater on Earth [START_REF] Oki | Global Hydrological Cycles and World Water Resources[END_REF]. It represents the largest reservoir of accessible freshwater for human water supply. Approximately 36% of the world population relies on groundwater for domestic purposes [START_REF] Döll | Impact of water withdrawals from groundwater and surface water on continental water storage variations[END_REF][START_REF] Taylor | Ground water and climate change[END_REF]). Nevertheless, the water renewal rates in groundwater flow systems are generally lower than in most surface water systems, and human perturbations on groundwater quantity and quality can result in long-term impacts [START_REF] Döll | Impact of water withdrawals from groundwater and surface water on continental water storage variations[END_REF][START_REF] Foster | Groundwater-a global focus on the 'local resource[END_REF]Famiglietti 2004;[START_REF] Burri | A review of threats to groundwater quality in the Anthropocene[END_REF][START_REF] De Graaf | Environmental flow limits to global groundwater pumping[END_REF]. Evaluating the extent of aquifers and the volume of water stored in regional groundwater flow systems is a challenging task because most conventional tools only allow for either point-specific investigations (e.g., wells and piezometers) or indirect measurements (e.g., geophysics). In addition, there is currently no remote sensing approach that provides precise water table depth or total groundwater storage estimates [START_REF] Mccabe | The future of Earth observation in hydrology[END_REF]. This further complicates the use of numerical modeling for simulating groundwater flow at the regional scale.

As a result, while groundwater protection measures are critically needed, their definition is particularly complex due to a lack of knowledge and data.

In vast watersheds of the Canadian Shield, water scarcity is generally not observed [START_REF] Oki | Global Hydrological Cycles and World Water Resources[END_REF]Kanae 2006, Döll 2009;[START_REF] Mekonnen | Four billion people facing severe water scarcity[END_REF]. However, major hydroelectric installations have resulted in modifications of surface water flows, especially in the James Bay area [START_REF] Dynesius | Fragmentation and Flow Regulation of River Systems in the Northern Third of the[END_REF][START_REF] Nilsson | Fragmentation and Flow Regulation of the World&#039[END_REF].

Inter-basin transfers led to reduced flow in some large rivers, therefore triggering saltwater intrusions within rivers, significant changes in erosion-sedimentation processes in rivers and degradation of water quality and associated impacts on freshwater ecosystems. Such changes are currently observed in large rivers of the James Bay area, such as the Rupert and Eastmain rivers [START_REF] Messier | Physical and Biological Modifications in Response to La Grande Hydroelectric Complex[END_REF][START_REF] Lepage | Salinity intrusion in the Eastmain River estuary following a major reduction of freshwater input[END_REF][START_REF] Schetagne | Suivi environnemental du complexe La Grande[END_REF][START_REF] Métivier | Monitoring saltwater intrusion in Rupert Bay, Québec, Canada, after the partial diversion of a major tributary[END_REF]. People of the Eeyou Istchee, the land of Cree First Nation in northwestern Quebec, experience these changes and are thus increasingly dependent on groundwater for drinking water supply.

In addition, given the ever-growing demand for natural resources, forestry and mining activities are expected to increase in Eeyou Istchee (MERN 2020a). This underlies the critical need to include groundwater protection as a criterion in land management.

For hundreds of years, Cree's have been aware of the groundwater «Moosh-june-behouk» in the regions they occupied. They possess key knowledge in regional hydrogeology, but this knowledge is not currently available in the scientific literature.

While part of the knowledge required for groundwater protection might stem from local-scale studies associated with specific problems, in the absence of a robust largescale spatial distribution and extent of aquifers, the sustainable development and establishment of efficient protection measures could be severely compromised. Among others, groundwater resources in southern Quebec are increasingly documented owing to the Groundwater Knowledge Acquisition Program, supported by the Quebec Ministry of the Environment [START_REF] Larocque | Results from the Quebec Groundwater Knowledge Acquisition Program[END_REF]. Global studies further allow for comparing the state of water resources in various regions of the world [START_REF] Fan | A simple hydrologic framework for simulating wetlands in climate and earth system models[END_REF][START_REF] De Graaf | A highresolution global-scale groundwater model[END_REF][START_REF] Gleeson | The global volume and distribution of modern groundwater[END_REF]). In such studies, however, the investigation mesh size is generally too coarse to draw conclusions about groundwater flow systems in specific watersheds of the Canadian Shield. In addition, the optimal approaches required for mapping groundwater resources must consider the differences between shields, platforms and orogens and the landscape features and unconsolidated deposits inherited from the most recent glaciation-deglaciation cycles. Therefore, new strategies aimed at evaluating groundwater resources in specific geologic provinces, such as shield environments of the northern hemisphere, could improve regional and global assessments of groundwater resources.

In this context, geographic information system (GIS)-based approaches stand as promising tools for improving groundwater protection during land management. The potential of GIS-based approaches as tools that can be used in MCA has already been proven for natural resources management [START_REF] Malczewski | GIS-based multicriteria decision analysis: a survey of the literature[END_REF][START_REF] Huang | Multi-criteria decision analysis in environmental sciences: Ten years of applications and trends[END_REF]) and specifically applied to groundwater [START_REF] Jha | Groundwater management and development by integrated remote sensing and geographic information systems: prospects and constraints[END_REF][START_REF] Hajkowicz | A Review of Multiple Criteria Analysis for Water Resource Planning and Management[END_REF][START_REF] Joerin | Using multicriteria decision analysis to assess the vulnerability of drinking water utilities[END_REF]Diaz-Alcaide and Martinez-Santos 2019). Nevertheless, the normalization and weighting of criteria included in MCA are common procedures in index calculations, and this could lead to limitations associated with subjectivity and collinearity between criteria. Such limitations can be partly compensated by using an analytical hierarchic process (AHP). This approach provides objective calculations to evaluate the weights of criteria based on their relative importance degree assigned in a pairwise comparison by an expert panel or individual user [START_REF] Saaty | The Analytic Hierarchy Process[END_REF]. The MCA-AHP approach has already been tested for the development of groundwater potential indices in fractured bedrock aquifers (Murthy and Mano 2009;[START_REF] Ta | Cartographie des eaux souterraines en milieu fissuré par analyse multicritère[END_REF][START_REF] Machiwal | Assessment of Groundwater Potential in a Semi-Arid Region of India Using Remote Sensing, GIS and MCDM Techniques[END_REF][START_REF] Adiat | Assessing the accuracy of GIS-based elementary multi criteria decision analysis as a spatial prediction tool -A case of predicting potential zones of sustainable groundwater resources[END_REF][START_REF] Agarwal | Delineation of groundwater potential zone: An AHP/ANP approach[END_REF][START_REF] Mallick | Geospatial and geostatistical approach for groundwater potential zone delineation[END_REF][START_REF] Yao | Cartographie du potentiel en eau souterraine du bassin versant de la Lobo (Centre-Ouest, Côte d'Ivoire): approche par analyse multicritère[END_REF][START_REF] Aladejana | Testing the ability of an empirical hydrological model to verify a knowledge-based groundwater potential zone mapping methodology[END_REF][START_REF] Akinlalu | Application of multi-criteria decision analysis in prediction of groundwater resources potential: A case of Oke-Ana, Ilesa Area Southwestern, Nigeria[END_REF][START_REF] Singh | Assessing the accuracy of GIS-based Multi-Criteria Decision Analysis approaches for mapping groundwater potential[END_REF][START_REF] Arulbalaji | GIS and AHP Techniques Based Delineation of Groundwater Potential Zones: a case study from Southern Western Ghats, India[END_REF][START_REF] Hr | Delineation of groundwater potential zones using remote sensing (RS), geographical information system (GIS) and analytic hierarchy process (AHP) techniques: a case study in the Leylia-Keynow watershed, southwest of Iran[END_REF][START_REF] Saranya | Groundwater potential zone mapping using analytical hierarchy process (AHP) and GIS for Kancheepuram District, Tamilnadu, India[END_REF][START_REF] Saravanan | Delineation of groundwater potential zone using analytical hierarchy process and GIS for Gundihalla watershed, Karnataka, India[END_REF]. With a view to continuing to improve these methods, the development of complementary and independent approaches allowing direct comparisons over large areas as well as the analysis of the sensitivity of MCA-AHP analyzes appear to be relevant avenues.

Fitting in the preestablished context, the objective of this study is to provide an evaluation of groundwater resources in vast watersheds of the James Bay area in Northern Quebec (Canada), within the Canadian Shield. The specific objectives are (1) to develop broadly applicable GIS-based approaches for evaluating groundwater resources in granular aquifers of shield areas, and (2) to use the aforementioned tools to propose recommendations for groundwater protection. A GIS-based MCA-AHP is developed to calculate an aquifer potential index that can be mapped at the regional scale. The results from the MCA-AHP calculations are compared to a geometrical approach that allows for estimating groundwater volumes in unconsolidated deposits as an alternative to the lack of validation data in large and heterogeneous watersheds.

These two approaches are used jointly to identify key areas for groundwater protection as a tool for land management. A sensitivity analysis also is further performed to assess impacts of subjectivity in the MCA-AHP calculations.

Study area

General context

The study area is in northwestern Quebec (Canada) and covers approximately 185,000 km 2 (Figure 2 The region is part of the Superior Geological Province, a component of the Canadian Shield (Figure 2.1b). It covers a large portion of the Abitibi greenstone belt, an assemblage of volcano-sedimentary and plutonic rocks known for their mining potential. The area is coveted for its natural resources (mainly forestry and mining), its hydroelectric potential and the presence of vast unconfined granular aquifers, which are recognized for groundwater quantity and quality [START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF][START_REF] Nadeau | A GIS-based approach for supporting groundwater protection in eskers: Application to sand and gravel extraction activities in Abitibi-Témiscamingue, Quebec, Canada[END_REF][START_REF] Melcc | Suivi hydrologique de différentes stations hydrométriques[END_REF]. The area is characterized by a cold and humid continental climate, with annual average temperatures near the freezing point and cumulative annual precipitation in the form of snow and rain reaching approximately 1000 mm/year.

Glacial history

The geomorphological features of the landscape in northwestern Quebec are largely inherited from the last Pleistocene-Holocene glaciation-deglaciation cycle (see Hardy 1982a;[START_REF] Veillette | New evidence for northwestward glacial ice flow, James Bay region, Quebec[END_REF][START_REF] Veillette | The spectacular cross-striated outcrops of James Bay, Quebec[END_REF][START_REF] Veillette | Ice-flow evolution of the Labrador Sector of the Laurentide Ice Sheet: a review, with new evidence from northern[END_REF][START_REF] Veillette | Uncovering the hidden part of a large ice stream of the Laurentide Ice Sheet, northern Ontario, Canada[END_REF][START_REF] Dubé-Loubert | Glacial and nonglacial events in the eastern James Bay lowlands, Canada[END_REF]. Glacial and postglacial events (Figure 2 Another striking feature is the Sakami moraine, extending over 630 km, and representing deposits that accumulated during a stagnation of the ice front as it receding to the east (Hardy 1982b). The landscape is also characterized by several vast glaciofluvial formations (eskers; Figure 2.2c). Their crests are generally aligned according to the orientation of the ice flow. Some of these formations are recognized for constituting vast aquifers that supply large quantities of high-quality groundwater [START_REF] Nadeau | A GIS-based approach for supporting groundwater protection in eskers: Application to sand and gravel extraction activities in Abitibi-Témiscamingue, Quebec, Canada[END_REF][START_REF] Nadeau | Stratigraphic sequence map for groundwater assessment and protection of unconsolidated aquifers: A case example in the Abitibi-Témiscamingue region, Québec, Canada[END_REF][START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF][START_REF] Rey | Using water stable isotopes for tracing surface and groundwater flow systems in the Barlow-Ojibway Clay Belt, Quebec, Canada[END_REF]. Melt waters from the receding glaciers constrained by topographic features and the ice margin allowed for the development of the proglacial Lake Barlow-Ojibway [START_REF] Veillette | Evolution and paleohydrology of glacial Lakes Barlow and Ojibway[END_REF]. The lake waters covered a large portion of the study area until their final drainage into the James/Hudson Bay outlets, approximately 7700 years ago [START_REF] Clarke | Paleohydraulics of the last outburst flood from glacial Lake Agassiz and the 8200BP cold event[END_REF]). Finegrained, deep-water glaciolacustrine sediments were deposited in this lake (Figure 2.2d). These silt and clay sediments now cover vast areas of the lowlands where they mask older unconfined sediments and the rugged bedrock surface. The Barlow- Ojibway clay belt is punctured by several bedrock hills and glaciofluvial formations found at higher elevations. Littoral sands were also deposited in Lake Barlow-Ojibway, and such units were reshaped by erosion following lake drainage (Figure 2.2e). The uppermost layer of unconsolidated deposits consists of organic material (peat) that has accumulated since the deglaciation period (Figure 2.2f). Organic deposits often form vast peatland complexes that are likely to play a key role in the regional surface and groundwater flow systems [START_REF] Rosa | A graphical approach for documenting peatland hydrodiversity and orienting land management strategies[END_REF]. During the glacial retreat, ice stream events associated with the Cochrane surge reshaped the surface of fine-grained, deepwater glaciolacustrine sediments in the western part of the region and deposited thin layers of clay glacial sediments (Hardy 1982a;[START_REF] Veillette | Uncovering the hidden part of a large ice stream of the Laurentide Ice Sheet, northern Ontario, Canada[END_REF]Figure 2.2g).

Following the drainage of Lake Barlow-Ojibway, the continent was strongly affected by glacial subsidence. This subsidence allowed for the Tyrrell Sea submersion, which reached present-day elevations ranging between 198 m along the Harricana River [START_REF] Hardy | La déglaciation et les épisodes lacustre et marin sur le versant québécois des basses terres de la baie de James[END_REF] and 274 m along the Eastmain River [START_REF] Brousseau | Caractérisation des dépôts meubles et reconstitution paléogéographique quaternaire de la région du réservoir Eastmain[END_REF]. Fine-grained, deep-water marine sediments of the Tyrrell Sea were deposited over a vast portion of the James Bay lowlands (Figure 2.2h). The isostatic rebound allowed for a gradual migration of the coastline to its current position. The coastal, deltaic and beach sediments deposited when the Tyrrell Sea waters receded now form terraces that are incised by rivers [START_REF] Hardy | La déglaciation et les épisodes lacustre et marin sur le versant québécois des basses terres de la baie de James[END_REF][START_REF] Hétu | La plaine alluviale de la rivière Eastmain, Jamésie (Québec): Stratigraphie, sédimentologie, paléoécologie[END_REF]Figure 2.2i).

Data sources

Digital elevation model (DEM)

The digital elevation model (DEM) used here consists of a 100 m x 100 m aggregation (6841 pixels × 5491 pixels) of the 10 m x 10 m DEM from Quebec's Ministry of Energy and Natural Resources (MERN 2020b). At the regional scale, this resolution provides a compromise between rapid computer calculations and intelligible visual representations. The precisely known elevations of 114 geodetic markers were used to evaluate the precision of the DEM. The selected DEM presents an average absolute error of 4.7 m with a standard deviation of 3.9 m. The precision of the DEM was further compared with other elevation models available for the entire study area, and the results are shown in the supplementary material provided in Appendix A-1.

Surficial deposits

There is currently no comprehensive, high-resolution map of surficial geology available at the scale of the study area. The data from the third ecoforestry inventory (MRNF 2009) and from the surficial deposits of Northern Quebec (MFFP 2016) were therefore used jointly. These inventories include data related to the surficial units at the 1: 20,000 scale. A simplified and homogenized surficial geology map was built using these data, after aggregation of surface units into seven main classes based on their texture (Table 2.1). The texture association for all surficial units included in the ecoforestry inventories are detailed in the supplementary material provided in Appendix A-2. It includes lakes, reservoirs and ponds; and the streams and rivers are structured and oriented as a function of water flow.

Boreholes and geology

As of 2020, the Geomining Information System (SIGEOM) from the Quebec Ministry of Energy and Natural Resources (MERN 2018) contained 80,045 diamond-drilling boreholes that are included within the study area. Among these, 73,019 boreholes provide information related to unconsolidated sediments thickness. Most of these exploration boreholes are located along major geological structures that present confirmed or potential mineralization. In addition, the Quebec Ministry of Energy and Natural Resources drilled 169 boreholes specifically aimed at evaluating unconsolidated sediments thickness and stratigraphy [START_REF] Allard | Géochimie de la fraction fine des tills et analyse des matériaux indicateurs de sediments glaciaires et fluvioglaciaires provenant de forages Rotasonic, région de la rivière Octave et Wawagosic (Abitibi) (Geochemistry of the fine fraction of tills and analysis of indicator materials of glacial and fluvioglacial sediments from Rotasonic boreholes, Octave River region and Wawagosic (Abitibi))[END_REF]. These boreholes likely provide the most reliable data within the study area. The SIGEOM database are also the main reliable source for major bedrock faults and bedrock geologic unit.

Methods

Two different GIS-based approaches were developed for evaluating groundwater resources (Figure 2. [START_REF] Saranya | Groundwater potential zone mapping using analytical hierarchy process (AHP) and GIS for Kancheepuram District, Tamilnadu, India[END_REF][START_REF] Saravanan | Delineation of groundwater potential zone using analytical hierarchy process and GIS for Gundihalla watershed, Karnataka, India[END_REF]. The second approach consists of geometric calculations that allow for spatially estimating the saturated thickness of granular aquifers. While the MCA-AHP and geometric approaches are independent, they both allow for assessing groundwater resources at the regional scale. They are ultimately used jointly to identify areas where groundwater protection should be prioritized. Data compilation, treatment, transformation and calculations are performed using ArcGIS 10.5.

The focus is set on aquifers associated with unconsolidated deposits. This choice is explained by the fact that the shallow granular aquifers of the Canadian Shield are generally characterized by a better water quality than fractured bedrock aquifers [START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF][START_REF] Bondu | An exploratory data analysis approach for assessing the sources and distribution of naturally occurring contaminants (F, Ba, Mn, As) in groundwater from southern Quebec (Canada)[END_REF]. Shallow granular aquifers further represent a more readily accessible water resource for local communities. In addition, these aquifers are known to have a high vulnerability and are likely to be associated with groundwater dependent ecosystems (GDE; [START_REF] Nadeau | A GIS-based approach for supporting groundwater protection in eskers: Application to sand and gravel extraction activities in Abitibi-Témiscamingue, Quebec, Canada[END_REF][START_REF] Rosa | A graphical approach for documenting peatland hydrodiversity and orienting land management strategies[END_REF].

Unconsolidated sediments thickness

While borehole data are substantial (over 73,000 boreholes), their spatial distribution is inadequate for interpolating the thickness of unconsolidated sediments at the scale of the study area using conventional interpolation methods as kriging or inverse distance weighting. The method developed by [START_REF] Nadeau | Stratigraphic sequence map for groundwater assessment and protection of unconsolidated aquifers: A case example in the Abitibi-Témiscamingue region, Québec, Canada[END_REF] was adapted to overcome this limitation. It allows for estimating unconsolidated sediments thickness based on the distance from the closest bedrock outcrop. The topography of the outcropping bedrock is also used as an analogue to predict the shape of the buried bedrock surface. In other words, the surface roughness and slope of the buried bedrock is likely to be similar to the surface roughness and slope of nearby bedrock outcrops.

The method used to evaluate the surface roughness of bedrock and till outcrops is presented in Figure 2.4. For each bedrock/till outcrop, the variations in ground elevation are calculated as a function of the distance from the limit of the outcrop with surrounding sediments (Figure 2.4a). Similarly, the thickness of unconsolidated sediments with respect to the distance from the closest bedrock/till outcrop can be evaluated using the available borehole data (Figure 2.4b). The trends evaluated using the two aforementioned approaches can be plotted in a graph presenting the average elevation range with respect to the distance from the bedrock/sediment limit

While this approach remains approximate, it allows for performing calculations that are coherent with field data and observations. The trend obtained by [START_REF] Nadeau | Stratigraphic sequence map for groundwater assessment and protection of unconsolidated aquifers: A case example in the Abitibi-Témiscamingue region, Québec, Canada[END_REF] for the Abitibi-Témiscamingue region (also within the Canadian Shield) is shown for comparison in Figure 2.4c.

Groundwater flow conditions: confined vs unconfined aquifers

The delimitation of unconfined vs confined aquifers relies on an evaluation of the extent of (1) fine-grained, deep-water sediments (either glaciolacustrine or marine) and

(2) clayey till deposits associated with the Cochrane surge. These sediments are considered regional aquitards, given the widely recognized low hydraulic conductivity of clay with respect to sand and gravel. Here, water-saturated sand and gravel units found beneath a clay layer are considered confined aquifers. Nevertheless, the aquitards are covered in places by younger deposits. There is a need, therefore, to evaluate the extent of clay deposits beyond the areas where they are observed at the surface. In the absence of precise stratigraphic data that could be interpolated at the scale of the study region, the extent of clay deposits was estimated based on geological observations and criteria based on the knowledge of the glacial history (section 2.2.1 General context): Based on these conditions, an evaluation of water depths in Lake Barlow-Ojibway and in the Tyrrell Sea can be used to identify areas where conditions were conductive to the accumulation of clayey sediments (using a 100 m x 100 m grid). The maximum elevation reached by Lake Barlow-Ojibway was evaluated by interpolating the isolines of highest washing limits previously identified by [START_REF] Veillette | Evolution and paleohydrology of glacial Lakes Barlow and Ojibway[END_REF], and restricted to the west side of the Sakami moraine. Similarly, the maximum elevation reached by the Tyrrell Sea was estimated using a linear interpolation of the elevation of ancient shorelines as identified by [START_REF] Hardy | La déglaciation et les épisodes lacustre et marin sur le versant québécois des basses terres de la baie de James[END_REF] and Brosseau (2008) along the Harricana and Eastmain Rivers, respectively. This approach allows accounting for the apparent tilt of the lake and sea paleo-surfaces, the latter resulting from the non-uniform isostatic rebound that occurred in the territory following deglaciation.

Multi-criteria analysis: Unconfined aquifer potential (AP)

The first approach used to calculate the AP for unconfined areas consists in applying a MCA on a 100 m x 100 m grid:

AP = ∑ 𝑤 𝑗 • 𝑥 𝑖𝑗 (2)
where 𝑤 𝑗 is the weight defined for criterion 𝑗 and 𝑥 𝑖𝑗 is the value assigned to component 𝑖 for criterion 𝑗 following standardization. The potential transmissivity (𝑇'), the surface slope (𝑆), the distance to a feature of the hydrographic network (DFHN), the drainage network density (DND), the distance to a major bedrock fault (DMBF) and the topographic wetness index (TWI) are the applicable criteria selected based on previous studies (Diaz-Alcaide and Martinez-Santos 2019). Each criterion is evaluated for each cell (100 m x 100 m) of the territory based on the available data (section 'Data sources'). These cell values are then standardized on a scale of 1 to 5 based on theoretical assumptions, higher values being attributed to conditions that are more conductive to the development of aquifers. The criteria and associated standardized values are provided in Table 2.2. When applied for the six criteria, Eq. 2 can be developed as Eq. 3 to express the AP index:

AP = 𝑤 𝑇′ • 𝑥 𝑇′ + 𝑤 𝑆 • 𝑥 𝑆 + 𝑤 DFHN • 𝑥 DFHN + 𝑤 DND • 𝑥 DND + 𝑤 DMBF • 𝑥 DMBF + 𝑤 TWI • 𝑥 TWI (3)
The AP is a dimensionless qualitative aggregation of physical parameters which provides information on storage and filling capacity. It was assumed that the criteria are independent but a bias could be introduced by a potential collinearity between them.

A sensitivity analysis was conducted to identify the impact of modifications in input parameters. The map of regional distribution of the AP is presented in section 2.5.3.1. The transmissivity (𝑇) [L 2 /T] of an aquifer corresponds to the product of its saturated thickness (b) [L] and hydraulic conductivity (𝐾) [L/T]:

𝑇 = 𝐾 • 𝑏 (4)
Here, in the absence of a precise knowledge of the hydraulic heads and hydraulic conductivities at the scale of the study area, a quantitative evaluation of 𝑇 is not feasible. To overcome this limitation, a semi-quantitative evaluation of the potential transmissivity (𝑇') is proposed:

𝑇′ = 𝑥 TEXT • 𝑥 UST (5)
Where 𝑥 UST is the value of the standardized UST and 𝑥 TEXT is the value of the standardized surficial units textures (TEXT) based on their expected hydraulic conductivity (Table 2.3). For a given TEXT, the areas characterized by higher UST are more conductive to the development of a significant aquifer in comparison to areas characterized by lower UST [START_REF] Nadeau | Stratigraphic sequence map for groundwater assessment and protection of unconsolidated aquifers: A case example in the Abitibi-Témiscamingue region, Québec, Canada[END_REF]. Then, typical K values measured in unconsolidated sediments from the Canadian Shield [START_REF] Cloutier | Projet d'acquisition de connaissances sur les eaux souterraines de l'Abitibi-Témiscamingue (partie 1) [Abitibi-Témiscamingue groundwater knowledge acquisition project[END_REF][START_REF] Cloutier | Projet d'acquisition de connaissances sur les eaux souterraines de l'Abitibi-Témiscamingue (partie 2) [Abitibi-Témiscamingue groundwater knowledge acquisition project[END_REF][START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF] [START_REF] Domenico | Physical and Chemical Hydrogeology[END_REF].

Soil surface slopes (S)

The slope of soil surface is likely to impart a significant control on the runoff vs infiltration ratio but their physical effect is still debated [START_REF] Chen | Green-Ampt infiltration model for sloping surfaces[END_REF]. Here, it is assumed that steeper slopes are likely to be more prone to surface runoff, while flatter grounds are likely to be more prone to infiltration and subsequent groundwater recharge [START_REF] Aller | DRASTIC : A standardized system for evaluating ground water pollution potential using hydrogeologic settings[END_REF][START_REF] Brito | Characterization of maximum infiltration areas using GIS tools[END_REF][START_REF] Pomeroy | Fuzzy logic-based assessment for mapping potential infiltration areas in lowgradient watersheds[END_REF]. The surface slopes (in degrees) were calculated from the DEM and their standardized classes are inspired for the slope parameter as define in the DRASTIC vulnerability index [START_REF] Aller | DRASTIC : A standardized system for evaluating ground water pollution potential using hydrogeologic settings[END_REF].

Hydrography and drainage density (DFHN and DND)

In areas of unconfined aquifers, the presence of lakes and ponds is considered an indication of a shallow water [START_REF] Rouleau | Aperçu de l'hydrogéologie en socle précambrien au Québec et des problématiques minières[END_REF], major fault zones are likely to be associated with an increased fracture density, a condition that is conducive to increased permeability and sustainability of water supply [START_REF] Sander | Lineaments in groundwater exploration: a review of applications and limitations[END_REF]. In contrast, a fault zone could act locally as a barrier to groundwater flow and reduce bedrock permeability (Gleeson and Novakowski 2009).

Topographic wetness index (TWI)

The topographic wetness index (TWI) allows for a relative assessment of the potential water saturation of soils [START_REF] Beven | A physically based, variable contributing area model of basin hydrology / Un modèle à base physique de zone d'appel variable de l'hydrologie du bassin versant[END_REF]Kirkby 1979, Sorensen et al. 2006). The TWI can yield mitigated results in some environments, but it provides a relevant tool for evaluating regional patterns (Tompson and Moore 1996;[START_REF] Rodhe | Wetland occurrence in relation to topography: a test of topographic indices as moisture indicators[END_REF][START_REF] Grabs | Modeling spatial patterns of saturated areas: A comparison of the topographic wetness index and a dynamic distributed model[END_REF]). The index is based on the hypothesis that soil moisture increases as the drainage area (𝛼) increases and surface slope (𝛽) decreases:

TWI = 𝑙𝑛 ( 𝛼 𝑡𝑎𝑛(𝛽) ) (6)
The TWI calculation is sensitive to the DEM used for calculating 𝛽 and to the method used for evaluating 𝛼 [START_REF] Mattivi | TWI computation: a comparison of different open source GISs[END_REF]. Here, the wetness index from Saga-GIS software (Bohner and Selige 2006) is used because it is better suited for areas where surface slopes are weak.

2.4.3.6 Criteria weighting using the analytical hierarchic process (AHP)

The MCA allows for the evaluation of AP according to weighted criteria. The criteria (Table 2.2) are calculated using an AHP. An importance degree is assigned to each pair of criteria (Table 2.4), either by an expert panel or by an individual user, in a pairwise comparison matrix (Table 2.5). The weight (𝑤) attributed to criterion (j) corresponds to the eigenvector of this matrix, as calculated using a procedure adapted from [START_REF] Malczewski | GIS and Multicriteria decision analysis[END_REF]:

Step 1 -Summation of values of each column of the pairwise comparison matrix;

Step 2 -Creation of a standardized table by dividing each component of the table by the corresponding summation (as evaluated in step 1, above);

Step 3 -Division of the sum of normalized values by the number of criteria for calculating the average of components in each line. One of the strengths of this approach is the ability to calculate a coherence index (CI)

to verify if degrees of importance attributed in the pairwise comparisons are coherent with each other:

CI = 𝜆-𝑛 𝑛-1 (7)
where n is the number of included factors and 𝜆 is the average of the coherence vector, which is evaluated by:

Step 1 -Multiplying the first column of the original comparison table by the weight of the first criteria (and so on);

Step 2 -Calculating the sum of each line;

Step 3 -Dividing each sum by the weight of the associated criteria.

For each CI, a consistency ratio (CR) can be evaluated using a random index (RI) associated with the number of included factors (for 𝑛 = 6, RI = 1.24; Saaty 1980):

CR = CI RI (8)
Consistency ratios ≤ 10 % indicate a coherence in the weights attributed to criteria, whereas 𝐶𝑅 values > 10% suggest that the criteria weighting procedure should be revised [START_REF] Saaty | The Analytic Hierarchy Process[END_REF]. The results of the AHP conducted here to assign weights and their coherence ratio are presented in Table 2.5. A specific flowchart of the main step of the MCA-AHP is provided in the Appendix A-3 and a spreadsheet detailing the calculations is provided in electronic supplementary materials. The second approach used here is based on geometric calculations that allow for estimating the saturated thickness of unconfined granular aquifers. The available wells and piezometers where hydraulic heads can be measured do not allow for interpolating the position of the water table at the scale of the study area, mainly because measurements are sparse and unevenly distributed. As an example, approximately 2240 groundwater wells with information on groundwater depth are included in Quebec's Hydrogeologic Information System for the study area. Nevertheless, the vast majority of these wells are highly concentrated in the municipalized area found at the southern limit of the study area. These data could therefore not be used to interpolate groundwater depth due to their poor spatial distribution. Nevertheless, in areas of unconfined granular aquifers, it can be postulated that the elevation of streams, rivers and lakes provide local information on the position of the water table. Under the hydroclimatic conditions prevailing in the study area, it is most likely that unconfined granular aquifers discharge into streams and lakes. This was indeed observed in other areas of the Canadian Shield located south of the study region [START_REF] Nadeau | A GIS-based approach for supporting groundwater protection in eskers: Application to sand and gravel extraction activities in Abitibi-Témiscamingue, Quebec, Canada[END_REF][START_REF] Nadeau | Stratigraphic sequence map for groundwater assessment and protection of unconsolidated aquifers: A case example in the Abitibi-Témiscamingue region, Québec, Canada[END_REF][START_REF] Rey | Using water stable isotopes for tracing surface and groundwater flow systems in the Barlow-Ojibway Clay Belt, Quebec, Canada[END_REF] and suggests that in the immediate vicinity of lakes and streams, groundwater heads in unconfined aquifers are likely to be equal (or higher) to the elevation of surface waters.

The water table elevation (WTE) was therefore estimated using a triangulated irregular network (TIN) created from points (vertices) of the hydrographic network. In unconfined aquifer areas, all lakes in contact with granular deposits were selected. To ensure the continuity and the representativeness of the WTE, lakes with areas > 10 km 2 were selected in sectors of confined aquifers. The streams and rivers with a Strahler order > 4 were also selected in confined aquifer areas, while all persistent streams were selected in unconfined aquifer areas. The elevations of the vertices defining the contour of water bodies (lakes and large rivers) and the lines of selected streams were extracted from the DEM and then used to produce a TIN. The approximate WTE was obtained by converting this TIN in a 100 m × 100 m grid that is coherent with the DEM. A subsequent correction was made to restrict WTE to values below the surface elevation.

The WTE allows for estimating groundwater depth (GD) (Eq. 9) and the saturated sediments thickness (SST) (Eq. 10) in areas of unconfined aquifers:

GD = DEM -WTE (9) SST = UST -GD ( 10
)
where DEM is the value of the digital elevation model (section 'Digital elevation model') and UST is the value of unconsolidated sediment thickness (section 'Unconsolidated sediments thickness'). All these parameters have units of length [L].

Eq. 10 allows for estimating the extent of unconfined granular aquifers along the z-axis, while the extent of these aquifers in the xy-plane can be estimated based on the data associated with groundwater flow conditions (section 'Groundwater flow conditions: confined vs unconfined aquifers' ). This process further provides evaluation of the volume of unconfined granular aquifers.

Confined aquifer delimitation

In confined aquifer areas, the AP is assumed to be proportional to the UST. Under the hydro-climatic conditions prevailing in the study area, it is most likely that the sand and gravel units found beneath aquitards are water-saturated [START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF][START_REF] Nadeau | Stratigraphic sequence map for groundwater assessment and protection of unconsolidated aquifers: A case example in the Abitibi-Témiscamingue region, Québec, Canada[END_REF]). In the absence of precise stratigraphic data to calculate the thickness of the aquitard layers and the precise extent of confined granular aquifers, the AP were attributed based on the assumption that the highest AP are associated with the highest UST values as observed in neighboring regions [START_REF] Nadeau | Stratigraphic sequence map for groundwater assessment and protection of unconsolidated aquifers: A case example in the Abitibi-Témiscamingue region, Québec, Canada[END_REF]. were further reported in this area by several authors [START_REF] Roy | Insights on the events surrounding the final drainage of Lake Ojibway based on James Bay stratigraphic sequences[END_REF][START_REF] Allard | Constraining the age of the last interglacial-glacial transition in the Hudson Bay lowlands (Canada) using U-Th dating of buried wood[END_REF][START_REF] Dubé-Loubert | Glacial and nonglacial events in the eastern James Bay lowlands, Canada[END_REF]. The UST calculated for the lowlands, near James Bay, presents a greater heterogeneity north of the Rupert River where thicker sediments are observed, mainly along the main river valleys.

Here, Eq. 1 was used to represent the increase in UST with respect to the distance to the closest bedrock outcrop. The approach has the advantage of allowing for the calculation of an average regional trend, as it accounts for all available borehole and bedrock outcrop data. Nevertheless, it is not suited for representing local features of the bedrock morphology, since it is based on average rather than local values. In order to provide a quantitative representation of this limitation, the UST compiled from the most reliable boreholes (n=169; Allard and Dubé-Loubert 2016) were plotted against the trend defined by Eq. 1 (Figure 2.7). The average difference between the predicted vs measured UST is 2.7 m and the mean absolute difference is 11.1 m-a significant value that reflects bedrock irregularity. For each of the distances reported on the x-axis, the predicted trend intersects the range defined by the average UST (measured from boreholes) ± standard deviation (σ). These observations suggest that the calculated UST are relevant when interpreted at the scale of the study area, but do not allow for precise local estimates. While this represents a limitation with respect to the evaluation of available water resources at the local scale, the approach remains relevant for identifying vast areas characterized by high AP. This information is critical for including groundwater protection in the context of land use planning since vast contiguous areas, rather than sparse and isolated small areas, are generally considered.

Groundwater flow conditions and aquifer delimitation

The extent of the fine-grained, deep-water glaciolacustrine and marine sediments (hereafter the "clay plain") is used to differentiate between unconfined and confined granular aquifers. Overall, the clay plain covers approximately 50% of the study area and is restricted to the western part of the study area. Along the Sakami moraine, the ice margin acted as a barrier preventing the submersion of the proglacial Lake Barlow- surrounding the unconfined aquifer are prone to the accumulation of groundwater within the unconfined, unconsolidated sediments. These characteristics are highly similar to that of the most productive unconfined glaciofluvial aquifers of the Abitibi-Témiscamingue region [START_REF] Nadeau | A GIS-based approach for supporting groundwater protection in eskers: Application to sand and gravel extraction activities in Abitibi-Témiscamingue, Quebec, Canada[END_REF][START_REF] Nadeau | Stratigraphic sequence map for groundwater assessment and protection of unconsolidated aquifers: A case example in the Abitibi-Témiscamingue region, Québec, Canada[END_REF][START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF].

The central part of the study region (Zone 2) consists of a transition zone between confined and unconfined aquifers located to the west and east, respectively. The unconfined granular aquifers of this zone are generally characterized by AP values ≤ 3.

The limited thickness of unconsolidated sediments in this zone and the lesser abundance of lakes suggest lower AP. Zone 2 extends towards the north to reach the Eastmain area, near James Bay. Towards the east, the Sakami moraine marks a clear limit between Zone 2 and Zone 3.

The eastern part of the study region is identified as Zone 3 in Figure 2.9. This zone is located outside of the limits of the clay plain. Most of the unconfined aquifers of this zone present AP < 3. The unconfined aquifers presenting the highest AP are largely associated with the Sakami moraine and with glaciofluvial formations located further east. Some wide areas of unconfined aquifers presenting AP = 5 are also found at the northeastern limit of Lake Mistassini (label B in Figure 2.9). In the southern part of Zone 3, the areas characterized by AP = 5 are mainly associated with littoral sands and glaciofluvial formations (label C in Figure 2.9).

Aquifer delimitation based on geometric calculations

The geometric calculations suggest thicknesses of unconfined granular aquifers reaching the maximum value of 26 m (see Eq.1; Figure 2.7), for regional average and standard deviation values of 8.5 m and 6.0 m, respectively. Based on the geometric calculations, it is estimated that 17% of the granular formations are not saturated with water. Generally, the highest aquifer thicknesses correspond to the areas where the highest AP were evaluated based on the MCA-AHP. This can be seen in Figure 2.10, where the saturated sediments thickness (SST) are reported according to the five classes of AP. The results from both approaches appear to be coherent, as the greatest SST (> 15 m) are mainly associated with the highest classes of AP (AP ≥ 4). Similarly, the average calculated SST decreases with decreasing AP. One of the main limitations of the geometric approach is the inability to account for the topography of eskers and moraines, as the elongated crests associated with these formations are not well represented by the regional interpolations. Owing to this limitation, the results suggest that the highest SST are found at the margin of eskers and moraines, with lower SST values along the crests. Nevertheless, based on several observations conducted in similar geologic formations in the Abitibi-Témiscamingue region, it is expected that the highest SST values should be observed in the core of eskers and moraines [START_REF] Nadeau | A GIS-based approach for supporting groundwater protection in eskers: Application to sand and gravel extraction activities in Abitibi-Témiscamingue, Quebec, Canada[END_REF]. The geometric calculations also allow for a first order estimate of the volume of groundwater stored in unconfined granular aquifers. Assuming an average total porosity of 30% for sand and gravel aquifers and uncertainties of approximately

Uncertainties and sensitivity analyses

Several factors can entail uncertainties in the MCA-AHP and the geometric calculations proposed here. One advantage of the MCA-AHP approach is that a sensitivity analysis can be performed, therefore allowing for quantifying the impact of different criteria on AP calculations. The MCA-AHP used for estimating AP includes six criteria that are weighted according to pairwise comparisons (section 'Multi-criteria analysis: unconfined aquifer potential (AP)'). While the approach relies on objective mathematical calculations, the choice of input parameters (mesh size, selected criteria, limits of classes used for standardization) and pairwise comparisons imply a subjective process. A sensitivity analysis was therefore conducted to identify the impact of modifications in the MCA-AHP input parameters on the calculated AP. Modifications were brought to the selected criteria, number and limits of classes, weights and spatial resolution (mesh size) according to 21 scenarios (Table 2.6).

The relative difference between the tested scenarios and the reference state (ref, used for calculations in this study) are evaluated by calculating, for each scenario, the total area associated with each AP class, and normalizing these values to the results obtained for the reference scenario (Figure 2.11). The data reveals that the relative difference between the tested scenarios and the reference state are greatest for AP classes 1 and 5. This was expected since these two classes of AP are the least abundant over the study area. The sensitivity analysis further allows for identifying scenarios that tend to overor under-estimate the abundance of AP within class #5, since this is the target class in this study. The data further support the interpretation that the number of classes used for partitioning input data and the approach used for evaluating the weight of criteria are the key parameters driving results of the MCA-AHP. In addition, the sensitivity analysis revealed that parameter T' has a major influence on calculations. This highlights the effect of the rugged bedrock surface causing zones with low UST and low AP near bedrock outcrops. User defined 5 Jenks 1 The limits of classes associated with distances (DFHN, DMBF) are defined by the user in every scenario. For the other criteria, the limits of classes can be either user defined as describe in Table 2.2 or calculated using natural breaks (Jenks). The scenario Ref corresponds to the MCA-AHP used in the current study.

Here, the reference state is considered most robust based on two main arguments:

1. A higher number of classes for partitioning the data associated with input parameters allows a more precise differentiation of AP values; 2. User-defined weights attributed using pairwise comparisons allow better accounting for the available knowledge on the studied groundwater flow system. Some additional factors, not included in the sensitivity analysis, can further influence AP calculations. Data precision in terms of x-y (longitude-latitude) geolocation stands as one potential issue. In vast areas, low resolution datasets are sometimes used for reducing the time required for computer calculations and the volume of data storage.

When multiple low-resolution data sources are used jointly, some of the errors can addup or cancel each other. Here, the uncertainty on AP values resulting from the potential errors in data sources was not quantified in the sensitivity analysis. These errors could occur if a feature of an input data source is not precisely geolocated, and if the x, y error on this feature causes it to be mapped in a grid cell that is adjacent to its exact position. This limitation should be considered when interpreting AP maps. In addition, Figure 2.11 Results of the sensitivity analysis. Each bar corresponds to a scenario, with associated color as identified in Table 2.6. The relative difference is presented in percentage of cells assigned for classes by the tested scenario (in color) compared to the reference state (in black).

the rugged bedrock surface and the uncertainties on the calculated water table elevation used in the geometric approach could yield local errors. In flat areas such as the James Bay lowlands, small errors in the UST and GD interpolations could cause errors in aquifer volume calculations over large areas. The 100 m grid resolution could further entail local scale errors, especially in areas where steep slopes are observed.

Nevertheless, the estimations proposed in this study are based on regional trends built using quantitative data and on hypotheses and assumptions that are with pre-existing scientific literature. Despite their inherent limitations, the approaches proposed here allow for improving the understanding of the regional groundwater flow system in an area where data are sparse and unevenly distributed.

Recommendations for groundwater protection and land management

Considering the limits of the approaches developed, their use to improve land use planning implies certain precautions. Here, it is proposed that the large zones associated with high aquifer potentials most likely represent areas where the groundwater resource is abundant. Such sectors should be given special attention in the context of land use planning. For example, more exhaustive local hydrogeological studies could be requested by governmental authorities if projects implying potential impacts on water resources are proposed in these areas. In the absence of local hydrogeological studies and field validations, it was agreed that the precautionary principle should prevail to protect the water resources.

Groundwater availability, quality, usage and vulnerability, along with the presence of groundwater dependent ecosystems (GDE), are key parameters for identifying aquifer areas where increased protection is required during land management. The results of the MCA-AHP and geometric approaches presented here allow for the identification of areas that present the highest AP and groundwater volumes in unconsolidated aquifers. Despite their inherent uncertainties, the results from the MCA-AHP and geometric approaches are used jointly to identify areas where increased awareness for groundwater protection is needed (Figure 2.12; see also electronic supplementary material). These areas are useful to guide land management since they provide a first large-scale overview of the resource allowing one to better prioritize field investigations for protection purposes. Priority-1 areas are defined as areas with AP = 5 and a SST > 10 m. When only one of these criteria is met, the priority is assumed to decrease. The lowest priority ranking (Priority 4) is assigned where none of the criteria are met (Figure 2.12).

One key feature of Zone 1 is the vast sector of unconfined granular aquifers associated with an outcrop of the Harricana moraine (Figure 2.12a). This zone likely represents a key target for groundwater protection; with some segments identified as Priority 1 in terms of groundwater protection needs. The proximity of the Abitibi greenstone belt, known for its high mining potential (MERN 2020a), is a main concern and requires special attention. Furthermore, this vast unconfined aquifer area likely represents the main recharge area for surrounding confined aquifers. Northwards, the road accesses to this area are extremely limited, suggesting that the current and near future groundwater usage will be very limited over most of the zone. If groundwater protection measures are applied therein, they should be based on the long-term management of the resource for future generations. In terms of short-term groundwater needs, the areas of thick sediments located in the immediate vicinity of access roads and hunting camps should be prioritized for protection. In addition, areas of minerotrophic peatlands (fens) could indicate the presence of GDE. Such areas, especially when spatially associated with high AP values, might correspond to highly esker and Harricana moraine [START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF]) are identified as Priority 1 areas

(Figure 2.12b).
In Zone 2, key areas of confined aquifers that should be prioritized for protection correspond to the main river valleys where large aquifers are found. Protecting such areas could further allow for an increased protection of rivers and associated ecosystems. In Zone 2 and Zone 3, the Sakami moraine stands as a unique geological feature where several areas of unconfined granular aquifers presenting Priority-1 areas are found (Figure 2.12c). Field observations further reveal that groundwater springs are spatially associated with this formation, supporting the hypothesis that GDE are highly likely to be found along the margin of the moraine. The moraine also marks the eastern limit beyond which roads and accesses to the land become extremely limited. Besides, its well-drained characteristics make it the ideal geologic environment for road construction (Seppäla 1999) and sand and gravel extraction [START_REF] Nadeau | A GIS-based approach for supporting groundwater protection in eskers: Application to sand and gravel extraction activities in Abitibi-Témiscamingue, Quebec, Canada[END_REF], and therefore a hotspot for human pressure. For these reasons, the Sakami moraine should stand as a critical target for groundwater protection. The Priority-1 areas located at the northeastern limit of Lake Mistassini also deserve attention (Figure 2.12d). It is noteworthy that this area corresponds to the southern limit of the Otish Mountains, an area with a strong potential for mining. Further environmental assessments in this area should consider groundwater protection as a key aspect of environmental studies, given that it could correspond to one of the greatest unconfined granular aquifers in the area and that it is part of a regional flow system feeding the largest natural lake in Quebec (Canada), making it a unique hydrological feature of Eeyou Istchee. At the regional scale, it is recommended to pay a special attention to Priority-1 areas located near local communities (Figure 2.12) and current house groupings and hunting camps. In such areas, the protection of shallow unconfined aquifers could allow for developing a network of community wells that would become available for a safe access to drinking water. The water quantity and quality commonly observed in shallow unconfined aquifers found in glaciofluvial formations [START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF] would likely offer the ideal conditions for installing and operating groundwater wells.

The definition of groundwater protection zones in such areas seems much needed given the fact that glaciofluvial formations are also targeted for sand and gravel extraction for road construction.

Larger-scale application of the MCA-AHP approach

Because different hydrological processes occur at different spatiotemporal scales, from local to regional groundwater flow systems, the approaches used to study groundwater must be adapted to the scales at which processes are assessed [START_REF] Blöschl | Scalue issues in hydrological modelling: a review[END_REF]. In large areas where data are limited and sparsely distributed, the use of large-scale GIS-based approaches can be useful to target areas where further local-scale data acquisition or increased environmental protection is required. Such approaches are also useful to compare different areas and provide a continental or even global perspective on the groundwater resource [START_REF] Fan | A simple hydrologic framework for simulating wetlands in climate and earth system models[END_REF][START_REF] Gleeson | The global volume and distribution of modern groundwater[END_REF][START_REF] Huscroft | Compiling and Mapping Global Permeability of the Unconsolidated and Consolidated Earth: GLobal HYdrogeology MaPS 2.0 (GLHYMPS 2.0)[END_REF]). Yet, the coherence of regional-scale GIS-based approaches can be affected by the large scale of the territory included in the analyses. Here, an attempt was made to apply the MCA-AHP approach beyond the limits of the study area. The MCA-AHP was thus applied to a 490,000 km 2 area, using a 500 m × 500 m grid and expanding the study area towards the north to include vast watersheds flowing into Ungava Bay (see Figure 2.1 vs Figure 2.13). The datasets used for evaluating the AP values are the same as those presented in section 'Data sources', except for the DEM that came from the Shuttle Radar Topography Mission [START_REF] Farr | The Shuttle Radar Topography Mission[END_REF]). The limit between confined and unconfined aquifers is based on an estimate of the extent of marine clays from post-glacial seas. The Tyrrell Sea limit is at 265 m above sea level (asl) [START_REF] Vincent | Le Quaternaire récent de la région du cours inférieur de la Grande rivière[END_REF], while the Iberville Sea

The dominance of glacial deposits punctuated by numerous bedrock outcrops as well as the rugged surface topography suggest that unconfined aquifers are shallow and of limited extent, therefore limiting AP values. Groundwater reserves in unconsolidated deposits appear to be particularly scarce in the Otish and Torngat Mountains and in the Koksoak River watershed. Dryer climatic conditions, along with increased permafrost continuity, could further reduce AP values north of the La Grande River. The sensitivity analysis suggested that the use of a coarser mesh does not make the approach unsuitable for a regional scale evaluation of aquifer potential. Therefore, it is proposed that the MCA-AHP has the potential for application over greater regions of the Canadian Shield using a coarser mesh.

Conclusions

This study focused on the development of a GIS-based MCA-AHP and a geometric approach for evaluating groundwater resources in a vast area of the Canadian Shield.

Both the MCA-AHP and the geometric approach yielded coherent results within the study region, suggesting that most of the parameters used to target aquifers are relevant.

These approaches allowed for the identification of priority areas for groundwater protection, altogether covering approximately 1500 km 2 . Accounting for these priority groundwater protection areas in land management is critical to ensure safe access to water for future generations and for protecting groundwater dependant ecosystems.

Such protection measures related to groundwater are critically needed, especially in Eeyou Istchee, where surface waters are extensively impacted by hydroelectric installations. To protect aquifers in the studied area, it is necessary to expand the gaze outside of hydrogeology and understand the course modern society is taking Eeyou istchee. By considering the political, business, and individual decisions that were, are and will be made, on how it will impact the land and its water resources. In the recent past, infrastructure developments such landfills, buildings, roads and sand and gravel extraction sites located in groundwater recharge areas and unconfined aquifers might threaten the groundwater resource. It seems hydrogeology considerations were often overlooked in the environmental assessments prior to the developed projects. The data provided here offer an opportunity to better protect aquifers for future generations.

The MCA-AHP seems to stand as a promising approach for creating bridges between coarser global approaches and local investigations. The tests realized here revealed that coherent results can be obtained from the MCA-AHP applied on a mesh size of 100 m vs 500 m, despite the fact that the coarser analysis relies on less precise elevation data.

Even if the granular aquifers are scattered and discontinuous, the observed regional trends provide insights about the spatial distribution of groundwater resources, information that is critically needed to improve land management. This study's findings suggest that the approach has the potential to be applied over other vast areas of the Canadian Shield as well as other shield areas of the northern hemisphere such as the Baltic Shield and Greenland.

CHAPITRE III BILAN HYDROLOGIQUE DE BASSINS VERSANTS EN MILIEU BORÉAL SUR LE BOUCLIER CANADIEN HYDROLOGICAL MASS BALANCE OF BOREAL WATERSHEDS OF CANADIAN SHIELD

Résumé

Cette étude vise à quantifier les composantes du bilan hydrologique de grands bassins versants du Bouclier canadien dans la région de la Baie-James au nord-ouest du Québec, Canada. Les travaux ont été menés sur six rivières du Bouclier canadien drainant au total plus de 185 000 km 2 des écozones du Bouclier boréal, de la Taïga du Bouclier et des plaines hudsoniennes du Canada. Les mesures de débit des rivières, la géochimie de l'eau superficiel (δ 2 H-δ 18 O et conductivité électrique de l'eau), des données de télédétection et des modèles SIG sont utilisés conjointement pour calculer les bilans hydrologiques. L'approche permet de répartir la contribution annuelle de l'eau de pluie, de la fonte des neiges, du ruissellement de surface, de l'évaporation, de la transpiration et du drainage des eaux souterraines sur le débit des grandes rivières. Les résultats suggèrent que le ruissellement des eaux de pluie (30 à 61% des précipitations totales) et la fonte des neiges (18 à 40% des précipitations totales) sont les principales contributions au débit des rivières, tandis que la contribution des eaux souterraines aux rivières représente < 12% des précipitations totales. Les bilans de masse des isotopes stables qui permettent d'estimer l'évaporation à l'échelle du bassin versant montrent des ratios par rapport aux précipitations totales entre 2% et 10% et suggèrent que la transpiration a une composition isotopique proche de l'eau de pluie estivale. De plus, les bilans hydrologiques suggèrent que la réserve totale d'eau stockée dans la partie active des bassins versants représente 10 à 20% des précipitations annuelles totales tandis que les flux d'eaux souterraines au-delà des bassins versants sont négligeables. Les bilans hydrologiques ont été davantage documentés pour fournir des informations sur la saisonnalité des flux des bassins versants. Les observations sont finalement utilisées pour proposer des recommandations pour le suivi des rivières en milieu boréal du Bouclier canadien et pour mieux orienter de futurs travaux de recherche.

Abstract

This study aims at evaluating the hydrological balance of large watersheds of the Canadian Shield in the James Bay area in Northwestern Quebec, Canada. The focus is set on six rivers of the Canadian Shield altogether draining more than 185,000 km 2 of the Boreal Shield, Taiga Shield and Hudson Plains ecozones of Canada. River discharge measurements, geochemical data (δ 2 H-δ 18 O and electrical conductivity of water), remote sensing, and GIS models are used jointly to calculate water balances. The approach allows for partitioning the influence of rainwater, snowmelt, surface runoff, evaporation, transpiration, and groundwater discharge to the hydrological balances of watersheds. On an annual basis, the results suggest that runoff from rainwater (30-61% of total precipitation) and snowmelt (18-40% of total precipitation) are the main contributions to river discharge, while the contribution of groundwater discharge to rivers represents <12% of the total precipitation. Over the study area, this contribution represents 2-5 km 3 of water. The stable isotope mass balances allow for estimating watershed-scale evaporation over inflow ratios ranging between 2 and 10% and confirm that transpiration has an isotopic composition close to summer rainwater. The hydrological balances further suggest that the total pool of water stored in the active portion of watersheds represents 10 to 20% of the total annual precipitation, while the export of groundwater beyond watersheds is negligible. The seasonal trends in the hydrological balances of monitored watersheds were further documented to provide insights into the sensitivity of watersheds as they face climate change. The observations are used to propose recommendations for monitoring of rivers in the Canadian Shield and to identify future research needs.

Introduction

A quantitative evaluation of the hydrological balance of watersheds is needed to better understand the impacts of climate change and human activities on water resources.

While field measurements and remote sensing allow for estimating precipitation, evaporation, transpiration and surface water discharge in watersheds, the contribution of groundwater in the hydrological balance of watersheds remains poorly quantified.

In addition, the water storage within the active portion of watersheds (i.e., the portion of a watershed where hydrological fluxes have a significant impact on the annual hydrological balance of the watershed as a whole) is rarely measured, partly because the volume of aquifers is often not well documented. Nevertheless, it is recognized that groundwater represents the largest freshwater reservoir at the global scale [START_REF] Oki | Global Hydrological Cycles and World Water Resources[END_REF][START_REF] Gleeson | The global volume and distribution of modern groundwater[END_REF]) and that the estimated contribution of groundwater to global river discharge reach ~10% [START_REF] Oki | Global Hydrological Cycles and World Water Resources[END_REF]. Since the direct evaluation of groundwater fluxes is complex, especially in vast and remote watersheds, groundwater-surface waters interactions are most often studied based on observations on surface waters. River monitoring is particularly relevant to the evaluation of hydrological and geochemical processes taking place at the watershed scale. As a leading example, geochemical balances based on the fluxes exported by large rivers allowed for quantifying chemical weathering processes at the global scale [START_REF] Gaillardet | Global silicate weathering and CO2 consumption rates deduced from the chemistry of large rivers[END_REF][START_REF] Dupré | Rivers, chemical weathering and Earth's climate[END_REF].

In cold and humid climates, the precipitation pattern, the intensity and duration of the evapotranspiration and freezing periods control the seasonality of the river flow regimes. Several studies have highlighted the benefits of this seasonal pattern for boreal ecosystems [START_REF] Schindler | Sustaining aquatic ecosystems in boreal regions[END_REF][START_REF] Nilsson | Boreal Riparian Vegetation Under Climate Change[END_REF][START_REF] Mustonen | Defining the natural flow regimes of boreal rivers: relationship with benthic macroinvertebrate communities[END_REF]. Some of these ecosystems depend on groundwater/surface-water hydrological exchanges [START_REF] Alley | Flow and Storage in Groundwater Systems[END_REF][START_REF] Hayashi | Effects of Ground Water Exchange on the Hydrology and Ecology of Surface Water[END_REF][START_REF] Bertrand | Review: From multi-scale conceptualization to a classification system for inland groundwaterdependent ecosystems[END_REF]Kløve et al. 2014a). These exchanges vary in time and space and depend on numerous meteorological, geological, hydrogeological, geomorphological, and cryogenic processes [START_REF] Zektser | Groundwater fluxes in the global hydrologic cycle: past, present and future[END_REF]. In the east coast area of James Bay (northwestern Quebec, Canada), the water balance of some watersheds is also extensively affected by large hydroelectric complexes within which waters are diverted and vast reservoirs are flooded. Such anthropogenic modifications have significantly altered the natural pattern of streamflow [START_REF] Woo | Streamflow hydrology in the boreal region under the influences of climate and human interference[END_REF][START_REF] Woo | Effects of Reservoirs on Streamflow in the Boreal Region[END_REF].

There is no comprehensive quantitative evaluation of the groundwater/surface-water interactions in large watersheds of the east coast of James Bay. This lack of knowledge imposes major limitations for land management (Kløve et al. 2014a). The precise delineation and quantification of groundwater/surface-water interactions is critically needed to develop strategies to protect the biodiversity associated to GDEs [START_REF] Kløve | Protection of groundwater dependent ecosystems: current policies and future management options[END_REF][START_REF] Sophocleous | Interactions between groundwater and surface water: the state of the science[END_REF]. Obtaining this knowledge is especially challenging in vast and remote areas where data acquisition is complex and costly.

To overcome this challenge, remote sensing data and spatial precipitation models are increasingly used to study the hydrological balance of watersheds [START_REF] Sahoo | Reconciling the global terrestrial water budget using satellite remote sensing[END_REF][START_REF] Pan | Multisource Estimation of Long-Term Terrestrial Water Budget for Major Global River Basins[END_REF][START_REF] Yao | Estimation of the terrestrial water budget over northern China by merging multiple datasets[END_REF]). Nevertheless, groundwater storage and fluxes are hardly detectable. The change in the water storage that can be obtained from GRACE satellites does not allow researchers to distinguish between shallow groundwater and surface water storage, and the low resolution of the data makes them more suitable for continental or global scale analyses. Field monitoring is therefore required, and several studies revealed the potential of δ 2 H-δ 18 O isotope mass balances for partitioning the sources contribution of streamflow in watersheds set in cold and humid climates [START_REF] St Amour | Isotopic time-series partitioning of streamflow components in wetlanddominated catchments, lower Liard River basin, Northwest Territories, Canada[END_REF][START_REF] Kalbus | Measuring methods for groundwater – surface water interactions: a review[END_REF]Gleeson et al. 2009;Orlova and Brandifeum 2014;[START_REF] Gibson | Isotope-based partitioning of streamflow in the oil sands region, northern Alberta: Towards a monitoring strategy for assessing flow sources and water quality controls[END_REF][START_REF] Rosa | Processes governing the stable isotope composition of water in the St. Lawrence river system, Canada[END_REF][START_REF] Bansah | Evaluating the Effects of Tracer Choice and End-Member Definitions on Hydrograph Separation Results Across Nested, Seasonally Cold Watersheds[END_REF][START_REF] James | The isotope hydrology of the Muskoka River Watershed, Ontario, Canada[END_REF]. Stable isotopes of water are also extensively used to quantify evaporation processes at various scales [START_REF] Craig | Deuterium and oxygen 18 variations in the ocean and marine atmosphere[END_REF][START_REF] Gibson | Estimating Evaporation Using Stable Isotopes: Quantitative Results and Sensitivity Analysis for Two Catchments in Northern Canada: Paper presented at the[END_REF][START_REF] Gibson | Regional water balance trends and evaporation-transpiration partitioning from a stable isotope survey of lakes in northern Canada[END_REF][START_REF] Skrzypek | Estimation of evaporative loss based on the stable isotope composition of water using Hydrocalculator[END_REF][START_REF] Rosa | Processes governing the stable isotope composition of water in the St. Lawrence river system, Canada[END_REF][START_REF] Rey | Using water stable isotopes for tracing surface and groundwater flow systems in the Barlow-Ojibway Clay Belt, Quebec, Canada[END_REF].

The current research builds on the aforementioned studies and aims to improve the knowledge of water balance of watersheds of the Canadian Shield, on the east coast of James Bay (northwestern Quebec, Canada). The approach focuses on the use of largescale datasets (remote sensing and spatial models) coupled with river discharge and isotopic (δ 2 H-δ 18 O of water) monitoring. The focus is on evaluating the groundwater hydrological functions in the water balance of watersheds and quantifying groundwater/surface-water interactions.

Study area

The study area covers 185,000 km 2 of the Canadian Shield in the Anischnabe (South)

and Cree (North) First Nations territories in northwestern Quebec, Canada The region is part of the Superior Geological Province, a component of the Canadian Shield (Figure 3.1b). Surficial deposits inherited from the last glaciation-deglaciation cycle cover most of the lowlands and bedrock depressions [START_REF] Hardy | La déglaciation et les épisodes lacustre et marin sur le versant québécois des basses terres de la baie de James[END_REF]. The extent and architecture of these surficial deposits and the characteristics of associated granular aquifers are described in section 2.5.2. Watersheds of the study area are part of the Hudson, James, and Ungava Bays basin, where measurable changes in discharge are to shallow groundwater (𝑆𝐺𝑊) recharge. The plants will intercept and transpire part of this water in such a way that water is returned directly to the atmosphere [𝑞 𝑇 ].

Part of the surface runoff from rain (qR) and snow (qS) will reach large rivers. The 𝑆𝐺𝑊 can either discharge in springs, streams, lakes, rivers, and in inland discharge areas

[qSGW] or reach a deeper groundwater flow system (𝐷𝐺𝑊) that could allow for groundwater flow out of the watershed [𝑞 𝐷𝐺𝑊 ], either through coastal groundwater discharge or inter-basin hydrological exchanges. Ultimately, the surface runoff and shallow groundwater that reach the hydrographic network will partly be evaporated within rivers and connected lakes [𝑞 𝐸 ] and exported out of the watershed as river discharge (𝑄 𝑅𝐼𝑉 ). The following expression can thus be used to describe the hydrological balance of rivers:

𝑄 𝑅𝐼𝑉 + 𝑞 𝐸 = 𝑞 𝑆 + 𝑞 𝑅 + 𝑞 𝑆𝐺𝑊 = 𝐼 -𝛽 (2)
Where 𝛽 represents the proportion of water that could not reach the hydrographic network:

𝛽 = 𝑞 𝑇 + 𝜌 = 𝐼 -𝑞 𝑆 -𝑞 𝑅 -𝑞 𝑆𝐺𝑊 (3)
In Equation 3, the 𝛽 term corresponds to the amount of water lost by watersheds due to transpiration (𝑞 𝑇 ), plus a 𝜌 term representing the groundwater balance including both the change in water storage within the watershed (±∆𝑠) and the deep groundwater flow not reaching rivers (𝑞 𝐷𝐺𝑊 ). This groundwater balance is expressed by the rearrangement of Equation 1as Equation 4where ET is the evapotranspiration

(ET = 𝑞 𝐸 + 𝑞 𝑇 ): 𝜌 = 𝑞 𝐷𝐺𝑊 ± ∆𝑠 = 𝐼 -𝐸𝑇 -𝑄 𝑅𝐼𝑉 (4)
The aim of the current research is to decipher the magnitude of these processes in a boreal subcontinental scale watershed.

The geochemical mass balance of watersheds

The hydrological balance (Equation 1) can further be coupled to geochemical balances using the isotopic composition (𝛿) and electrical conductivity (𝜎) of water: undergoes a gradual heavy isotopes enrichment during the evaporation process. When plotted in a 𝛿 2 𝐻 -𝛿 18 𝑂 graph, residual liquid waters will tend to define local evaporation lines (LEL; basal plane in Figure 3.3). Evaporation also contributes to increase the electrical conductivity of surface waters, as dissolved solids tend to remain in the liquid phase during evaporation while the evaporated moisture has few dissolved solids (𝜎 𝐸 ≈ 0). The 𝛿 𝐸 term can be estimated using a simplified version of the Craig and Cordon (1965) model [START_REF] Gibson | Estimating Evaporation Using Stable Isotopes: Quantitative Results and Sensitivity Analysis for Two Catchments in Northern Canada: Paper presented at the[END_REF]):

𝐼𝛿 𝐼 =
𝛿 𝐸 = 𝛿 𝑅𝐼𝑉 -ℎ𝛿 A -𝜀 1-ℎ (7)
where 𝛿 𝑅𝐼𝑉 is the flux-weighted isotopic composition of river, ℎ the relative humidity, 𝛿 𝐴 is the moisture composition of ambient air (assumed to be in equilibrium with the flux-weighted average isotopic composition of local precipitation), and 𝜀 is the total isotope fractionation [START_REF] Gibson | Estimating Evaporation Using Stable Isotopes: Quantitative Results and Sensitivity Analysis for Two Catchments in Northern Canada: Paper presented at the[END_REF][START_REF] Gibson | Regional water balance trends and evaporation-transpiration partitioning from a stable isotope survey of lakes in northern Canada[END_REF][START_REF] Skrzypek | Estimation of evaporative loss based on the stable isotope composition of water using Hydrocalculator[END_REF]. The 𝛿 2 𝐻 -𝛿 18 𝑂 of atmospheric moisture (𝛿 𝐴 ) during the iceoff period can be estimated based on the assumption that it is in equilibrium with the stable isotope composition of precipitation and using an average air temperature (≈15°C), representing the prevailing conditions during the ice-off season. The liquidvapour isotopic fractionation factor (𝛼 𝑉 𝐿 ) can be estimated based on the [START_REF] Horita | Liquid-vapor fractionation of oxygen and hydrogen isotopes of water from the freezing to the critical temperature[END_REF] equations.

It is assumed that the isotopic composition of shallow and deep groundwater [𝛿 𝑆𝐺𝑊 , 𝛿 𝐷𝐺𝑊 ] is similar to the amount weighted average isotopic composition of precipitation, consistent with previous observations conducted in the Abitibi-Témiscamingue region (south of the study area) [START_REF] Rey | Using water stable isotopes for tracing surface and groundwater flow systems in the Barlow-Ojibway Clay Belt, Quebec, Canada[END_REF]). The groundwater sampled in the latter study systematically plots on the LMWL, which suggests that there is no significant evaporation prior to infiltration or that the evaporation occurring directly on soils and vegetation is complete (i.e., 100% of the pool of water on soils and vegetation available for evaporation is returned to the atmosphere).

As groundwater flows, its dissolved solids concentration and electrical conductivity (𝜎 𝑆𝐺𝑊 ) increase as a function transport time and distance due to increased water-rock interactions. This concept is well illustrated by the conceptual cross-sections proposed by Tòth (1999). Pre-existing data on groundwater in the southern part of the study area revealed that groundwater from unconfined granular aquifers, unconfined bedrock aquifers, and confined aquifers present distinct electrical conductivities (𝜎 𝑈𝐺𝐴 ; 𝜎 𝑈𝐵𝐴 ; 𝜎 𝐶𝐴 , respectively, in Figure 3.2a). Overall, in this area, unconfined aquifers tend to present lower 𝜎 than confined aquifers, and granular aquifers tend to present lower 𝜎 than bedrock aquifers [START_REF] Rey | Using water stable isotopes for tracing surface and groundwater flow systems in the Barlow-Ojibway Clay Belt, Quebec, Canada[END_REF][START_REF] Cloutier | Projet d'acquisition de connaissances sur les eaux souterraines de l'Abitibi-Témiscamingue (partie 1) [Abitibi-Témiscamingue groundwater knowledge acquisition project[END_REF]2015;2016).

Here, from a soil water perspective, transpiration is considered a non-fractionating process as water consumption by plants does not significantly influence the 𝛿 2 𝐻 -𝛿 18 𝑂 and electrical conductivity of the remaining soil water. The isotopic composition (𝛿 𝑅𝐼𝑉 ) and electrical conductivity (𝜎 𝑅𝐼𝑉 ) of river water will represent an integrated composition inherited from hydrogeochemical processes that occurred at the scale of the watershed. The rain (𝑅) and snow (𝑆) endmembers present lower 𝜎 than the 𝑆𝐺𝑊 endmember.

Sources of river discharge

The seasonal changes in 𝛿 𝑅𝐼𝑉 follow a cycle with an enrichment in heavy isotopes in response to a greater contribution from summer precipitation and evaporation during the ice-off period, an increased 𝜎 during winter baseflow due to a greater groundwater contribution, and a marked decrease in 𝜎 and 𝛿 values during the spring freshet due to dilution from melting snow. The loop extends beyond the triangle defined by these three endmembers in the 𝜎 vs 𝛿 plots during the summer period because of evaporation, which causes heavy isotope enrichment and a slight increase in 𝜎 in surface waters.

The star symbol in the Figure 3.3 represents the flux-weighted average value for 𝛿 𝑅𝐼𝑉 .

As presented in the conceptual model described above (Equation 2), the 𝑄 𝑅𝐼𝑉 term represents a mixture of 𝑞 𝑅 , 𝑞 𝑆 , and 𝑞 𝑆𝐺𝑊 that is affected by losses caused by evaporation (𝑞 𝐸 ). The following expression can thus be used to describe the conceptual model of hydrological fluxes in the watersheds and to partition streamflow components:

𝑄 𝑅𝐼𝑉 = 𝑞 𝑆 + 𝑞 𝑅 + 𝑞 𝑆𝐺𝑊 -𝑞 𝐸 = 1 (8)
The 𝛿 and 𝜎 can further be used:

𝛿 𝑅𝐼𝑉 = q 𝑆 𝛿 s + q R 𝛿 R + q SGW 𝛿 SGW -𝑞 𝐸 𝛿 E (with 𝛿 2 𝐻 and 𝛿 18 𝑂) (9)
σ 𝑅𝐼𝑉 = q 𝑆 σ s + q R σ R + q SGW 𝜎 SGW -𝑞 𝐸 σ E (10)

Equations 8 to 10 can be solved for 𝑞 𝑅 , 𝑞 𝑆 , 𝑞 𝑆𝐺𝑊 , and 𝑞 𝐸 using matrix calculations. The matrix product (𝐴𝑋 = 𝐵; Equation 11) is solved for the unknowns using inverse matrix (𝐴 -1 𝐵 = 𝑋; Equation 12) where 𝐴 is the matrix of determinants (𝑎), 𝑋 is the matrix of unknowns (𝑥), and 𝐵 is the matrix of results (𝑏):

( 𝑎 1,1 𝑎 1,2 … 𝑎 1,𝑛 𝑎 2,1 𝑎 2,2 … 𝑎 2,𝑛 ⋮ ⋮ ⋱ ⋮ 𝑎 𝑛,1 𝑎 𝑛,2 ⋯ 𝑎 𝑛,𝑛 ) × ( 𝑥 1 𝑥 2 ⋮ 𝑥 𝑛 ) = ( 𝑏 1 𝑏 2 ⋮ 𝑏 𝑛 ) (11) ( 1 1
𝛿 𝑇 = 𝛽𝛿 𝛽 -𝜌𝛿 𝐼 𝛽-𝜌 (15)
The 𝛿 𝑇 values can provide insights on the main sources of water used by plants.

Riverine process during the ice-on period

After the high waters of the fall period, river discharge is expected to present a gradual decrease during the ice-on period, until reaching winter baseflow conditions (Figure 3.2b). Under ice-on conditions, the contribution from runoff decreases and large rivers mainly drain large lakes and groundwater. During this period, the water in rivers freezes and the isotopic evolution of river water over time can potentially be deduced from ice cores if the effective ice-water isotopic separation (𝜀 𝑒𝑓𝑓 ), representing the difference between the isotopic composition of the basal ice (𝛿 𝑏𝑎𝑠𝑎𝑙 𝑖𝑐𝑒 ) and that of the sub-ice water (𝛿 𝑠𝑢𝑏-𝑖𝑐𝑒 𝑤𝑎𝑡𝑒𝑟 ) is known:

𝜀 𝑒𝑓𝑓 = 𝛿 𝑏𝑎𝑠𝑎𝑙 𝑖𝑐𝑒 -𝛿 𝑠𝑢𝑏-𝑖𝑐𝑒 𝑤𝑎𝑡𝑒𝑟 (16)
If 𝜀 𝑒𝑓𝑓 can be estimated, ice cores collected at the end of the ice-on period could provide key insights on the isotopic composition of river water during the winter period, a time at which river monitoring and water samples collection is difficult.

Methods

The approach used here focuses on atmospheric fluxes (𝐼 and 𝐸𝑇), river discharge measurements (𝑄 𝑅𝐼𝑉 ), and the isotopic composition (𝛿 2 𝐻 -𝛿 18 𝑂) and electrical conductivity (𝜎) of water. The data are used to solve the equations included in the theoretical and conceptual models, which allow to estimate atmosphere-surface watergroundwater hydrogeochemical exchanges.

Q 𝐵𝑓 = 4.5 × 10 -9 (m 𝑠 ⁄ ) × A (17)
where Q 𝐵𝑓 is the river baseflow (m 3 /s) and 𝐴 is the catchment area (m 2 ). This empirical relation is based on a compilation of historical baseflow measurements from 10 hydrometric stations located in study area (MELCC 2020).

Geochemical

Results

The complete dataset from this study is provided in Appendices C to G. A summary of hydrological (𝑅, 𝑆, 𝐸𝑇, 𝑄) and geochemical (𝛿 2 𝐻 -𝛿 18 𝑂 and σ) results is presented below.

Hydrological results

The average annual atmospheric inflows (I) for entire region ranged between 880 ± 11% mm/yr in 2018 and 849 ± 10% mm/yr in 2019 (Table 3 The results obtained from discharge measurements are shown in Figure 3. period (November-April). The results suggest that peak discharge over winter baseflow ratios rang between 3.7 and 16.5.

Geochemical results

The complete geochemical dataset (𝛿 2 𝐻 -𝛿 18 𝑂; 𝜎; T; DO; pH; ORP) from this study is provided in Appendices D to G, while Table 3.3 and Figure 3.6 provide a summary of available results. The 𝛿 2 𝐻 and 𝛿 18 𝑂 results obtained for river waters ranged from -128.3‰ to -75.7‰

and -17.7‰ to -9.9‰, respectively (Figure 3.7b-c). The lowest and highest 𝛿 2 𝐻 -𝛿 18 𝑂 values are recorded during spring snowmelt and during summer baseflow, respectively. The 𝛿 2 𝐻 and 𝛿 18 𝑂 of winter river baseflow range from -107.8‰ to -85.2‰ and -14.5‰ to -11.62‰, respectively. The 𝜎 measured in rivers suggest values ranging from 10.4 μs/cm to 159.8 μs/cm, with an average of 47.5 μs/cm (Table 3.3).

The 𝜎 values increased slightly during winter and reach a maximum at the end of the ice-on period, followed by a marked decrease during spring flood.

The 𝛿 2 𝐻 and 𝛿 18 𝑂 results obtained for 192 river ice samples allowed for constructing 20 vertical profiles of river ice. Overall, the length of the ice cores ranged between approximately 0.3 m and 1.0 m and the isotopic composition 𝛿 2 𝐻 -𝛿 18 𝑂 of the ice samples ranged from -66.4‰ to -213.1‰ and -8.6‰ to -29.0‰, respectively. The complete dataset are reported in Appendix G.

Discussion

The following discussion focuses on evaluating the hydrological flux presented in the conceptual and theoretical models (section 3.3). The extent of the calculations depends on the availability of data for each sub-watershed (see Figure 3.1). Detailed calculations are proposed for rivers that have been subject to intensive monitoring, while qualitative interpretations are proposed for the sub-basins where the data are fragmentary. The advantages and limits of the methods used and developed herein are then discussed, along with recommendations pertaining to future research needs.

Atmospheric fluxes and latitudinal gradients

The main data used for documenting the atmospheric inflows (𝐼, 𝑅, and 𝑆) to watersheds are the IMERG and CaPA datasets (for hydrological fluxes from rain and snow) and the RCWIP (for the isotopic composition of precipitation). The 𝐼 calculated from the IMERG dataset are 20% lower than that calculated from the CaPA model.

The models available for estimating the 𝐸𝑇 term also show discrepancies, as the MODIS dataset yields estimates of 𝐸𝑇 that are 29% lower than those of the SSEBop dataset.

The IMERG and CaPA models both suggest dryer conditions for the northernmost watersheds in comparison to the southernmost watersheds (Table 3.2). This pattern is consistent with previously reported climate data that suggests colder and dryer conditions in more northern areas within the study region [START_REF] Phillips | Les climats du Canada[END_REF][START_REF] Berteaux | Portrait du climat du nord du Québec et du Labrador pour la période 1981-2010. Rapport présenté au Consortium Ouranos sur la climatologie régionale et les changements climatiques[END_REF]. The MODIS and SSEBop data further suggest that the ET fluxes decrease northwards, likely due to a shorter ice-off period characterized by lower temperatures. The 𝐼 and 𝐸𝑇 gradients (see section 3.5.1) do not completely balance each other and, as a result, the watershed-area-normalized river discharges (expressed in mm/yr) show a gradient of -25 mm/° φ. This finding implies that the area-normalized hydrological fluxes from the study area decrease northwards. Climate change could entail changes in these gradients. It is expected that hydrological fluxes, including precipitation, evapotranspiration, and river discharge will increase by 2050 in the study area [START_REF] Hennigs | État des connaissances relatives aux changements climatiques et à l'adaptation dans le territoire Eeyou Istchee Baie-James. Rapport présenté au Comité consultatif pour l'environnement de la Baie James[END_REF]. Theses changes in hydrological fluxes could generate proportionally more significant impacts in the northeast portion of the study area,

where current climatic conditions are dryer and colder.

Within the study region, the RCWIP data suggest -5.2‰/°φ and -0.63‰/°φ gradients in the 𝛿 2 𝐻 and 𝛿 18 𝑂 of precipitation, respectively (based on annual flux-weighted values of watershed). A similar gradient is also observed in the isotopic composition of rivers, where it is observed with -5.4‰/°φ and -0.67‰/°φ for 𝛿 2 𝐻 and 𝛿 18 𝑂, respectively. The river water samples revealed that 𝛿 𝑅𝐼𝑉 tend to plot below the LMWL in 𝛿 2 𝐻 -𝛿 18 𝑂 cross plots. The flux-weighted average river discharge is enriched by approximately 1 to 9‰ and 0.7 to 2.0‰ for 𝛿 2 𝐻 and 𝛿 18 𝑂, respectively, compared to the modelled flux-weighted average isotopic composition of precipitation. These observations likely reflect the heavy isotope enrichment inherited from evaporation processes occurring in watersheds. The annual flux-weighted river composition at site CC presents the greatest heavy isotope enrichment (15‰ and 3.0‰ for 𝛿 2 𝐻 and 𝛿 18 𝑂, respectively), but the peak discharge associated with snowmelt was likely missed at this site.

The hydrological balance of watersheds

The hydrological balance of watersheds is assessed by evaluating the terms included in Equation 4. The available results and data allow for quantifying the 𝐼, 𝐸𝑇, and 𝑄 𝑅𝐼𝑉 terms in 2018 for 9 sub-watersheds (Table 3.2). The term 𝜌 is calculated annually (i.e., from one winter baseflow to the other) and represents the residual between 𝐼, 𝐸𝑇, and 𝑄 𝑅𝐼𝑉 (Equation 4). The estimated 𝜌 values range between -208 mm/yr and 131 mm/yr. These values correspond to proportions of -20% and 15% of average 𝐼. One key feature is that the calculated 𝜌 terms are smaller in sub-watersheds where hydrometric stations provide daily measurements river discharge (𝑄 𝑅𝐼𝑉 in Table 3.2).

In these cases, the 𝜌 values are less than 7% of mean 𝐼.

The 𝜌 term is more variable in watersheds where there are no hydrometric stations. In these watersheds, the river discharges are interpolated between a limited number of field measurements, thus yielding a higher uncertainty than for rivers where a hydrometric station allows for calculating daily values of discharge. Overall, the most precise and reliable discharge data suggest that 𝐼 -𝐸𝑇 ≈ 𝑄 𝑅𝐼𝑉 . While the annual discharge is not precisely known for all watersheds, the calculations suggest that the water transit within aquifers could act as a buffer to hydrological variations in rivers, without modifying the total annual discharge from rivers. It is expected that the hydrological response of rivers will be more dynamic in watersheds where the direct runoff/recharge ratio is greater and dampened in watersheds where the direct runoff/recharge ratio is smaller.

𝑄
While the ∆𝑆 term is likely to be close to zero on an annual basis, significant changes are observed on monthly timescales. Here, the data suggest that the water storage in watersheds reach maximum values from autumn to the end of the ice-on period, due to humid conditions at the end of the ice-off period followed by low river discharge and snow accumulation during the ice-on period. The water storage in watersheds reach minimum values during the early summer period due to high discharge during freshet and increased 𝐸𝑇 during early ice-off period.

At the watershed scale, the data suggest a storage of water representing an average of 10 to 20% of the annual 𝐼. This transient storage presents a cycle with an amplitude of 10 to 18%. Based on this observation, it can be postulated that the reduction of the duration of ice-on period due to global warming will reduce the maximum storage in watersheds in response to a decrease in snow accumulation and an earlier peak in 𝐸𝑇 at the onset of the ice-off period. 3.4).

These scenarios allowed for quantifying the sensitivity of calculations. suggest that the runoff from rain reaching rivers represents 30 to 61% of 𝐼, and that the runoff of the melting of snow reaching rivers represents 18 to 46% of 𝐼. The shallow groundwater flux reaching rivers is estimated to range between 0 and 12% of 𝐼. It is further estimated that 2 to 10% of 𝐼 is returned to the atmosphere by evaporation. This quantitative estimation of evaporation is similar to the evaporation over inflow ratios (𝐸/𝐼) of watersheds estimated based on the [START_REF] Craig | Deuterium and oxygen 18 variations in the ocean and marine atmosphere[END_REF] model (Table 3.6), following the calculation procedures presented in [START_REF] Gibson | Estimating Evaporation Using Stable Isotopes: Quantitative Results and Sensitivity Analysis for Two Catchments in Northern Canada: Paper presented at the[END_REF], [START_REF] Gibson | Regional water balance trends and evaporation-transpiration partitioning from a stable isotope survey of lakes in northern Canada[END_REF], [START_REF] Telmer | Isotopic constraints on the transpiration, evaporation, energy, and gross primary production Budgets of a large boreal watershed: Ottawa River Basin, Canada[END_REF], and [START_REF] Skrzypek | Estimation of evaporative loss based on the stable isotope composition of water using Hydrocalculator[END_REF]:

𝐸 = 𝑖 (𝛿 𝑅𝐼𝑉 -𝛿 𝑖 ) 𝑚(𝛿 * -𝛿 𝑅𝐼𝑉 ) ( 18
)
where 𝐸 is the total evaporation within the watershed [L 3 /T], 𝑖 (here, 𝑖 = 𝐼) is the total watershed inflow [L 3 /T] (also used in Equations 1-5), 𝛿 𝑅𝐼𝑉 is the flux-weighted average isotopic composition of river discharge (‰), 𝛿 𝐼 is the flux weighted average isotopic composition of watershed inflow (‰), 𝛿 * is the limiting isotopic composition (‰) [START_REF] Gat | Isotope hydrology of inland sabkhas in the Bardawil area[END_REF][START_REF] Gat | Properties of the isotopic species of water: the ''isotope effect[END_REF], and 𝑚 is the enrichment slope (see [START_REF] Skrzypek | Estimation of evaporative loss based on the stable isotope composition of water using Hydrocalculator[END_REF] for details on these parameters). *From 2018 average precipitation data with 𝜌 ≠ 0, **From 2018 average precipitation data with 𝜌 = 0, ***From 2019 average precipitation data with 𝜌 = 0

The calculations from the scenarios suggest that the 𝜎 𝑆𝐺𝑊 values used for partitioning streamflow components strongly influence the results. Site HA is used to discuss the sensitivity of calculations to the 𝜎 𝑆𝐺𝑊 values. At this site, the calculated annual contribution of groundwater to river discharge ranges between 8 and 25% of (𝐼) for the lowest and highest 𝜎 𝑆𝐺𝑊 endmember values, respectively. For large-scale watersheds, a composite value integrating multiple aquifer water sources is more likely to represent regional flow pattern. In that sens, for the monitored watersheds, the river discharge water from groundwater is constrained between 1.7 to 4.9 km 3 /yr. This regional contribution reprensents 3% of precipitation. This groundwater discharge compares well to other calculations realized in watersheds of the southerm Canadian Shield (Gleeson et al. 2009).

The seasonality of river water sources

The data and calculations allow for evaluating the seasonal cycle in the sources of river the εeff of 2.84 ± 0.34‰ and 17.9 ± 3.9 ‰ found by [START_REF] Gibson | Stable isotopes in river ice: identifying primary over-winter streamflow signals and their hydrological significance[END_REF] in northwestern Canada.

The isotopic patterns measured in zones with congelation ice are considered to represent the conditions illustrated on the left side of the conceptual model (Figure 3.9).

The stable conditions observed during the ice-on period are attributed to the constant supply of water from large lakes and groundwaters. In these zones, river water freezing must exceed ice erosion. Seven of the 20 ice cores samples did not show clean congelation ice. In these samples, the comparison between the isotopic composition of sub-ice water and that of ice reveals negative 𝜀 𝑒𝑓𝑓 values. These values are most likely caused by the presence of a basal ice that is depleted in heavy isotopes due to snow-ice mixing, as illustrated in the central part of the conceptual model (Figure 3.9). [START_REF] Gibson | Stable isotopes in river ice: identifying primary over-winter streamflow signals and their hydrological significance[END_REF] suggested that the proximity of rapids and the presence of layers of slush or liquid water are non-ideal conditions for congelation ice formation.

Another plausible explanation for the absence of congelation ice at these sites is that ice erosion exceeds ice formation. Such conditions were observed in areas where the river section is reduced (i.e., shallower or narrower), causing higher flowrates under the ice.

Here, it was also observed that some ice cores presented a heterogeneous structure with variable 𝛿 2 𝐻 and 𝛿 18 𝑂 compositions. These cores are likely associated with the conditions illustrated on the right side of the conceptual model (Figure 3.9). In these zones, the results and observations suggest that the ice cover consists of an accumulation of ice debris in zones where the morphology of the riverbed entails conditions that are prone to ice jam. Overall, the observations made here suggest that ice cores from rivers should be collected in areas where rivers are narrow and deep if the research objective is to document the isotopic composition of river water during the ice-on period. These conditions would likely be conductive to the development of In the Harricana River in Amos (site HA), the presence of rapids significantly increases the streamflow velocity, leaving the water surface free of ice all year long. It was therefore possible to monitor the isotopic composition of water at this site during the ice-on period. The results suggested a progressive depletion in heavy isotopes in river water during the ice-on period. Nevertheless, the results revealed that river water during the ice-on baseflow is enriched in heavy isotopes compared to groundwater, suggesting a non-negligible contribution from the rain endmember throughout the winter.

Moreover, the increasing value of the 𝐷 𝑥𝑠 observed during winter at this site suggests that with ice cover, evaporated water at the surface of lakes is mixed with underlying water or freeze in small tributary river and further released into the spring flood.

The isotopic composition of transpired moisture

Evaluating the isotopic composition of the moisture transpired by plants based on field measurements is a challenging task, especially in vast and remote areas. In this research, the isotopic composition of the evaporated moisture is estimated based on Equation 15 (section 3.4). In the monitored watersheds where sufficient data are available for deciphering the sources of discharge, the calculated 𝛿 𝑇 values range from -107.7‰ to -55.1‰ for 𝛿 2 𝐻 and from -14.81‰ to -8.14‰ for 𝛿 18 𝑂 (Table 3.7). The data presented in Table 3.6 reveal that 𝛿 𝑇 are generally enriched in heavy isotopes compared to the flux-weighted average isotopic composition of precipitation (𝛿 𝐼 ). This finding suggests that the water transpired by plants mainly corresponds to summer rains. While this might seem obvious, it raises an important question on groundwater recharge processes: if transpiration preferentially removes water from summer rains, why is the groundwater isotopic composition very similar to the flux-weighted average isotopic composition of precipitation (𝛿 𝐼 )? One plausible explanation is that the infiltration process could be less efficient during the snowmelt period due to frozen soils, therefore limiting the recharge of waters that are significantly depleted in heavy isotopes and balancing the effect of transpiration (which limits the recharge of waters that are significantly enriched in heavy isotopes). This explanation would imply that within the study area, groundwater recharge mainly occurs shortly after the snowmelt period, when transpiration by plants is still limited, and towards the end of the ice-off season, when colder conditions and less daylight reduce vegetation growth. The seasonality in river water sources revealed by the data highlights the need for river water monitoring. The results suggest that time series should focus on three main periods:

-The end of the ice-on period: Samples collected during this period are critical to evaluate the maximum groundwater contribution. Theses samples can further provide information on the σ_SGW term, which is needed for hydrograph separation. Since temporal changes are unlikely to occur during winter baseflow, only a limited number of samples is required for this period.

-The spring freshet: The sampling frequency must be increased during this period to calculate the sources of river waters during the peak discharge event.

However, water sampling during the spring flood is difficult due to safety issues. During this period, unstable ice conditions are often observed along rivers shorelines, which limits access for sampling. In addition, large blocks of ice are transported in the rivers during snowmelt, which limits the possibility of using small boats for sampling. Further research aimed at developing safe and efficient protocols for the collection of river samples during snowmelt are, therefore, much needed.

-The end of the ice-off period: Samples collected at the end of the ice-off season provide key information on evaporation processes, as the results suggest that the greatest heavy isotope enrichments are observed during this period.

Sampling this period is also important to capture the contribution of rainwater runoff during the high discharge event observed before the onset of the ice-on period.

These recommendations are mainly suited to river monitoring in vast and remote areas characterized by cold and humid climates. Monitoring programs implemented on rivers in such areas can be costly and logistically complex. Therefore, focusing on key periods of rivers' hydrographs could help improve knowledge of hydrological processes occurring in watersheds.

Limits, challenges, and advantages of the proposed approaches

The monitoring approach used in this study allows for documenting seasonal changes in hydrological processes occurring within watersheds. Nevertheless, sampling during the snowmelt period is especially complex. In vast areas such as the study region, the snowmelt can be diachronous due to differences in latitude and altitude, among other aspects. In addition, the spring freshet can be very short if weather conditions such as warm temperatures and heavy rains contribute to snowmelt. Such events, however, can only be predicted on short notice. Therefore, it is challenging for a research team to conduct intense monitoring during snowmelt. One of the most promising approaches to this problem would be to develop research partnerships with local communities.

Such partnerships would encourage knowledge sharing. Focusing on sites located near discharge monitoring stations would help reduce the required field procedures, which would also facilitate collaborations with local communities.

The results from this study suggest that the isotopic compositions of precipitation from RCWIP provide reliable values that could be used for future large-scale studies. Being able to use RCWIP is an advantage because monthly monitoring of precipitation can be challenging in remote areas. One shortcoming pertaining to precipitation is that, in the present study, the potential effect of sublimation on the snowpack was not considered, while previous studies suggested that the loss of snow by sublimation can reach 40% in boreal regions [START_REF] Pomeroy | Fuzzy logic-based assessment for mapping potential infiltration areas in lowgradient watersheds[END_REF]). In addition, the effects of isotopic fractionation during snowmelt were not considered here [START_REF] Taylor | Isotopic evolution of a seasonal snowpack and its melt[END_REF]. This lack could lead to uncertainties in the evaluation of rivers water budgets. Future research focusing on river monitoring would benefit from evaluating these procedures through a monitoring of the snowpack.

Another shortcoming of the present study is that the electrical conductivity of groundwater remains poorly documented. This lack of knowledge increases the uncertainty associated with the quantification of groundwater discharge to rivers.

While a comprehensive assessment of groundwater quality in the northern part of the Canadian Shield is not yet feasible, mainly due to the lack of groundwater wells, one plausible solution would be to monitor groundwater electrical conductivity in groundwater springs. Such springs provide a direct access to groundwater where automated probes could be installed for a low-cost monitoring of groundwater electrical conductivity. Finally, a better knowledge of the residence time of runoff water over large watersheds would be needed to improve water budgets.

Conclusions

This study provides new knowledge on the hydrological balance of large watersheds of the Canadian Shield in the James Bay area, in northwestern Quebec, Canada. Such knowledge is critically needed in vast and remote areas such as the northern part of the Canadian Shield, where significant impacts of climate change are expected and where human impacts are continuously increasing. This will allow for better predicting the potential impacts of changes in the hydrologic balance on local communities and GDEs. In addition, there is currently a critical lack of knowledge about groundwater resources in northern Canada, which is partially explained by the lack of groundwater wells and infrastructures that can be used for studying aquifers.

The approach used in this study made it possible to combine remote sensing, field monitoring, and spatial model data. The study provides novel insights on the role of groundwater in the hydrological budgets of watersheds. The results reveal that interbasin exports of groundwater and groundwater discharge to the ocean domain are insignificant in the study area. The data further suggest that the groundwater fluxes account for less than 12% of the rivers' discharge on an annual basis. The calculations also suggest that the total amount of water stored in the active layer of watersheds is equivalent to 10 to 20% of the total precipitation reaching the watersheds. Moreover, the results allowed for documenting atmospheric fluxes from watersheds. Evaporation estimates from remote sensing data and isotopic mass balances are in agreement; both suggest evaporation over inflow ratios in the range of 2 to 10%. While the isotopic composition of transpired moisture was not measured, it was estimated based on isotope mass balance calculations. The results suggest that the transpired moisture is plotting close to the LMWL and that it is enriched in heavy isotopes compared to the flux-weighted average isotopic composition of precipitation. These findings confirms that summer rains are preferentially transpired by plants and that the pool of soil water from the infiltration that occurred during snowmelt is not the main source of water for plants. The results further provided insights on the riverine processes occurring during the ice-on period. Ice cores collected at the end of the ice-on period and containing a sufficient thickness of congelation ice revealed that the isotopic composition of river water is invariant during the ice formation period. This suggests that the frequency of river monitoring during winter could be reduced. Overall, the data lead to the conclusion that river monitoring is critically important during three main periods: (1) at the end of the ice-on period to capture the groundwater contribution during winter baseflow, (2) during the spring flood to capture the peak discharge, and (3) at the end of the ice-off period to estimate evaporation over inflow ratios and capture the high waters of the fall.

Overall, the approaches developed and applied here could be transferred to other areas of the Canadian Shield to improve land management by allowing a better consideration of the water resources and associated ecosystems. 2) En période estivale, l'inertie thermique de l'eau permet de maintenir une température du sol plus basse dans les secteurs alimentés par l'eau souterraine;

3) L'humidité du sol est plus stable dans les secteurs alimentés par l'eau souterraine en comparaison aux secteurs où la nappe est plus profonde.

Sur la base de ces hypothèses, le POEDES est exprimé selon l'équation 5 :

𝑃𝑂𝐸𝐷𝐸𝑆 = 1 - (𝑇𝑉𝐷𝐼 ̅̅̅̅̅̅̅ + 𝐷𝜎 10
) 2

(5)

Où, POEDES est adimensionnel et prend des valeurs sur une échelle allant de 0 à 1, où 4 9,9 32,2 8,9 11,6 3,9 15_25 16 -12,4 5,4 29,7 6,3 12,1 3,8 15_26 17 -14,4 4,2 33,4 6,3 12,1 3,7 16_24 16 -17,9 6,2 32,4 6,3 11,7 3,6 16_25 14 -18,4 4,9 34,3 4,5 12,6 3,0 16_26 18 -16,6 11,4 31,5 10,1 11,2 4,0 17_24 15 -18,2 6,9 31,8 9,8 11,9 4,5 17_25 17 -15,6 5,7 31,0 6,6 11,7 3,2 17_26 22 -14,4 5,8 31,2 5,9 12,4 3,4 18_24 13 -15,7 5,4 30,6 5,9 12,7 2,4 18_25 19 -12,2 7,2 27,9 7,1 12,4 4,0 18_26 18 -15,5 4,4 27,6 10,3 9,7 6,6 19_24 15 -16,9 6,3 31,0 6,0 12,3 3,1 19_25 19 -15,9 5,8 32,2 6,7 13,6 3,0 19_26 21 -14,8 3,2 31,6 6,8 12,9 4,3 20_24 23 -14,1 6,1 28,5 7,7 13,7 4,4 20_25 20 -14,2 4,8 30,2 6,5 14,2 moyenne sur 36% des bassins versants avec un écart-type de 16%. Les lacs couvrent en moyenne 13% du territoire avec un écart-type 6%. Sur la base des échantillons d'eau prélevés dans les rivières (section 3.5.2), la contribution moyenne d'eau souterraine aux cours d'eau correspond en moyenne à 37 mm/an/L 2 avec un écart-type de 38 mm/an/L 2 . Les valeurs varient entre 3 mm/an/L 2 et 128 mm/an/L 2 . Les taux d'évaporation sur les bassins versants sont en moyenne de 45 mm/an/L 2 avec un écarttype de 29 mm/an/L 2 , alors que la transpiration atteint en moyenne 262 mm/an/L 2 avec un écart-type de 54 mm/an/L 2 . Les taux de ruissellement pour les précipitations sont en moyenne de 57% avec un écart-type de 14% pour la pluie, et de 93% avec un écarttype de 23% pour la fonte des neiges. Les valeurs négatives de certains flux et les taux de ruissellement supérieur à 100% sont probablement inhérentes à la logistique des travaux de terrain et aux conditions climatiques qui impliquent un échantillonnage asynchrone dans les différents bassins versants. L'utilisation des données sous la forme centrée réduite permet de pallier cette contrainte et s'appuie sur l'hypothèse voulant que la variance entre les valeurs des flux attribuées aux sous-bassins versants est représentative de la variance des différents contextes d'écoulement en présence. La densité moyenne du réseau de drainage est de 0,5 m -1 dans un rayon de 1000 m avec un écart-type de 0,1 m -1 . La distance moyenne à un cours d'eau est de 833 m avec un écart-type de 235 m. Les proportions de milieux humides mis en place au-dessus d'un aquifère à nappe captive (sur une couche argileuse) et sur un aquifère à nappe libre sont en moyenne de 17% et 5% avec des écarts-types respectifs de 8% et 4%. L'indice topographique d'humidité TWI moyen des bassins versants atteint 9,7 avec un écarttype de 0,4. Les indices spatiaux moyen de NDVI et TVDI sont de 0,32 avec des écartstypes de 0,01 et 0,34 avec des écarts-types de 0,06. L'indice moyen des aquifères potentiels est de 2,9 avec un écart-type de 0,2. L'épaisseur de sédiments saturés est 11,0 m avec un écart-type de 2,4 m. Le niveau de priorisation tel que défini à la section 2.5.5 est en moyenne de 2,5 avec des écarts-types de 0,2. Une valeur basse est révélatrice d'un niveau de priorisation élevé. Finalement, sur les bassins versants caractérisés, le POEDES moyen est de 0,58 avec un écart-type de 0,04. Variable CP1 CP2 CP3 CP4 GW -0,34 0,67 -0,04 -0,54 E 0,70 0,44 0,38 0,09 DFHN -0,82 0,33 0,21 0,33 Pente.moy 0,94 -0,05 -0,12 0,23 Confine -0,982 -0,04 0,06 0,11 Gra 0,69 0,56 0,32 0,13 Roc 0,91 -0,32 -0,15 0,03 TWI.moy -0,93 0,10 0,25 -0,04 POEDES -0,01 0,91 0,07 0,25 Zone_priorit. 0,06 -0,824 0,50 0,01 ESS -0,12 0,82 -0,31 0,20 AP 0,57 0,60 0,25 -0,37 Des contraintes logistiques ont cependant réduit la portée des résultats dans la mesure où un suivi annuel n'a pas pu être réalisé sur l'ensemble des sites d'échantillonnage.

Ces contraintes ont néanmoins mis en évidence la grande valeur des échantillons prélevés et la complexité de réaliser des suivis hydrogéochimiques sur d'aussi vastes territoires. Au niveau des sites les moins accessibles, un nombre limité d'échantillons a pu être prélevé, ce qui empêche de calculer des flux annuels représentatifs.

Cependant, l'ensemble des données obtenues sont rendues disponibles et pourront servir de référence pour de futurs travaux dans le secteur.

L'étude d'un territoire nordique du Bouclier canadien implique la prise en compte du couvert de glace et des conditions de gel pendant une période de plus de cinq mois.

Dans ces conditions, le suivi temporel de la composition géochimique des cours d'eau s'avère complexe. Pour obtenir des données représentatives des processus hydrologiques ayant lieu lors de cette période, des travaux de carottage de glace ont été réalisés juste avant la fonte. Les échantillons de glace ainsi prélevés peuvent être utilisés afin de suivre l'évolution de la composition isotopique de l'eau des rivières tout au long de l'hiver, en s'appuyant sur le facteur de séparation isotopique effectif entre la glace et l'eau. L'analyse des profils verticaux des carottes de glace a révélé que les tronçons de rivières étroite et profonde seraient plus favorables à la formation d'une glace de congélation idéale pour documenter la composition isotopique de l'eau s'écoulant sous le couvert de glace. L'identification de ces caractéristiques morphologiques fluviales permettra de mieux cibler les sites pour de futurs travaux analogues.

Les indices spatiaux permettant de cibler de potentiels écosystèmes dépendant de l'eau souterraine sur de vastes territoires sont de plus en plus employés au sein de la Les bilans isotopiques appliqués à l'échelle des bassins versants tendent à confirmer que les plantes utilisent davantage l'eau issue des précipitations plutôt que les eaux de fonte des neiges. La baisse anticipée des apports d'eau issue de la fonte des neiges dans le contexte de réchauffement climatique modifierait ainsi peu la provenance de l'eau du bilan hydrique des plantes. Selon les bilans hydrologiques et isotopiques réalisés, les conditions hydro-climatiques des bassins versants étudiés soutiennent une réserve d'eau active de l'ordre de 10 à 20 % des précipitations annuelles. Ces réserves stockées au sein du couvert de neige, dans les lacs ou dans les aquifères présentent un emmagasinement total qui fluctue selon une amplitude de 10 à 18 % au gré des saisons.

Les carottes de glace prélevées dans les conditions favorables ont montré une évolution progressive de la composition isotopique de l'eau de la rivière, avec un appauvrissement progressif en isotopes lourds en période de gel, suggérant une contribution accrue des pluies estivales au début de la période de gel et une augmentation de la contribution relative de l'eau souterraine lors de l'hiver. Les résultats suggèrent aussi que la formation du couvert de glace n'influence pas significativement la composition isotopique de l'eau des rivières. En ce qui concerne le régime hydrogéochimique saisonnier des rivières, les travaux ont permis d'identifier trois périodes spécifiques pour le suivi de l'eau des rivières, soit ( 1 A1). The precisely known elevations of 114 geodetic markers were used to evaluate the precision of these DEMs. The data provided in Table A1 reveal that the DEM from the Quebec's Ministry of Energy and Natural Resources (MERN), with a resolution of 10m, provides the greatest precision with absolute error of 3.9 m.

Nevertheless, given the scale of the study region, this DEM was modified to a 100m mesh to allow for faster computer calculations. The latter, identified as MERN-B in Table A1, was used for further GIS-based calculations, while with less precision. For comparable resolution, the DEM from Shuttle Radar Topography Mission (STRM) of the NASA looks more accurate with an absolute error of 4.0 m and a standard deviation of 4.2 m. This DEM is constantly higher than the marker with generally positive error, a trend not observed in other products where more variability existed. 2. Texture association for all of the surficial units included in the ecoforestry inventories

The comprehensive surficial geologic map used for this study is from the merge of two distinct data sources as mention in section 'Data sources' in the main article.

This merge based on their respective texture ensure simplicity of product and ultimately their integration into a multicriteria analysis. In the Table A2, all units mapped and their assign class for the MCA are listed. Polygons delimiting the extent of units at the scale of 1:20,000 have been rasterized in a grid of 100 m. This step allows eliminating geometric artefacts between the polygon and the border of maps. Over the area, other products are available with field validation but their extent is too limited into the study are. The system uses ultrasonic pulses to measure the velocity of suspended sediments transported in water based on the frequency difference between the emitted and received signal [START_REF] Simpson | Discharge Measurements Using a Broad-Band Acoustic Doppler Current Profiler[END_REF]. This velocity is interpreted as an analog of the water velocity. The profiler further simultaneously calculates the depth of the riverbed and its own speed using a bottom-tracking technology [START_REF] Le Coz | Emplois de l'ADCP en rivière : une revue de synthèse[END_REF]). The ADCP profiler thus allows for measuring the total discharge within a river section when data are collected from one shoreline to the other. Here, the ADCP was attached to a boat (zodiac, canoe, or trimaran) and displacement was ensured using a small motor, paddles or ropes, depending on field conditions. A hitch was designed to attach the ADCP device to the boats, to ensure the stability of the profiler, and to keep the sensors underwater while conducting measurements. Profiles were repeated along the same sections with the objective of obtaining a reproducibility within ± 5% between four measurements. When this reproducibility target was not met, more measurements were collected. The locations of transects were selected according to safety criteria, flow conditions and to avoid human infrastructures that could significantly alter measurements.

Precipitation sampling method

Rain and snow were collected using narrow mouth 1 L HDPE bottles with 4.8 cm funnels equipped with plastic screens at their base to block debris and insects. Paraffin oil was added to the bottles to prevent evaporation. Local collaborators were responsible for collecting precipitation and sending the tightly sealed bottles by mail to the laboratory each month. Snow samples were collected in 10.1 cm inside diameter

All Weather P-2000 standard rain gauges (without cap) by the local collaborators.

Snow samples were left to melt overnight and subsequently stored in 1 L HDPE bottles with paraffin oil to prevent evaporation. Daily samples were mixed to obtain monthly composite samples. All composite precipitation samples were subsequently subsampled in 60 mL HDPE bottles and kept at 4°C until stable isotope analysis.

The five composite samples from the snowpack were collected using 1 m HDPE tubes.

Samples were left to melt at room temperature. The melted snowpack samples were then immediately used to measure the electrical conductivity (𝜎) and subsampled in 60 mL HDPE bottles that were kept at 4°C until further 𝛿 2 𝐻 -𝛿 18 𝑂 analyses.

River Water sampling method

River water samples were collected either using a small boat or directly from the bank, depending on field conditions. A rigid pole was used to collect samples approximatively 1 m below the water surface directly into HDPE bottles. Sampling bottles were rinsed three times with river water before sampling. The in situ parameters (pH, dissolved oxygen, conductivity, redox potential and temperature) were systematically measured using a YSI 556 MPS multi-parameter probe.

River ice coring method

The snow cover on the ice surface was removed before coring. Once extracted, the ice cores were removed from the ice corer and tightly wrapped in clean plastic bags using annotations to identify the depth and polarity of samples. During transport, ice cores were kept at temperatures below zero and then stored in a freezer until further subsampling procedures. After cleaning ice debris from the hole in the ice, river water under the ice cover was collected using a specially designed pole to ensure that only water from the river located approximately 1 m under the ice was collected. The pole design allows for submerging the HDPE bottle with its mouth facing downwards until the desired sampling depth is reached, and then reversing the polarity of the bottle to allow it to fill with river water. This precaution is conducted to avoid potential contamination of the river water sample by ice crystals from the coring procedure. The samples were taken directly from 60 ml HDPE bottles. The samples were stored in a refrigerator maintained at 4°C until stable isotope analyses.

The samples were collected to allow for a comparison with the corresponding bottom ice samples to evaluate the water-ice fractionation factor. These subsampling steps required for the ice cores to be conducted under freezing temperatures to avoid melting, either in a freezing room of the UQAT or outside when the air temperature was significantly below the freezing point. The ice cores were systematically cut in half lengthwise with a band saw and sub-sampled into 5 cm long half cylinders. These half cylinders were cut in layers to allow the collection of one subsample at each 5 cm interval. The external part of the subsamples was then removed using a band saw. These procedures were performed to allow for collecting small subsamples from the center of the core, while the external part of the core (potentially contaminated by under-ice water flowing into the corer during the field sampling procedure) was discarded. The extra half of the core (cut lengthwise) was kept in a freezer as a backup. The ice subsamples were stored in clean, tightly sealed bags and left to melt at 4°C in a refrigerator. As soon as the ice melted completely, the water was transferred to 30 ml HDPE bottles, tightly sealed, and kept at 4°C until stable isotope analyses.

Stable isotope analysis

Each water sample was injected (1 microliter) and measured 10 times. The first two injections of each sample were rejected to limit memory effects. ‰ ‰ ‰ ‰ WM WM-H18-G-1-A1 1 -121,9 0,4 -16,3 0,1 -8,6 0,0 8,1 WM WM-H18-G-1-B1 2 -104,6 0,3 -14,3 0,1 -7,4 0,1 9,5 WM WM-H18-G-1-B2 3 -113,5 0,1 -15,7 0,0 -8,0 0,1 11,9 WM WM-H18-G-1-B3 4 -123,1 0,3 -17,0 0,1 -8,5 0,0 12,8 WM WM-H18-1 5 -91,43 0,34 -12,87 0,08 -6,77 0,11 subice-water WM WM-H18-G-2-A1 1 -134,5 0,2 -17,9 0,0 -9,6 0,1 8,4 WM WM-H18-G-2-A2 2 -90,5 0,2 -12,2 0,0 -6,4 0,0 7,0 WM WM-H18-G-2-B1 3 -112,9 0,6 -15,1 0,1 -7,9 0,1 7,8 WM WM-H18-G-2-C1 4 -106,5 0,1 -14,1 0,0 -7,3 0,1 6,4 WM WM-H18-G-2-C2 5 -120,6 0,2 -15,8 0,1 -8,0 0,0 5,9 WM WM-H18-2 6 -91,27 0,11 -12,68 0,03 -6,63 0,06 subice-water HJ HJ-H18-G-1-A1 1 -134,8 0,2 -17,6 0,1 -9,5 0,1 6,0 HJ HJ-H18-G-1-B1 2 -130,2 0,1 -16,9 0,1 -8,9 0,0 4,7 HJ HJ-H18-G-1-C1 3 -90,3 0,3 -10,8 0,1 -5,9 0,0 -3,7 HJ HJ-H18-G-1-C2 4 -119,2 0,2 -15,7 0,2 -8,5 0,0 6,6 HJ HJ-H18-G-1-D1 5 -107,5 0,1 -14,6 0,0 -7,9 0,1 8,9 HJ HJ-H18-G-1-D2 6 -67,6 0,2 -8,6 0,1 -4,7 0,0 1,0 HJ HJ-H18-G-1-D3 7 -114,6 0,2 -15,3 0,0 -8,1 0,0 8,1 HJ HJ-H18-E-1 8 -87,33 0,10 -11,85 0,01 -6,35 0,08 subice-water HJ HJ-H18-G-2-A1 1 -87,9 0,1 -11,1 0,1 -6,0 0,0 0,9 HJ HJ-H18-G-2-B1 2 -133,3 0,4 -16,9 0,1 -9,1 0,0 2,2 HJ HJ-H18-G-2-C1 3 -200,7 0,1 -26,2 0,0 -14,2 0,1 9,3 HJ HJ-H18-G-2-C2 4 -152,7 0,5 -20,5 0,1 -10,9 0,1 11,4 HJ HJ-H18-G-2-C3 5 -99,1 0,3 -13,2 0,1 -7,0 0,0 6,1 HJ HJ-H18-G-2-D1 6 -92,1 0,6 -11,9 0,1 -6,4 0,1 3,3 HJ HJ-H18-G-2-D2 7 -123,6 1,0 -16,8 0,1 -9,0 0,0 10,7 HJ HJ-H18-G-2-D3 8 manquant HJ HJ-H18-G-2-D4 9 -69,7 0,8 -9,5 0,1 -5,1 0,1 6,0 HJ HJ-H18-G-2-E1 10 -72,6 0,4 -9,9 0,0 -5,3 0,1 6,7 HJ HJ-H18-G-2-E2 11 -70,7 1,2 -9,6 0,2 -5,1 0,1 6,1 HJ HJ-H18-E-2 12 -87,15 0,15 -11,84 0,03 -6,36 0,07 subice-water BL BL-H18-G-1-A1 1 -213,1 0,6 -29,0 0,1 -15,4 0,1 19,1 BL BL-H18-G-1-A2 2 -177,6 0,2 -23,4 0,2 -12,5 0,0 9,5 BL BL-H18-G-1-B1 3 -68,8 1,2 -9,2 0,1 -4,7 0,1 4,9 BL BL-H18-G-1-C1 4 -69,7 1,3 -9,5 0,2 -4,8 0,1 6,3 BL BL-H18-G-1-D1 5 trop mince BL BL-H18-G-1-E1 6 -67,2 1,1 -8,9 0,0 -4,6 0,1 3,8 BL BL-H18-G-1-F1 7 -68,1 0,3 -8,6 0,1 -4,7 0,1 0,8 BL BL-H18-G-1-G1 8 -68,9 0,6 -8,9 0,1 -4,8 0,2 2,1 BL BL-H18-G-1-H1 9 -68,0 1,4 -8,6 0,1 -4,8 0,1 1,1 BL BL-H18-1-surf 10 -88,64 0,27 -12,02 0,02 -6,33 0,05 subice-water NBJ NBJ-H18-G-2-A1 1 -81,5 1,1 -10,5 0,1 -5,9 0,1 2,9 NBJ NBJ-H18-G-2-A2 2 -81,1 1,8 -10,9 0,3 -5,8 0,1 5,8 NBJ NBJ-H18-G-2-A3 3 -80,5 1,1 -10,8 0,3 -5,8 0,1 5,6 NBJ NBJ-H18-G-2-A4 4 -80,6 0,3 -11,0 0,1 -5,7 0,1 7,4

Annexe G (suite) : Compilation des résultats d'analyse isotopique des échantillons de glace.

Site Échantillon_code ordre δ 2 H stdev δ 18 O stdev δ 17 O stdev Note Dxs ‰ ‰ ‰ ‰ ‰ ‰ NBJ NBJ-H18-G-2-A5 5 -78,0 0,6 -10,4 0,1 -5,5 0,0 5,1 NBJ NBJ-H18-G-2-B1 6 -77,8 0,2 -10,5 0,1 -5,4 0,1 6,4 NBJ NBJ-H18-G-2-B2 7 -77,2 0,7 -10,1 0,1 -5,8 0,1 4,0 NBJ NBJ-H18-G-2-B3 8 -73,8 1,4 -9,9 0,2 -5,3 0,2 5,6 NBJ NBJ-H18-G-2-B4 9 -72,2 0,1 -9,5 0,1 -5,0 0,1 4,0 NBJ NBJ-H18-G-2-B5 10 -73,3 0,5 -9,7 0,0 -5,0 0,0 4,4 NBJ NBJ-H18-G-2-B6 11 -72,7 0,5 -9,5 0,0 -5,0 0,1 3,6 NBJ NBJ-H18-G-2-B7 12 -71,1 0,3 -9,5 0,1 -4,7 0,0 4,9 NBJ NBJ-H18-E 13 -90,25 0,44 -12,13 0,06 -6,60 0,11 subice-water HBJ HBJ-H18-G-2-A1 1 -70,0 0,4 -9,6 0,0 -5,2 0,1 7,0 HBJ HBJ-H18-G-2-A2 2 -68,7 0,2 -9,5 0,1 -5,0 0,1 7,1 HBJ HBJ-H18-G-2-B1 3 -66,4 0,3 -9,1 0,0 -4,8 0,2 6,5 HBJ HBJ-H18-G-2-B2 4 -67,0 1,4 -9,1 0,1 -4,7 0,2 6,0 HBJ HBJ-H18-G-2-B3 5 -70,1 0,5 -9,7 0,1 -4,9 0,1 7,4 HBJ HBJ-H18-G-2-B4 6 -70,6 0,5 -9,8 0,1 -5,2 0,1 7,7 HBJ HBJ-H18-G-2-B5 7 -72,1 0,2 -9,9 0,1 -5,3 0,1 7,0 HBJ HBJ-H18-G-2-B6 8 -72,2 0,5 -10,0 0,0 -5,3 0,1 8,1 HBJ HBJ-H18-G-2-B7 9 -71,6 0,5 -9,9 0,1 -5,3 0,1 7,4 HBJ HBJ-H18-G-2-B8 10 -72,4 0,5 -9,9 0,1 -5,4 0,1 6,9 HBJ HBJ-H18-G-2-B9 11 -73,1 0,7 -10,0 0,1 -5,3 0,1 6,5 HBJ HBJ-H18-G-2-C1 12 -73,7 0,4 -10,0 0,1 -5,2 0,2 6,3 HBJ HBJ-H18-G-2-C2 13 -74,1 0,6 -10,2 0,1 -5,2 0,1 7,8 HBJ HBJ-H18-G-2-C3 14 -72,5 0,6 -9,9 0,1 -5,2 0,2 6,5 HBJ HBJ-H18-G-2-C4 15 -71,9 0,7 -9,7 0,1 -5,1 0,1 6,1 HBJ HBJ-H18-G-2-C5 16 -70,8 0,3 -9,5 0,1 -5,2 0,0 5,5 HBJ HBJ-H18-G-2-C6 17 -71,5 0,5 -9,9 0,1 -5,3 0,0 7,7 HBJ HBJ-H18-G-2-C7 18 -71,0 0,1 -9,5 0,1 -5,2 0,1 4,9 HBJ HBJ-H18-G-2-C8 19 -71,4 0,6 -9,6 0,1 -5,3 0,1 5,4 HBJ HBJ-H18-E-1 20 -90,39 0,42 -12,40 0,07 -6,64 0,02 subice-water BS BS-H18-G-1-A1 1 -136,7 0,3 -17,6 0,1 -9,5 0,0 4,2 BS BS-H18-G-1-B1 2 -98,7 0,3 -13,2 0,1 -7,1 0,0 7,3 BS BS-H18-G-1-B2 3 -91,4 0,1 -11,8 0,1 -6,5 0,0 3,4 BS BS-H18-G-1-B3 4 -104,0 0,3 -14,0 0,1 -7,6 0,0 8,0 BS BS-H18-G-1-B4 5 -93,9 0,7 -12,1 0,2 -6,6 0,0 3,0 BS BS-H18-G-1-B5 6 -93,3 0,4 -11,8 0,2 -6,2 0,0 1,1 BS BS-H18-E 7 -85,20 0,07 -11,53 0,02 -6,14 0,09 subice-water HA HA-H18-G-1-A1 1 -101,1 0,5 -13,8 0,1 -7,3 0,1 9,1 HA HA-H18-G-1-A2 2 -114,0 0,1 -15,3 0,1 -8,2 0,1 8,1 HA HA-H18-G-1-A3 3 -91,1 0,4 -12,3 0,0 -6,6 0,0 7,0 HA HA-H18-G-1-B1 4 -100,3 0,4 -13,5 0,1 -7,2 0,0 7,5 HA HA-H18-1 5 -85,95 0,26 -11,93 0,03 -6,17 0,07 subice-water HA HA-H18-G-2-A1 1 -92,6 0,2 -12,5 0,1 -6,8 0,0 7,3 HA HA-H18-2 2 -85,68 0,12 -12,04 0,10 -6,31 0,07 subice-water HV HV-H18-G-1-A1 1 -108,3 0,4 -14,2 0,2 -7,7 0,1 5,1 HV HV-H18-G-1-A2 2 -106,6 0,2 -13,9 0,1 -7,7 0,1 4,7

Annexe G (suite) : Compilation des résultats d'analyse isotopique des échantillons de glace.

Site Échantillon_code ordre δ 2 H stdev δ 18 O stdev δ 17 O stdev Note Dxs ‰ ‰ ‰ ‰ ‰ ‰ HV HV-H18-G-1-B1 3 -71,4 0,3 -9,0 0,1 -4,9 0,0 0,8 HV HV-H18-E 4 -85,92 0,24 -11,62 0,07 -6,48 0,06 subice-water NBJ NBJ-H18-G-1-A1 1 -75,2 0,4 -9,8 0,2 -5,4 0,1 2,9 NBJ NBJ-H18-G-1-A2 2 manquant NBJ NBJ-H18-G-1-A3 3 -75,9 0,4 -10,1 0,2 -5,5 0,0 4,6 NBJ NBJ-H18-G-1-A4 4 -74,9 0,6 -10,0 0,1 -5,4 0,0 4,9 NBJ NBJ-H18-G-1-A4 5 -76,2 0,3 -9,9 0,1 -5,4 0,0 doublon 3,2 NBJ NBJ-H18-G-1-B1 6 -75,5 0,1 -9,7 0,1 -5,4 0,1 2,2 NBJ NBJ-H18-G-1-B2 7 -75,7 0,8 -10,0 0,2 -5,4 0,1 4,5 NBJ NBJ-H18-G-1-B3 8 -76,9 0,6 -10,5 0,2 -5,3 0,1 7,0 NBJ NBJ-H18-G-1-B4 9 -74,1 0,2 -10,3 0,1 -5,1 0,0 8,3 NBJ NBJ-H18-G-1-C1 10 -76,3 0,3 -10,3 0,3 -5,3 0,0 6,5 NBJ NBJ-H18-G-1-C2 11 -73,7 0,3 -10,3 0,1 -5,2 0,1 8,7 NBJ NBJ-H18-G-1-C3 12 -73,5 0,5 -9,7 0,3 -5,1 0,1 4,2 NBJ NBJ-H18-G-1-C4 13 -73,3 0,5 -9,7 0,2 -5,2 0,1 4,6 NBJ NBJ-H18-G-1-D1 14 -75,2 0,6 -10,1 0,2 -5,4 0,1 5,4 NBJ NBJ-H18-G-1-D2 15 -76,9 0,5 -10,4 0,1 -5,5 0,0 6,2 NBJ NBJ-H18-G-1-E1 16 -78,4 0,4 -10,7 0,1 -5,5 0,1 7,5 NBJ NBJ-H18-G-1-E2 17 -77,1 0,4 -10,2 0,1 -5,4 0,0 4,4 NBJ NBJ-H18-G-1-F1 18 -74,7 0,6 -9,7 0,1 -5,2 0,1 2,5 NBJ NBJ-H18-G-1-F2 19 -75,0 0,4 -9,5 0,2 -5,1 0,0 0,7 NBJ NBJ-H18-G-1-G1 20 -77,0 0,4 -10,2 0,1 -5,4 0,0 4,8 NBJ NBJ-H18-G-1-G2 21 -77,3 0,2 -10,3 0,2 -5,4 0,0 5,4 NBJ NBJ-H18-G-1-H1 22 -75,3 0,7 -10,1 0,1 -5,2 0,1 5,5 NBJ NBJ-H18-G-1-H1 23 -71,7 0,4 -10,4 0,1 -5,3 0,0 doublon 11,8 NBJ NBJ-H18-G-1-H2 24 -75,5 0,4 -10,2 0,1 -5,4 0,0 6,4 NBJ NBJ-H18-E 25 -90,25 0,44 -12,13 0,06 -6,60 0,11 subice-water HBJ HBJ-H18-G-1-A1 1 manquant HBJ HBJ-H18-G-1-A2 2 -69,9 0,1 -10,0 0,1 -5,2 0,0 10,0 HBJ HBJ-H18-G-1-A3 3 -68,4 0,6 -9,7 0,1 -5,1 0,0 9,0 HBJ HBJ-H18-G-1-A4 4 -69,6 0,3 -9,6 0,1 -5,0 0,1 7,2 HBJ HBJ-H18-G-1-B1 5 -70,6 0,5 -9,9 0,3 -5,1 0,1 8,6 HBJ HBJ-H18-G-1-B2 6 manquant HBJ HBJ-H18-G-1-B3 7 manquant HBJ HBJ-H18-G-1-B4 8 -72,7 1,0 -10,3 0,0 -5,2 0,0 10,1 HBJ HBJ-H18-G-1-C1 9 -73,0 0,6 -9,8 0,2 -5,2 0,2 5,2 HBJ HBJ-H18-G-1-C2 10 -74,1 0,4 -10,0 0,3 -5,2 0,1 5,7 HBJ HBJ-H18-G-1-C3 11 -71,7 0,5 -9,6 0,2 -5,3 0,2 5,3 HBJ HBJ-H18-G-1-C2 12 -72,9 0,4 -10,0 0,1 -5,5 0,0 doublon 6,9 HBJ HBJ-H18-G-1-C3 13 -75,9 0,2 -10,8 0,2 -5,6 0,0 doublon 10,4 HBJ HBJ-H18-G-1-cul 14 -71,6 0,1 -9,5 0,1 -5,1 0,3 4,7 HBJ HBJ-H18-E-2 15 -90,04 0,23 -12,69 0,02 -6,63 0,08 subice-water ORBJ PACES-CNE-156 1 -118,9 0,6 -15,9 0,1 0 à 3.6 cm 8,0 ORBJ PACES-CNE-157 2 -112,6 0,2 -14,9 0,1 3.6 à 7.2 cm 6,6 ORBJ PACES-CNE-158 3 -107,6 0,4 -14,2 0,0 7.2 à 10.8 cm 5,9 Quantifier la ressource en eau stockée dans les aquifères peu profonds et déterminer leurs interactions avec l'eau de surface est essentiel pour mieux évaluer la sensibilité de l'eau souterraine et des écosystèmes qui en dépendent découlant du changement climatique et des impacts de l'activité humaine. Cependant, l'évaluation de la ressource en eau est complexe sur des territoires vastes éloignés où les données préexistantes sont rares et limitées. Cette étude se concentre sur l'utilisation d'approches SIG à grande échelle pour produire de nouvelles informations sur les ressources en eaux souterraines peu profondes. La zone d'étude couvre 185 000 km 2 englobant de vastes bassins versants boréaux du bassin de la Baie-James (Canada). Une analyse spatiale multicritères (AMC) basée sur des données géologiques et hydrogéologiques a été couplée à des calculs géométriques afin d'évaluer le volume des aquifères granulaires peu profonds et d'identifier les zones où la protection des eaux souterraines doit être prioritaire. Des données satellitaires et des modèles spatiaux combinés à des données in situ, telles que des mesures de débit et la composition géochimique de l'eau (isotopes stables [δ 2 H-δ 18 O] et conductivité électrique), ont ensuite été utilisés pour calculer le bilan hydrologique terrestre et estimer les composantes du cycle de l'eau. Enfin, un nouvel indice basé sur la stabilité de l'indice de sécheresse Température-Végétation, calculé à partir de 290 images du satellite Landsat 8, a révélé l'emplacement de potentiels occurrences d'écosystèmes dépendant des eaux souterraines (POEDES).

Les résultats ont révélé que les aquifères potentiels (AP) les plus importants sont mal distribués et principalement associés à des formations fluvioglaciaires. Ils contiennent en moyenne 8,5 m de sédiments saturés pour une évaluation globale d'environ 40 km 3 ± 4 km 3 d'eau sur 7% du territoire. La contribution annuelle des eaux souterraines au débit des rivières est généralement inférieure à 12% des précipitations d'un bassin versant atteignant 2,3 km 3 sur le territoire. La contribution de l'eau provenant des pluies et de la fonte des neiges est de 30 à 61% et de 18 à 40%, respectivement, avec des ratios évaporation sur les apports de 2 à 10%. Le débit présente des cycles saisonniers avec une contribution plus élevée d'eau de fonte des neiges pendant la crue printanière et de l'eau de pluie à l'automne, une composition évoluant vers celle des eaux souterraines pendant la période de glace et une contribution stable de l'eau souterraine tout au long de l'année. De plus, le bilan de masse isotopique suggère que les plantes transpirent principalement de l'eau issue des précipitations. Les terres proches des cours d'eau dans les vallées principales présentent les plus haut POEDES. Une analyse en composantes principales des caractéristiques des bassins versants révèle des corrélations significatives entre les paramètres des eaux souterraines estimés au cours de cette étude. Une topologie de bassins versants est proposée à partir des deux premières composantes qui expliquent 77% de la variance. A partir de celle-ci, les bassins versants au sud du territoire dépendent plus des écoulements souterrains que ceux situés au nord qui montrent une surface plus imperméable. Quatre regroupements de bassins versants identifient les secteurs où une vigilance accrue des eaux souterraines est nécessaire pour préserver la ressource en eau et la biodiversité qui dépend des eaux souterraines.

  isotopic composition of precipitation and observed isotopic composition of average river discharge. b) 3D graphical representation of the isotopic composition and electrical conductivity of water. The positions of endmembers needed for calculations are also shown. c) 𝛿 and 𝜎 times-series in monitored river……………………………………………………………………3.8Annual partitioning of monitored watershed based on the endmembers further express in proportion of inflow in the Table 4. b) Schematic representation of annual natural pattern of streamflow sources contribution. The thickness of the blue cycle refers to a qualitative representation of the discharge pattern. c) Hydrograph the Harricana River in Amos with partition of the streamflow components. zone d'étude en milieu boréal. b) Identification des sous-bassins versants caractérisés. c) Identification des images Landsat 8 traitées…………………….……………………………….. l'indice Temperature-vegetation Dryness index ………………………………………………...………………... 4.5 Répartition du nombre d'images satellites Landsat 8 superposées répondant aux critères de sélection………………………………….. 4.6 Distribution spatiale du Temperature-vegetation Dryness index moyen pour l'ensemble du territoire d'étude et distribution spatiale du décile de l'écart-type des Temperature-vegetation Dryness index calculés sur chaque pixel du territoire d'étude.…………………………………….. 4.7 Distribution spatiale de l'indice du potentiel d'occurrence d'écosystème dépendant de l'eau souterraine (POEDES) sur le territoire d'étude. ……………………………………………………... 4.8 Représentation graphique des deux premières CP réalisées sur l'ensemble des variables intégrées dans l'analyse statistique détaillées dans le Tableau 4.1………………………………………………..….. 4.9 Représentation graphique des deux premières composantes principales (CP) réalisées sur les variables suite aux itérations et représentation graphique de la 2 e et 3 e CP réalisées sur les variables suite aux itérations. ……………………........................................................ ...... 4.10 Étendue de la moyenne des indices du potentiel d'occurrence d'écosystème dépendant de l'eau souterraine (POEDES) calculés pour les codes de végétation identifiés dans la cartographie écoforestière réalisée sur le territoire……………………….……………………….. 4.11 Typologie des bassins versants et positionnement des sous-bassins versants…………………………………………………………….….

Figure 1 .

 1 Figure 1.1 a) Distribution spatiale par continent du volume d'eau souterraine à l'échelle du globe (valeur exprimée en km 3 ). b) Distribution des aquifères de l'Amérique du Nord et recharge associée (Données tirées du World-wide Hydrogeological Mapping and Assessment Programme [WHYMAP]).

1. 3

 3 Le contexte de la Baie-James En milieu boréal, les grands espaces, la faible densité de population et la présence de grands fleuves ont favorisé l'implantation de grands complexes hydroélectriques. Ces complexes ont modifié l'écoulement de l'eau de surface notamment dans le grand bassin hydrologique des baies de James, d'Hudson et d'Ungava. Le transfert d'eau interbassins, qui consiste à dévier une quantité significative d'eau d'un bassin versant dans le bassin versant voisin pour augmenter la capacité d'un réservoir, fait augmenter le débit des certains cours d'eau au détriment d'autres. Il en résulte notamment des changements dans la salinisation des estuaires, dans le processus d'érosion, une dégradation de la qualité de l'eau disponible et une modification des interactions eau

  Figure 1.2 Représentation conceptuelle des éléments traités dans la présente thèse. N = Neige, P = Pluie, E = Évaporation, T = Transpiration, Q = Débit des rivières, PEDES = potentiels écosystèmes dépendant de l'eau souterraine.

  Figure 2.1 Study Area. a) Included watersheds. Legend: 1. Harricana River 2. Nottaway River 3. Broadback River 4. Rupert River 5. Pontax River 6. Eastmain River. b) Geology of Superior Province and position of the Continental water divide; c) Synthesis of main components of the last glacial retreat and post-glacial events. The extent of the Tyrell sea transgression, limit of Cochrane surge and Barlow-Ojibway submersion are based on previous work by Dionne (1972), Dubé-Loubert et al. (2012) and Veillette (1994), respectively.

  .1c) resulted in the deposition of unconsolidated deposits, presenting a complex architecture that dictates the geometry of shallow groundwater flow systems. The oldest unconsolidated unit of the region covering bedrock (Figure2.2a) consists of a compact coarse till(Hardy 1982a; Figure 2.2b). Till outcrops are commonly observed in high grounds located near bedrock outcrops, at elevations exceeding the maximum elevation reached by younger sediments. During the last deglaciation, northwestern Quebec was located at the boundary between two large ice masses-namely, the Nouveau-Québec glacier, receding to the northeast, and the Hudson glacier, receding to the northwest(Hardy 1982a). One of the most striking features resulting from this setting is the Harricana interlobate moraine, which extends more than 500 km (Figure2.1c).

Figure 2 .

 2 Figure 2.2 South-North composite stratigraphic column of the study area based on Hardy (1982a). The thickness of sediment layers shown in the figure is unknown and not representative of actual thickness. Legend: a. Bedrock, b. Glacial till, c. Glaciofluvial sands and gravels, d. Fine-grained, deep-water glaciolacustrine sediments, e. Glaciolacustrine, aeolian and alluvial coastal sands, f. Organic sediments, g. Cochrane till, h. Marine clay, i. Coastal and alluvial sands.

  Figure 2.3 Flowchart of the main methodological steps of this study.

Figure 2 . 4

 24 Figure 2.4 Evaluation of unconsolidated sediments thickness. a) Method for evaluating the average increase in bedrock/till surface elevation (ΔEA) as a function of the distance to the closest contact with sediments; b) Strategy for evaluating the average increase in unconsolidated sediments thickness (ΔEB) as a function of the distance to the closest bedrock/till outcrop; c) Estimated trends representing average surface roughness (∆E ̅̅̅̅ ) as a function of the distance to the contact between sediments and bedrock/till outcrops. This trend is a regional analog of the unconsolidated sediments thickness (UST).

1 -

 1 The surficial units mapped as fine-grained, deep-water sediments (either glaciolacustrine or marine) and Cochrane till deposits are part of the clay extent; 2-Clay deposits could be buried in areas where younger unconsolidated sediments are outcropping; 3-The Sakami moraine represents the eastern boundary of the glaciolacustrine Barlow-Ojibway submersion extent; 4-In Lake Barlow-Ojibway, clay deposits are restricted to areas where the depth of water was ≥ 50 m (Veillette 1983); 5-The maximal elevation of clay deposits compared to the shoreline elevation reached by the Tyrrell Sea submersion suggests that marine clay deposits are restricted to areas where the depth of the Tyrrell Sea water was ≥ 30 m.

Figure 2 . 5

 25 Figure 2.5 Proposed ranges of hydraulic conductivity for various sediments and for the fractured bedrock (The bold line represents the trend of hydraulic conductivity used for standardisation). Legend: (A) Cloutier et al. (2016); (B) Freeze and Cherry (1979); (C) Domenico and Schwartz (1990).

2. 5

 5 Figure 2.6 Unconsolidated sediments thicknesses used for identifying groundwater flow conditions (Unconfined aquifers are represented in green and confined aquifers are in brown).

Figure 2 . 7 1

 271 Figure 2.7 Comparison of unconsolidated sediments thickness (UST) from high quality boreholes (Allard and Dubé-Loubert 2016) with the predicted trend of regional UST (see Eq. 1).

Figure 2 . 8

 28 Figure 2.8 Maps of standardized values of criteria needed for Aquifer Potential calculations. a) T': Potential transmissivity; b) DFHN: Distance from hydrographic network; c) S: Slope; d) DND: Drainage network density; e) DMBF: Distance to a major bedrock fault; f) TWI: Topographic wetness index.

Figure 2 .

 2 Figure 2.10 Calculated area of saturated sediments thickness (SST) for the five classes of aquifer potential (AP) in unconfined aquifers. UST: Unconfined sediment thickness.

(

  Figure 3.1a). It includes the watersheds of the Harricana, Nottaway, Rupert, Broadback, Pontax, and Eastmain Rivers of the James Bay Basin within the Boreal Shield, Taiga Shield, and Hudson Plains ecozones[START_REF] Marshall | A National Ecological Framework for Canada: Attribute data. Ecosystem Stratification Working Group[END_REF]. Under natural conditions, the rivers of the area generally present a bimodal hydrograph with lowest discharge at the end of the ice-on period (winter baseflow), highest discharge during the spring snowmelt followed by summer baseflow, and slightly increased discharge in the fall season. The northern part of the study area is characterized by extensive hydroelectric installations that modify the hydrology of large watercourses such as the Rupert, Eastmain, Opinaca, and La Grande Rivers. In these rivers, temporal changes in discharge are mainly a function of water management by Hydro-Quebec. The Harricana, Nottaway, Broadback, and Pontax Rivers are among the last naturally flowing, non-fragmented large boreal rivers of the James Bay area[START_REF] Dynesius | Fragmentation and Flow Regulation of River Systems in the Northern Third of the[END_REF][START_REF] Nilsson | Fragmentation and Flow Regulation of the World&#039[END_REF][START_REF] Grill | Mapping the world's free-flowing rivers[END_REF]).

Figure 3 .

 3 Figure 3.1 a) Monitored watersheds and outlets, data acquisition and extent of the aquifer type over study area. The site HA is the intensively monitored outlet. H = Harricana river. N = Nottaway river. B = Broadback river. R = Rupert river. P = Pontax river. E = Eastmain river. ADCP = Acoustic Doppler Current Profiler. UBA = Unconfined bedrock aquifer. UGA = Unconfined granular aquifer. CA = Confined aquifer. b) Location and bedrock geology of the study area.

Figure 3 . 3

 33 Figure 3.3 Conceptual representation of the annual geochemical cycle of river discharge (loop) (modified from St Amour et al. 2005). The streamflow isotopic composition (𝛿 𝑅𝐼𝑉 ) results from the mixing of water from rain (𝛿 𝑅 ), snowmelt (𝛿 𝑆 ), and groundwater (𝛿 𝑆𝐺𝑊 ). Evaporation causes heavy isotope enrichment along the LEL and explains the offset. The blue drawing is a loop projected in a 3D space. The drawings in red are projections on the rear faces of the 3D graph. The drawings in black are the projection on the basal plan of the 3D graph. Black dots are used to identify the rain (𝑅), snow (𝑆), shallow groundwater (𝑆𝐺𝑊), and evaporation (E) end-members. The star symbol represents the flux-weighted average composition of streamflow. The dashed lines represent the theoretical mixing triangle formed by the 𝑆, 𝑅, and 𝑆𝐺𝑊 end-members, before evaporation.

Figure 3 .

 3 Figure 3.3 provides a conceptual representation of the annual geochemical pattern of river discharge. This graph shows the influence of different water fluxes on the isotopic composition and electrical conductivity of rivers. The basal plane consists of a 𝛿 2 𝐻 -𝛿 18 𝑂 graph on which the LMWL and 𝐿𝐸𝐿 are represented. The 𝛿 𝑅 and 𝛿 𝑆 both plot on the 𝐿𝑀𝑊𝐿, while 𝛿 𝑆𝐺𝑊 plots at the intersection between the 𝐿𝐸𝐿 and 𝐿𝑀𝑊𝐿. The 𝛿 𝐸 plots to the left of the 𝐿𝑀𝑊𝐿 and is aligned with the LEL. The two planes constituting the rear faces of the 3D block illustrate the relationship between the 𝛿 and 𝜎 of water.

  3.4.1 Hydrological Measurements and Data Sources 3.4.1.1 Atmospheric fluxes from remote sensing and models Precipitation (rain and snow) and evapotranspiration (ET) can be estimated based on remote sensing datasets and spatial models. Details of the data used are provided in Table 3.1. These datasets are aggregated to obtain mean and range of the cumulative values of monthly water inflows and outflows for each of the studied watersheds. The precipitation fluxes from the beginning of November to the end of March are assumed to represent the snow water equivalent of the snowpack.

  flow are provided in Appendix B-1. Overall, 53 stream flows were measured according to 236 profiles realized at 22 sites in 2018-2019 (Figure 3.1a; Appendix C). In addition, three hydrometric stations operated by the Government of Quebec (MELCC 2020) are located in the study area (HA, BM, WM; Figure 3.1a). A total of 27 ADCP profiles were carried out near these hydrometric stations to allow for comparisons. Discharge measurements were scheduled to collect data during the spring freshet, the summer baseflow, and the autumn high waters. The baseflow at the end of winter (ice-on period) is assumed to be proportional to the watershed area according to:

  Figure 3.4 Monthly cumulative water fluxes estimated from remote sensing data and spatial models. The solid lines present average values calculated for entire study area, while the envelopes illustrate the range of average values calculated for monitored subwatersheds.

Figure 3 . 5

 35 Figure 3.5 Temporal variations in river discharge. The solid black line represents the average river flow, estimated based on data from hydrometric stations. The solid green line is interpolated river flow based on field measurements (black points). The baseflow values (red points) are calculated from Equation1and are assigned at the end of the ice-on period (mid-April).

Figure 3 .

 3 Figure 3.7 a) Calculated and observed isotopic composition of precipitation and observed isotopic composition of average river discharge. b) 3D graphical representation of the isotopic composition and electrical conductivity of water. The positions of endmembers needed for calculations (rain (𝑅), snowmelt (𝑆), groundwater (𝐺), and evaporated moisture (𝐸)) are also shown. c) 𝛿 and 𝜎 times-series in monitored river. The orange line is associated with site HA, which is the site with most samples, while black lines are used for all other sites.

  𝐷𝐺𝑊 and 𝛥𝑠 terms of Equation 4 likely have a limited impact on the annual hydrological balance of watersheds. The annual discharge could be assessed by 𝐼 -𝐸𝑇 = 𝑄 𝑅𝐼𝑉 , and the groundwater fluxes directly reaching the James Bay through the 𝑄 𝐷𝐺𝑊 are negligible in comparison to freshwater fluxes from rivers. These interpretations are consistent with previous work by Luijendijk et al. (2020), who found that the groundwater discharge into the James Bay from coastal watersheds of study area reached a maximum of 0.14 km 3 /yr, corresponding to <0.1% of 𝐼. These findings further suggest that the role of groundwater fluxes in the hydrological balance of watersheds is mainly associated with processes where water transfers from recharge areas to rivers through shallow groundwater flow systems. Under such conditions,

  3.6.3 Deciphering the sources of river dischargeBased on the assumption that 𝜌 is negligible, the cumulative discharge from the studied watersheds can be estimated from river discharge, with values reaching 109 km 3 in 2018 and 100 km 3 in 2019. These fluxes correspond to 65% and 62% of the total precipitation (𝐼) for the same years. The remaining fraction is attributed to evapotranspiration. In comparison,[START_REF] Déry | Interannual variability and interdecadal trends in Hudson Bay streamflow[END_REF] estimated an annual mean discharge of 96.5 km 3 for the same rivers between 1964 and 2008. The sources of river runoff are evaluated based on Equation 12. The 𝛿 𝑅 and 𝛿 𝑆 are estimated for each watershed based on flux-weighted data from the RCWIP[START_REF] Terzer | Global isoscapes for δ 18 O and δ 2 H in precipitation: improved prediction using regionalized climatic regression models[END_REF]. The 𝜎 𝑅 and 𝜎 𝑆 are estimated based on field data and assumed to be homogeneous in space and constant in time. The isotopic composition of the evaporating moisture is estimated based on Equation8. The 𝛿 𝐼 is used as an analog of the 𝛿 𝑆𝐺𝑊 , consistent with observations realized near the study area[START_REF] Castelli | Hydrogéochimie des sources associées aux eskers de l'Abitibi[END_REF][START_REF] Rey | Using water stable isotopes for tracing surface and groundwater flow systems in the Barlow-Ojibway Clay Belt, Quebec, Canada[END_REF]) and with assumptions used in previous studies[START_REF] St Amour | Isotopic time-series partitioning of streamflow components in wetlanddominated catchments, lower Liard River basin, Northwest Territories, Canada[END_REF][START_REF] Gibson | Isotope-based partitioning of streamflow in the oil sands region, northern Alberta: Towards a monitoring strategy for assessing flow sources and water quality controls[END_REF].The 𝜎 𝑆𝐺𝑊 is estimated based on the data pertaining to the 𝜎 of groundwater in the main types of aquifers of the Canadian Shield in Abitibi-Témiscamingue. These aquifers include unconfined granular aquifers (UGA), unconfined bedrock aquifers (UBA), and confined aquifers (CA), with respective approximated 𝜎 values of ~150 μs/cm, ~400 μs/cm, and ~500 μs/cm (based on results from[START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF][START_REF] Rey | Using water stable isotopes for tracing surface and groundwater flow systems in the Barlow-Ojibway Clay Belt, Quebec, Canada[END_REF]. For each watershed, the proportion of each type of aquifer (in terms of % coverage in the 𝑥𝑦 plane) are estimated based on the groundwater flow conditions data provided in section 2.5.1. For calculating these proportions, only the shallowest aquifer is considered in every point of the territory (Figure3.1a). For example, the bedrock aquifers are only considered if no granular aquifer is found above. The value assigned to 𝜎 𝑆𝐺𝑊 is expected to represent the composite electrical conductivity of groundwater in unconfined aquifers. This seems realistic as groundwater/surface-water hydrological exchanges are likely to occur preferentially in shallow unconfined aquifer areas, whereas the contribution of deeper confined aquifers to surface waters might be restricted to deep groundwater discharge zones found in large lakes and rivers. Previous studies also suggested that the groundwater contribution from unconfined granular aquifers and peatlands largely contribute to river discharge in James Bay lowlands[START_REF] Orlova | Surface Water and Groundwater Contributions to Streamflow in the James Bay Lowland, Canada. Arctic[END_REF]. It thus seems realistic to propose that the shallower groundwater flow systems found in unconfined aquifers are more dynamic than the deeper groundwater flow systems found in confined aquifers. This interpretation is also in agreement with a previous study suggesting that the average residence time of groundwater in confined aquifers of the Abitibi-Témiscamingue region can reach thousands of years[START_REF] Boucher | Groundwater age investigation of eskers in the Amos region, Quebec, Canada[END_REF], while most geochemical results reveal that groundwaters found in the unconfined aquifers of the study area show low dissolved solids contents, suggesting short water residence times[START_REF] Cloutier | Atlas hydrogéologique de l'Abitibi-Témiscamingue[END_REF]). While it is not feasible to precisely constrain the 𝜎 𝑆𝐺𝑊 value, three scenarios are tested for solving Equation 12 based on different endmember values for the 𝜎 𝑆𝐺𝑊 (Table

  waters(Figure 3.8a). The cycle is illustrated conceptually in Figure3.8b, where the iceon baseflow conditions are characterized by the highest groundwater contributions. Following the ice-on baseflow, the relative contribution of meltwaters can reach up to 80% of the discharge during the spring freshet (Figure3.8b-2-3).During the ice-off season, the discharge gradually decreases to reach summer baseflow conditions. During this period, the relative contribution of the rain endmember gradually increases and a slight increase in the shallow groundwater contribution is also observed (Figure3.8b-4). The onset of the fall highwaters following the ice-off baseflow conditions then reveals a period during which the contribution from rainwaters can reach 80% (Figure3.8b-5). The hydrograph separation for site HA provides a summary of the observations (Figure3.8c). The hydrograph separation reveals that the groundwater contribution to river discharge is stable during the year compared to the contributions from the rain and snow endmembers. The contribution from runoff originating from snowmelt is observed through the entire ice-off period, while that of runoff originating from rainwater is bimodal, with an increased relative contribution in the fall period. The evaporation flux is observed through the entire iceoff period and reaches a maximum value in the fall, when waters that have undergone summer evaporation are drained by rivers.

Figure 3 .

 3 Figure 3.8 a) Annual partitioning of monitored watershed based on the endmembers (Rain (R), Snowmelt (S), Shallow groundwater (SGW) further express in proportion of total inflow (𝐼) in Table 5. b) Schematic representation of annual natural pattern of streamflow sources contribution. The thickness of the blue cycle refers to a qualitative representation of the discharge pattern. The blue gradient refers to the evolution of the shallow groundwater proportion. c) Hydrograph the Harricana River in Amos (site HA) with partition of the streamflow components (endmembers). d) Relative contribution of runoff sources of all sample. Empty circle refers to partitioning streamflow for uncompleted annual geochemical data that not allowing calculation of 𝛿 𝐸 . Colors refer to season. Purple = Spring, Green = Summer, Orange = Autumn, Blue = Winter. The black line connects the sample of the Harricana River at the Amos (site HA).

  congelation ice, as illustrated at the left of the conceptual model. The zones where the rivers are wide should be avoided as there might be a channel of active flow and zones of stagnant water below the ice. The zones corresponding to the patterns depicted on the right side of the conceptual model should also be avoided. The ice does not provide much information on the isotopic composition of the river water during the ice-on season.

Figure 3 . 9

 39 Figure 3.9 Schematic representation of ice formation and associated isotopic composition. Despite the stable isotopic composition of baseflow (𝛿 𝐴 = 𝛿 𝐶 = 𝛿 𝐸 ), only congelation ice could be use to reconstruct isotopic composition of river during ice-on. The 𝛿 𝐵 , 𝛿 𝐷 , 𝑎𝑛𝑑 𝛿 𝐹 refer to the isotopic composition of the ice core.

1 -

 1 canadien. Plus spécifiquement, l'étude vise à: Développer un indice représentatif du potentiel d'occurrence d'écosystèmes dépendant de l'eau souterraine (POEDES) sur la base d'un suivi temporel de l'indice TVDI; 2-Identifier les corrélations entre les variables hydro-géo-topographiques des sous-bassins versants montrant les tendances favorables aux interactions eau de surfaceeau souterraine; 3-Évaluer et cartographier la diversité hydro-géo-topographique des bassins versants dans un optique de préservation de la ressource en eau souterraine. Une analyse en composantes principales des caractéristiques des bassins versants, qui permet la classification de bassins versant à grande échelle (Park et al. 2021), servira de base à l'élaboration d'une typologie permettant de définir cette diversité. Le présent chapitre est complémentaire à l'identification des aquifères potentiels (identifiés au chapitre 2; Nadeau et al. 2021) et de la contribution de l'eau souterraine au réseau hydrographique (calculée au Chapitre 3). Il est nécessaire pour évaluer le niveau de dépendance des bassins versants et des écosytèmes à l'eau souterraine. Ces connaissances seront utilisées afin d'émettre des recommandations relatives à la gestion du territoire en vue de limiter les impacts des changements climatiques et du développement anthropique sur la biodiversité dépendante de l'eau souterraine. Elles serviront aussi d'outils pour orienter les futures démarches visant la préservation de EDES pour l'ensemble du territoire boréal. 4.2 Les EDES en milieu boréal Les EDES se définissent comme étant les milieux dont la composition, la structure et la fonction relève de l'approvisionnement en eau souterraine. Cela correspond à l'appellation anglo-saxonne largement employée Groundwater dependent ecosystems (GDE; Klove et al. 2011a). La Figure 4.1 conceptualise les différents milieux susceptibles d'héberger des EDES au sein d'un bassin versant. Ces EDES peuvent se développer dans les aquifères, profiter de la résurgence d'eau souterraine à la surface du sol ou s'approvisionner de l'eau souterraine peu profonde par le développement de leur système racinaire[START_REF] Eamus | A functional methodology for determining the groundwater regime needed to maintain the health of groundwater-dependent vegetation[END_REF]). L'intensité des interactions et des connections entre l'eau souterraine et l'eau de surface varie dans le temps, dans l'espace[START_REF] Amoros | Les concepts d'hydrosystème et de secteur fonctionnel dans l'analyse des systèmes fluviaux à l'échelle des écocomplexes[END_REF]) et à différentes échelles[START_REF] Griebler | Ecological assessment of groundwater ecosystems -Vision or illusion?[END_REF]

Figure 4 . 1

 41 Figure 4.1 Schéma conceptuel des milieux susceptibles d'héberger des écosystèmes dépendant de l'eau souterraine. (1) Résurgence le long des systèmes de failles. (2) Résurgence gravitaire ponctuel et (3) diffuse dans les dépôts meubles. (4) Nappe accessible par le système racinaire des plantes. (5) Exfiltration d'eau dans les rivières et (6) les lacs. (7) Apport régional d'eau souterraine dans la dynamique hydrologique des estuaires. (8) Résurgence diffuse le long du littoral.

Figure 4 .

 4 Figure 4.2 a) Localisation de la zone d'étude en milieu boréal. b) Identification des sous-bassins versants caractérisés. Ils s'insèrent dans les grands bassins versants de la rivière Harricana (orange), de la Rivière Nottaway (vert), de la Rivière Broadback (bleu), de la rivière Rupert (rouge), de la rivière Pontax (jaune) et de la rivière Eastmain (violet). c) Identification des images Landsat 8 traitées.

Figure 4 .

 4 Figure 4.2 a) Localisation de la zone d'étude en milieu boréal. c) Identification des images Landsat 8 traitées. c) Identification des sous-bassins versants caractérisés. Ils s'insèrent dans les grands bassins versants de la rivière Harricana (orange), de la Rivière Nottaway (vert), de la Rivière Broadback (bleu), de la rivière Rupert (rouge), de la rivière Pontax (jaune) et de la rivière Eastmain (violet).

Figure 4 . 3

 43 Figure 4.3 Organigramme illustrant la méthode adoptée.

Figure 4 ( 1 )

 41 Figure 4.4 : Représentation graphique l'indice TVDI dans l'espace de la température du sol (LST) en fonction de la végétation (NDVI). Le gradué de couleur réfère à l'interprétation du niveau de sécheresse du sol où banc est indicateur de sécheresse et bleu foncé indicateur d'un sol humide. (Modifiée de Standholt et al. 2002)Figure 4.3 Organigramme illustrant la méthode adoptée.

L

  'analyse en composantes principales (ACP) est utilisée pour révéler l'existence de corrélations entre plusieurs variables quantitatives calculées pour plusieurs entités. Les variables représentent ici les différentes caractéristiques hydro-géo-topographiques (Tableau 4.1) tandis que les entités font références aux sous-bassins versants. Les valeurs des variables calculées pour les sous-bassins versants sont centrées réduites pour uniformiser leur unité tout en conservant leur variance. L'ACP permet d'extraire l'information contenue dans le jeu de données et de la représenter de façon synthétique sous la forme de composantes principales. L'analyse permet d'illustrer la variance expliquée par un éventail de paramètres en seulement quelques composantes. La représentation graphique de ces composantes permet de décrire les interrelations entre les variables par l'utilisation de vecteurs. L'orientation et la longueur de ces vecteurs définissent le niveau de corrélation entre le plan des composantes principales et la variable du jeu de données. À terme, la démarche vise à identifier les corrélations entre le POEDES et les autres variables. L'analyse permet de positionner les bassins versants par rapport aux axes des composantes principales montrant les particularités qui les caractérisent. Des similitudes et différences émergent pour former une typologie de bassins versants basée sur l'interprétation de la signification des CP. 4.5 Résultats 4.5.1 Les indices spatiaux Au total, 290 images satellites Landsat 8 ont été traitées pour obtenir le NDVI, la LST et le TVDI à l'échelle de la zone d'étude. La Figure 4.5 montre la superposition des images répondant aux critères de sélection pour la période de référence sur l'ensemble du territoire. Certains secteurs ont fait l'objet d'une superposition de 73 images. La totalité du territoire terrestre a pu être couverte par un minimum de 6 images superposées. Le NDVI varie entre 0 et 0,6 sur l'ensemble du territoire et les LST calculées atteignent 40°C avec une valeur moyenne de 15°C. Le Tableau 4.2 fournit les résultats des composantes du calcul de TVDI détaillées par les équations 1, 2 et 3. La pente (a) de la droite de sécheresse permettant de définir la LSTmax selon le NDVI montre une variabilité significative au cours de l'été. La pente au mois de juin atteint une valeur moyenne de -19,4 avec un écart-type de de 7,3. L'ordonnée à l'origine (b) y atteint une température maximale moyenne de 35,6°C. Dans les mois subséquents, la pente diminue graduellement atteignant en moyenne -15,4, -14,3 et -10,6 pour les mois de juillet, août et septembre. Les températures maximales moyennes diminuent aussi progressivement avec des valeurs de 34,4°C, 30,0°C et 22,2°C. Les températures minimales fluctuent en moyenne entre 9,6 et 14,2°C. Le mois de juillet est celui qui montre les températures minimales les plus élevées.

Figure 4 . 5

 45 Figure 4.5 Répartition du nombre d'images satellites Landsat 8 superposées répondant aux critères de sélection.

Figure 4 . 5

 45 Figure 4.5 Répartition du nombre d'images satellites Landsat 8 superposées répondant aux critères de sélection.

  Figure 4.6a représente la distribution du TVDI moyen calculé sur l'ensemble du territoire d'étude. Les valeurs se situent sur une échelle allant de 0 à 1 avec une valeur moyenne de 0,36 et un écart-type de 0,12. Les valeurs les plus élevées sont essentiellement observées dans les basses terres, au pourtour de le Baie-James et sur certains affleurements rocheux, tandis que les valeurs plus faibles se concentrent généralement autour des plans d'eau (grands lacs et littoral) et dans les basses terres du plateau de l'Abitibi situées au sud du territoire d'étude. Les écart-types varient de 0 à 0,45 avec une moyenne et un écart-type de 0,10 et 0,03, respectivement. La Figure 4.6b montre la répartition spatiale des déciles de l'écart-type. L'écart-type est plus élevé dans les basses terres sur le pourtour de la Baie-James et sur certains affleurements rocheux. Les écart-types plus faibles montrent une répartition spatiale moins systématique. Une quantité moindre de données semble favoriser un écart-type plus faible comme en témoigne les bandes mauves d'orientation nord-sud sur la Figure 4.6b. 4.5.2 Le POEDES Les valeurs calculées pour le POEDES se situent entre 0 et 1. Sur l'ensemble du territoire, le POEDES moyen est de 0,54, avec un écart-type de 0,18. Le patron de dispersion présenté à la Figure 4.7 montre des POEDES plus faibles dans les basses terres de la Baie-James ainsi que sur certains hauts topographiques correspondant à des affleurements rocheux. Les dépressions, et plus spécifiquement celles sillonnées par le réseau hydrographique de surface montrent des valeurs de POEDES plus élevées.

Figure 4 .

 4 Figure 4.6 a) Distribution spatiale du Temperature-vegetation Dryness index (TVDI) moyen pour l'ensemble du territoire d'étude. b) Distribution spatiale du décile de l'écart-type des Temperature-vegetation Dryness index calculés sur chaque pixel du territoire d'étude.

Figure 4 .

 4 Figure 4.6 a) Distribution spatiale du TVDI moyen pour l'ensemble du territoire d'étude. b) Distribution spatiale du décile de l'écart-type des TVDI calculés sur chaque pixel du territoire d'étude.

Figure 4 . 7

 47 Figure 4.7 Distribution spatiale de l'indice du potentiel d'occurrence d'écosystème dépendant de l'eau souterraine (POEDES) sur le territoire d'étude.

  Tableau 4.3 Sommaire des variables hydro-géo-topographiques des bassins versants (Variables détaillées au Tableau 4.1).

Figure 4 . 8

 48 Figure 4.8 Représentation graphique des deux premières composantes principales (CP) réalisées sur l'ensemble des variables intégrées dans l'analyse statistique détaillées dans le Tableau 4.1. Les vecteurs rouges réfèrent aux variables jugées les moins bien projetés sur le plan ou redondants.

4. 6

 6 Figure 4.10 présente le POEDES moyen calculé selon le type de végétation identifié par les données de la cartographie écoforestière réalisée sur le territoire (MRNF 2009; 2012). L'indice moyen s'avère « élevé » (>0,75) pour les bétulaies et certaines sapinières. En contrepartie, l'indice est généralement faible (<0,5) pour les milieux non forestiers tels que les tourbières et certaines landes ainsi que pour certaines pessières noires. La quantification de l'apport d'eau souterraine pour contribuer aux besoins écologiques des milieux où le POEDES dépasse la portée de ce travail de thèse. Une certaine réserve doit donc être respectée dans l'interprétation de la dépendance de ces milieux à la présence d'eau souterraine. Des travaux supplémentaires seraient nécessaires pour mieux comprendre la relation qui existe entre le POEDES et la végétation. Néanmoins, une attention particulière devrait être portée à ces couverts forestiers dans une perspective où une dégradation ou une modification des conditions hydrologiques changerait les composantes du bilan hydrologique des bassins versants. De plus, tel que précisé sur la Figure 4.10, l'efficacité du POEDES semble être illustrée par le niveau trophique des tourbières. Sur le territoire, les tourbières ombrotrophes,

Figure 4 .

 4 Figure 4.11 Typologie des bassins versants et positionnement des sous-bassins versants. RR : Rugueux-Ruissellement. GL : Granulaires libre. BTMH : Basses-Terres -Milieux humides. PDES : potentiellement dépendant de l'eau souterraine. NDES : non dépendant de l'eau souterraine.

Figure 4 . 1

 41 Figure4.12 Représentation cartographique de la projection des sous-bassins versants sur l'axe de la composante principale 2 comparant leur potentiel dépendance à la ressource en eau souterraine (PDES).

  communauté scientifique mondiale. Cependant, peu de travaux avaient été menés en milieu boréal et au sein du Bouclier canadien, où le manque de données hydrogéologiques représente une entrave majeure à l'évaluation de la ressource en eau souterraine. Dans cette thèse, le suivi temporel réalisé sur 290 images thermiques du satellite Landsat 8, a permis de développer un indice de potentiel d'occurrence d'écosystèmes dépendants de l'eau souterraine (POEDES). Cet indice reflète la stabilité du niveau d'humidité du sol au cours de la période estivale, ce paramètre étant déduit de la température du sol. La démarche s'appuie sur l'hypothèse voulant que l'inertie thermique de l'eau contribue à maintenir la température des sols humides plus basse lors de la période estivale, en comparaison aux sols secs. Ainsi, les secteurs alimentés par l'eau souterraine montrent des températures plus basses et plus stables tout au long de la période estivale en comparaison aux secteurs où la nappe d'eau souterraine est profonde. L'indice ainsi développé permet de localiser, à grande échelle, les secteurs de résurgence d'eau souterraine ou de nappe affleurante. Ces connaissances sont essentielles pour localiser d'éventuels écosystèmes dépendants de l'eau souterraine et, ultimement, tenir compte de la biodiversité dépendante de l'eau souterraine dans les politiques de gestion du territoire. 5.1.2 Résultats Les résultats ont montré une épaisseur moyenne de sédiments augmentant progressivement depuis l'affleurement du socle rocheux pour atteindre une épaisseur de 26 m au-delà de 1150 m de distance. L'analyse géométrique a, quant à elle, permis d'estimer le volume d'eau souterraine s'écoulant dans des conditions de nappe libre à 40 km 3 ± 4 km 3 réparti sur 7 % du territoire. Limités spatialement aux formations fluvioglaciaires, les aquifères granulaires à nappe libre auraient un volume d'eau moyen estimé à 3,1 m/L 2 . Cette donnée complète l'information disponible selon laquelle les aquifères du Bouclier canadien sont restreints à des systèmes d'écoulements locaux et peu profonds (Figure 1.1). Les travaux réalisés ont permis de documenter les fluctuations saisonnières dans la composition isotopique de l'eau de six grandes rivières sur 21 sites distincts. Les données acquises révèlent par leur hétérogénéité, une diversité de contextes hydrologiques au sein des bassins versants drainés par les rivières. La contribution calculée de l'eau souterraine au débit annuel des cours d'eau a montré des proportions variables selon les rivières, avec des valeurs correspondant à <12 % des précipitations annuelles. Ces flux représentent un volume d'eau de 2,3 km 3 et montrent une faible saisonnalité. Ils peuvent toutefois représenter une proportion significative du débit des rivières, notamment à la fin de la période du couvert de glace (≈60 % calculé sur la rivière Harricana à Amos [HA]). Les eaux superficielles des bassins versants localisés aux plus hautes latitudes de la région ont des compositions isotopiques et des conductivités électriques se rapprochant davantage de celles des précipitations comparativement aux eaux superficielles des bassins versants situés plus au sud où la contribution d'eau souterraine aux débits des rivières est plus significative.

1 .

 1 ) la fin de la période du couvert de glace pour évaluer la proportion maximale atteinte par la contribution de l'eau souterraine au débit des rivières, (2) la crue printanière pour quantifier la source et les flux d'eau lors de crue maximale annuelle et (3) la fin de la période libre de glace pour estimer l'évaporation, l'apport des pluies à la crue printanière et le flux annuel d'eau souterraine vers les rivières. Procéder à l'échantillonnage des rivières lors de ces périodes permettra de limiter les contraintes logistiques et les coûts de futurs programmes de surveillance tout en assurant un suivi relativement complet des processus hydrologiques. Les valeurs de POEDES calculées sur les bassins versants ont montré une corrélation avec les flux d'eau souterraine contribuant au réseau hydrographique et aux secteurs montrant de fortes probabilités d'aquifère potentiel. La corrélation existant entre ces trois approches calculées de façon indépendante tend à indiquer que le POEDES fournit des résultats cohérents et qu'il pourrait être employé comme un outil supplémentaire pour localiser des écosystèmes dépendant de l'eau souterraine (EDES). Il pourrait également être un outil exploitable à l'échelle régionale afin de cibler les secteurs prioritaires pour l'acquisition de données terrain sur les écosystèmes dépendant de l'eau souterraine. 5.1.3 Perspectives Les changements climatiques anticipés vont occasionner des modifications dans l'ampleur et la saisonnalité des flux hydrologiques. L'augmentation de la température va également modifier la composition isotopique des eaux naturelles. Il est donc d'intérêt de poursuivre à long terme le suivi de la composition isotopique des précipitations (pluie et neige) et des cours d'eau. Une telle surveillance permettra de suivre l'évolution des flux d'eau souterraine alimentant les rivières des grands bassins versants de la Baie-James. Ceci sera indispensable dans la mesure où ces flux d'eau souterraine soutiennent une biodiversité. Cependant, les travaux menés sur le terrain au cours de ce projet ont demandé la mise en place d'une logistique complexe et la collaboration de plusieurs personnes pour assurer un suivi rigoureux des paramètres étudiés. En ce sens, le développement d'un réseau de surveillance intégrant la collaboration des universités, du gouvernement et des communautés locales sera nécessaire pour assurer le partage de connaissances et pour maintenir un suivi à long terme. Le potentiel d'occurrence d'écosystèmes dépendant de l'eau souterraine (POEDES) devrait aussi être développé et suivi au cours des années ou au cours d'une même année. Ce suivi permettrait d'établir la quantité minimale d'images à analyser pour établir un POEDES représentatif. Cela s'avère important car les résultats ont montré que le nombre d'images traitées et superposées influence la stabilité de l'indice d'humidité TVDI nécessaire au calcul du POEDES. De plus, le suivi intra-annuel permettra d'observer l'existence d'une saisonnalité dans la dépendance des milieux aux eaux souterraines. Des travaux ciblés le long de résurgences ponctuelles et diffuses d'eau souterraine seront aussi requis afin de mieux connaitre le niveau d'influence de la température des résurgences d'eau souterraine sur les températures de la surface du sol calculées depuis des capteurs thermiques satellitaires. De plus, l'utilisation d'images en infrarouge thermique à haute résolution obtenues par imagerie aérienne serait complémentaire à l'utilisation d'images satellitaires. Le développement d'un algorithme permettant d'automatiser le traitement d'images faciliterait aussi le déploiement du POEDES sur d'autres vastes territoires et sur une plus longue période. 5.2 Volet sociétal Au sein du territoire d'étude, l'écoulement naturel des eaux de surface a été significativement modifié par l'implantation de grands complexes hydroélectriques ayant eu pour effet de réduire l'accès à l'eau potable pour certaines communautés locales. L'identification de secteurs stratégiques caractérisés par un potentiel aquifère élevé et où des mesures de protection sont jugées nécessaires fournira les connaissances requises afin de déployer des mesures visant à assurer un accès à l'eau pour les générations futures. Le territoire comprend une quantité non négligeable d'aquifères granulaires à nappe libre et ces derniers sont inégalement dispersés et généralement peu accessibles. Des précautions particulières devraient être prises afin de protéger ces aquifères, notamment ceux identifiés prioritaires, advenant le développement de projets pouvant impacter la ressource en eau. Une attention particulière devra aussi être portée sur les secteurs à proximité des communautés locales et des campements. Les aquifères peu profonds à nappe libre pourraient constituer la meilleure source pour l'approvisionnement en eau potable. A titre d'exemple, un réseau de puits communautaires développé au sein des aquifères granulaires à nappe libre situés près des routes pourrait offrir un accès sécuritaire à l'eau potable et par le fait même accroître la qualité de vie de la population. Les infrastructures routières et les sites d'extraction de ressources minérales devraient davantage intégrer la disponibilité de l'eau souterraine dans leur planification pour limiter leurs impacts sur la ressource. Dans une perspective de gestion durable du territoire boréal, des mesures devront aussi être mises en place pour protéger la ressource en eau souterraine et ainsi contribuer au maintien des flux d'eau souterraine vers les rivières et ce, spécialement dans les secteurs où ils sont les plus élevés. Une détérioration des conditions d'écoulement de l'eau souterraine pourrait y avoir des conséquences nuisibles sur les écosystèmes des rivières et affecter les populations qui dépendent des ressources naturelles. Les travaux réalisés dans cette thèse ont aussi permis de hiérarchiser différentes zones du territoire selon le degré d'importance de l'eau souterraine dans les bassins versants. La cartographie, intégrant diverses connaissances sur l'eau souterraine : volume des réserves, flux d'échange dans les rivières et présence potentielle d'écosystèmes dépendants de l'eau souterraine, fournit aux gestionnaires du territoire un outil adapté au contexte régional pour intégrer des modulations en matière d'aménagement durable du territoire. De plus, l'épaisseur de dépôts meubles estimée pour évaluer le potentiel aquifère représente une information nécessaire à l'évaluation des coûts et contraintes associés aux campagnes de forage (prospection minière) et au développement d'infrastructures (routes, barrages et bâtiments). Cette information constitue une connaissance complémentaire pour identifier différents scénarios d'exploitabilité cohérents avec la disponibilité de ressource en eau souterraine et les écosystèmes qui en dépendent. Le développement de projets pourra être mieux encadré pour assurer la préservation du milieu naturel si de telles informations sont prises en compte. 5.3 Reproductibilité et ouvertures Les méthodes et approches élaborées dans cette thèse ont été utilisées dans le contexte du milieu boréal du nord-ouest du Québec, au sein du Bouclier canadien. Le déploiement de l'analyse multicritère pour définir le potentiel aquifère sur l'ensemble du nord du Québec a révélé une reproductibilité satisfaisante de l'approche à grande échelle, et ce malgré une diminution de la résolution des données. Le tableur Excel sur le processus analytique hiérarchique (AHP), fourni en matériel électronique supplémentaire fourni dans le cadre du Chapitre II, pourra aider à établir les poids aux critères sur d'autres territoires. Les critères d'applicabilité des méthodes ne sont pas exclusifs au territoire québécois. Ils sont présents dans d'autres secteurs de l'hémisphère nord, tels que l'ensemble du Bouclier canadien ainsi que sur le Bouclier scandinave et du Groenland. À notre connaissance, ces régions nordiques sont rarement ciblées dans les recherches à portée mondiale concernant les impacts de l'activité humaine sur les eaux souterraines. Cela pourrait s'expliquer par le fait que les conditions démographiques et hydrologiques qui prévalent dans la plupart de ces zones ne sont pas alarmantes. Néanmoins, ces zones représentent une superficie significative de la surface de la Terre. L'amélioration des connaissances sur la ressource en eau souterraine dans ces régions est donc indispensable pour mieux évaluer et représenter la disponibilité des réserves d'eau douce à l'échelle planétaire. La demande accrue de ressources naturelles, notamment d'énergie (hydroélectricité) et de matières ligneuses et minérales, pourrait générer des impacts humains importants dans un proche avenir. L'acquisition de connaissances sur la répartition de la ressource en eaux souterraines sur ces territoires, telles que produites dans cette thèse, pourrait donc fournir des données indispensables pour assurer la protection des aquifères nordiques. Finalement, les travaux réalisés dans le cadre de cette thèse montrent bien la nécessité d'utiliser des approches d'analyses spatiales pour aborder le rôle de l'eau souterraine dans les bassins versants du milieu boréal. Des efforts supplémentaires seraient nécessaires pour confronter les différentes approches développées au cours de cette thèse à différentes échelles spatiales. Un secteur bien documenté et instrumenté tel que le site d'Amos sur la rivière Harricana (site HA) pourrait être idéal pour investiguer plus en détails la fiabilité des approches proposées. Néanmois, les connaissances complémentaires acquises par l'intégration de données multidisciplinaires sur ces vastes territoires peu accessibles participeront à améliorer la gestion durable de la ressource en eau souterraine et la protection de la biodiversité associée. À l'échelle de la planète, il est dans l'intérêt de chacun d'agir pour préserver cette biodiversité. Notre qualité de vie, celle de ceux qui nous entourent et surtout celle de nos enfants Precision of digital elevation models (DEM) available for entire study area Four different DEMs are available as open access and cover the entire study area (Table

  47 NEIGE AMOS H18-1 -142,31 1,15 -19,10 0,14 HA Neige-Amos-H18-2 4,40 48 NEIGE AMOS H18-2 -164,77 1,64 -22,57 0,28 HA Neige-Amos-H18-3 7,54 49 NEIGE AMOS H18-3 -149,28 0,61 -20,80 0,07 HBJ HBJ-H18-N

  L'HYDROGÉOLOGIE EN MILIEU BORÉAL QUÉBÉCOIS: BILAN HYDROLOGIQUE ET GESTION DU TERRITOIRE Spécialité : Sciences et Génie de l'Environnement. Mots clefs : Réserves et écoulements de l'eau souterraine, potentiel aquifère, bilans géochimiques, potentiel occurrence d'écosystème dépendant de l'eau souterraine, SIG.

  conceptuels, des approches d'analyses spatiales et des méthodes adaptées à l'étude des grands espaces. Elle utilise notamment les SIG pour gérer, représenter, analyser et

	importants étaient épars et principalement associés à des formations fluvio-glaciaires
	traiter des données géologiques, hydrogéologiques et géochimiques à l'ensemble du telles que les moraines Harricana et Sakami. Les calculs géométriques ont permis
	territoire. d'estimer que le volume total d'eau souterraine stockée dans les aquifères granulaires
	libres atteignait environ 40 km 3 . Utilisées conjointement, les deux approches ont révèlé
	Spécifiquement, cette thèse vise trois objectifs spécifiques : que les aquifères granulaires libres peu profonds qui nécessitaient une protection accrue
	1. Coupler des méthodes d'analyse spatiale multicritère et géométrique pour des eaux souterraines représentaient environ 5 % du territoire. Dans les zones
	délimiter de potentiels aquifères granulaires peu profonds; d'aquifères granulaires captifs, les aquifères potentiels étaient situés dans les vallées

2. Contraindre les paramètres du bilan hydrologique sur la base de bilans isotopiques et d'équations de mélange pour estimer la contribution de l'eau souterraine au réseau hydrographique de surface; 3. Développer un indice spatial basé sur les variations thermiques de la surface du sol localisant les potentiels écosystèmes dépendant de l'eau souterraine et déterminer si l'occurrence de ce type de milieu est corrélée aux conditions hydrologiques et morphologiques des bassins versants. 1.6 Organisation de la thèse La thèse est structurée de façon à répondre aux objectifs spécifiques au sein de chapitres distincts. Ainsi, les chapitres 2, 3 et 4 sont construits selon un format d'article scientifique avec une mise en contexte, une description du territoire d'étude sensiblement la même pour chaque article, les méthodes utilisées, une description détaillée des résultats, une discussion et une conclusion. Une conclusion générale (Chapitre 5) synthétisant les principales retombées de l'étude et les perspectives complète la thèse. CHAPITRE II APPROCHES D'ANALYSE SPATIALE POUR L'ÉVALUATION ET LA PROTECTION DE LA RESSOURCE EN EAU SOUTERRAINE DANS LES GRANDS BASSINS VERSANTS DU BOUCLIER CANADIEN Résumé Cette étude s'est concentré sur le développement de deux approches d'analyse spatiale utilisées conjointement pour évaluer les ressources en eau souterraine associées aux aquifères granulaires du Bouclier canadien. La première approche était une analyse multicritères (MCA) utilisant un processus hiérarchique analytique (AHP) basée sur des données géologiques et hydrogéologiques pour classer la probabilité de trouver des réserves en eau souterraine facilement disponibles dans un territoire spécifique. La seconde approche reposait sur des calculs géométriques réalisés au sein d'un SIG pour évaluer l'étendue et le volume des aquifères. Les approches ont été appliquées sur un maillage de 100 × 100 m dans une zone de 185 000 km 2 (≈ 18,5 millions de mailles) correspondant aux bassins versants de la région de la Baie-James au nord-ouest du Québec. L'approche MCA-AHP a révélé que les aquifères avec les potentiels les plus fluviales et les basses terres. Une analyse de sensibilité menée sur l'approche MCA-AHP a révélé que la taille du maillage de la grille n'affectait pas significativement les résultats. Par conséquent, l'approche a été étendue vers le nord à un territoire de 490 000 km 2 (selon un maillage de 500 × 500 m) atteignant la région de la baie d'Ungava. La méthode proposée pourrait être adaptée et appliquée à d'autres zones de bouclier.

Table 2 .

 2 1 Classification of main surface units based on their dominant texture.

		Classification
	Surface units	(surface unit texture)
	Till cover (> 1 m), till veneer (< 1 m)	Glacial till
	Cochrane till, fine-grained deep-water sediments (marine and lacustrine)	Clay
	Glaciolacustrine and alluvial deposits	Silty sand
	Proglacial sediments, proglacial outwash, shallow glaciolacustrine and shallow marine sediments	Sand
	Frontal and interlobate moraine, ice-contact deposits	Gravel
	Organic	Peat
	Bedrock outcrop	Bedrock

Table 2

 2

	.2 Standardization of the thematic layer included in the multi-criteria approach
	on a scale of 1 to 5.
	Layers	Dimensions MCA	Classes of	Standardization
			Symbol	values	(dimensionless)
	Potential	[ ]	T'	15-25
	transmissivity			9-14
				5-8
				3-4
				1-2
	Soil Surface	°	S	0-2
	slope			2-6
				6-12
				12-18
				> 18
	Topographic	[ ]	TWI	> 11,5
	wetness index			9,5-11,5
				8,5-9,5
				7,5-8,5
				< 7,5
	Distance to a	[L]	DMBF	< 200
	major bedrock			200-400
	fault			400-600
				600-1000
				> 1000
	Distance to a	[L]	DFHN	0-500
	feature of the			500-1000
	hydrographic			1000-1500
	network			1500-2000
				> 2000
	Drainage	[L/L 2 ]	DND	0-0,5
	network density			0,5-1,0
				1,0-1,5
				1,5-2,0
				> 2,0

table where

 where

		2.4.3.4 Distance to a major bedrock fault (DMBF)
		Bedrock faults could cause linear depressions in the bedrock surface due to a lower
		resistance to erosion. These depressions are filled with greater thicknesses of
		unconsolidated sediments, thus locally increase the aquifer volume. While the
		Precambrian bedrock is likely to present a generally low hydraulic conductivity
				surface-and groundwater interactions are
	likely	to	occur	(Ala-aho et al. 2013;	Barthel and Banzhaf 2016,
	Larocque and Broda 2016; Arnoux et al. 2017). Here, ratings are evaluated based on
	the distance from the closest lake or pond, up to a distance of 2000 m. The density of
	the surface drainage network is also used as an indication of the ability of surface soils
	to route precipitation into surface runoff. Well-developed and denser drainage systems
	are interpreted as an indication of surface runoff dominance over infiltration and
	groundwater recharge (Mardhel et al. 2021). Drainage density was calculated with the
	cumulative length of stream segments within a 1000 m radius of each cell over the
	entire study area.

Table 2

 2

	.4 Degrees of importance attributed to criteria during pairwise comparisons
	(modified from Saaty 1980).
	Less important Extreme Very high High Average	Equal	More important Average High Very high Extreme
	1/9	1/7	1/5	1/3	1	3	5	7	9

Table 2

 2

		.5 Pairwise comparisons realized by authors and weights of criteria included in
	the multicriteria analysis (see spreadsheet provided in the electronic supplementary
	material for details).
		DND	S	DMBF	TWI	DFHN	T'	Weight (w)
	DND 1	3	5	1/3	1/5	1/5	0.100
	S	1/3	1	1	1/3	1/5	1/5	0.048
	DMBF 1/5	1	1	1/5	1/7	1/7	0.038
	TWI 3	3	5	1	1/3	1/5	0.139
	DFHN 5	5	7	3	1	1/3	0.260
	T'	5	5	7	5	3	1	0.416
	Σ	14.53	18.00	26.00	9.87	4.88	2.08	1
								CR= 0.087

Table 2 .

 2 6 Sensitivity analysis of the multicriteria analysis.

	Scenario	ID color	Selected	Criteria	Number of	classes	Limits of	classes 1	Weights	Mesh
	Ref		6		5		User defined	User defined	100 m
					2		Jenks
					3 5		Jenks User defined	Equal values
			6		5		Jenks
					2		Jenks
					3		Jenks		User defined
					5		Jenks			100 m
			5 (without DND)
			5 (without S)
			5 (without DMBF) 5 (without TWI)	5		User defined	User defined
			5 (without DFHN)
			5 (without T')
					2		Jenks
					3 5		Jenks User defined	Equal values
			6		5 2		Jenks Jenks			500 m
					3 5		Jenks		User defined

  𝑅𝛿 𝑅 + 𝑆𝛿 𝑆 = 𝑄 𝑅𝐼𝑉 𝛿 𝑅𝐼𝑉 + 𝑞 𝐷𝐺𝑊 𝛿 𝐷𝐺𝑊 + 𝑞 𝐸 𝛿 𝐸 + 𝑞 𝑇 𝛿 𝑇 ± ∆𝑆𝛿 ∆𝑆 (5) 𝐼𝜎 𝐼 = 𝑅𝜎 𝑅 + 𝑆𝜎 𝑆 = 𝑄 𝑅𝐼𝑉 𝜎 𝑅𝐼𝑉 + 𝑞 𝐷𝐺𝑊 𝜎 𝐷𝐺𝑊 + 𝑞 𝐸 𝜎 𝐸 + 𝑞 𝑇 𝜎 𝑇 ± ∆𝑆𝜎 𝛿 𝑅 is enriched in heavy isotope with respect to snow 𝛿 𝑆 . The electrical conductivity of atmospheric inflows (𝜎 𝐼 ), including rain (𝜎 𝑅 ) and snow (𝜎 𝑠 ), are generally low due to low total dissolved solids in precipitation. The evaporating vapour (𝛿 𝐸 ) will be depleted in heavy isotopes as the residual liquid water (surface water)

∆𝑆 (6)

where 𝛿 is the isotopic composition of water (either 𝛿 2 𝐻 or 𝛿 18 𝑂, expressed in ‰ against Vienna Standard Mean Ocean Water

[VSMOW]

) and 𝜎 is their electrical conductivity of water (in µ𝑆/𝑐𝑚), with subscripts corresponding to the terms defined in Equation 1 (Figure

3

.2a). The stable isotope compositions of atmospheric inflows (𝛿 𝐼 ), including rain (𝛿 𝑅 ) and snow (𝛿 𝑠 ), vary according to seasonality, latitude, and altitude. They can be interpolated over vast regions based on quantitative models

(Bowen and Ravenaugh 2003;[START_REF] Terzer | Global isoscapes for δ 18 O and δ 2 H in precipitation: improved prediction using regionalized climatic regression models[END_REF]

. When plotted in a 𝛿 2 𝐻 -𝛿 18 𝑂 graph (basal plane in Figure

3

.3), the stable isotope composition of precipitation from mid-to high latitude regions generally define local meteoric water lines (LMWL) where

  1 -1 𝛿 18 𝑂 𝑠 𝛿 18 𝑂 𝑅 𝛿 18 𝑂 𝑆𝐺𝑊 -𝛿 18 𝑂 𝐸 𝛿 2 𝐻 𝑠 𝛿 2 𝐻 𝑅 𝛿 2 𝐻 𝑆𝐺𝑊 -𝛿 2 𝐻 𝐸The isotopic composition (𝛿 2 𝐻 -𝛿 18 𝑂) and electrical conductivity (𝜎) of water are used as conservative tracers in Equation12. The values that must be attributed to endmembers can represent different periods of time, depending on the water residence time in the catchment. As illustrated in Figure3.2b, precipitation falling between November and March accumulates in the snowpack, to be further released during the snowmelt period in late April. A single flux-weighted average value can be used for 𝛿 𝑆 because the spring freshet generally consists of one major event in the hydrograph of large boreal rivers. Here, the isotope fractionation occurring during the melting process

	𝛿 𝛽 =	𝐼𝛿 𝐼 -𝑄 𝑅𝐼𝑉 𝛿 𝑅𝐼𝑉 -𝑞 𝐸 𝛿 𝐸 𝛽	(13)
	with:
	𝛽𝛿 𝛽 = 𝑞 𝑇 𝛿 𝑇 + 𝜌𝛿 𝜌	(14)
	3.3.4 The stable isotope composition of transpired moisture
	The water budget conceptualized in Figure 3.2 and represented by Equation 2 includes
	a 𝛽 term that corresponds to the proportion of water not reaching the hydrographic
	𝜎 𝑆 network. This portion of water corresponds to transpiration (𝑞 𝑇 ), changes in water 𝜎 𝑅 𝜎 𝑆𝐺𝑊 0 ) -1 × ( 1 𝛿 18 𝑂 𝑅𝐼𝑉 𝛿 2 𝐻 𝑅𝐼𝑉 𝜎 𝑅𝐼𝑉 ) = ( 𝑞 𝑆 𝑞 𝑅 𝑞 𝑆𝐺𝑊 𝑞 𝐸 storage (±∆𝑆), and deep groundwater flow (𝑞 𝐷𝐺𝑊 ; Equation 3). Equation 5 can be ) (12) reorganized for isolating 𝛿 𝛽 :

[START_REF] Taylor | Isotopic evolution of a seasonal snowpack and its melt[END_REF]

) is assumed to impart a negligible effect on the isotopic compositions of large rivers during snowmelt. The 𝛿 𝐸 term can be estimated based on Equation 9, while 𝛿 𝑆𝐺𝑊 are assumed to be stable over time. The 𝛿 𝑅 term can be more variable and cause uncertainties in calculations due to the unknown residence time of water in the watersheds. Here, a single, flux-weighted average 𝛿 𝑅 value is used to solve Equation

12

. At the scale of the studied watersheds, it is assumed that the distances and periods required for the transfer of rainwater to large rivers is sufficient to buffer temporal changes in 𝛿 𝑅 . The terms 𝜎 𝑆 and 𝜎 𝑅 can be assumed to be close to zero, while the σ 𝑆𝐺𝑊 term, if not known a priori, must be estimated based on assumptions related to the hydrogeological context of watersheds.

Under the conditions prevailing in the study area, it can be assumed that 𝛿 𝜌 ≈ 𝛿 SGW ≈ 𝛿 𝐼 and the stable isotope composition of the transpired moisture [𝛿 𝑇 ] can be estimated:

Table 3 .

 3 1 Remote sensing datasets used in this study.

	Water type	Dataset	Spatial resolution	Temporal resolution	Reference
	Precipitation	IMERG CaPA	0,1° (7 × 11 km) Monthly Huffman et al. (2009) 10 km Daily Fortin et al. (2015)
	ET	MODIS/terra 1 km SSEBop 1 km	8 days Monthly Senay et al. (2013) Running et al. (2019)
	3.4.1.2 Discharge measurements
	In situ river discharges were measured using a Sontek RiverSurveyor M9 acoustic

Doppler current profiler (ADCP). Specific details of the method used to measure the

  Figure 3.1a) from May 2018 to October 2019. In addition, five composite samples were also collected to estimate the average 𝛿 2 𝐻 -𝛿 18 𝑂 and 𝜎 of the snowpack (Figure 3.1a). Samples were collected at the end of the ice-on period,

	Measurements and Data Sources 3.4.2.3 River ice sampling
	3.4.2.1 Precipitation
	River ice cores were collected at 11 outlets (Figure 3.1a) at the end of winter in 2018
	During this study, 61 monthly composite samples of precipitation were collected at and 2019 using a Kovacs Mark II corer. This corer allows for retrieving 9 cm diameter
	ice cores in 1 m long sections. Specific details about the ice coring method are provided in Appendix B-4. three monitoring stations (3.4.2.2 River water sampling 3.4.2.4 Stable isotope analyses (δ 2 H-δ 18 O)
	River waters were collected downstream of rapids or near the center of the main
	channel of rivers to obtain well-mixed samples. Specific details of this sampling
	method are provided in Appendix B-3. Overall, 222 samples were collected within 21
	sub-watersheds (Figure 3.1a). Sampling frequency depended on the sites' accessibility
	and is mainly synchronized with discharge measurements. The Harricana River
	(sampled in the city of Amos; site HA in Figure 3.1a) is used as the main reference
	station because it was monitored at a greater frequency (daily monitoring during spring

just before snowmelt, to obtain composite samples of the snow accumulated during the entire winter. Specific details about the sampling method are provided in Appendix B-2. The isotopic composition (𝛿 2 𝐻 -𝛿 18 𝑂) of precipitation (including rain and snow) can further be interpolated over vast regions based on quantitative models

(Bowen and Ravenaugh 2003)

. Here, in situ isotopic composition of precipitation is used for comparisons with the regionalized cluster-based water isotope prediction model (RCWIP,

[START_REF] Terzer | Global isoscapes for δ 18 O and δ 2 H in precipitation: improved prediction using regionalized climatic regression models[END_REF]

.

freshet for a total of 75 samples collected over 18 months). The number of samples collected at each site is provided in Appendix D. All samples were analyzed at the Geotop laboratory (Université du Québec in Montréal) with a LGR (Los Gatos Research) model T-LWIA-45-EP Off-Axis Integrated Cavity Output Spectroscopy (OA-ICOS) device. Specific details of stable isotope analyses are provided in Appendix A-5. Results are given in delta units (δ) in ‰ vs VSMOW. Overall, 491 samples containing waters from rivers, rain, snow, snowpack, and ice cores were analyzed for stable isotopes of the water molecule (δ 2 Hδ 18 O).

Table 3

 3 River discharge (𝑄), and groundwater balance (𝜌 = 𝐼 -𝐸𝑇 -𝑄)). SA : Entire Study Area. Italic 𝑄 refers to data from hydrometric stations.

	.2). The precipitation

watersheds and approximately 700 mm/yr in the northern watersheds. Late summer and early autumn periods are characterized by the highest humidity, with precipitation reaching an average of 150 mm/month in September and October 2018. During the iceon period, conditions are dryer, with average precipitation (mainly snow) < 50 mm/month. The precipitation accumulated in the snowpack at the end of the ice-on period represented a water equivalent of 29% and 32% of total inflow in 2018 and 2019, respectively. The month of April marks a transition, with both liquid and solid precipitation and measurable ET in some areas. The ET decrease with increasing latitudes (≈ 35 and 20 mm/°φ in 2018 and 2019, respectively). From north to south, the ET reached 230 to 390 mm/yr and 275 to 370 mm/yr in 2018 and 2019, respectively, with an average of 305 ± 17% mm/yr and 321 ± 28% mm/yr for the entire study area.

In spring and early summer, the ET flux was equal to rainfall in 2018 and 2019.

Maximum ET rates are observed in late spring and early summer, reaching approximately 100 mm/month. The ET fluxes decrease from August to November. Table

3

.2 Summary of hydrological result (Total inflow (𝐼), Rain (𝑅), Snow (𝑆), Evapotranspiration (𝐸𝑇), .2 (suite) Summary of hydrological result (Total inflow (𝐼), Rain (𝑅), Snow (𝑆), Evapotranspiration (𝐸𝑇), River discharge (𝑄), and groundwater balance (𝜌 = 𝐼 -𝐸𝑇 -𝑄)).

Table 3 .

 3 3 Statistical summary of geochemical results for each water type.

			River			Precipitation		Snowpack			Ice core
		δ	2 H	δ	18	O	σ	δ	2 H	δ	18 O	σ	δ	2 H	δ	18 O	σ	δ	2 H	δ	18 O
		‰	‰	μs/cm	‰	‰	μs/cm	‰	‰ μs/cm	‰	‰
	Min	-128.27 -17.70 0.50 -200.4 -27.07 3.50 -196.64 -25.04 4.35 -213.14 -29.03
	Max	-75.67 -9.90 159.80 -50.3 -6.85 126.00 -142.31 -19.10 13.54 -66.40 -8.57
	Mean -100.1 -13.5 47.46 -97.0 -13.3 18.17 -168.38 -21.92 6.77 -91.34 -12.13
	Sample		222				211		61			54		13			13		192

Table 3 .

 3 4 Description of values assigned to electrical conductivity (σ) of groundwater. UGA = Unconfined granular aquifers. UBA = unconfined bedrock aquifers. CA = Confined aquifers.

	Scenario	Origin aquifer of	Assignated 𝜎 𝑆𝐺𝑊
		groundwater
	Reference	Mainly from unconfined	Weighted average of the approximated
	value	aquifers	𝜎 of the UGA and UBA
	1	A composite of main aquifers	Weighted average of the approximated 𝜎 of the UGA, UBA and CA
	2	Mainly from unconfined granular aquifers	Approximated 𝜎 of the UGA
	3	Mainly from unconfined bedrock aquifers	Approximated 𝜎 of the UBA

The estimated values of reference for 𝜎 𝑆𝐺𝑊 range between 289 and 393 μs/cm (Table

3

.5). These endmembers are shown in Figure

3

.7b to allow for a comparison with river waters. The results obtained from Equation

12

are shown in Table

3

.6. The results are expressed as proportions of the annual 𝐼 for each sub-basin. The results

Table 3 .

 3 5 End-members values of watersheds. 𝑅: Rain, 𝑆: Snowmelt, 𝑆𝐺𝑊: Shallow groundwater, 𝐸: Evaporation, S1-S5: scenarios, UGA: Unconfined granular aquifer, UBA: Unconfined bedrock aquifer, CA: Confined Aquifer. The provided isotopic compositions correspond to the weighted average composition of precipitation from Regionalized cluster-based water isotope prediction for months of April to October (Rain), for months of November to March (Snowmelt) and for all months (Groundwater). The electrical conductivity of water (𝜎) is assumed to be invariant over the year. The 𝜎 values used as reference (RV) and for additional scenarios (S1-S3) are based on aquifers type proportion.

	Catchments	𝑅	𝛿 2 𝐻 (‰) 𝑆 𝑆𝐺𝑊	𝑬	𝑅	𝛿 18 𝑂 (‰) 𝑆 𝑆𝐺𝑊 𝑬	𝑅 𝑆	𝜎 (μs/cm) 𝑆𝐺𝑊 RV S1 S2 S3	Aquifer (%) UGA UBA CA
	BL 1	-76.6 -144.1 -94.1 -235.5 -10.9 -19.6 -13.2 -38.4			316 351 150 400 27 53 19
	BM 1 -77.4 -144.7 -95.0 -250.0 -11.0 -19.7 -13.3 -40.8			317 399 150 400 18 37 45
	BS 1	-75.5 -142.5 -92.6 -222.0 -10.8 -19.4 -13.0 -36.1			289 367 150 400 28 35 37
	CC 1 -83.4 -158.6 -103.9 -223.1 -11.8 -21.4 -14.4 -35.2			355 361 150 400 17 79 4
	CW 1 -83.3 -157.7 -104.8 -257.0 -11.8 -21.3 -14.5 -40.7			332 381 150 400 19 52 29
	HA 1 -75.4 -141.0 -92.9 -245.9 -10.7 -19.2 -13.0 -39.9			309 429 150 400 14 24 63
	HJ 1 HV 1 -75.2 -140.8 -92.7 -244.8 -10.7 -19.2 -13.0 -39.9 -77.9 -144.5 -95.4 -255.4 -11.1 -19.6 -13.3 -41.9 NBJ 1 -80.0 -148.9 -99.1 -11.3 -20.2 -13.8 WM 1 -81.1 -150.9 -101.0 -260.4 -11.5 -20.4 -14.0 -40.6	19.50 ± 13.60	6.77± 1.73	324 464 150 400 6 294 421 150 400 16 22 62 14 80 325 410 150 400 15 36 48 348 398 150 400 14 53 33
	BRBJ 2 -83.6 -155.9 -107.1 -280.4 -11.7 -21.0 -14.7 -44.3			377 436 150 400 5	47 48
	BRN 2 -85.3 -161.2 -108.4		-12.0 -21.7 -15.0				13 74 13
	Cold 2 -89.7 -158.6 -117.7		-12.4 -21.3 -16.0				1	40 60
	ERBJ 2 -88.8 -164.7 -114.2 -276.8 -12.3 -22.0 -15.6 -43.6			381 454 150 400 3	36 61
	E13k 2 -88.3 -156.1 -108.4		-12.3 -21.1 -14.9				3	34 63

Table 3 .

 3 5 (suite) End-members values of watersheds. 𝑅: Rain, 𝑆: Snowmelt, 𝑆𝐺𝑊: Shallow groundwater, 𝐸: Evaporation, S1-S5: scenarios, UGA: Unconfined granular aquifer, UBA: Unconfined bedrock aquifer, CA: Confined Aquifer. The provided isotopic compositions correspond to the weighted average composition of precipitation from Regionalized cluster-based water isotope prediction for months of April to October (Rain), for months of November to March (Snowmelt) and for all months (Groundwater). The electrical conductivity of water (𝜎) is assumed to be invariant over the year. The 𝜎 values used as reference (RV) and for additional scenarios (S1-S3) are based on aquifers type proportion.

	Catchments	𝑅	𝛿 2 𝐻 (‰) 𝑆 𝑆𝐺𝑊	𝑬	𝑅	𝛿 18 𝑂 (‰) 𝑆 𝑆𝐺𝑊 𝑬	𝑅 𝑆	𝜎 (μs/cm) 𝑆𝐺𝑊 RV S1 S2 S3	Aquifer (%) UGA UBA CA
	E26k 2 -89.9 -165.1 -115.7		-12.5 -22.1 -15.8				3	35 62
	HBJ 1 -81.1 -150.2 -99.8		-11.4 -20.3 -13.8				4	7 89
	OM 1 -78.1 -146.4 -95.8 -254.7 -11.1 -19.9 -13.4 -40.9			356 401 150 400 12 57 31
	Joli 2 -89.5 -164.2 -116.0		-12.4 -22.0 -15.8				0	27 73
	ORBJ 2 -90.4 -166.6 -115.9 -280.3 -12.5 -22.3 -15.8 -44.0 Fish 2 -90.5 -162.1 -116.3 -12.5 -21.7 -15.8 PRBJ 2 -87.6 -162.8 -112.6 -290.3 -12.2 -21.8 -15.4 -45.7 RRBJ 2 -86.7 -162.5 -111.2 -273.2 -12.1 -21.9 -15.3 -43.7	19.50 ± 13.60	6.77± 1.73	381 437 150 400 4 393 470 150 400 1 364 433 150 400 7 357 365 150 400 16 78 6 49 47 27 73 43 50
	RRN 2 -88.0 -165.7 -111.8		-12.3 -22.3 -15.4			352 352 150 400 19 80 1
	TH 1 -82.1 -152.4 -101.3		-11.5 -20.6 -14.0			361 488 150 400 1	7 92
	HT 1 -78.8 -146.0 -96.5		-11.2 -19.8 -13.4			393 488 150 400 0	11 89

Table 3 .

 3 6 Quantitative estimation of terrestrial fluxes where flux-weighted estimation of river discharge (hydrological and geochemical values) are available. The results are expressed in proportion of total inflow (𝐼). 𝑞 𝑅 : Rain runoff, 𝑞 𝑆 : Snowmelt runoff, 𝑞 𝑆𝐺𝑊 : Shallow groundwater contribution, 𝑞 𝐸 : Evaporation, 𝑞 𝑇 : Transpiration (𝑞 𝑇 = ET-𝑞 𝐸 ), ET : Evapotranspiration, β : Water not reach hydrographic network (β = 𝐼-𝑞 𝑆 -𝑞 𝑅 -𝑞 𝑆𝐺𝑊 ), 𝜌 : Groundwater balance, E/I : Evaporation ratio from[START_REF] Craig | Deuterium and oxygen 18 variations in the ocean and marine atmosphere[END_REF]

	Catchments	mm/yr	I	± %	𝑞 𝑅 %	± %	%	𝑞 𝑆	± %	𝑞 𝑆𝐺𝑊 % ± %	𝑞 𝐸 %	± %	%	β	± %	𝑞 𝑇 %	± %	δ 2 H	E/I	δ 18 O
	BS*	1085 5 50 2 25 3	2 2.5 8	1	23 5	26 6	9% 8%
	BL*	1044 10 61 1 30 2	1 6.4 7	1	8 9	28 6	5% 6%
	BM*	1017 9 49 2 28 2	2 1.8 4	2	21 9	32 6	4% 4%
	CC*	898 11 34 1 18 3	7 0.1 9	2	41 11 27 10 15% 13%
	CW*	897 14 48 2 34 4	9 0.4 8	2	10 13 27 7	5% 6%
	WM*	920 11 35 0 40 1	3 0.1 10 2	23 11 28 8	2% 5%
	HV*	1046 5 35 3 20 4 12 0.4 5	2	33 5	33 7	3% 3%
	HA*	1026 5 30 3 22 3 12 0.3 4	1	37 5	34 6	2% 3%
	HJ*	978 6 38 2 20 4	9 0.4 2	2	34 6	38 6	1% 2%
	OM**	951 10 38		27			1	5		34			34		2% 3%
	ERBJ***	658		39		22			1	3		37			37		4% 4%
	RRBJ***	873		47		23			1	3		29			29		4% 3%
	ORBJ***	654		37		22			2	3		39			39		4% 4%
	BRBJ***	848		42		35			1	5		22			22		1% 2%
	PRBJ***	735		45		34			0	4		20			20		1% 2%

Table 3 .

 3 7 Stable isotopic composition of transpired water compared to isotopic composition of inflow.

		𝑇		𝐼		Mean enrichment
		δ 2 H	δ 18 O	δ 2 H	δ 18 O	δ 2 H	δ 18 O
	HV*	-78.2 -11.10 -92.70 -13.00 -14.52	-1.90
	HA*	-77.5 -11.02 -92.90 -13.00 -15.36	-1.98
	HJ*	-86.9 -12.21 -95.40 -13.30 -8.49	-1.09
	BS*	-78.0 -11.09 -92.60 -13.00 -14.64	-1.91
	BL*	-78.3 -11.15 -94.10 -13.20 -15.82	-2.05
	BM*	-78.9 -11.21 -95.00 -13.30 -16.13	-2.09
	CC*	-92.2 -12.91 -103.90 -14.40 -11.73	-1.49
	CW*	-77.7 -11.06 -104.80 -14.50 -27.12	-3.44
	WM*	-55.1 -8.17 -101.00 -14.00 -45.86	-5.83
	OM**	-78.0 -11.10 -95.80 -13.40 -17.84	-2.30
	BRBJ*** -85.3 -11.95 -107.10 -14.70 -21.80	-2.75
	PRBJ*** -95.4 -13.20 -112.60 -15.40 -17.19	-2.20
	ERBJ*** -106.5 -14.60 -114.20 -15.60 -7.71	-1.00
	ORBJ*** -107.7 -14.74 -115.90 -15.80 -8.21	-1.06
	RRBJ*** -107.6 -14.81 -111.20 -15.30 -3.58	-0.49

*From 2018 average precipitation data with 𝜌 ≠ 0, **From 2018 average precipitation data with 𝜌 = 0, ***From 2019 average precipitation data with 𝜌 = 0 3.6.7 Recommendations pertaining to river monitoring

  1 représente le plus haut potentiel, 𝑇𝑉𝐷𝐼 ̅̅̅̅̅̅̅ est la valeur moyenne des TVDI calculés pour chaque pixel du territoire pour la période de référence (mois estivaux) et 𝐷𝜎 est le décile attribué à l'écart-type des TVDI au sein du jeu de données regroupant Tableau 4.1 Description des caractéristiques des sous-bassins versants considérées. T : Topographie. G : Géologie. H : Hydrologie. Hg : Hydrogéologie. IS : Indice spatiaux.

	Tableau 4.1 (suite) Description des caractéristiques des sous-bassins versants 4.4.3 Analyse en composantes principales
	considérées. T : Topographie. G : Géologie. H : Hydrologie. Hg : Hydrogéologie. IS :
	Indice spatiaux.
	Caractéristique Description MNA moy. Altitude moyenne du bassin versant Caractéristique Description	Catégorie T Catégorie	Chapitre de la thèse II Chapitre de la thèse
	minimale observées sur le bassin versant MNA d'aquifère granulaire. Pente moyenne (°) calculée à partir du définissant la présence potentiel Zone priorit. Pente_moy parcelles du bassin versant MNA écart Écart entre l'altitude maximale et Valeur de priorisation moyenne des	T T	Hg		II II	II
			Proportion du bassin versant où la nappe Estimation de la contribution de	G Hg		II	III
	Confiné GW	phréatique régionale est potentiellement confinée sous une couche d'argile l'eau souterraine (mm/an) au débit du cours d'eau principal sur le
			glaciolacustres ou glaciomarines. bassin versant
		T	Proportion du bassin versant recouvert Estimation de la transpiration sur le	G	H		II	III
	Gra		de dépôts meubles sans couverture bassin versant (mm/an)
		qr/R	argileuse Ratio de la contribution du			H		III
	Roc		Proportion du bassin versant où le socle ruissellement d'eau de pluie au rocheux affleure débit du cours d'eau principal par	G			II
	Lac		Proportion du bassin versant occupé par rapport aux précipitations un plan d'eau pluvieuses calculées sur le bassin	H			II
	MH	TVDI et	Proportion du bassin versant occupé par versant un milieu humide Indice TVDI moyen du bassin	H	IS		II	IV
			Proportion du bassin versant occupé par versant	H			II
	MH confiné TVDI moy. un milieu humide sus-jacent à une Indice TVDI moyen du bassin			IS		IV
			couche argileuse versant
		Proportion du bassin versant occupé par TVDI et dec Décile moyen des écart-types de	H	IS		II	IV
	MH Libre	un milieu humide en condition l'indice TVDI sur le bassin versant
		POEDES	d'aquifère à nappe libre Indice du potentiel d'occurrence			IS		IV
			Indice de saturation du sol moyen du d'écosystèmes dépendant de l'eau	H			II
	TWI moy.	bassin versant. (Topographic Wetness souterraine moyen du bassin
			Index) versant
	DND	Long.	Densité moyenne du réseau de drainage Longitude du centroïde du bassin sur le bassin versant versant.	H Variable supplémentaire	II
	Lat. DFHN	Distance moyenne par rapport au réseau Latitude du centroïde du bassin hydrographique versant.	H Variable supplémentaire	II
	DMBF	géologique majeure Distance moyenme d'une faille	G			II
	AP		bassin versant Indice moyen d'aquifère potentiel sur le	Hg		II
	ESS		sédiments saturés du bassin versant Estimation de l'épaisseur moyenne de	G			II

l'ensemble des écarts-types calculés sur chaque pixel du territoire. L'utilisation du décile permet de transformer et d'étaler les écarts-types calculés des indices TVDI uniformément sur une échelle de 0 à 1 selon une précision décimale, où 1 représente la plus grande variance.

4.4.2 Les variables hydro-géo-topographiques

La caractérisation hydro-géo-topographiques des sous-bassins versants présentés à la Figure 4.2b s'effectue sur la base des données compilées et produites au cours de cette thèse. Les caractéristiques considérées pour l'analyse statistique des bassins versants sont rapportées au Tableau 4.1. Leur source, leurs incertitudes et les traitements réalisés sont précisés dans les chapitres correspondants.

Table A1

 A1 Quebec's ministry of natural resources (MRN); 4 MERN DEM after transformation into a coarser grid. The errors are systematically calculated by subtracting the DEM elevation value from the geodetic marker's known elevation. The values are based on 114 DEM vs markers comparisons.

	Evaluation of the precision / accuracy of available DEMs
		STRM 1 NED 2 MERN 3 MERN-B 4
	Grid resolution (m)	90	50	10	100
	Average error (m)	3.8	2.7	1.6	2.6
	Average absolute error (m)	4.0	4.6	3.9	4.7
	Standard deviation (m)	4.2	5.4	5.0	5.5
	Absolute standard deviation (m) 3.9	3.8	3.4	3.9
	Greatest negative error (m)	-1.9	-9.4	-8.9	-8.8
	Greatest positive error (m)	23.6 19.4	20.7	20.7
	Sum of errors (m)	427.9 305.4 178.9	299.5
	Sum of absolute errors (m)	457.2 528.2 448.9	541.2
	1 NASA Shuttle Radar Topography Mission (SRTM); 2 USGS National Elevation
	Dataset (NED); 3

Table A2

 A2 Classification of surficial geologic units mapped in the used product, based on associated texture as presented in Table2.1 of the main article.

	3. The main steps of the MCA-AHP
	ANNEXE B
	Class	Class
	MCA 2 Class MCA DÉTAILS SPÉCIFIQUES DES MÉTHODES UTILISÉES DANS LE CADRE DU 1 1BA Ablation till (average thickness > 1 m) 1BC Rogen moraine (ribbed) 1BC Rogen moraine (ribbed) 1BD Drumlins 1BD 1BG De Geer moraine 1BG De Geer moraine 1BI Interlobate moraine 5 1BIM Interlobate moraine (thin) 1BI2A Interlobate moraine between 0.5 -1 m) 5AM Marine sediments (deep water) -Thin 5AR Deep-water marine deposits 3AN ancient 3AN Ancient fluvial deposits 1BIY Interlobate moraine (average thickness 5A Marine sediments (deep water) 5A Deep-water marine deposits 3AE recent 3AE Recent fluvial deposits glaciofluvial deposits 1 m 3 3AC Current with water) average thickness between 0.5 -3A Alluvial deposits 1BI Interlobate moraine 4GSY Glaciolacustrine sediments (shallow 3 Fluvial deposits 3 Fluvial deposits 1A Glacial deposits without particular morphology / Undifferentiated till 1A Undifferentiated till (average thickness > 1 m) 1AB Glacial block field 1AB Glacial block field (4 ha) 1AD Washed till 1AD Washed till 1AS Washed till of sedimentary rock R1 R1A 1B Undifferentiated glacial deposit (form) MCA Class MCA 1AAY Cochrane till (2 nd advance) average 4GAY Glaciolacustrine sediments (deep water) 7E Thick organic deposits (> 1m) 7BR Organic deposits: rippled bog CHAPITRE III thickness between 0.5 -1 m 1BF Frontal moraine average thickness between 0.5 -1 m 7R Very thin organic deposits 7BS Organic deposits: structured bog 1BF Frontal moraine 5 1BFY Frontal moraine (average thickness between 0.5 -1 m) 4GS Glaciolacustrine sediments (shallow 4GS Glaciolacustrine sediments; 7T Thin organic deposits 7BU Organic deposits: uniform bog 1BF1A Frontal moraine with till R1BF 1BF2A Frontal moraine with ice contact 1BF2B D deltaïc sediments water) thickness between 0.5 -1 m Frontal moraine with glaciofluvial 4 4GSM Thin glaciolacustrine sediments (shallow 7TY Thin organic deposits -average deposits water) shallow water (thickness < 0.25m) 4P Beach 4P Beach R4GS 7TM Thin Thin organic deposits 7FR Dépôts organiques: rippled fen 1. Discharge Measurement from Sontek RiverSurveyor M9 acoustic Doppler current 7FS Dépôts organiques: structured fen 7FU Dépôts organiques: uniform fen profiler (ADCP)
	1BN Corrugated moraine 1BP Loose till moraine 1P Loose till moraine 1BT Debris drag mounds 1AM Glacial deposits without particular morphology / Undifferentiated till (thin) 1AR Glacial deposits without particular morphology / Undifferentiated till (very thin) 1AY Glacial deposits without particular morphology / Undifferentiated till (average thickness between 0.5 -1 m) 1BDY Drumlins (average thickness between 0.5 -1 m) 1BPY Loose till moraine (average thickness between 0.5 -1 m) 1BR Very thin undifferentiated glacial deposit (form) 1BY Undifferentiated glacial deposit (form) average thickness between 0.5 -1 m 1M Thin 1Y Average thickness between 0.5 -1 m 1AA Thick Cochrane till (2 nd advance) R1AA 1AAM Thin Cochrane till (2 nd advance) 1AAR Very thin Cochrane till (2 nd advance) Glaciofluvial deposits 2 8A Weathering materials 2 8 8C Colluvium 1P Loose till 1AR Thin undifferentiated till 1BAR Thin ablation till thickness between 0.25 -1 m) (average 1AR_R Thin undifferentiated till (average thickness between 0.25 -1 m) with very frequent outcrops (<16 ha) 1AR_7 Thin undifferentiated till (average thickness between 0.25 -1 m) with very frequent organic deposits (<16 ha) 1BGR De Geer moraine (average thickness between 0.25 -1 m) 1PR Thin loose till moraine (average thickness between 0.25 -1 m) 1AA Clay till 1AAR Clay till (average thickness between 0.25 -1 m) (average thickness between 0.25 -Slope and weathering deposits Glaciofluvial deposits 4 1 m) 2B Proglacial deposits 5AY Marine sediments (deep water) average 5 8E Rocky scree 2BR Proglacial deposits (very thin) thickness between 0.5 -1 m 8P Skin sliding 2BP Deltaic esker 2BE Outwash R5A 8G landslide 2B Proglacial deposits 4 2BEY Outwash (average thickness between 0.5 5S Marine sediments (shallow water) 5S Shallow-water marine deposits 9 Eolian deposits 9 Eolian deposits 2BE Outwash 5SM Marine sediments (shallow water) -Thin 5SR Shallow-water marine deposits 4 9A Eolian deposits -active dune 9A Eolian deposits -active dune -1 m) (average thickness between 0.25 -9S Eolian deposits -stabilized dune 9S Eolian deposits -stabilized dune R2BE 2BD Glaciofluvial delta 2BD Glaciofluvial delta 1 m) R Rock R Rock (thin layer of deposits) 4 2A Ice contact deposits 5SY Marine sediments (shallow water) -5S_R Shallow-water marine deposits RS Sedimentary rock R_5A Rock (thin layer of deposits) with 2A Ice contact deposits 5 2AM Ice contact deposits (thin) average thickness between 0.5 -1 m with very frequent outcrops (<16 very frequent deep water deposits 2AR Ice contact deposits 2AY Ice contact deposits (average thickness ha) 1 (<16 ha) 2AE Esker (4 ha) R5S 6 Marine littoral deposits R_6 Rock (thin layer of deposits) with between 0.5 -1 m) 2AE Esker R6S 6_7 Marine littoral deposits with very very frequent littoral marine 2AK Kame, Kame terrace frequent organic deposits (<16 ha) deposits (<16 ha) 2AK Kame 6A Marine littoral deposits (current beach) 6R Thin marine littoral deposits R_7 Rock (thin layer of deposits) with 2AT Kame terrace 6S Marine littoral deposits (high beach) 6S_R Thin marine littoral deposits with very frequent organic deposits R2A 4A Lacustrine plain 4 3 4AY Lacustrine plain (average thickness -water) frequent outcrops (<16 ha) (deep water) R7T 7_R Organic deposits with very 4GAR Very thin glaciolacustrine sediments (<16 ha) frequent marine littoral deposits 4GAM Thin glaciolacustrine sediments (deep R7 7_6 Organic deposits with very R4GA 7 Organic deposits 7 Organic deposits water (thickness > 1m) 1 m 4GA Glaciolacustrine sediments (deep water) 4GA Glaciolacustrine sediments; deep beach) -average thickness between 0.5 -R4A 6SY Thin Marine littoral deposits (high between 0.5 -1 m) between 0.25 -1 m) 4A Lacustrine plain beach) over rock (average thickness Lacustrine deposits 6SM Thin Marine littoral deposits (high 6SR Thin sandy marine littoral deposits R2AE very frequent outcrops (<16 ha) (<16 ha)	2 --4 5

  Annexe G (suite) : Compilation des résultats d'analyse isotopique des échantillons de glace. Site Échantillon_code ordre δ 2 H stdev δ 18 O stdev δ 17 O stdev Note Dxs ‰ ‰

				0,43 -21,05 0,05
	Mean value
	HA	6,26	-152,12	-20,82
	HBJ	6,94	-178,65	-24,14
	CW	5,08	-166,27	-21,95
	BM	5,53	-180,81	-21,99
	NBJ	11,75	-167,06	-21,22
		Average σ 6,77	-168,38	-21,92
		Mean range σ 1,73
	standard deviation σ 2,55

μS/cm, while samples from the snowpack yielded 𝜎 values ranging between 4.4 μS/cm and 13.5 μS/cm, for an average of 6.8 μS/cm.

: average flux-weighted in

2018, 2 : average flux-weighted in 2019.

rd eco-forestry inventory Surficial deposit of Northern Quebec

(MELCC: project no. 2058049), la Fondation de l'UQAT (FUQAT: project no. 2058036) et le programme de développement des ressources humaines cries. d'Eastmain et du conseil de band,

Code

Description of units Code Description of units Code Description of units Code Description of units Code Description of units Code Description of units Code Description of units Code Description of units

valuable GDE where increased land protection is needed. Southwards, near Amos and Val-d'Or, the well-known granular aquifers associated with the Saint-Mathieu/Berry (surrounding Ungava Bay) limit is at 195 m asl [START_REF] Gray | Dynamics of the Late Wisconsin Ice Sheet in the Ungava Peninsula Interpreted from Geomorphological Evidence[END_REF]. The results suggest that AP values ≥ 4 only represent 9.3% of this territory (with 3.5% as confined aquifers and 5.8% as unconfined aquifers). The high AP values associated with unconfined aquifers are restricted to small glaciofluvial formations that are scattered in the territory. The latter are hardly noticeable in Figure 2.13 due to their small extent.

Figure 2.13 Aquifer potential values calculated based on the multicriteria analysis applied on a 490,000 km 2 territory (Québec, Canada).

Understanding riverine processes during the ice-on period

River monitoring during the ice-on season is challenging due to (1) safety issues during sampling on ice, especially at the beginning and at the end of the season, (2) difficulties related to discharge measurements under the ice, (3) water sampling constraints due to quick freezing in sampling tubes and bottles when the air temperature is significantly below zero, and (4) constraints associated with ice formation on probes used for in situ physicochemical measurements. Due to these constraints, among others, river monitoring is rarely performed during the winter. Here, an attempt was made to collect ice cores from rivers near the end of the ice-on period and to use the isotopic composition of the ice to document hydrological processes in rivers. Twenty vertical profiles of river ice were constructed from the available ice samples (𝑛 = 192; Appendix G). The collected ice cores are generally < 1 m in length and tend to present two main types of ice (Figure 3.9). A significant layer of snow-ice, which is a mixture of precipitation and surface water, generally overlies the congelation ice, which is formed by the freezing of streamflow water. This pattern was reported in previous studies [START_REF] Ferrick | Diffusion model validation and interpretation of stable isotopes in river and lake ice[END_REF][START_REF] Gibson | Isotopic characteristics of ice cover in a large northern river basin[END_REF][START_REF] Sophocleous | Interactions between groundwater and surface water: the state of the science[END_REF]. While congelation ice was observed in most ice cores (13 out of 20), only four of the cores had sufficient thicknesses to allow for reconstructing a profile based on 5 cm samples. In these four ice cores, the results revealed constant 𝛿 2 𝐻 and 𝛿 18 𝑂 values, suggesting a stable isotopic composition of river water during the ice-on period. The results did not show a progressive depletion in 𝛿 2 𝐻 and 𝛿 18 𝑂 in congelation ice, and Rayleigh-type conditions are thus unlikely in the studied watersheds [START_REF] Gibson | Isotopic characteristics of ice cover in a large northern river basin[END_REF]. The isotopic composition of sub-ice waters was compared to the isotopic composition of the basal congelation ice. This comparison allowed the calculation of the effective ice-water isotopic separation (𝜀 𝑒𝑓𝑓 ), with values of 2.66 ± 0.5 ‰ and 17.59 ± 2.29 ‰ for δ 18 O and δ 2 H, respectively. These results are similar to

8 -88,54 0,14 -12,37 0,04 -6,61 0,01 subice-water

Annexe G (suite) : Compilation des résultats d'analyse isotopique des échantillons de glace.

Où 𝑀 𝐿 correspond au facteur multiplicatif de redimensionnement spécifique à la bande 10 (3,3420×10 -4 ), 𝑄 𝑐𝑎𝑙 sont les valeurs quantitatives et calibrées des pixels de la bande 10, 𝐴 𝐿 correspond au facteur additionnel de redimensionnement spécifique à la bande 10 correspondant à 0,1 et 𝑂 𝑖 représente le facteur de correction de la bande 10 correspondant à 0,29 [START_REF] Wang | An Improved Mono-Window Algorithm for Land Surface Temperature Retrieval from Landsat 8 Thermal Infrared Sensor Data[END_REF], Ensuite, LY est converti en température (°C)

de luminosité (Brightness temperature : BT; Équation 2) :

Où 𝐾 1 et 𝐾 2 sont des constantes thermiques de conversion spécifique à la bande 10 (degré Kelvin) soit 774,8853 et 1321,0789 respectivement et 𝐿𝑌 correspond à la radiance spectrale du haut de l'atmosphère obtenue de l'équation 1, Ensuite, l'émissivité de la surface du sol est calculée sur la base du NDVIm, L'indice NDVI est utilisé pour estimer la proportion de végétation (PV) selon l'Équation 3 suivant la méthode des seuils [START_REF] Sobrino | Land surface temperature retrieval from LANDSAT TM 5[END_REF][START_REF] Wang | An Improved Mono-Window Algorithm for Land Surface Temperature Retrieval from Landsat 8 Thermal Infrared Sensor Data[END_REF][START_REF] Ermida | Google Earth Engine Open-Source Code for Land Surface Temperature Estimation from the Landsat Series[END_REF]; 4) :

Finalement, le LST est donnée par l'Équation 5 (Statoupoulou et Cartalis 2007) :

Où BT correspond à la température (°C) de luminosité (Équation 2), 𝜆 est la moyenne de longueur d'onde de la radiance émise et spécifique à la bande 10 dont la valeur est de 10,8, E est l'émissivité calculée à l'équation 4, et 𝜌 est la constante donnée par l'Équation 6 : Quantifying the water resources stored in shallow aquifers and their hydrological exchanges with surface waters is critical to better evaluate the sensitivity of groundwater resources and groundwater dependent ecosystems (GDEs) deriving from climate change and human impacts. The evaluation of water resources is still challenging in vast and remote areas, where pre-existing data are sparse and limited. This study focuses on the use of large-scale GIS-based approaches to produce new insights into shallow groundwater resources. The study area covers 185,000 km 2 encompassing vast boreal watersheds of the James Bay basin (Canada). A spatial multi-criteria analysis (MCA) based on geological and hydrogeological data was coupled with geometric calculations in order to evaluate the volume of shallow granular aquifers and to identify areas where groundwater protection should be prioritized. Satellite data and spatial models combined with in situ data, such as discharge and geochemical data (stable isotopes [δ 2 H-δ 18 O] and electrical conductivity of water), were then used to calculate terrestrial water budget and estimate water cycle components. Finally, a new index based on the stability of the Temperature-Vegetation dryness index, calculated from 290 images of the satellite Landsat 8, revealed location of potential occurrence of groundwater dependent ecosystem (POGDE).

Results revealed that the highest Aquifer Potential are sparsely distributed and mostly associated with glaciofluvial formations. They contain an average of 8.5 m of saturated sediment for approximately 40 km 3 ± 10% of water over 7% of the territory. The groundwater contribution to river discharge is generally below 12% of watershed precipitation and reach 2.3 km 3 over study area. The rainfall and snowmelt contributed is in proportion of 30-61% and 18-40%, respectively with evaporation over inflow ratios of 2-10%. The drained water presents seasonal patterns with highest snowmelt contribution during the spring flood, rainwater dominating in fall, a composition evolving towards that of groundwater during ice-on period and groundwater providing a stable baseflow contribution throughout the year. Furthermore, isotopic mass balance suggest that plant transpire mostly water from rainfall. The land near streams in main valleys present highest POGDE. A Principal Component Analysis of watershed characteristics revealed significant correlations between groundwater parameters estimated during this study. A watershed topology is proposed from the two first components which explain 77% of the variance. From that one, the southern watershed mainly depends on groundwater flows compared to the northern one, which had more impermeable surface. Four clusters of watersheds locate where increased groundwater vigilance and preoccupation are required to preserve water resources and the current biodiversity that depends on groundwater.