Figure 1: Athabasca Glacier Research Basin (AGRB, left) and Peyto Glacier Research Basin (PGRB, right) study areas in the Canadian Rockies, Alberta. Glacier ice and snow Hydrological Response Units (HRUs) used as examples for result discussions are highlighted in red.

2.2. Study Period and Environmental Conditions

The study was conducted between July 2017 and September 2021, and streamflow evaluation was performed for the four complete WYs inside this period (2018, 2019, 2020, and 2021). The 2018 year was the worst interior British Columbia wildfire season to date (Parisien et al. , 2023) – as these fires were generally upwind of the study basins they impacted the sites with smoke and soot deposition (Aubrey-Wakeet al. , 2022). The 2019 year did not have any major wildfires; however, field inspections revealed that albedo remained low in the Athabasca Glacier from soot and algae growth feeding off soot deposited in 2017 and 2018 (Aubry-Wake et al. , 2022a; Bertoncini et al. , 2022; Esser et al. , accepted manuscript). The 2020 year had no major wildfires and can be considered a control year with greatly reduced soot and algae. The 2021 year was characterized by intense heat and an unprecedented heat dome that dominated the region for many days in late June and early July (Lin et al. , 2022). Although there was also a considerable Western Canada wildfire season in 2021, only light smoke was observed in the study basins; therefore, this year was not considered a high-activity wildfire season. The 2017 year was also a high-activity wildfire season. However, DA evaluation was not performed in that year because Sentinel-2B satellite images were only available from July 2017 onwards and thus not spanning the whole glacier ablation period. The Sentinel-2B satellite launch conferred high revisit rates every 2 to 5 days in the region.

2.3. Cloud-computing Remotely Sensed Albedo Framework Implementation

High-resolution remotely sensed albedo was estimated using a framework developed by Shuai et al. (2011) for Landsat images, updated by Li et al. (2018) for Sentinel-2 images, and applied by Bertonciniet al. (2022) in the Columbia Icefield to assess the impacts of wildfires on albedo and net SW radiation. To extend the application for use in hydrological models, Bertoncini et al. (2022)’s framework was implemented in the Google Earth Engine (GEE) cloud-computing platform. The algorithm was slightly modified to run in GEE. The main differences from Bertoncini et al. (2022) include the following: the use of both MODIS Aqua and Terra platforms to retrieve BRDF, which was made to simplify the quality control (QC) steps with more observations; the BRDF inversion method allowed negative coefficients, but the inversion was run once more without the negative coefficient to alleviate the issue since there is no non-negative least squares option in GEE; and also in order to make the framework widespread applicable instead of using station-measured incoming SW radiation, ERA5-Land 9-km incoming SW radiation was utilized. In addition, Sentinel-2 reflectance atmospheric correction was performed using the 6S model (Vermoteet al. , 1997) through its Python implementation (Py6S) instead of Sen2Cor.

In summary, the framework utilizes BRDF information from MODIS to account for differences in spectral albedo due to sensor-solar angular variability. This BRDF information is then downscaled to Sentinel-2 20-m resolution, and the high spatial resolution spectral albedo is converted to a broadband albedo using Li et al. (2018) conversion equations for Sentinel-2. The algorithm can be applied worldwide because it was implemented in GEE. The algorithm was run for both study area basins, and the mean snow and ice 20-m albedo within each HRU was extracted for DA. The station-measured albedo at Athabasca Ice AWS was utilized to evaluate the remotely sensed albedo estimates. This evaluation standard error was used as the satellite albedo measurement error for DA in both basins since the pixel that Peyto Main AWS falls within is not representative of snow and ice albedo due to infrastructure and bare soil contamination.

2.4. CRHM Hydrological Modelling

CRHM was used to predict streamflow in both basins without (CTRL) and with albedo data assimilation (DA). The CRHM configuration employed was similar to that of Pradhananga et al. (2022), with an updated firn and ice HRU designation from the two 2021-08-12 Sentinel-2 images in Figure 1. CRHM is a modular physically based hydrological model with a glacier energy and mass balance and routing modelling modules (Pomeroyet al. , 2022). The models were set up with 90 and 65 HRUs for AGRB and PGRB, respectively. The CTRL run used off-the-shelf constant albedos for ice (0.30) and firn (0.55). For spin-up purposes, the model was run from 2015-10-01 to 2021-09-30. Both models were forced with hourly station-measured air temperature, relative humidity, wind speed, incoming short- and long-wave radiation, and precipitation using the Athabasca Moraine and Peyto Main AWSs (Figure 1). These forcing variables were quality-controlled using the Fang et al. (2019) methodology. No calibration was performed in the CTRL or DA runs. Another difference in model configuration from Pradhananga and Pomeroy (2022) is that this study used a glacier and firn melt module employing a katabatic calculation of turbulent energy fluxes based on Grisogono and Oerlemans (2001) and tested in Peyto Glacier by Munro (2004) and Aubry-Wake et al. (2022b). This katabatic module represents the contribution of sub-daily glacier katabatic winds to the overall melt of ice and firn, advancing upon the less physically based configuration utilized by Pradhananga and Pomeroy (2022).

Snow albedo was simulated using the same algorithm employed in the Canadian Land Surface Scheme (CLASS) (Verseghy, 2012). The CRHM version used in this study has four other albedo algorithms capable of simulating snow albedo, including from observations and simply using a constant albedo. The model can calculate albedo using Gray and Landine (1987)’s method, which accounts for snow-covered area depletion in shallow snowpacks. The model can also simulate albedo utilizing Bakeret al. (1990)’s method, which is based on a decay function that refreshes when snowfall occurs. Verseghy’s algorithm evolves upon Barker’s by accounting for differences between dry and wet snow and also taking into consideration daily mean temperatures, which is crucial in the context of heatwaves. Verseghy’s algorithm was chosen to be used in this study given the widespread use of CLASS inside land surface and hydrology prediction systems in cold regions, e.g., in the Modélisation Environmentale Communautaire (MEC) - Surface and Hydrology (MESH) model (Pietroniro et al. , 2007; Wheater et al. , 2022). The chosen albedo algorithm feeds albedo information to the utilized energy balance snowmelt model SNOBAL (Marks et al. , 1998); therefore, it is expected that albedo assimilation will also have a large impact on snow depth and SWE.

2.5. Ensemble Kalman Filter Assimilation Framework

The DA period started in July 2017 when Sentinel-2 images became available at a higher revisit frequency, with the addition of Sentinel-2B. One individual Sentinel-2 scene was sufficient to cover each basin. Sentinel-2 images with less than 30% cloud cover were selected for albedo generation and as a DA date. If clouds or shadows completely covered an HRU, DA was not performed, but assimilation was executed for the other HRUs on the same date. DA was conducted using an EnKF method, following Clark et al. (2006) and Lv and Pomeroy (2020). EnKF was run with 20 ensemble members. Station-measured variables were perturbed with standard deviations displayed in Table 1. The standard deviation for less uncertain forcings (air temperature, relative humidity, incoming short- and long-wave radiation) was reduced to half of that used by Lv and Pomeroy (2020) since there is evidence that previously employed perturbation standard deviations were disproportionately large for these variables. For instance, Tanget al. (2023) have shown that hydrological modelling uncertainty due to temperature forcing is closer to 2 ºC instead of the commonly used 5 ºC. Wind speed and precipitation remain highly uncertain variables, and their standard deviations were kept the same as Lv and Pomeroy (2020).

Table 1: Station-measured forcings perturbed by prescribed standard deviations. The type of perturbation is also displayed.