Glacial Melt and Potential Impacts on Water Resources in the Canadian Rocky Mountains

As a result of global climate change, glacial melt occurs worldwide. Major impacts are expected on the dynamics of aquifers and rivers in and downstream of mountain ranges. This study aims at quantifying the melt water input fluxes into the watersheds draining the Canadian Rocky Mountains and improving our knowledge about the fate of meltwater within the hydrological cycle. To this end, we use (1) time‐variable gravity data from GRACE satellites that are decomposed into water storage compartments; (2) an ensemble of glacier information: in situ observations, geodetic measurements, and a mass balance model; and (3) in situ surface water and groundwater level observations. The glacier mass balance model estimates a total ice mass change of ~43 Gt for the period 2002–2015, corresponding to an average of −3,056 (±2,275) MCM/yr (million cubic meters per year). 78% of the meltwater total flows west of the continental divide (to the Pacific Ocean), while 22% flows east of the continental divide (to the Arctic Ocean and Hudson Bay). However, the GRACE‐derived total water storage increases, suggesting that groundwater storage compensates for the glacial melt with an increase of 3,976 (±2,819) MCM/yr. A plausible explanation is that meltwater is not immediately flowing down in rivers but rather stored locally in aquifers. This hypothesis is supported by in situ river base flow observations, showing base flow increase in basins draining the ice melt, mostly west of the continental divide. Direct in situ evidences such as well water level time series are not sufficiently available to fully support this hypothesis.


Introduction
Over the last decades, it has been observed and reported that glaciers are melting and losing mass throughout the World, for example, in the Himalayas (Bolch et al., 2012;Farinotti et al., 2015), the Alps (Bauder et al., 2007;Haeberli et al., 2007;Rabatel et al., 2018), the Andes (Kozhikkodan Veettil & de Souza, 2017), and the Rocky Mountains (Clarke et al., 2015;Demuth et al., 2008). Fluctuations in mass of the World's glaciers and ice caps have significant impacts on global sea level (Cazenave et al., 2018;Nerem et al., 2018;WCRP Global Sea Level Budget Group, 2018), regional water cycles (Radic & Hock, 2011), and infrastructure (Carrivick & Tweed, 2016). Since the end of the Little Ice Age (~1850), glaciers and ice caps have been losing mass steadily with rates of mass loss over the past two decades being historically unprecedented (Zemp et al., 2015). Recent assessments of global glacial melt indicates that mountain glaciers and ice caps have lost 259 ± 28 Gt/yr (equivalent to 0.72 ± 0.08 mm/yr sea level equivalent) during -2009(Gardner et al., 2013Zemp et al., 2019). The continental ice sheets of Greenland and Antarctic, which contain~70 m of potential sea level rise in water equivalent ice mass, have been shrinking at rates of 344 ± 48 Gt/yr (equivalent to 0.99 ± 0.13 mm/yr sea level equivalent) during the 2005-2010 period (Shepherd et al., 2012). Contribution from ice sheets is therefore proportionately low relative to the rates measured for smaller glaciers and ice caps. Updated measurements of glacier change to 2014 however indicate that contributions to global sea level rise from ice sheets have increased to 1.26 mm/yr sea level equivalent, whereas contributions from mountain glaciers have remained constant. Regardless, mountain glaciers and ice caps will continue to be significant contributors to global sea level rise (IPCC, 2014(IPCC, , 2019 as they are expected to lose 43% to 74% of their mass between 2010 and 2100 (Huss & Hock, 2018), thereby imposing additional strain on water availability for human and natural systems functioning. Rates of mass loss from glaciers in western Canada have increased fourfold since the mid-2000s (Menounos et al., 2018), with projections indicating that up to 90% of current glacier mass in the Canadian Rockies and Interior Ranges will be gone by the end of the century (Clarke et al., 2015).
The mass balance of individual glaciers has been monitored through an internationally coordinated network of in situ observations with some records extending back in time over 100 years (Zemp et al., 2015). Precise monitoring of the individual components affecting the surface mass balance of a glacier (i.e., surface melt, runoff, internal refreezing, and snow accumulation) provides a robust indication of climate change and is critical for calibration and validation of measurements derived from remote sensing and computer models. While a large proportion of the Earth's glaciers is located in remote regions, satellite technology arises as a cost-effective solution to monitor global glacier change. Repeat satellite-based surface elevation measurements of an ice mass can be used to estimate ice thickness change, and with proper assumptions of density, changes in ice mass can be inferred (e.g., Foresta et al., 2018). Gravity field time series from the Gravity Recovery And Climate Experiment (GRACE; Tapley et al., 2004) provides low-resolution observations containing glacier mass changes since 2002.
GRACE data are an unprecedented source of information for quantitative hydrosciences (e.g., Castellazzi et al., 2016Castellazzi et al., , 2018Farinotti et al., 2015;Huang et al., 2012;Longuevergne et al., 2010). The system actively orbited Earth from April 2002 to November 2017, gathering more than 15 years of time-variable gravity measurements. Raw data (Level 1) contain the core information, that is, range rate between the two satellites, which can be directly linked to gravity changes from the mission start in April 2002. GRACE Level 2 solutions provide spectral information (Stokes coefficients), which can be converted into surface mass changes (Wahr et al., 1998). Numerous "ready to use" Level 3 data solutions were created (e.g., Bruinsma et al., 2010;Huang et al., 2012;Save et al., 2016). Most recent Level 3 solutions are "mascons" (e.g., Save et al., 2016); they mark a step-up in GRACE data processing, providing unfiltered, stripe-free, and high-resolution data to end users. Despite the number of Level 3 solutions available, Castellazzi et al. (2018) points out that no guidelines exist regarding the best solution to use for any given application.
Several authors designed methods to better extract information from GRACE, for example, restoring signal amplitude loss due to filtering using scaling factors (Landerer & Swenson, 2012;Long et al., 2015) or delimiting the area of influence of an expected signal to optimize its extraction from other signal contributors (Longuevergne et al., 2010). Regardless of the strategy adopted, the introduction of a priori knowledge on mass change distribution (Castellazzi et al., 2018;Farinotti et al., 2015;Long et al., 2016) leads to significantly better results than simple spatial averaging over an arbitrary influence footprint or a watershed. The spatial a priori can be obtained from models (Long et al., 2016) or inferred using proxies (Castellazzi et al., 2018;Farinotti et al., 2015).
In western Canada, as in most cold regions, receding glaciers and ice fields are affecting surface water and groundwater resources downstream and in general the dynamics of the water cycles. Huss and Hock (2018) defined globally the expected date of the run-off "peak water," which marks a shift from increasing to receding melt season river streamflow. Although they consider that all glacier melt is transported through rivers, it is also reasonable to expect changes in aquifer recharge and storage, given the extended relationships between river streamflow and aquifer dynamics. Indeed, it is widely accepted at catchment scale that aquifer recharge near or within mountainous regions originates directly from infiltration of liquid and/or solid precipitation or indirectly from melting glaciers through streamflow (Hood et al., 2006;Roy & Hayashi, 2008). However, the complex nature of surface-water/groundwater interactions limits our ability to forecast changes to groundwater systems downstream from melting glaciers.
Tóth (2009) demonstrated regional effects of various landforms on groundwater flow and occurrence of basinal groundwater fluxes. While mountains can be seen as impermeable landforms with steep slopes, aquifer systems can be well developed within the fractured network as exemplified in the Himalaya (Andermann et al., 2012) and the Rocky Mountains (Caine & Tomusiak, 2003). In these conditions, glacier meltwater can flow both downstream in the rivers and downwards in aquifers. Recently, researchers have studied a dense network of GPS stations in the Sierra Nevada Mountains (USA) to demonstrate that mountain aquifers could store significant amount of water and contribute to vertical surface displacement (Argus et al., 2017). They also show that parching of deep soil moisture or/and groundwater (unaccounted for in land surface models [LSMs]) occurs during droughts, which suggests a general underestimation of resilient water masses locally stored in mountain ranges.
Two recent studies explored the groundwater storage change (ΔGWS) in the prairies of Alberta, east and close to the Canadian Rocky Mountains (CRM). Huang et al. (2016) decomposed the GRACE signal to infer an increasing GWS trend gradient, higher in Eastern Alberta than close to the CRM, suggesting that glacier mass loss is not the only driver of this observed increase. Bhanja et al. (2018) studied the relations between in situ ΔGWS observations in Alberta, precipitation, and GRACE-derived ΔGWS. They observed poor degrees of correlation for most watersheds connected to the CRM, suggesting that recharge processes other than direct vertical infiltration from precipitation are probably occurring.
The dominant aquifer recharge processes occurring at the foothills remain uncertain: focused or diffused near-surface infiltration, surficial horizontal groundwater exchanges through shallow groundwater systems, or basal groundwater flux through deep-seated basin fractures. Conversely, little is known about the changes associated to glacial melt on downstream flow systems, infiltration into the ground and impacts on surfacewater and groundwater-dependent communities located downstream. This is a particularly worrying situation given that glacial meltwater inputs into downstream hydrological systems are expected to change in the upcoming decades. As such, filling this science gap is crucial to support water managers and policy makers in adopting appropriate water management plans and foster adaptation to expected upcoming changes of hydrologic regimes.
The objective of this study is to quantify changes in glacial meltwater inputs into the hydrological system, and to assess the related impacts on the different water storage components. For that purpose, we use an empirical glacier mass balance model to estimate changes in meltwater inputs to all basins draining the mountains range. We infer the total water storage change in and around the study area by applying a spatial constraint to GRACE data. We decompose the total water storage change into water storage components using the glacial melt modeling results, other models, and in situ data. Finally, we explore available in situ ΔGWS observations to support our interpretation of the geodetic measurements from GRACE, including base flow change analysis and water levels measured in observation wells.

Study Area
The study region covers an area of~200,000 km 2 in western Canada. It is centered over the CRM and includes parts of the Southern Interior Ranges, the Front Ranges, the Foothills, and the Plains. The area is referred to as CRM throughout the article. The glacierized region spans the jurisdictional boundary between the Canadian provinces of Alberta and British Columbia (Figure 1), which includes >3,000 km 2 of glacier cover. The CRM glacial meltwater contributes to summer flow of many major river systems that drain the eastern slopes, including, from north to south, the Athabasca, North Saskatchewan, Bow, and South Saskatchewan Rivers (Figure 1c). On the western side of the CRM, there are two main river systems: the Fraser River and the Columbia River ( Figure 1c). The area includes important regional-and local-scale bedrock aquifers mostly of clastic, carbonate, and shale interbedded lithology, as well as some unconsolidated aquifers located in numerous buried valleys in the Plains (Figure 1d; Rivera et al., 2018).
The water fluxes in the CRM experience high degrees of variation in space and time. Precipitation, for instance, has a very steep gradient on the eastern slopes of the CRM with an average value of~1,900 Figure 1. Location of the study area on a continental-scale map (a) and a regional-scale map (b). Major river systems and mean annual precipitation rates (c) for 1970-2000 (Fick & Hijmans, 2017). Glaciers, mountains ranges, and major aquifers (d) across the study area. Aquifers 40, 17, and 10 are potentially impacted by glacial melt. Glaciers discussed in this study are presented in green dots (Wo: Woolsey, Co: Columbia Icefield, Py: Peyto, Rr: Ram River, Hg: Haig Glacier). mm/yr, which drastically contrasts with the~430 mm/yr in the city of Calgary, only 100 km to the east ( Figure 1c). The regional gradient of precipitation is not well known due to limited high-elevation weather records. However, some values can be inferred from recent regional-scale studies. Precipitation and potential evapotranspiration in the plains part of the region ranges respectively from 400 to 500 and from 300 to 400 mm/yr (Allen et al., 2014;Wang et al., 2013). Recharge varies from 350 to 200 mm/yr from the eastern slopes to the foothills, down to 100 mm/yr or less in the Plains (Allen et al., 2014;Wang et al., 2011). On the eastern side of the CRM, the average river runoff varies remarkably among the main river systems (Wang et al., 2014); the mean annual stream flow is~220 mm/yr for the Peace River,~150 mm/yr for the Athabasca, 60 mm/yr for the North Saskatchewan, and~30 mm/yr for Bow River. River discharge represents approximately 44% of precipitation in the northern parts to 7% of precipitation in the southern parts, reflecting the effects of decreased precipitation and increased evapotranspiration. On the western side, there are two main river systems: the Fraser and the Columbia systems. The mean annual runoff of the upper Fraser varies from 800 to 1,200 mm/yr, including a large annual volume of water originating from the continental divide, where the glaciers are located.
Glacial mass balance in the CRM is governed primarily by climatic net surface mass balance, that is, accumulation gained primarily from snowfall minus mass lost due to melt runoff. As in other regions, glaciers of the CRM typically experience net mass loss at the lowest elevations where melt runoff exceeds accumulation (ablation zone), and net gain at higher elevations where accumulation exceeds losses due to melting (accumulation zone). Key to estimating annual mass changes across an entire glacier surface (Bn) for any given year is knowledge of the end-of-melt season equilibrium line altitude (ELA), which is the theoretical elevation that separates the zone of ablation from the zone of accumulation, and the mass balance gradient (db/ dz) which describes mass change with elevation.

Glacier Mass Balance Modeling 3.1.1. Model
We present an empirical model to estimate annual mass balance for 3,449 km 2 of glacier cover within the CRM for the period 2002-2015. The annual climatic net surface mass balance Bn is calculated as the average point net mass balance value (m WTE [water thickness equivalent]) across the study area according to equation (1): where n is the number of cells in the study area, H is the glacier surface elevation, and CRMa is entire area of the CRM.
The glacier surface elevation is obtained from the Global Digital Elevation Model (GDEMv2; https://asterweb.jpl.nasa.gov/gdem.asp). It is derived from optical stereo image pairs collected by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) multispectral sensor. It has a nominal vertical accuracy of ±17 m (95% confidence level) with a horizontal resolution of 75 m (ASTER Validation Team, 2011). It is resampled to a 100 m resolution raster and extracted to the extents of Randolf Glacier Inventory (RGIv6, RGI Consortium, 2017) over the study area. Absolute accuracy of the RGIv6.0 polygons is estimated to be~±2% for regions of comparable size as in this study, that is,~3,500 km 2 (Pfeffer et al., 2014). Glacier area and topography are considered constant throughout 2002-2015. Due to persistently negative mass balance documented for glaciers within the study region (Marshall, 2014;Zemp et al., 2015; WGMS website: https://wgms.ch), the area represented by the RGI polygons, digitized in 2004 and 2006, is therefore likely to be slightly overestimated for the end of the study period.
Regional topography across the CRM influences large scale climate patterns characterized primarily as a transition from the relatively continental (colder and drier) conditions of the eastern slopes of the Rockies, to the inner montane conditions (wetter and warmer) of the southwestern interior ranges (Østrem, 1966, 2006). As a result, glaciers along the southwest margin of the study area receive almost twice the snow accumulation as the eastern sector (Huss & Hock, 2015), where median glacier heights average ~400 m lower than along the continental divide (RGIv6.0). Large variations in topography and precipitation, combined with a paucity of data from monitored glaciers across the CRM study area, pose significant challenges to estimating regional scale mass balance based on key glaciological parameters such as db/dz and ELA.
Large variations in median glacier height across the study area suggest a similar magnitude of ELA variability. Failing to account for these variations in the modeled ELA will thus result in large systematic errors in modeled estimates of Bn. In this study, the spatial pattern of regional ELA is modeled as a function of median glacier height, which has been shown to exhibit a strong relationship with the "balanced budget" (0Bn) ELA (Braithwaite & Raper, 2009;Cuffey & Paterson, 2010). Differences between median glacier height and the observed 0Bn ELA, of −17 and −65 m for Peyto and Woolsey Glaciers (1966Glaciers ( -1975, respectively, indicate a similarly close relationship exists for representative glaciers in the CRM. Mapping the regional ELA as a function of median glacier height therefore establishes a reliable reference surface to which adjustments in height of annual ELA can be applied. The regional ELA across the CRM is mapped from linear trend lines through median height of individual glaciers (RGIv6.0) along five major transects running parallel and perpendicular to the continental divide CRM ( Figure S5, supporting information). Initially, a decreasing trend extending southwest from the continental divide based on median glacier height change along the "E-W_Central" transect is applied to the full study area. Next, a decreasing northwest (or increasing southeast) slope parallel to the continental divide is applied to the data set east of the continental divide, and for the western section north of Woolsey Glacier. Final adjustment is implemented to account for increases in median glacier height along the Selkirk and Monashee mountain ranges region south of Woolsey Glacier. The resulting three-dimensional pattern of 0Bn provides a regional reference to which annual ELA fluctuations measured at Peyto Glacier are applied, and from which mass balance at each 100 m grid cell across the CRM is calculated.
The mass balance model developed in this study is calibrated with in situ glacier mass balance measurements from the Peyto and Woolsey Glaciers located in these end-member climatic regions of the study area ( Figure 1d). Similar to Marshall et al. (2011), we use a single db/dz value estimated at 5.2 ± 1.2 mm/m·yr WTE, which is considered representative of the CRM and the 2002-2015 time-period. We calculate db/dz for individual years as the difference in Bn at the highest and lowest elevation of the glacier, divided by the total elevation range between upper and lower values of bn. Values of db/dz for Peyto (4.9, stdev = 1.2 mm/m·yr WTE) and Woolsey (5.4, stdev = 1.2 mm/m·yr WTE) Glaciers are then calculated as the average of all individual years over the period of measurement overlap between the two sites (1966 to 1975). Our derived value of db/dz agrees to within error with that from Peyto Glacier (i.e., 5.6, stdev = 1.6 mm/m·yr WTE.) for the time period corresponding to this study (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), and to the value (6.4 mm/m·yr WTE) used for the eastern slopes glaciers in Marshall et al. (2011). The value of db/dz of 9.8 mm/m·yr WTE reported for the Haig Glacier (2001-2007) is not considered by Marshall et al. (2011) as representative of the broader CRM.
Model validation is performed by comparing our results with independent estimates of glacier mass balance derived from geodetic measurements and previously published models. Results from this validation exercise provide an indication of the ability for our model to estimate mass balance at the ice field, drainage basin, and regional spatial scales. To assess the total meltwater contributed by each drainage basin, the meltwater volume map derived from this model is cut according to the delineation of the basins.

LiDAR/DEM Validation Data
Validation of the surface mass balance model output derived in this study is performed through comparisons with geodetically measured ice thickness changes over the Columbia Icefield. With the exception of solid ice mass discharge from the terminus (~800 m wide, ice thickness at the terminus unknown) of the Columbia Glacier that flows~60 m/yr into a proglacial lake (van Wychen et al., 2018), the Columbia Icefield as a whole loses mass almost entirely by surface melt run-off . Airborne Laser Terrain Mapper instrument was mounted in a DHC-6 Twin Otter aircraft guided by a flight management system that assisted terrain following so as to approximate a constant flying height of~2,100 m above ground. A scan angle of 25°and a pulse repetition frequency of 33 kHz were utilized, resulting in a raw planimetric resolution of~1 m over a nominal swath width of~2 km. Individual laser points have a nominal height measurement accuracy of ±10 cm (1 sigma), confirmed by preacquisition and postacquisition validation surveys over nearby infrastructure (for additional information, see Demuth, 2013;Hopkinson & Demuth, 2006;Ussyshkin et al., 2006). A swath width normalization was performed and set to 1,500 m centered immediately below the flight line. The LiDAR transect data were corrected from ellipsoidal elevations to orthometric (CGVD28) using the NRCan HTv2.0 height transformation tool (GPS.H). The LiDAR survey covers~40% (or 66 km2) of the surface area of the Columbia Icefield including the full elevation range (1,700-3,400 m above sea level [asl]) and all its major flow lines including Athabasca and Saskatchewan Glaciers.
Cloud-free optical stereo image pairs were collected over the entire Columbia Icefield in August 2016 by the Pleiades Satellite constellation. The Pleiades Satellite mission is a French Space Agency (Centre National d'Etudes Spatial) initiative in which twin satellites (1A and 1B) collect high-resolution multispectral stereo images along a polar orbit track (98.2°inclination) to provide global coverage over a 26 day cycle. Pleiades imagery is optimized for digital elevation model (DEM) generation over glacierized surfaces through high radiometric range (12-bit data) optical image acquisitions, which enhance contrast and minimize saturation over highly reflective ice and snow surfaces (Berthier et al., 2014). Stereo pairs collected over the Columbia Icefield are processed to 2 m horizontal resolution DEMs with the AMES Stereo Pipeline software (Shean et al., 2016) using automatic settings and without ground control (Marti et al., 2016). With the exception of a small data gap (250 m × 100 m) near the summit of Mount Snow Dome, the resulting DEM provides complete coverage over the Columbia Icefield.
Based on an average difference in surface height over stable bedrock terrain between the Pleiades DEM and LiDAR data, a vertical offset of 16 m is applied to the Pleiades DEM. A Vertical shifting of the Pleiades DEM results in an average difference of <1 ± 1.5 m between the two elevation data sets as measured over bedrock terrain at several locations surrounding the Columbia Icefield. Both elevation data sets were acquired in late summer; therefore, uncertainties in height due to variations in the snowpack between geodetic surveys are assumed to be minimal. Absolute accuracy of the height change differences between the 2010 LiDAR surveys and the 2016 Pleiades DEM is conservatively estimated at ±2 m.
Geodetically measured thickness changes covering nearly the full elevation range of the Columbia Icefield are binned into 100 m elevation bands, density corrected, and converted to a hypsometrically averaged value of thickness change (m WTE) across the entire area of the Icefield. The upper and lower extents of the geodetically measured thickness changes, that is, 3,436 m asl and 1,777 m asl, respectively, are extrapolated to the true upper and lower elevation limits of the Icefield to account for elevation bands lacking survey coverage. Elevation bands not sampled by the geodetic measurements account for less than 0.6% of the total Columbia Icefield area. Conversion of geodetically measured thickness changes across the Columbia Icefield to WTE values is based on assumptions about the near-surface ice and firn density. We convert thickness change into water equivalent using a firn density of 550 kg/m 3 across the accumulation zone (i.e., above 2,900 m asl; Schiefer et al., 2007), and an ice density of 900 kg/m 3 below 2,900 m asl.

In Situ and Independent Model Validation Data
Comparisons between our model, in situ data, and independent models are used to quantitatively assess the accuracy of our model. The Global Glacier Evaluation Model (GloGEM; Huss & Hock, 2015) models climatic surface net balance, and frontal ablation for tidewater glaciers, to estimate recent changes in glacier mass and predict evolution of the Earth's glaciers and ice caps to 2100 (excluding ice sheets). The model has been tuned to a variety of key parameters (i.e., accumulation, melt, refreezing, frontal ablation, and geodetic glacier change) and validated through comparisons with in situ glacier-wide annual surface mass balance measurements (RMSE = 0.73 m WTE; n = 3,148) obtained from the World Glacier Monitoring Service (WGMS, 2012). Validation of our modeled surface mass balance output at the regional scale is performed through comparisons with GloGEM across the CRM as a whole, and separately east and west of the continental divide. Validation at the glacier-wide scale is performed through comparisons between GloGEM and in situ data from the calibration sites Peyto Glacier (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) and Woolsey Glacier (1966Glacier ( -1975, and noncalibration sites of Ram River (2002)(2003)(2004)(2005)(2006)(2007)(2008) and Haig Glacier (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). Comparisons between modeled estimates of glacier thinning (GloGEM and this study) with previously published geodetic thickness changes measured over glaciers and icefields in the South Interior (S.I.) and Southern Rockies (S.R.) Ranges (Menounos et al., 2018;Schiefer et al., 2007) (Cheng et al., 2011(Cheng et al., , 2013. To explore the compromise between spatial resolution and noise reduction, we consider two filtering options, based on DDK8 (lighter) and DDK5 filters (stronger). The fourth GRACE solution is the stabilized solution from the Space Geodesy Research Group (GRGS RL04http://grgs.obs-mip.fr/grace) complete up to spherical harmonics degree and order 80. Throughout the article, these solutions are referred to as CSR MASCON, T96DDK8, T96DDK5, and GRGS, respectively. They comprise 159, 156, 156, and 150 near-monthly measurements from April 2002 to early 2017, representing the time-variable gravity field converted to WTE. GRACE Level 1B data processing (T96DDK solutions) and corresponding forward models are performed using a modified version of the GRACE Matlab Toolbox (GRAMAT; Feng, 2018). Trends in mass change over the time period considered are calculated for all time series after removal of seasonal variations by fitting and removing a stationary sine curve ( Figure S8, supporting information).
GRACE total mass storage signal is considered to be solely influenced by components of the water storage and the GIA (equation (2)).
where ΔTM corresponds to the total mass change as measured directly by the GRACE system, ΔGIA is the influence of the GIA, ΔSWS is the surface water storage change, ΔSS is the snow cover storage change, ΔSMS is the soil moisture storage change, ΔIS represents the change in ice storage (IS) of glaciers, and ΔGWS the change in GWS. In this study, we estimate ΔGWS as follows:

Signal Concentration and Decomposition
Glacier mass changes are spatially concentrated, unlike, for instance, soil moisture mass changes. Mass changes are therefore best recovered considering the specific signal smoothing inherent to GRACE data. Spatial concentration functions, as proposed by Longuevergne et al. (2010), are used to restore the amplitude of the signal dampened by truncation and/or filtering. This is best operated when the distribution of masses is known (Longuevergne et al., 2013). The spatial functions (one for each GRACE solution) are multiplied by their corresponding GRACE monthly map, producing concentrated monthly values of mass change for the study area. The same spatial functions are applied to the auxiliary data used for GRACE signal decomposition.
The spatial functions are based on the GRACE forward model of a glacier mass distribution map. A uniform and synthetic mass is allocated over the glacier area and the resulting map is truncated/filtered in the same way to each GRACE solution ( Figure 2). To simplify processing, the functions are set to 0 where the truncated/filtered signal represents less than 5% of the original signal amplitude. The resulting spatial function is uniformly weighted to fully retrieve the initial synthetic mass injected ( Figure 2). Truncation at Degree/Order 200 and a Gaussian filter of 50 km radius is used for the CSR RL05 MASCON solution (Figures S1-S3, supporting information). It is considered as an acceptable approximation of its resolution according to empirical tests and personal communication with H. Save (November 2016).
Several components are subtracted from GRACE-derived TM signal (equation (2)) to isolate glacier and groundwater contributions. GIA is subtracted by considering four models (Geruo et al., 2013;Lambert et al., 2013;Robin et al., 2017;Stuhne & Peltier, 2015). Water storage change in lakes is assessed by gathering data from the Historical Hydrometric Data portal (https://wateroffice.ec.gc.ca). The portal provides water level data for most of the lakes over the study area. Given that they are surrounded by significant slopes, volume changes are assumed to be a linear function of the lake area multiplied by level fluctuations. Surface water data are considered complete when available for at least 120 out of the 153 months considered in the GRGS solution (~80%). Time series are resampled according to GRACE temporal resolution, and the trend is calculated similarly to GRACE data. A 0.25°grid (matching the grid resolution used for GRACE processing) is created using trend values from each lake storage time series. The map is then truncated and filtered to make it compatible with the four GRACE solutions considered. ΔSMS data are derived from GLDASv1 and GLDASv2.1 (Rodell et al., 2004). Snow Storage data are derived from the snow water equivalent estimates derived from the NOAH model included into GLDASv2.1 (Rodell et al., 2004). All decomposition data, including GIA estimates, are filtered and truncated according to each of the GRACE solutions they are subtracted from.

Error Estimates
Error propagated into the GRACE signal decomposition mainly arises from uncertainties related to the choice of (1) the GRACE processing strategy and its related assumptions, (2) the GIA correction model, and (3) the model used for ΔSMS influence removal. GIA and ΔSMS estimation errors are derived by considering an ensemble of various representative sources/models. To conservatively estimate these errors, the variability range resulting from different models is considered, and only end-members (minimum and maximum of all models) are injected into GRACE signal decomposition. The error is considered as the distance between the two end-members and their median. Errors from GIA and ΔSMS estimates are considered uncorrelated; hence, their total error is computed using (err1 2 + err2 2 ) 1/2 . This error is propagated into the decomposition separately for each GRACE solution. To take into account the errors linked to the choice of GRACE data processing strategy, we consider the entire domain defined by extreme values in an ensemble of four decomposed GRACE solutions. SWS and SS removal are considered to propagate insignificant errors into the decomposition. Surface water influence is estimated through a comprehensive in situ data (A2-D2) East-west transect of the spatial function used to concentrate the glacier signal for each GRACE solution. A 0.25°p ixel map of the glaciers (blue lines) is truncated/filtered according to each solution (red lines), then a weighting function (black lines) is calculated to restore the signal loss. Assuming a uniform mass distribution within the glacierized pixels, it allows retrieving 100% of the mass initially injected. In other words, the 3-D integrals of the functions 1 and 4 are equivalent. Note that even though these functions optimize glacial mass recovery, GRACE signal is still influenced by other water storage compartments (SMS, SWS, SS) and GIA.

10.1029/2018WR024295
Water Resources Research analysis. Snow accumulation is expected to have a relatively small trend over decadal time-scales (compared to other components), since end of summer is generally snow-free across the study area. We consider its related error insignificant while calculating the total error.

In Situ ΔGWS Observations
Groundwater storage change is observed through river base flow analysis and by direct observation of water levels in wells. River base flow (Qb) is the portion of flow mostly related to groundwater-surface water interactions, and not impacted by excess runoff from precipitation. For this analysis, we use a simple graphical method based on the graph lower envelope (Kliner & Knezek, 1974;Linsley et al., 1975). Low-flow annual extremes are identified and linked to form a lower curve envelope representing Qb contribution to total stream flow. The average Qb flow during the GRACE era is compared to the average of the 30 previous years (blue lines, Figure 3). Twenty-nine monthly time series of river flow are taken into account; two examples are given on Figure 3. They were downloaded from the Historical Hydrometric Data portal (https://wateroffice.ec.gc.ca), their coordinates and names are provided in Table S2 (supporting information). Direct observations of groundwater storage change are based on well water level time series. Well level time series were retrieved from the Alberta Groundwater Observation Well Network (http://groundwater.alberta.ca/ WaterWells/) and from Groundwater Information Network (http://gin.gw-info.net/service/api_ngwds: gin2/en/gin.html). We selected 14 observation wells with long-term water level measurements across the eastern part of the study area. No data could be retrieved for the western side of the study area.

Sensitivity Analysis
Reliability of our model output is dependent on the range of realistic uncertainties established for key parameters db/dz and ELA. As established in section 3, we use end member values of db/dz ranging from 4 to 6.4 mm/m·yr WTE, with a midrange value of 5.2 mm/m·yr WTE. Uncertainties in regional ELA are established based on the following: i Representativeness of median glacier height as a proxy for the 0Bn elevation. Differences between median glacier height and the 0Bn elevation for Peyto and Woolsey Glaciers are measured at −17 and −67 m, respectively. Braithwaite and Raper (2009) report agreement between median glacier height and 0Bn to ±82 m based on analysis of data for 94 glaciers obtained from the World Glacier Monitoring Service database. ii Net difference between the average of all median glacier heights in the CRM and modeled 0Bn elevation.
Results indicate an average difference of +25 m between median glacier heights from RGIv6.0 and 0Bn values modeled along transects presented in this study ( Figure S6, supporting information). From the above, we estimate ±100 m as a relatively conservative range of uncertainty in regional ELA values for the CRM. Based on a strong relationship between db/dz and ELA over Peyto Glacier from 2002 to 2015 indicating a decrease in db/dz with increasing ELA (r = 0.7), we run our sensitivity analysis for ELA +100 m using the lower db/dz value of 4 mm/m·yr WTE, and ELA −100 m for the higher db/dz value of 6.4 mm/m·yr WTE.
Results from this sensitivity analysis indicate that use of a range of values for db/dz leads to convergence (divergence) of modeled values of Bn for a decreasing (increasing) ELA. Hence, the range of Bn for the upper and lower limits of db/dz decreases to a minimum of ±0.02 under less negative mass balance conditions (i.e., ELA minus 100 m), and increases to a maximum of ±0.13 m/yr WTE under more negative mass balance conditions (i.e., ELA plus 100 m). Based on model simulations using the 2002-2015 average regional ELA, the total range of uncertainty in Bn of 1.27 m/yr WTE spans from a maximum Bn of −0.18 m/yr WTE for a db/dz of 4 mm/m·yr WTE and regional ELA −100, to a minimum Bn of −1.45 m/yr WTE for a db/dz of 6.4 mm/m·yr WTE and regional ELA +100. We therefore estimate the total model error to be ±0.

Validation With In Situ and Independent Model Data
Comparisons between in situ measurements of glacier surface level change and results from modeling developed in this study and from GloGEM, provide validation to our modeling results at the regional and glacierwide scales (Table 1). Comparisons between GloGEM and mass balance modeled in this study over the period 2002-2015 indicate near-zero differences for the entire CRM and the region west of the continental divide. A larger discrepancy over the eastern CRM indicates a slightly higher rate of thinning estimated by GloGEM of 0.27 ± 0.66 m/yr WTE relative to our model. Over individual glacier basins, comparisons between modeled (GloGEM and ours) and in situ measurements indicate agreement to within 0.18 m/yr WTE. The largest discrepancies (−0.27 m) over individual glacier basins, however, occur between GloGEM and our model, but are still well within the uncertainty of both models, that is, ±0.73 m/yr WTE and ± 0.64 m/yr WTE respectively. Overall, the average standard deviation of all comparisons (i.e., ±0.62 m/yr WTE) is lower than the nominal uncertainties as obtained from the sensitivity analysis. This demonstrates confidence in the ability of our model to estimate mass change to within an acceptable margin of error over a variety of glacierized surfaces within the CRM.
Next, we compare modeled estimates of glacier thinning with previously published geodetic thickness changes measured over glaciers and icefields in the S.I. and S.R. regions (Menounos et al., 2018;Schiefer et al., 2007; see Figure 1 for locations). Geodetic thickness change measurements performed at decadal scale time intervals between 1985 and 2018 are compared to averaged model domain results and in situ mass balance measurements from Peyto Glacier (Table 2). Despite agreement to within error between all estimates of glacier thickness change, there exist relatively large discrepancies between some estimates. A maximum discrepancy of 0.87 m/yr WTE for all three periods occurred between the two model results (GloGEM and our model) for thickness changes over the S.I. during 1985-1999, while geodetically measured changes showed values midrange between the two modeled estimates for both S.I. and S.R. Modeled estimates of thickness change agreed closely for the 2000-2009 period (i.e., differences of 0.28 m/yr WTE and 0.09 m/yr WTE for S.I. and S.R. ranges, respectively) however, geodetically measured changes during this period for both the S.I. and S.R. ranges are 4 times lower than the average modeled values. Closest agreement between geodetic and modeled thickness changes occurred during the 2009-2018 time interval when all data sets showed maximum thinning rates for the entire 1985-2018 period.
A relatively low rate of thinning is estimated by our model for 1985-1999 in the S.I. ranges and by the 2000-2009 geodetic mass balance values from Menounos et al. (2018) for the S.I. and S.R., all showing anomalously low rates of thinning relative to the in situ mass balance measured at Peyto Glacier. The low rates of geodetically measured thinning may have been due to exclusion of 1°× 1°ASTER DEM tiles containing less than 5 km 2 of glacier ice. This population of glaciers typically experiences accelerated melting and shrinkage due to lower albedo and fragmentation (Demuth et al., 2008). The decadal scale averaging over which the geodetic changes are measured also distorts the melt characteristics during this time interval, particularly for 2003 in  (Bash & Marshall, 2014;Marshall, 2014;Shea & Marshall, 2007), the large deviations mentioned above likely reflect conservative thinning rates at their most negative margin of error.
Overall, close correspondence (r 2 = 0.70) between Bn for the entire CRM and in situ measurements of Bn for Peyto Glacier from 2002-2015 provides support to our modeled glacier mass balance estimations over the multiannual time scale ( Figure S7, supporting information). Specifically, the positive modeled mass balance year of 2010 corresponds to the second lowest melt year in the Peyto record since 2002, however the highest modeled mass loss in 2003 is only the fifth highest melt year of the in situ series. This discrepancy highlights a potential limitation of the model, which may not be accurately capturing variability in the ELA versus mass balance relationship due to anomalies in snow precipitation and/or temperature lapse rates that may occur in specific years. We note however that 2003 represents an extreme variation in the ELA versus mass balance relationship, and does not significantly impact the annual modeled glacier loss when averaged over the 14 year period of study.

Discussion
While glacier change across the CRM has been previously estimated through process-based glacier models (Huss & Hock, 2015;Marzeion et al., 2012;Radic & Hock, 2014) and geodetic mass balance measurements (Menounos et al., 2018;Schiefer et al., 2007), results from these studies lack either the spatiotemporal resolution, or proper validation specific to the region. First, the spatial extent of the geodetic change measurements do not conform to hydrological boundaries and therefore do not allow for glacier melt generated from this region to be partitioned into specific drainage basins and watersheds. Second, by integrating the available in situ glacier mass balance measurements, our empirical model provides an independent approach to the more sophisticated process-based models which rely on downscaling coarse resolution reanalysis data to estimate the primary mass balance parameters. Downscaling has been shown as fairly reliable for estimating summer melt over the complex terrain of Western Canada (Shea & Marshall, 2007). It does, however, induce significant uncertainty in precipitation estimates (Anderson, 2017;Tang et al., 2018;Trubilowicz et al., 2016), which is an important component of glacier mass balance in this region (Demuth et al., 2008). By integrating in situ glacier mass balance measurements (some of which are not used in any of the global models mentioned above) and large-scale trends in topography across the CRM, our model estimates glacier change at appropriate temporal and spatial resolutions specifically for the purposes of this study.
Our high-resolution glacier melt model allows us to estimate how the meltwater input to hydrological flow systems is partitioned into the four draining watersheds (Figure 6). Glacier meltwater is an additional load of 1717, 667, 384, and 288 MCM/yr to the most upstream sections of the Columbia, Fraser, Mackenzie, and Nelson watersheds, respectively. In addition to the natural snow/rain inputs, meltwater represents a surplus input flux. Most of the meltwater volume drains westward into the Columbia River (56%) and Fraser River (22%) watersheds. Eastward draining meltwater represents~22% of the total meltwater of the CRM. It is The empirical glacier model developed in this study represents the first successful application of a mass balance gradient model over the entire CRM (i.e., S.I. and S.R.). While a thorough validation with independent model data (i.e., GloGEM) has proven the robustness of our model for estimating glacier mass balance at the icefield and regional spatial scales, reducing systematic errors is necessary for increasing model reliability at the glacier-wide spatial scale; in particular, where local differences between median glacier height and ELA may vary significantly from the regional mean. Similarly, the tendency for mass balance gradient (db/dz) to decrease (increase) under increasingly (decreasingly) negative mass balance conditions may also introduce systematic biases depending on the degree of interannual climatic variability across the region. Reducing biases and uncertainties in modeled estimates of Bn over subregional spatial scales requires more robust knowledge of the variability of ELA and db/dz. Broad-scale mapping of end-of-season snowline (i.e., ELA) from satellite and airborne remote sensing methods (Bakke & Nesje, 2011;Shea et al., 2013) is needed to provide better control on medium to small scale variability thereby improving model performance.
Additionally, integrating basin-wide climate parameters such as snow precipitation and temperature lapse rates (Bash & Marshall, 2014), as well as accounting for ongoing glacier shrinkage during the modeling period would significantly improve interannual estimates of Bn. Reducing model biases also requires an increase in the number of carefully selected in situ monitoring observation sites that sample a representative range of elevations to capture the full range of climatic and glaciological conditions across the CRM. Modeling glacier mass balance using well calibrated key glaciological model parameters will provide a simplistic and robust approach to estimating ongoing changes in glacier mass at all spatial scales in the CRM as an alternative to the more complex process-driven climatic net balance models.  (Figure 7) show annual amplitude of~20 cm WTE and no obvious trend can be visually observed. The four time series resulting from the different processing strategies are consistent. It is observed that differences between filtering levels of unconstrained solutions (T96DDK5 and T96DDK8) do not imply important differences in noise level after spatial concentration.

Influence of Nonwater GIA
GIA is the ongoing viscoelastic relaxation of the Earth in response to the past presence of large ice masses. In Canada, it can potentially have a great influence on the mass changes detected by GRACE; for example, it is equivalent to a WTE change of about 30-60 mm/yr around Hudson Bay and~10 mm/yr in the Great Lakes region.
The GIA mass redistribution effect can be determined by either GIA modeling, or geodetic observation or a combination of the two. GIA modeling can be based on a de-glacial history models (such as ICE-5G) from Last Glacial Maximum, relative sea level changes and mathematical and physical modeling, and can be constrained to present-day geodetic observations . The GIA correction can also be based on direct GPS and absolute gravity observations (Lambert et al., 2013).  et al. 2015) and Innu/Laur16 (Simon et al., 2015(Simon et al., , 2016. ICE-6G_C is used below N48°and Innu/Laur16 above N52°. A Gaussian transfer function was applied for a smooth transition from one model to the other between N48°and N52°. This combined model constrains the two GIA models at the sparse GPS observation. The study region is close to the zero-GIA uplift line (the so-called "hinge line") in North America, where a vertical crustal motion velocity of 1-2 mm/yr is observed/modeled, comparing to ≥10 mm/yr around the Hudson Bay. It also lies within a tectonic transition zone (between the Cordillera and the Interior Plains; Figure 1; Flück et al., 2003), which makes the GIA effect challenging to estimate. GIA is detected by GRACE as a major mass increase in the CRM (Table 3). While the four GIA models presented are recent and state-of-the-art, we note that results vary greatly from model to model (Table 3 and supporting information Table S1 and Figure S4) The maximum difference between GIA models can attain nearly 100% of the minimum GIA effect, transferring significant uncertainty into GRACE signal decomposition.

Surface Water (ΔSWS) and Snow Storage (ΔSS) Contributions
After analyzing storage variations of the main lakes and reservoirs located across the study area, it is observed that eight lakes are showing significant water storage trends (above 2 MCM/yr; Figure 8a). Trend rate for each lake is computed from storage change time series (Figure 8c), then spread over a 0.25°grid and transferred into GRACE resolution according to the characteristics of each solution (Figure 8b). Kinbasket Lake is the main surface water body influencing GRACE signal, due to its large volume, significant storage trend, and its central position within GRACE's sensing footprint. The resulting mass change map in WTE has an amplitude of [0-0.5] cm/yr.

Water Resources Research
The Snow Storage (SS) contribution is expected to be low, as the snow is almost completely melted at the end of each summer season (September). According to GLDASv2.1 NOAH snow water equivalent estimates, it varies within the range [80 120] MCM/yr depending on the footprint and resolution considered.  Figure 8. Assessing the influence of storage change in lakes (SWS) over GRACE signal (GRGS solution is presented here as example, i.e., truncation at Deg./Ord. 80 and no filtering). Lakes and dams with significant storage trend (>2 MCM/yr) are shown in red in (a). Their influence over GRACE ΔTM signal depends on their position within the influence footprint. SWS time series (c) are truncated/filtered to GRACE resolution and spatially weighted by the spatial functions ( Figure 2). An example of the resulting GRGS-compatible SWS trend map is shown in (b). It is observed (b) that Lake Kinbasket has the largest influence on GRACE signal, due to its position at the center of the footprint (a, b) and its large storage change during the study period (c).

Soil Moisture (ΔSMS) Contribution
Five LSM-derived Soil Moisture Storage (ΔSMS) estimates are derived over each GRACE mass concentration function ( Figure 9). All models show a negative trend. Table 4 presents the ΔSMS mass trend observed within the different GRACE solutions and over their respective spatial footprints. Based on five LSMs, a minimum and maximum ΔSMS influence is computed for each GRACE solution (last two rows in Table 4 ). Modeling land-surface processes is complex and although all models presented are physically based, assumptions occur while parametrizing and representing these processes. Each model might, therefore, be better suited for specific soil conditions and climate. While the NOAH model is considered as the most realistic globally, several papers show that an ensemble of LSMs provide a comprehensive estimate of the uncertainty in estimating total ΔSMS (Kato et al., 2007). 4.2.5. GRACE ΔTM Signal Decomposition GRACE total mass (ΔTM) trend over the CRM is slightly positive, with mass change rates in the range [795, 1,109] MCM/yr (Table 5 and Figure 7). The decomposition of GRACE ΔTM signal allows estimation of the combined ΔIS (glaciers) and ΔGWS trends at 1,074 (±1,388) MCM/yr (Table 5 and section 3). By subtracting the result from the glacier melt modeling and propagating the corresponding error, the ΔGWS trend is estimated at 3,976 (±2,819) MCM/yr (Table 5 and section 4). As expected, the ΔSMS trend removal (Table 5 and section 2) injects large uncertainties into the decomposition workflow, with standard deviations of ΔSMS trend estimates corresponding to~50% of the mean value (Tables 4 and 5). Since the spatial concentration functions (Figure 2) are designed to fully recover the glacier mass change from GRACE data, truncation and filtering do not need to be applied to the glacier melt estimation before being incorporated into the decomposition workflow.

Discussion
We present a state-of-the-art and comprehensive decomposition of GRACE data, with particular attention given to concentrated mass recovery applicable to hydrological features such as glaciers and large lakes. We include four GRACE solutions, five SMS change estimates, four GIA models, and an in situ analysis of SWS change. Error estimation relies on comparing different models and so assumes that discrepancies between models are a reliable estimation of errors. This assumption is particularly challenging to verify in  mountainous terrain where hydrological heterogeneity is high and poorly represented in global models and where auxiliary data are for the most part unavailable (e.g., GPS stations, SMS probes, gauging station in small streams, and well level monitoring data).
The two main sources of uncertainties in the GRACE signal decomposition arise from IS, GIA, and IS contribution removal. IS estimation errors are discussed in section 4. GIA correction can be subject to a high degree of uncertainty (Caron et al., 2018). Our analysis over the CRM shows that the GIA correction is subject to uncertainty of~50%. The mean GIA effect is estimated at +5.5 mm/yr and standard deviation at 2.7 mm/yr in WTE ( Figure S4, supporting information). It is important to note that these predictions tend to overestimate the GIA effect due to the following reasons: (1) GPS observation are subject to the elastic rebound due to glacial melt within the CRM; (2) The ratio being used to convert the GPS vertical velocity is derived from the absolute gravity and GPS stations east of the CRM, which is likely greater than within the CRM; (3) ICE-6G (VM5a) does not model lateral viscosity variation of the Earth's mantle, which is relatively higher in the Interior Plains than in the Cordillera (see Figure 1). Therefore, the magnitude of the GIA correction from the four models is considered as the upper limit. ΔSMS removal using global models is also an important source of error. In mountainous regions, it is typically limited by the spatial density of precipitation and temperature data, thwarting the representation of high spatial frequencies in ΔSMS. ΔSMS models do not account for the particularly shallow and highly water-conductive soil layers occurring over the region. In addition, frozen conditions at high altitudes leads to stable ΔSMS. Considering only end-members from an ensemble of models may only partially account for these sources of uncertainty. It is reasonable to think, however, that ΔSMS trend estimates are unrealistically negative, which may lead to a slight overestimation of the GWS increase (Table 5).

Observable Impacts on Groundwater Stocks 4.3.1. River Base Flow Changes
Glaciers are known to make significant contributions to stream flow (e.g., Schaner et al., 2012). They also contribute to base flow through direct drainage of late-summer melt and through direct and indirect aquifer recharge (with/without groundwater-surface water interactions; Immerzeel et al., 2012). Late-summer stream flow usually increases during glacial melt but could also decrease if cumulative glacier volume changes result in a decrease in excess run-off (Demuth & Pietroniro, 2002;Stahl & Moore, 2006). At the scale of large glaciers or assemblages of icefields and their outlet glaciers, however, there is evidence that their relative contribution to streamflow is still increasing as a result of global warming (Pradhananga et al., 2017). For a limited number of subregions of the CRM, it has been reported that because of severe reductions in glacier area and retreat to higher elevations, glacier contributions to stream flow have  . Positive numbers show increasing river base flow (and groundwater storage) during GRACE era, and vice versa. River base flow flowing eastward into Alberta has been stable or slightly increasing (a). River base flow flowing westward into British Columbia has been slightly increasing along the Fraser system and strongly increasing along the Columbia system (b). The Kootenay branch (including the Duncan River, with an increase of 1,109 MCM/yr) is probably a major contributor to this increase (b). already begun to decline (Demuth & Pietroniro, 2002;Huss & Hock, 2018;Stahl & Moore, 2006). Less is known about the impacts of glacial melt on river base flow. As the total annual water inputs to draining rivers are increasing due to glacial melt, increases in base flow and aquifer storage around streams are expected (Immerzeel et al., 2012). The late-summer melt also directly contributes to higher base flow regimes in upstream sections, as the melt compensates for low precipitation and high evaporation rates during that period. Downstream, the excess base flow is expected to decrease the late-summer discharge of surficial aquifers into rivers and, as a consequence, to increase groundwater storage.
The eastern slopes of the CRM drain 22% of glacial meltwater through its three main river systems (Figure 10a). Stable or slightly increasing base flow volumes are observed within these systems. There are important regional-scale aquifer systems in that area; hence, a slight increase in base flow can potentially reflect a large groundwater storage increase. The western slopes drain 78% of glacial meltwater through two main river systems ( Figure 10b). Both systems are showing a significant increase in base flow.
A large increase is measured in one of the tributaries (i.e., Duncan River) of the Columbia system, potentially playing a role in the surplus measured downstream.

Hydraulic Head Changes
Bedrock aquifers are dominant in the eastern part of the study region (Alberta; Figure 11), where 84% of extraction wells draw water from bedrock aquifers and 16% from shallow unconsolidated deposits (Alberta Environment, 2007). Most of the aquifers on the eastern side of the study area are of clastic and carbonate nature where groundwater flows through fractures and/or double-porosity media. The western part of the study area (British Columbia), where~78% of the total meltwater flows ( Figure 6), is relatively poor in permeable rocks, and aquifers are mostly composed of shale and volcanic formations, that is, lavas and pyroclastic materials sometimes overlain by sedimentary rocks and/or shale ( Figure 11; Rivera et al., 2018). Unfortunately, publicly available records of hydraulic head change only cover the east side of the CRM, which only drains 22% (671 MCM/yr) of the total annual meltwater volume ( Figure 6).

. Discussion
The most likely mechanisms for groundwater generated in the immediate vicinity of glaciers are (1) direct infiltration of precipitation as rain or snow; (2) infiltration of water from the melting of glaciers. Those mechanisms provide recharge to local-, intermediate-, and, eventually, regional-scale aquifers. The two sides (British Colombia to the West, Alberta to the East) of the CRM are different in terms of climate, geography, physiography, geology, and hydrogeology; hence, repartition of surplus glacial water into water storage compartments can occur in different proportions.
Groundwater recharge in the local groundwater flow systems may quickly contribute to the rivers' base flow to the west. Contrastingly, groundwater fluxes in larger (and more resilient) flow systems may become stagnant in zones associated with broad valleys in the foothills to the East. Groundwater may also infiltrate into deeper systems and remain stored for very large periods within the mountains range, before eventually discharging in the Foreland nearest to the Foothills to the east. The study area is isolated and mountainous; data are lacking to clearly identify those mechanisms. Through river base flow analysis and direct groundwater level observation, this section aims at providing clues to better understand the impact of glacial melt on the regional groundwater flow systems.
River base flow increased during 2002-2015 for each river system draining the CRM. This occurred in relative proportions of the amount of glacier melt inputs to each corresponding watershed (Figures 6 and 10). However, the aforementioned change in the base flow regime in the Columbia system could, however, be attributed not only to melting glaciers, but also to river flow regulation on account of the significant dam operations in that system. Most major dams (e.g., Rivelstoke, Arrow, Kinibasket, and Kootenay Lakes) are present along this flow system, implying important direct anthropogenic influence on river flow, which can lead to groundwater storage and base flow increase. Regardless of its driving factor, increased groundwater storage influences GRACE signal.
Direct, in situ records of groundwater head change supports conclusions from base flow analysis that groundwater storage has increased during the GRACE era, particularly in the southeastern portion of the study area. Most of these wells are located relatively far from the CRM glaciers; hence, regional groundwater storage increase might not be solely related to recent (scale of few decades) increase in glacial melt. This is supported by increasing trends in soil moisture and snow accumulation (Table 5), suggesting increased precipitation and groundwater recharge during 2002-2015. Melt-driven increase in groundwater storage could occur in relation to increased streamflow and surface-groundwater interactions. Direct mass transfers from glaciers to groundwater at the time scale of our study are unlikely due to the important transport times related to these regional groundwater flow systems (Tóth, 2009). The relationship between water from glacier melt and groundwater in the area of this study needs further investigations. Observation of base flow regimes and well levels suggests that further work in that direction would be beneficial, particularly using hydrogeochemical tools whereby hydrological fluxes can be more fully characterized.

Conclusion: Glacial Melt Drainage and Fate of Meltwater
We present an approach to rigorously estimate mass changes in the CRM and to understand the fate of glacial meltwater in this region. Specifically, we present (1) an empirical model for describing the glacial melt at high resolution and quantifying surplus of water inputs to each major drainage basin of the CRM; (2) a spatially constrained GRACE data decomposition with particular attention given to concentrated mass recovery applicable to significant hydrological features such as glaciers and large lakes; (3) a complete surface water storage and glacial isostatic adjustment analysis for the region; (4) an analysis of the mass balance for every water storage compartment for 2002-2015; and (5) an in situ data analysis attempting to confirm the groundwater storage increase, as deduced through the geodetic-based mass balance.
Our results show that glacial mass change in the CRM occurs at an average rate of −3,056 (±2,275) MCM/yr for the period 2002-2015. Glacier meltwater is an additional load of 1,717, 667, 384, and 288 MCM/yr to the most upstream sections of the Columbia, Fraser, Mackenzie, and Nelson watersheds, respectively. Most of the meltwater flows into the Columbia (56%) and Fraser (22%) drainage basins.
The GRACE ΔTM trend is surprisingly stable given the modeled net mass loss from glacial melt. The decomposition workflow points out that positive ΔGIA and ΔGWS are largely responsible for compensating the glacial mass loss in the GRACE signal ( Figure 12). The GRACE signal decomposition workflow allows observing the glacier and groundwater storage changes, which shows a combined water mass change of 1,074 (±1388) MCM/yr over 2002-2015. By subtracting the glacier melt contribution using the melt model results, we show that groundwater storage has increased by 3,976 (±2,819) MCM/yr over 2002-2015 ( Figure 12). This suggests compensation between ice mass loss and groundwater storage increase, which is in part confirmed by exploring in situ river base flow change observations. A clear spatial relationship appears between the meltwater drained by each watershed measured through melt modeling and the increased base flow detected in the in situ data analysis. This supports the hypothesis that glacial melt plays a role in the groundwater storage increase in the CRM.
ΔGWS increases can arise in two ways: (1) natural replenishment due to increased recharge from solid or liquid precipitation; (2) direct or indirect (through surface water flow) recharge from glacial melt. It is important to keep in mind the spatial sensitivity of GRACE for masses other than glaciers while interpreting our results. For example, a large ΔGWS increase at the border of the signal concentration functions (Figure 2) may have a similar impact on the concentrated ΔTM signal to a smaller ΔGWS increase in the center. Importantly, given the aquifer settings in the CRM, the former is more likely than the latter. Finally, it is important to keep in mind that other factors such as dam/reservoir operations, or changes in climate conditions can play a significant role in groundwater storage changes. Figure 12 summarizes the results from this study. It also illustrates that the ΔGWS increase potentially occurs in aquifers with contrasting response times and resilience. A surplus of input into surficial aquifers will translate into storage increase faster than in deeper, more resilient, aquifers. Nevertheless, the total storage of deep aquifers is important, with potentially large storage accumulations after sustained abovenormal input rates. Our study shows that groundwater flow could delay glacier meltwater transfers to rivers by tens to hundreds of years, which (1) would have important consequences on the dynamics of sea level rise and (2) should be taken into account while interpreting time-variable gravity measurements over melting land ice. More studies are needed to understand storage changes in the different aquifers systems of the