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.