Topographie data

The main topographic data used in this study is a I-meter resolution digital elevation model (DEM) derived from two WorldView-2 (WV2) satellite stereo images acquired on July 31 , 2010, covering the lower glacier and September 18, 2010, covering the upper glacier. The merged DEM covers the whole glacier and surrounding topography. It was further extended with tiles from the Canadian Digital Surface Model (CD SM) (20-meter resolution) in order to include ail topography susceptible to cast shadow on the glacier. The merged DEM was resampled to 100 meters resolution using nearest neighbour interpolation to allow faster ca1culation with the mass balance mode!. The resulting grid is comprised of355 columns and 300 rows, and has approximately the same extent as the NARR gridpoint used for model forcing (Fig. 1 b). The glacier extent of August 2009 from Tennant and Menounos (2013) was used as boundary of the 2010 WV2 DEM. The slope, aspect and sky-view factors were derived from this DEM with ArcGIS toolboxes in order to be used as inputs in the mass balance mode!.

The fim area was delimited by a mean snowline delineated from LANDSAT satellite images from year 1986 to 2013 (Fig. lc). 18 cloud-free images were chosen near the end of the hydrological season (September 30) and used to map the mean snowline position at the end of summer. Image dates ranged between August 22 and October 2, which was necessary to find cloud-free images close to the end of the ablation period to provide the best estimate of the equilibrium line altitude (ELA) as possible (Fig. lc). To take into account historical glacier recession in mass balance simulations, multi-temporal DEMs and glacier boundaries from Tennant and Menounos (2013) were used to update the glacier geometry over time in the mass balance mode!. Only the area inside the glacier extent were updated from the main DEM for every year available (1979; 1986; 1993; 2001) except 2009. We used spline spatial interpolation to convert the elevation point of Tennant and Menounos (2013) into a grid. We used the glacier extents to extract and merge the interpolated elevations on the glacier with the offglacier topography of the main DEM (2010). We discarded the elevation points ofTennant and Menounos (2013) for 2009 because the points in the accumulation area were not dense enough to perform the spatial interpolation. The WV2 2010 DEM was used instead in the dynamical mass balance simulation. Two static mass balance simulations were performed using the 1979 and 2009 geometries and compared with the dynamical simulation in which the glacier geometry was adjusted when DEMs were available (1979; 1986; 1993; 2001; 2009), in order to see the impact of glacier recession on mass balance. The sensitivity analysis was performed using the GLIMS extent of 2011 before we had access to elevation points and glacier extents of Tennant and Menounos (2013). An ArcGIS model builder with every manipulation and parametrization of the DEMs is available in the supplementary material.

Mass balance observations

Mass balance is measured by the GSC since 2012 using the glaciological method (Ednie et al., 2017) and the observations available were used to validate the modelled mass balance. End-of-winter mass balance observations (bw) include snow depth soundings at, and between ablation stakes along the glacier centerline (Fig. lc). Snow depth were converted to snow water equivalent (SWE) using snow density measured at a network of reference snowpits dug in the accumulation zone, near the ELA and in the ablation zone, and complemented by snow cores (Ednie et al., 2017). End-of-summer ablation (bs) was measured at a network of 13 stakes along the glacier centerline (Fig. lc). The number of bs observations varied between years due to sorne stakes emerging completely from the ice surface and falling off before field visits, or because stakes in the upper part of the glacier could sometimes not be reached by foot during field visits. bw observations are more numerous because the upper glacier was accessed by helicopter at the end of winter and the surveys do ne on skis, and because additional snow soundings were made between ablation stakes.

The annual mass balance (ba) was calculated by summing the winter and summer balance data. The data obtained over these 5 years were used to validate the mass balance of the model. A second, independent validation was perfonned by comparing the mass balance reconstructed by the model with cumulative geodetic mass changes from 1979 to 2009 estimated by Tennant and Menounos (2013). Tennant and Menounos (2013) created multi-temporal DEMs of the entire Columbia Icefield using Interprovincial Boundary Commission Survey (IBCS) maps for glacier extent and photogrammetry of aerial photographs from 1979, 1986 and 1993. For year 1999, they used a 30 m resolution Landsat 5 Thematic Mapper (TM) Image of September 2001 for glacier extent and the Shuttle Radar Topography Mission (SRTM) DEM of February 2000, which they atlributed to best represent the glacier surface at the end of the 1999 summer ablation season, due to the penetration of the radar wave in the following year’s winter snowpack. Finally, they used a Satellite Pour l’Observation de la Terre 5 (SPOT 5) images with a resolution of 2.5 m for the year 2009. Points matched on stereoscopic image pairs were gridded to a 100 m resolution in the ablation area and to 200 m in the accumulation are a where low contrasts resulted in a smaller number of elevation points. The point clouds were interpolated to a continuous DEM surface with a resolution of 100 m to allow their comparison with the simulated mass balance. The uncertainties on elevation changes were detennined by analyzing elevation changes on stable terrain around the icefield (Tennant and Menounos, 2013).

Reanalysis data

While the precision of the on-glacier AWS data is useful to characterize the glacier microclimate, its short and discontinuous record is not adequate to drive a physically based, distributed glacier mass balance model over a longer historical period. Therefore, meteorological reanalysis data were used to force the mass balance model over the 1979-2016 period, and the A WS data used to apply a first-order bias correction to the reanalysis data. Data from the NARR were chosen for this study because of its higher temporal (3 h) and spatial (32 km) resolution compared to other commonly used products, such as ERA interim and NCEP reanalyses. NARR data for the gridpoint closest to Saskatchewan glacier (see Fig. 1) were acquired from the National Center for Environmental Prediction (NCEP) at the National Centers for Atmospheric Research (NCAR). The NARR gridpoint elevation is 2430 m a.s.\., i.e. 237 m higher than the elevation of the glacier A WS used to downscale the data. The following variables were used: instantaneous values of air tempe rature and relative humidity values at 2 m above the surface (TMP2m-ANL, RH2m-ANL), wind speed vectors at 10 m above the model surface (U and Y wind components: UGRDIOm-ANL, YGRD10m-ANL), surface 3-hourly accumulated precipitation (APCPsfc-ACC), 3-hourly averaged surface downward shortwave radiation fluxes (DSWRFsfc-A VE) and instantaneous values of total cloud cover (CC: TCDCclm-3hr).

The 3-hourly data were interpolated to an hourly interval using linear interpolation and compared with the glacier A WS observations from 2014-2016. 3-hourly NARR variables were interpolated to the center of the hourly averaging interval used by the A WS datalogger, but attributed an end-of-the hour time tag as used by the datalogger. For instantaneous variables (ANL) the concurrent time tag was used for the interpolation while for averages (AVE) the time at the center of the averaging interval was used. The 3-hourly accumulated (ACC) precipitation totals were disaggregated to hourly values by dividing the 3-hour totals into three exact quantities. Both the incoming solar radiation and air temperature have strong diurnal cycles. Over the year, solar noon varies between 12 h 41 to 12 h 56 and sunshine duration between 7 h 45 to 16 h 45. The 3-hourly NARR data cou Id thus underestimate the daily peaks in solar radiation and air temperature, especially sin ce midday NARR 3 hourly average value spreads between Il h 00 and 14 h 00. However, given that solar noon occurs near the middle of this interval the NARR midday solar radiation average may in fact weil approximate the peak mid-day value, while the 14 h 00 instantaneous temperature value is close to the time of maximum daily temperature. Still, to reduce the probability of the diurnal cycle being attenuated in the interpolated NARR data, a shape-preserving piecewise cubic interpolation was used to interpolate these variables to an hourly interval.

The NARR meteorological variables were downscaled to the A WS site using a simple bias correction procedure (Cannon et al., 2015 ; Wetterhall et al., 2012) using Meteo Lab (ML) Matlab® Toolbox. As precipitation was not measured at the glacier A WS, the homogenized historical precipitation record was used to downscale the NARR precipitation data. Three simple bias correction methods were tested (scaling, delta and empirical quantile mapping – EQM) and applied to the temperature, relative humidity, wind speed and precipitation data. Because errors in incoming solar radiation can originate from improper representation of the atmospheric transmissivity and cloud cover in NARR, but also from differences in shading between the NARR smoothed topography and the real topography, a time-varying scaling bias correction method was used to correct the NARR shortwave radiation data.

A mean diurnal multiplicative correction factor was calculated by scaling the mean observed diurnal SWin cycle with that of the houri yinterpolated NARR. A separate diurnal correction factor was calculated for each month of the year, to account for the seasonality in sun angle and related errors between NARR and observations. The bias correction methods were evaluated against the glacier A WS data using split sample cross-validation, and compared with the baseline performance, i.e. without corrections to the NARR variables. The A WS data was split into two one-year sub-periods; downscaling methods were calibrated on one period and validated on the other, and then both sub-periods were inverted, and the mean validation statistics calculated. For precipitation the whole historical record was used, so validation sub-periods are longer th an for other variables. The Pearson correlation coefficient (R), mean error (bias), root mean square error (RMSE) and mean absolute error (MAE) were used for performance assessment. The downscale performance was evaluated at both hourly and daily time intervals.