Evaluation of surface mass-balance records using geodetic data and physically-based modelling, Place and Peyto glaciers, western Canada

Abstract

Reliable, long-term records of glacier mass change are invaluable to the glaciological and climate-change communities and used to assess the importance of glacier wastage on streamflow. Here we evaluate the in-situ observations of glacier mass change for Place (1982–2020) and Peyto glaciers (1983–2020) in western Canada. We use geodetic mass balance to calibrate a physically-based mass-balance model coupled with an ice dynamics routine. We find large discrepancies between the glaciological and geodetic records for the periods 1987–1993 (Place) and 2001–2006 (Peyto). Over the period of observations, the exclusion of ice dynamics in the model increased simulated cumulative mass change by ~10.6 (24%) and 7.1 (21%) m w.e. for Place and Peyto glacier, respectively. Cumulative mass loss using geodetic, modelled and glaciological approaches are respectively − 30.5 ± 4.5, − 32.0 ± 3.6, − 29.7 ± 3.6 m w.e. for Peyto Glacier (1982–2017) and − 45.9 ± 5.2, − 43.1 ± 3.1, − 38.4 ± 5.1 m w.e. for Place Glacier (1981–2019). Based on discrepancies noted in the mass-balance records for certain decades (e.g. 1990s), we caution the community if these data are to be used for hydrological model development.

1. Introduction

Glaciers and ice sheets are the largest source of fresh water on earth. Glacier area/volume adjustment in response to climate changes has significant impacts on water resources and downstream population (Demuth and Pietroniro, 2003; Demuth and Keller, 2006; Comeau and others, 2009; Intsiful and Ambinakudige, 2021). Western Canadian glaciers lost considerable mass over the last several decades with increasing rates in the period since the mid-1990s (Larsen and others, 2007; Arendt and others, 2009; Zemp and others, 2015; Demuth, 2018; Hugonnet and others, 2021). Menounos and others (2019) observed a fourfold increase in mass budget between 2000–2009 and 2009–2018 for the western North American glaciers in response to changes in temperature, precipitation and storm tracks (Walters and Meier, 1989; Demuth and Keller, 2006; Shea and Marshall, 2007; Demuth and others, 2008; Marshall and others, 2011; Menounos and others, 2019). Glacier mass loss in western Canada is projected to increase in the coming decades due to increases in annual and seasonal mean temperature, and reduced snowfall amounts (Schiefer and others, 2007; Clarke and others, 2015; Vincent and others, 2018; Hock and others, 2019b).

Glacier mass balance is the sum of accumulation and ablation over a specific period of time, typically one year (Cogley and others, 2011). Glacier surface mass balance can be directly measured in the field using the glaciological method (e.g. Young, 1981; Cogley and others, 2011). Another approach, the geodetic method, measures the elevation change between two points in time and multiplies this change by the estimated ice/snow/firn density of a given surface (Sapiano and others, 1998; Beedle and others, 2014). Unlike the glaciological method, geodetic methods implicitly include changes in mass due to basal melt or subglacial erosion, and internal changes in mass due to refreezing or melt (e.g. Zemp and others, 2013; Andreassen and others, 2016). Finally, mass balance can also be estimated using physically-based models that differ in their level of complexity (e.g. Hock, 1999; Huss and others, 2009; Radić and others, 2015; Clarke and others, 2015; Mortezapour and others, 2020).

Mass-balance records may contain error due to biased or incomplete sampling, and changes in measurement programmes (e.g. changes in sampling methods, sampling frequency or logistical challenges that lead to sparse sampling) (O'Neel and others, 2019). The uncertainties of glaciological measurements range between ± 0.1 and ± 0.4 m a⁻¹ (Braithwaite, 1986; Cogley and Adams, 1998; Cogley, 2009; Beedle and others, 2014; Andreassen and others, 2016; Pelto and others, 2019). The uncertainties of geodetic mass-balance measurements have often been reported to be less than those of glaciological measurements (e.g. Andreassen and others, 2016; Klug and others; Pelto and others, 2019), but to achieve low uncertainties the geodetic approach requires accurate elevation models and boundary delineation, adequate temporal separation between consecutive surveys and reliable assumptions for snow and firn density (Fischer, 2011; Huss, 2013; Pelto and others, 2019). The geodetic method calculates the true or hydrological mass balance. However, it is also subject to the effects of dynamic readjustment from a previous climate forcing episode and may be out-of-phase with the reaction/volume adjustment of the glacier (Demuth and Ednie, 2016).

Long-term glacier mass change records are needed to assess the sensitivity and long-term response of glaciers to recent climate change. These records can also be used to improve physically-based models of glacier runoff. Canada is home to over 25 000 individual glaciers (Arendt, 2017), yet records of glacier mass balance that exceed 40 years exist for only five glaciers and three ice caps (WGMS, 2021). Three of these glaciers (Peyto, Place and Helm) are located in western Canada. None of these records has, to our knowledge, been evaluated for long-term biases that may arise from gross error, change in measurement protocol or other factors. Reanalysis of glacier mass-balance data improves the data quality and uncertainty estimates (Zemp and others, 2013; Thomson and others, 2017). Homogenization of glaciological mass-balance series is important to understand the long-term climate changes (Huss and others, 2009).

In this study, however, we lack all the measurements needed to carry out a complete homogenization of glaciological mass-balance records. We, instead, evaluate the glaciological records against geodetic observations, and compare annual mass balances calculated with a calibrated mass-balance model against the glaciological observations to identify years where the reported observations may be biased. We use historical air photos and Light Detection and Ranging (LiDAR) laser altimetry data collected for Peyto and Place glaciers to generate elevation models and calculate geodetic mass balance. The specific objectives of this study are to:

1. 1. Calculate geodetic mass balance at the two long-term observation sites in western Canada using historical air photos and LiDAR data;

2. 2. Use a physically distributed mass-balance model forced with climate reanalysis data to simulate the mass-balance series of Peyto and Place glaciers;

3. 3. Assess the importance of updating glacier elevation and extent on surface mass balance of Place and Peyto glaciers; and

4. 4. Identify possible reasons for the differences in glaciological, geodetic and modelled glacier mass balances.

2. Study area

Peyto Glacier (51.66° N, 116.56° W) is located in the Canadian Rockies in Banff National Park, Alberta, Canada (Fig. 1). The glacier has a northeast aspect having complex contributing basins in its upper reaches converging to a typical valley glacier snout configuration, with elevations ranging from 2100 to 3200 m above sea level (a.s.l). Temperature records available from a climate station located on the glacier (Pradhananga and others, 2021) for the period 2013–2018 show that the daily average temperature varied between + 15 and − 30° C. Based on the nearest meteorological station at Bow Summit (Fig. 1), total winter precipitation varies between 400 and 800 mm and in summer between 200 and 500 mm.

Fig. 1. Locations of the climate stations and extents of Place and Peyto glaciers in the beginning (1981/1982) and end (2020) of study period. Coordinate system used: EPSG 3857.

Place Glacier (50.42° N, 122.6° W) is located in the southern Coast Mountains, British Columbia, ~450 km west of Peyto Glacier (Fig. 1). The glacier faces northwest in its lower reaches almost orthogonal to the aspect of the accumulation region, and ranges in elevations from 1800 to 2650 m a.s.l. Based on temperature records available for Pemberton meteorological station, the monthly average temperature varied between + 22 and − 6° C during the last 10 years (Vincent and others, 2020). Annual mass-balance measurements on Peyto and Place Glaciers were started in 1966 and 1965, respectively, as part of Canada's contributions to the International Hydrological Decade (Østrem, 2006; Demuth and Keller, 2006).

3. Data and methods

Our approach uses (1) digital elevation data derived from historical air photos and recent airborne LiDAR surveys to calculate geodetic glacier mass balances, (2) climate reanalysis data to force a mass-balance model, (3) glaciological mass-balance data collected and assessed by Environment Canada and Natural Resources Canada and subsequently reported by the World Glacier Monitoring Service (WGMS, 2018), and (4) homogenized climate station data (Pemberton) and automatic weather station data (Bow Summit) to check the reliability of the reanalysis data (Fig. 2).

Fig. 2. Methodology to analyse and validate the mass-balance records. Geodetic mass balance is calculated from the aerial photos and LiDAR data. Mass-balance model consists of SnowModel and ice dynamics model, glaciological mass balance is available from WGMS records.

3.1 Data

3.1.1 Aerial photographs

Aerial negatives from the Canadian National Air Photo Library, GeoBC and the Air Photo Distribution branch of the Alberta government were photogrammetrically scanned at various resolutions (Table 1). The scale of the air photographs ranged from 1:15 000 to 1:90 000. We used MicMac, an open source photogrammetry suite, to generate elevation models from the stereo photos (Rupnik and others, 2017). MicMac uses Structure from Motion (SfM) to solve the camera pose and reconstruct the 3D geometry from the 2D scenes using bundle adjustment on a collection of overlapping images and creates a sparse point cloud (Snavely and others, 2006; Westoby and others, 2012; Pierrot-Deseilligny and others, 2014). A multiview stereo image matching (MVS; dense image matching) technique is then applied to the scaled and georeferenced sparse point cloud to increase the density of the points and generate surface models at multiple resolutions (Hirschmuller, 2007; Jaud and others, 2019). Relative orientation (determination of the position and orientation of one image to another) failed for four of the 13 photos acquired in 1997 for Peyto Glacier due to fresh snow (90% of glacier surface area was snow-covered) and poor dynamic range of the scanned negatives. We thus excluded these photos for further analysis. Ground control points (GCPs) from stable terrain surrounding each glacier provide absolute orientation (Table 1). We used a LiDAR hillshade image to select the GCPs. However, LiDAR data often do not cover the entire overlapping regions of the aerial photos and therefore, we also used 3 m PlanetScope imagery and SRTM DEM where LiDAR coverage does not exist to select uniformly distributed GCPs from the overlapping regions. We used SIFT (scale invariant feature transform) feature matching algorithm with a template window of 2000 by 2000 pixels, and we retained any template regions that exceeded the MicMac default correlation threshold (0.2). Fresh snow, poor contrast and missing coverage often hindered accurate elevation point measurements and resulted in data voids in the higher accumulation zones of Peyto and Place glaciers (Table 1, Supplementary Fig. S1).

Table 1. Descriptions of geodetic data used for mass-balance measurements

AA: Accumulation Area; UNBC: University of Northern British Columbia

3.1.2 LiDAR data

We completed multiple airborne LiDAR surveys for both glaciers (Table 1) using a Riegl Q-780 scanner with dedicated inertial measurement unit and global navigation positioning system (IMU-GNSS). Processing steps for the LiDAR data are similar to those described by Menounos and others (2019) and Pelto and others (2019). Post processing of the flight trajectories using the PosPac Mobile Mapping Suite (Applanix), with Trimble CenterPoint RTX yielded horizontal and vertical positional accuracy better than ± 0.15 m (Pelto and others, 2019). These georeferenced point clouds were then processed using RiPROCESS and exported to LAS (LiDAR data exchange file) format, and LAStools (https://rapidlasso.com/lastools/) were used to convert the point data to a DEM of 5 m resolution.

3.1.3 Climate data

We obtained meteorological fields used to drive the mass-balance model (described later) from ERA5 Land, a global dataset that is dynamically downscaled from the ERA5 reanalysis product, produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) (Hersbach and others, 2020; Muñoz-Sabater and others, 2021). ERA5 Land has a horizontal resolution of 0.1° (9 km) and is available from January 1950 at an hourly time step. Surface observations from a range of sources were used to evaluate the ERA5 Land climate fields. Homogenized climate data for climatological stations in Canada (monthly data, having a lot of gaps) are available from Environment Canada (Mekis and Vincent, 2011; Vincent and others, 2020), and we use Pemberton Station which is closest to Place Glacier (Fig. 1) (50.31° N, − 122.73° E, 204 m a.s.l, data availability: 1913–1991 for monthly precipitation, 1913–2020 for monthly temperature). In addition, we also analysed precipitation and temperature data recorded at Bow Summit meteorological station (Fig. 1) (51.70° N, − 116.47° E, 2080 m a.s.l, data availability: November 2008 to present), from Alberta Climate Information Service (ACIS, http://agriculture.alberta.ca/acis/) located within 10 km from the Peyto glacier terminus (Pradhananga and others, 2021).

3.1.4 Glaciological mass-balance data

In our study, we lack original field books for the years before 1990 for both glaciers and 1998–2014 for Peyto Glacier and therefore the point mass-balance data of individual stakes for many years are not available. In addition, the stakes were often removed, renamed or moved to new locations which has not been well documented and therefore reduced our confidence to use the point mass-balance data. Hence, we used Peyto and Place glaciological mass-balance data reported to the WGMS (https://wgms.ch/latest-glacier-mass-balance-data/). Annual surface mass-balance observations exist for both Place (1965–2020) and Peyto (1966–2018) glaciers, with some breaks in between. Seasonal and annual mass-balance values for 100 m elevation bands are available for 1965-1974, 1981-1989 and 1994-1995 for Place Glacier, and for 1966–1990, and 1993–1995 for Peyto Glacier. There are no mass-balance observations reported for Place Glacier for the year 2014. Glaciological mass balance of 2019 and 2020 for Place Glacier are available from WGMS (2021) and https://wgms.ch/latest-glacier-mass-balance-data/. Mass-balance observations for Peyto Glacier are not reported for the years 1992, 2019 and 2020. The reporting of glacier mass-balance data for Place and Peyto has not been consistent over the last 50 years. The elevation band mass-balance data, survey periods, glacier accumulation and ablation areas along with the uncertainties are incomplete for these glaciers (e.g. WGMS, 2021). Reporting in some years, for example, consists of mass balance per elevation band, whereas such data are not available for other years. This reporting contrasts with the long-term point raw glaciological mass-balance data available for all USGS benchmark glaciers (McNeil and others, 2021).

3.1.5 Geodetic mass-balance data

Prior to our calculations of geodetic mass balance, DEM pairs were co-registered following the methods described in Nuth and Kääb (2011) and Dehecq and others (2021). For both glaciers, the 2019 LiDAR DEMs were used as master DEMs and all other DEMs were co-registered to these surfaces. To calculate the geodetic mass balance, elevation differences between two co-registered DEMs were clipped to the glacier extent of the earliest year. As fresh snow and poor contrast often hindered accurate elevation point measurement from the aerial images, we first excluded erroneous elevation difference (dh) values that exceeded three standard deviations of the average dh calculated for 50 m elevation bins (extracted from the 2019 LiDAR DEM). We then fit a third-order polynomial to the average dh of the 50 m bins (McNabb and others, 2019) and integrated the fitted curve over the glacier hypsometry to estimate total volume change (dV). LiDAR-based elevation difference data contained few outliers (<3%). Glacier volume change for these datasets are the product of mean elevation change and the area of each elevation band. To compare to the glaciological approach (Zemp and others, 2013), we calculate the geodetic mass balance (B) in metres of water equivalent (m w.e.) by dividing the volume change by the average surface area (A) of the glacier extents and assuming an ice density conversion factor (ρ _i) of 850 kg m⁻³ and density of water (ρ _w) as 1000 kg m⁻³ (Huss, 2013).

(1)$$B = {{\rm d}V\over A}\cdot{\rho_{\rm i}\over \rho_{\rm w}}$$

To calculate geodetic mass balance of 2019–2020 using LiDAR elevation models, we followed the method described in Pelto and others (2019) and applied separate density values to ice (910 kg  m⁻³) and snow (590 kg m⁻³) to calculate the mass balance (little identifiable firn existed on either glacier at the end of the 2019 and 2020 ablation season).

Uncertainty in geodetic mass balance (σ _B) measurements were calculated following Brun and others (2017) using Eqn (2–4). Uncertainty arises from (1) the uncertainty on the rate of elevation change (σ _Δh), (2) the uncertainty on glacierized area (σ _A) and (3) the uncertainty on volume to mass conversion (σ _Δb) as follows.

(2)$$\sigma_{\Delta {\rm h}} = \sigma_{{\rm dh}}\sqrt{{A_{{\rm corr}}\over 5A}}$$

Where A is the glacier area, and A _corr = π.L ², where L is the decorrelation length (~500 m for western North American glaciers, Menounos and others, 2019) and σ _dh is the standard deviation of the rate of elevation change on the stable terrain. For more details, please refer Brun and others (2017) and Menounos and others (2019).

Uncertainty on volume change (σ _ΔV) was then calculated using Eqn (3)

(3)$$\sigma_{\Delta {\rm V}} = \sqrt{( \sigma_{\Delta {\rm h}}( p_{\rm A} + 5( 1-p_{\rm A}) ) A) ^2 + ( \sigma_{\rm A} {\rm d}h) ^2}$$

Where dh is the mean elevation change rate, p _A is the fraction of surveyed area after outlier removal, and σ _A = 0.1A, considering 10% error in glacier area (Kääb and others, 2012).

Finally, the random uncertainty on geodetic mass balance was obtained by using Eqn (4)

(4)$$\sigma_{\Delta {\rm B}} = \sqrt{( \sigma_{\Delta {\rm V}} \rho_{\rm i}) ^2 + ( \sigma_{\rho_{\rm i}} \Delta V) ^2}$$

Where, ρ _i = 850 kg m⁻³, is the volume to mass conversion factor used here, $\sigma _{\rho _{\rm i}} = 60$ kg m⁻³, is the uncertainty of ρ _i, and ΔV is the volume change.

The random uncertainty on geodetic mass balance using 2019 and 2020 LiDAR elevation models was obtained by using Eqn (5)

(5)$$\sigma_{\Delta {\rm B}} = \sqrt{( \sigma_{\Delta {\rm V}} \rho_{{\rm iL}}) ^2 + ( \sigma_{\Delta {\rm V}} \rho_{{\rm sL}}) ^2 + ( \sigma_{\rho_{{\rm iL}}} \Delta Vi) ^2 + ( \sigma_{\rho_{{\rm sL}}} \Delta Vs) ^2}$$

Where, ρ _iL = 910 kg m⁻³, assumed ice density; ρ _sL = 590 kg m⁻³, assumed snow density; $\sigma _{\rho _{{\rm iL}}} = 10\, {\rm kg}\, {\rm m}^{-3}$, uncertainty for ice density; $\sigma _{\rho _{{\rm sL}}} = 90$ kg m⁻³, uncertainty for snow density; ΔVi is the volume change of ice; ΔVs is the volume change of snow.

3.2 Methods

We use a physically based surface mass-balance model (SnowModel) with gridded elevation and glacier extent as inputs to calculate glacier mass balance. We also couple an ice dynamics model to SnowModel to update the glacier surface and extent after every year and calculate mass balance based on the updated surface and extent. In one experiment (ice dynamics off), we feed the first available glacier surface and extent (1981 for Place, 1982 for Peyto) to SnowModel and use yearly meteorological inputs to calculate surface mass balance for each year. In this first experiment, the glacier surface elevation and extent remains constant through the period of simulation. In the second experiment (ice dynamics on), we feed the first available surface and extent to the SnowModel, calculate mass balance for one year, feed this mass balance to the ice dynamics model, update the glacier surface elevation and extent for the next year and use these new data for the next year's simulation (Fig. 3). We herein use the term ‘ice dynamics on’ to refer to the simulation. We use the results with ‘ice dynamics off’ only to highlight the effects of glacier ice dynamics on mass balance. These two sets of experiments help gauge the potential biases imposed if driving a mass-balance model within a treatment of dynamical changes of the glacier that arise from both shrinkage (marginal retreat) and changes due to evolution of the glacier surface.

Fig. 3. Methodology for calculating mass balance (a) ice dynamics off, (b) ice dynamics on. The red boxes and solid black lines indicate these inputs are used only once during initial set up of the model. The blue boxes and blue dashed lines mean these inputs are changing every year.

3.2.1 SnowModel

We used the physically-based SnowModel (Liston and Elder, 2006b), to estimate the seasonal and annual mass balance of both glaciers. The spatially distributed surface energy balance model uses the interactions of spatially and temporally varying atmospheric forcing conditions with the local topography and vegetation to generate snow water equivalent distribution over glaciers and non-glacierized surfaces (Liston and Elder, 2006b). SnowModel requires inputs of meteorological forcing data (10min - 1 day time step) and spatially distributed fields of topography and vegetation (1–200 m resolution) (Liston and Elder, 2006a). We used precipitation, wind speed and direction, air temperature and relative humidity from ERA5 Land as meteorological inputs. The glacier topography for different years is available from the DEMs and the land cover class data for the first year of the model run is generated using the corresponding geodetic glacier extents.

SnowModel assumes that snow cover distribution and evolution for each point in space and time is independent of the lateral energy exchanges from neighbouring points and governed by the mass and energy linked to the melt and sublimation of snow. SnowModel consists of four submodels: (1) MicroMet, (2) EnBal, (3) SnowPack and (4) SnowTran-3D. MicroMet uses the meteorological inputs and interpolates and distributes those across the simulation domain using Barnes objective analysis scheme (Barnes, 1964, 1973; Koch and others, 1983). Using default or user defined monthly lapse rates of precipitation and temperature, the model adjusts those to the elevations provided by the topography data (Kunkel, 1989; Thornton and others, 1997). The MicroMet module generates fields of incoming and outgoing shortwave and longwave radiation, and surface pressure, considering a constant albedo for the glacier (Liston and Elder, 2006b). EnBal calculates the surface energy exchanges based on the meteorological conditions provided by MicroMet, and SnowPack simulates the snow depth and water equivalent evolution based on the precipitation and melt fluxes available from MicroMet. The SnowTran-3D module redistributes snow based on wind fields. Given the challenges of downscaling wind fields we have not used the SnowTran-3D submodel in our analysis. SnowPack calculates 16 output variables during the simulation period, and we use total runoff (R), solid precipitation (P) and sublimation (L) to calculate the mass balance (B = P-L-R). SnowModel does not consider the effects of mass accumulation by rain and subsequent refreezing in the accumulation area, and errors introduced by disregarding this mass-balance component are unknown. In our study, snowpack evolution and snow and ice melt are simulated with an hourly time step at 30 m spatial resolution. We used default values for the lapse rates in the MicroMet module (temperature: 4.4–8.8° C km⁻¹, precipitation: 0.20–0.35 km⁻¹).

3.2.2 Ice dynamics model

A glacier responds to long-term climate change by adjusting its geometry. We simulate annual changes in ice thickness and extent for both Peyto and Place glaciers with a vertically-integrated ice dynamics model based on the mass continuity equation (Clarke and others, 2015). This model assumes non-sliding (considering negligible effect of sliding on mass balance; Clarke and others, 2015) isothermal shallow ice approximation. For any grid cell (ice column), elevation change per unit time represents mass added or removed from its surface and the rate of ice leaving it (flux divergence):

(6)$${\partial H\over \partial t} = B - \nabla\cdot\vec{q}$$

(7)$$\nabla\cdot\vec{q} = ( \Gamma H^{n + 2}\mid{\nabla S}\mid^{n-1}\nabla S) $$

where $\nabla \cdot\vec{q}$ is the ice flux divergence, B is mass balance, S is the glacier surface elevation, $\nabla S$ is the surface slope, H is the glacier thickness, Γ is a positive constant depending on the values of Glen's flow law coefficient and exponent, ice density, gravity acceleration, sliding law coefficient and sliding law exponent, and n is the exponent in Glen's flow law (more details can be found in Clarke and others, 2015). Equation (6) suggests that the thickness of the glacier may change due to both the climate-induced mass balance (B) and ice dynamics ($\nabla\cdot\vec{q}$). Based on annual mass balances calculated with SnowModel, a derived glacier bed topography, and a DEM, the ice dynamics model calculates the thickness of the glacier for a given time step in the model. Glacier bed topography is obtained using the ice thickness grids provided by Farinotti and others (2019), and subtracting this from the SRTM DEM. The spin up time for the ice dynamics model was 40 years for both Place and Peyto glaciers. Glacier surface elevations were updated with monthly time-steps and grid cells were classified as ice when the ice thickness was >10 m.

3.2.3 Model calibration

Geodetic mass balances derived from co-registered LiDAR DEMs of 2006 and 2017 for Peyto Glacier and those of 2006 and 2016 for Place Glacier were used to calibrate SnowModel. The calibration procedure involved adjustments to the prescribed elevation of the meteorological input data. We initialized SnowModel with the ERA5 station elevation obtained by dividing the geopotential height (m² s⁻²) by gravity for the closest grid point and adjusted the elevation until the cumulative mass balance simulated with SnowModel matched the geodetic mass balance (Fig. 4). We used default values for all other parameters in SnowModel.

Fig. 4. Methodology for model calibration against the geodetic mass balance.

3.2.4 Initial calibration: without ice dynamics

First, we used the LiDAR data of 2006, 2016 and 2017 for an initial calibration of SnowModel. Using the 2006 DEM and landcover and meteorological data of 2006–2016 for Place and 2006–2017 for Peyto Glacier, SnowModel calculates the mass balance for this period. We compared this with the geodetic mass balance of the same period and if the difference is greater than the uncertainty of the geodetic mass balance (± 0.1 m w.e.) for this period, we updated the station elevation as stated above. Thus, we obtained an initial tuned elevation for the forcing data (Fig. 4).

3.2.5 Second calibration: with ice dynamics

SnowModel was also initialized with the forcing data elevations obtained in the first calibration step, and the first DEM and landcover classification. The mass balance after one year was calculated with SnowModel, and glacier surface elevations and extents were updated with the ice dynamics model. Updated surface and extents are then fed to SnowModel for the next year run. We repeated this process from 1981 to 2020 for Place and 1982 to 2020 for Peyto and the cumulative mass balance for the calibration period is calculated and compared with the corresponding geodetic balance. The station elevation for the forcing data is updated, and the steps are repeated until the difference of geodetic and modelled mass balance is <0.1 m w.e (Fig. 4). We use fixed dates (1st October − 30th September) for each hydrological year (e.g. Huss and others, 2009).

3.3 Uncertainty assessment of glaciological and modelled mass balance

The glaciological observations reported to WGMS do not include uncertainty estimates. The standard errors of glaciological observations have been estimated using Eqn (8) (Dyurgerov and others, 2009; Zemp and others, 2013):

(8)$$\sigma = {\sigma_{{\rm PoM}}\over \sqrt{N}}$$

where σ _PoM m w.e. is the standard deviation of the glaciological observations over the period of measurement, and N is the number of years.

ERA5 Land data do not provide uncertainty estimates for its parameters and it is recommended to use the uncertainties of the equivalent ERA5 fields (Muñoz-Sabater and others, 2021). Based on two years of ensemble spread of ERA5 temperature and precipitation fields when the modelled mass balance was either a maximum or a minimum and calculating the change in model mass balance by changing glacier albedo (as glacier albedo changes with time due to changes in the the distribution of snow, ice and firn) by ± 0.05 (default albedo for ice in SnowModel = 0.40), we consider an average uncertainty of ± 0.5 m w.e a⁻¹. We acknowledge that uncertainties exist in glacier thickness, that, for individual glaciers range between − 20% and $+ 60\percnt$ (Fig. S1 in Farinotti and others, 2019). Other studies also report considerable differences in ice thickness estimated by employing different models (Farinotti and others, 2017, 2021; Otto and others, 2022). These uncertainties are expected to influence ice flux most and likely have minimal impact on modelled surface mass balance, except in areas of strong vertical motion (zones of submergence and emergence).

4. Results

In 1947, Place Glacier was 4.37 km² (Menounos and Schiefer, 2008). By 2018, the glacier area was reduced to 2.68 km² (Fig. 1). A proglacial lake (0.08 km²) near the glacier terminus formed prior to 1981 and with further glacier retreat another lake (0.04 km²) formed around 2005. Peyto Glacier had an area of 13.98 km² in 1948 (Mlynowski and Menounos, 2012). By 2006, Peyto Glacier disintegrated into two separate glaciers that cover an area of 11.16 km² (Fig. 1) in 2018. A new proglacial lake has also formed at the terminus of the glacier in 2009. The lake has since increased in size and presently occupies ~0.16 km².

LiDAR DEM difference maps over the calibration period (2006–2016 for Place Glacier, 2006–2017 for Peyto Glacier) show the magnitude of recent elevation change for both glaciers (Fig. 5). The glacier-wide mass change for Place and Peyto glaciers, respectively, is − 1.5 ± 0.1 m w.e. a⁻¹ (− 0.04 ± 0.01 Gt) and − 0.90 ± 0.1 m w.e. a⁻¹ (− 0.11 ± 0.01 Gt). Geopotential heights of the nearest ERA5 Land grid points for Place and Peyto glaciers are 1490 and 2310 m a.s.l respectively. After the initial calibration, the adjusted elevations for the grid points were 1510 and 2230 m and after the second calibration with ice dynamics, adjusted elevations were 1560 and 2270 m a.s.l respectively. Geodetic mass budgets were mostly negative for both glaciers (Table 2), although for Place Glacier, the mass loss is comparatively more negative than Peyto Glacier (Fig. 6), especially in recent decades.

Fig. 5. LiDAR surface elevation change (m) maps for Place (2016–2006) and Peyto (2017–2006) glaciers (outlines of 2006 for Place and Peyto glaciers in blue).

Fig. 6. Time series of geodetic, modelled and glaciological mass balances for Place and Peyto glaciers. The shaded grey, blue and red zones denote the uncertainties of geodetic, modelled and glaciological mass balance respectively. Place Glacier annual glaciological mass balance is missing for 2014 and Peyto Glacier annual glaciological mass balances are missing for 1992, 2019 and 2020 from WGMS records.

Table 2. Geodetic mass balance for Place and Peyto glaciers

4.1 Comparison of annual mass-balance results per epoch

Our analysis compares modelled and glaciological annual and seasonal (summer and winter) mass balance for years where seasonal values are available. In addition, we also compare the glacier-wide mass balance from modelled and glaciological mass-balance series with specific mass balance based on geodetic data for each period. The modelled and glaciological annual mass balances for both glaciers show substantial variability between epochs (Fig. 6, Supplementary Table S2). We observed that glaciological mass balances are sometimes substantially more negative (1994, 1995) or more positive (2011, 2012) than modelled or geodetic mass balance for Place Glacier. Prior to 2000, only two years are modelled to have positive mass balances (1996 and 1999), and the model shows no positive mass balance years after 2000. For Peyto Glacier, the model simulates nearly balanced mass budgets in 1996 and 1999, whereas glaciological records show positive budgets in 1996 and 2000. The maximum modelled and observed mass loss at Peyto Glacier occurred in 1998. For Place Glacier, the greatest modelled mass loss occurs in 2014, when no glaciological records exist. The mean bias of the differences of glaciological from modelled annual mass balances is − 0.05 m w.e for Place Glacier and 0.01 m w.e for Peyto Glacier. Values of R² and RMSE between glaciological and modelled series for the whole period and periods of geodetic measurements for both Place and Peyto glaciers are provided in Supplementary Table S1. Averaged over the periods of geodetic mass-balance calculations and corrected for the dates of the geodetic measurements, we note that modelled mass balances mostly agree with both geodetic and glaciological records (Fig. 7), considering the range of uncertainty in all of these data sources. Differences between the geodetic and glaciological series are greatest for the periods 1987–1993 (Place) and 2001–2006 (Peyto), (Fig. 7). The modelled and geodetic mass balance closely match with each other for the period 2006–2016 for Place Glacier, whereas the glaciological data for years 2010, 2011 and 2012 denote years of mass gain which are not reflected in the simulated data (Fig. 6).

Fig. 7. Comparison of specific geodetic mass balance for different periods with glaciological and modelled series.

4.2 Cumulative mass-balance results

Over the period 1981–2019, cumulative geodetic, glaciological and modelled mass balances for Place Glacier are − 45.9 ± 5.2 m w.e. (− 0.15 ± 0.02 Gt), − 38.4 ± 5.1 m w.e. (− 0.10 ± 0.01 Gt) and − 43.1 ± 3.1 m w.e. (− 0.11 ± 0.01 Gt), respectively (Table 3; Fig. 8). Between 1982 and 2017, we observe cumulative geodetic, glaciological and modelled mass balances of − 30.5 ± 4.5 m w.e. (− 0.35 ± 0.05 Gt), − 29.7 ± 3.6 m w.e. (− 0.31 ± 0.04 Gt) and − 32.0 ± 3.6 m w.e.(− 0.34 ± 0.03 Gt), respectively, for Peyto Glacier (Table 3; Fig. 8). We record an average mass loss for Place [1981–2019] and Peyto [1982–2017] glaciers of − 1.2 ± 0.1 and − 0.9 ± 0.1 m w.e. a⁻¹, respectively. Based on glaciological records, these values are − 1.1 ± 0.1 and − 0.9 ± 0.1 m w.e.  a⁻¹, respectively. Year to year glaciological versus modelled cumulative mass balance of both glaciers match closely (Peyto: R² = 0.99, RMSE = 2.45 m w.e.; Place: R² = 0.98, RMSE = 1.82 m w.e.). Inspection of cumulative modelled, glaciologicaland geodetic mass balances reveals times where the modelled and observed records diverge (Fig. 8). At Place Glacier, for example, positive mass balance occurs for both the modelled and geodetic balances in the mid-1990s whereas the glaciological records display strong mass loss (Table 3). For other periods, modelled mass balances generally agree with the geodetic and glaciological results, though glaciological mass balances are slightly less negative than modelled or geodetic results in the period after 2000. At Peyto Glacier, all three approaches produce consistent results, though glaciological mass balances are primarily less negative than either modelled or geodetic approaches (Fig. 8). Neglecting ice dynamics (ice dynamics off), modelled cumulative mass loss is overestimated by 10.6 m w.e. at Place Glacier over 39 years and 7.1 m w.e. for Peyto Glacier over 38 years (Fig. 8).

Fig. 8. Modelled (with ice dynamics on and off) (1981–2020), glaciological (1981–2019) and geodetic mass balance (differencing with 1981 DEM) for Place Glacier, and modelled (with ice dynamics on and off) (1982–2020), glaciological (1982–2018) and geodetic (differencing with 1982 DEM) mass balance for Peyto Glacier.

Table 3. Geodetic mass balance with respect to reference years (1981 for Place, 1982 for Peyto), corresponding cumulative glaciological and modelled mass balance, and glacier area

4.3 Comparison of annual and seasonal mass balance by elevation

Further insights into the performance of SnowModel can be obtained by examining the annual (Fig. 9) and seasonal mass balances (Fig. 10) as a function of elevation. For Place Glacier, we compare annual mass-balance grids produced from the 1981–1993 DEM differencing (as 1989 DEM for Place is not available; see Table 1) and the glaciological and modelled annual mass-balance grids from 1981 to 1989 (Fig. 9). The area-altitude distribution of Place Glacier (using the DEM of 1993) shows that most of the glacier area lies below 2400 m elevation (Supplementary Fig. S2). The 1993 DEM of Place has no data gaps, and we assume that hypsometry did not change significantly from 1981 to 1993. At lower elevations, all three mass balance series accord. Above the median elevation (2020 m) of the glacier, the geodetic mass balance is more negative than modelled or glaciological records. The change in modelled mass balance is always positive with increasing elevation, however for both glaciological and geodetic approaches the change is either negative or zero at the higher elevations.

Fig. 9. Box plots for elevation band wise annual glaciological (1981–1989), modelled (1981–1989) and geodetic (1981–1993) mass balance grids for Place Glacier and annual glaciological (1982–1989), modelled (1982–1989) and geodetic (1982–1986, 1986–1991) mass balance grids for Peyto Glacier.

Fig. 10. Box plots for elevation band wise glaciological and modelled (Place: 1981–1989, 1994–1995, Peyto: 1982–1990,1993–1995) winter and summer mass balance.

We also compare the elevation binned estimates of modelled and observed mass balance for Peyto Glacier for the period 1982–1991 (Fig. 9) and include geodetic estimates of mass change simply to highlight the importance of ice dynamics as a function of elevation. The hypsometry plot of Peyto Glacier (using the DEM of 1991, assuming no significant change of hypsometry from 1982 to 1991) shows that most of the glacier area lies below 2900 m elevation (Supplementary Fig. S2). At lower elevations (median elevation = 2580 m), both modelled and glaciological estimates are comparable and generally mirror geodetic mass change (Fig. 9). Glaciological observations tend to be more negative at lower elevations compared to the geodetic mass balance. The geodetic mass balance is comparatively more negative than glaciological or modelled mass balance in some of the higher elevation bins for both Place and Peyto glaciers.

We further compare modelled and observed seasonal mass balances by elevation bands (Fig. 10). For Place Glacier, central tendency and variability of glaciological and modelled seasonal mass balance per elevation band are comparable (Fig. 10). Summer glaciological balance is more variable than modelled mass balance, however. For Peyto Glacier, the model and glaciological observations diverge, notable where glaciological observations record higher mass loss during summer at lower elevations and greater mass gain during winter, especially at highest elevations. However, ice dynamics also play significant role to control the mass balance and we discuss the role of ice dynamics and their effect on modelled and observed mass balances in greater detail below.

4.4 Mass balance for 2020

Due to COVID19 travel restrictions, glaciological mass-balance observations were not collected in 2020 for Peyto Glacier. However, we are able to calculate the annual balance using LiDAR surveys conducted in fall 2020 (Table 1) for both glaciers, and compare these with the modelled estimates. The 2019–2020 geodetic mass budget of Place Glacier is − 1.73 ± 0.33 m w.e., while the glaciological and modelled mass balance are − 1.65 ± 0.13 and − 1.38 ± 0.50 m w.e. respectively. Geodetic calculations indicate that Peyto Glacier was nearly balanced in 2019–2020 (− 0.10 ± 0.29 m w.e) while the modelled mass balance for the same period is − 0.31 ± 0.50 m w.e.

5. Discussion

Below, we discuss factors that contribute to differences between glaciological observations and geodetic measurements and our approach to identify possible errors in the modelled and glaciological records.

5.1 Glaciological mass balance

Glaciological mass-balance measurements involve end-of-winter accumulation measurements (snow pits or cores to measure density and snow probing to measure snow depths), and end-of-summer probing and ablation stakes to measure snow and ice melt. Measurements are taken in both accumulation and ablation regions of the glacier and interpolated between the measured points using a range of methods to estimate glacier-wide mass loss or gain (Østrem and Stanley, 1966). Glacier run-off measured using glaciological mass-balance methods is often an input to hydrological models in basins where glacier melt is an important source of streamflow (e.g. Schaefli and others, 2005; Konz and Seibert, 2010). Errors in mass-balance measurements will thus impact the model calibration and subsequent model predictions.

The glaciological mass-balance values used here are those reported to WGMS (2018). We note that while these measurements are often checked for gross errors, they do not represent a homogenized data product. Others (Huss and others, 2009; Andreassen and others, 2016; Klug and others; O'Neel and others, 2019) adjust historical records of glaciological mass balance using original field notes and standardizing methods of data interpolation and density conversion, and adjusting for glacier geometry over the period of measurement. Huss and others (2009), for example, homogenized the 50 years of glaciological records of winter and annual balance of two glaciers from Switzerland using a mass-balance model and geodetic mass balance and observed a good agreement of the homogenized series with the conventional records for Gries Glacier, but found a significant bias for Silvretta Glacier. Klug and others homogenized the glaciological and geodetic mass balances of an Austrian Glacier (Hintereisferner) for 2001–2011 and reported three years where the two have significant differences and concluded that the glaciological observation network might be inadequate for those years. O'Neel and others (2019) homogenized the point mass-balance records of five benchmark glaciers of the US Geological Survey by calibrating those against geodetic mass balance and a mass-balance model. In four of the five glaciers they studied (Gulkana Glacier, South Cascade Glacier, Wolverine Glacier and Lemon Creek Glacier), geodetic mass balance was more negative compared to the corresponding glaciological values. The record from South Cascade Glacier, which is comparable in size to Place Glacier, showed a bias of ~5 m w.e. over 32 years. We also observed more mass loss using geodetic method than glaciological mass balance for both glaciers (e.g. mass loss in Gt in Table 3), though glaciological mass loss has been slightly underestimated due to missing data for one of the years (Supplementary Table S2). O'Neel and others (2019) mentioned that one of the reasons of this difference could be an underestimation of ablation and/or overestimation of accumulation in the glaciological measurements. We have also observed much higher glaciological winter mass gain than modelled in the upper elevations of Peyto Glacier.

To determine the glaciological mass balance for the glaciers in our study, a second-order polynomial is applied to point observations to estimate the summer and winter mass balances in 100 m elevation bins. At elevations that lack point observations, mass-balance values are altered using expert knowledge or past mass-balance gradients. The use of expert knowledge normally occurs at high elevations where the mass-balance values derived from the polynomial equation starts to deviate from an expected reality. When lacking a clear equilibrium line altitude (ELA) value in the stake data, late summer Landsat imagery is used to estimate the annual transient snow line and used as the ELA. This may partly explain why we see a large difference in modelled vs. glaciological winter mass balance at higher elevations for Peyto Glacier. Snow redistribution by wind or gravity could also be a factor responsible for this difference; both processes are not simulated in the present study.

Unfortunately, robust homogenization is not possible in our study as we lack point mass-balance data of the individual stakes for many years. A less ideal but feasible option for us is to compare the records based on elevation bins. We then compare the geodetic and glaciological observations against the modelled series. The dates of the glaciological surveys vary between years, so this potentially introduces an additional source of error in the comparisons with modelled mass balances. The dates for glaciological observations are not available from the WGMS (2018) records for most years and therefore we cannot correct this error. Demuth and Keller (2006) stated that while calculating the glaciological mass balance of Peyto Glacier, the glacier hypsometry was kept constant in some periods, whereas adjusted in some other periods. Moreover, glaciological mass balance measures only the surface mass balance, whereas the geodetic mass balances include basal and internal mass balances. In addition, refreezing meltwater could affect mass balance at higher elevations.

Years without measurement in the glaciological record is also an issue. There was no mass-balance measurements of Place Glacier in 2014, and for Peyto Glacier in 1992. According to Demuth and Keller (2006) and Fluctuations of Glaciers 1990–1995 (Haeberli and others, 1995), there are no mass-balance observations for Peyto Glacier in either 1991 or 1992. They also stated that mass-balance programme started with new personnel and focus in 1993, before which it went through a very uncertain period (Haeberli and others, 1995). The large mass loss of Place Glacier in 1994–1995 compared to the model and geodetic observations may arise from uncertainties in the glaciological observations. In addition, for Peyto Glacier, the 1991 mass-balance measurement is identical as the one in the WGMS database for the year 1993, but it is unclear where this value comes from.

The interannual variability of the glaciological mass-balance series is considerably higher than the modelled series. This may be due in part to inadequate stake data, particularly in inaccessible accumulation regions, and the spatial interpolation required to estimate glacier-wide mass balances (Beedle and others, 2014; Pelto and others, 2019). Errors in ice and snow density measurements, varying glaciological survey dates and climate variability will also contribute to the interannual variability in glaciological mass balances (Zemp and others, 2013). However, our modelled mass balances may also not capture the true variability of complex systems. SnowModel requires model calibration, but the model parameters may not perfectly represent the reality. Here we used default values of all the parameters for SnowModel and forced the model using reanalysis data. A range of tunable model parameters, such as glacier albedo, temperature, precipitation and their default monthly lapse rates made it a difficult choice about which of these should be tuned for calibration. We used station elevation to calibrate the model as SnowModel then uses its MicroMet module to interpolate the station precipitation and temperature values to the whole domain of analysis and applies default monthly lapse rates to consider for elevation changes.

5.2 Geodetic mass balance

The mass balance of the glaciers in western Canada is primarily controlled by winter precipitation and summer temperature and mass balance for Place Glacier has been observed to be predominantly negative after 1976 because of the shift of the Pacific Decadal Oscillation to a warm phase (Moore and Demuth, 2001). In our analysis we observed a predominantly negative geodetic mass balance for Place Glacier, with an increasing negative trend in the most recent decades because of climate warming. Such increased mass loss in the recent decades has also been observed for some other reference glaciers in this region, for example Columbia Glacier, South Cascade Glacier, Rainbow Glacier, among others (WGMS, 2021).

DEMs used to calculate geodetic mass balance prior to 2006 were derived from aerial photographs, and those from 2006 onward were based on LiDAR data. The accuracy of glacier surface elevation measurements derived from air photos is limited by low dynamic range and poor contrast in the imagery, and in particular in accumulation areas, where fresh snow often limits the number of matching (tie) points. Lack of suitable tie point resulted in unrealistic elevation differences in the accumulation region and voids (20–40%) in the elevation difference grids and introduced high uncertainties in the mass-balance measurements. For example, fresh snow appears on the Place Glacier photos for 1987, leading to a data void of almost 40% in the corresponding DEM (Table 1, Supplementary Fig. S1). To obtain the elevation band wise annual mass-balance plots for Place Glacier, we used photos from 11 September 1981 and 31 August 1993, both of which had much less snow cover (Table 1, Supplementary Fig. S1). The DEMs generated using these photos had similar elevation ranges leading to less interquartile range in the higher elevations of the elevation band wise plots (Fig. 9). Air photos from 30 July 1982 were used to generate the 1982 DEM for Peyto Glacier (Table 1). The extensive snow cover on this date affected the quality of the DEM. This explains why we see much more variation in geodetic mass balance in the higher elevation bands of Peyto Glacier (Fig. 9).

5.3 Validation of modelled results

Substantial differences (> 1 m w.e.) in the modelled and glaciological annual mass balances at Place Glacier occur in the years 1995 and 2011 (Fig. 6). To assess these departures, we examine the relation between glaciological mass balances and estimates of ELA. We use optical imagery to assess the transient snowline, which can approximate the ELA if captured at the end of the ablation season. We used Landsat data for 1995 and RapidEye data for 2011, both imaged on 11^th September, at the end of ablation season.

The 1995 ELA was reported to be 2602 m a.s.l, i.e. close to the highest glacier elevation (~2650 m), such that the AAR was assumed 0. Reported and modelled mass balances for 1995 were − 2.486 ± 0.13 and − 1.18 ± 0.5 m w.e. (Supplementary Table S2), respectively. Analysis of the 1995 Place Glacier imagery (Fig. 11) gives an AAR of ~0.41 and an ELA of 2132 m, which should correspond to a mass balance of $-0.33\, {\rm m}\ {\rm w.e}$. according to the ELA vs. mass-balance relationship for Place Glacier based on WGMS records (Supplementary Fig. S3).

Fig. 11. Place Glacier as observed on 11-09-1995 (Landsat 5 TM), 11-09-2011 (RapidEye). Outlines are generated by digitizing the glacier boundary. Contour lines are extracted using 2019 LiDAR DEM.

For 2011, the reported ELA is 2010 m a.s.l and AAR is 0.58, while the optical imagery suggests an AAR of 0.61 (Fig. 11) and ELA of 1985 m. Based on the ELA–mass balance relation and the optical imagery, the expected mass balance for 2011 would be 0.27 m w.e. and the reported mass balance is 0.355 ± 0.13 m w.e. , however the model gives a mass balance of − 0.96 ± 0.5 m w.e.. The difference in 1995 and 2011 modelled mass balance is significantly lower (0.12 m w.e.) than that of glaciological mass balance ( 2.8 m w.e.). ERA5 monthly winter precipitation (40 and 44 mm respectively) and summer temperature (9.88 and 9.7°C respectively) in 1995 and 2011 reveal that the values were very close to each other, resulting in close modelled mass balance in these years.

We also compared the ELAs from modelled mass-balance grids with the transient snow lines (TSL) available from PlanetScope optical images. We identified few clear images from July to September 2018 and 2019 for both glaciers. Overall, we observe strong correlation (r = 0.97 for Place and r = 0.94 for Peyto) between the observed TSL and modelled ELA (Supplementary Fig. S13). The results indicate that though our model is able to capture the increase of ELA from lower to higher elevation as dates change from July to September (Supplementray Fig. S13), but the ELA is higher (100–200 m) in the model than the observed TSL. Factors that could explain this discrepancy include: (i) negative bias in modelled winter balance; and (ii) positive bias in modelled summer balance. We suspect that some of this bias could arise from neglected physics in our model simulation, namely the redistribution of snow by wind. The lack of snow redistribution may also explain the model's winter negative bias in mass balance observed for some of the elevation bins (Fig. 10).

For Peyto Glacier, the cumulative mass balance for different periods for which geodetic data are available match closely with those from glaciological and modelled series. The difference in specific mass balance of model and glaciological series for the period 2001–2006 with the geodetic mass balance may be a result of interpolation of the significant voids (40%), especially in the higher accumulation regions of the 2001–2006 DEM difference, which was the result of excessive snow cover on the 2001 photographs and data void in the corresponding DEM (Table 1, Supplementary Fig. S1).

Over the last 40 years, we observe a modelled glacier mass (area) loss of 59% (41%) and 34% (21%), respectively, for Place and Peyto glaciers. The projected (2015–2100) percentage changes in mass of the glaciers of western Canada and US reported by Hock and others (2019a) indicates that around 20–60% mass loss occurs in the first half of the period. The mean and maximum projected mass loss over 85 years is around 60% and 75% for a low emission (RCP2.6) scenario and around 75% and 95% for a high emission (RCP8.5) scenario, respectively (Hock and others, 2019a). Based on these values, the observed mass loss of Place Glacier over ~40 years (this study) is higher than the low emission scenario. The reason could be the smaller size of Place glacier as Hock and others (2019a) also mention that the small glaciers lose relatively more than the bigger ones, and often disappear completely. Our area change estimates are much less than those observed by Hock and others (2019a). This is perhaps expected as glacier length/area change is a delayed response to climate than glacier mass changes (Haeberli, 1995). Hock and others (2019a) also mention that similar values of regional glacier area loss and mass loss projected by the models are unexpected and need further analysis.

5.4 Elevation dependence of mass balance

Differences in the modelled and glaciological mass balance by elevation (Figs 9, 10) may be the result of a number of factors that include prescribed albedo, precipitation lapse rates and wind redistribution, among others. SnowModel assumes a glacier wide constant albedo, whereas the albedo decreases at lower elevations in summer due to exposed ice. This partly explains why the glaciological mass loss at lower elevations is more negative than modelled values. Another aspect could be that the default monthly precipitation lapse rates used in the model are insufficient to simulate correct precipitation at different elevations for Peyto Glacier. As we calibrated our model using the annual mass balance, the underestimation of mass loss at lower elevation in summer is compensated by the underestimation of mass gain at the higher elevations in winter. Mass redistribution by wind in winter could also play an important role at Peyto Glacier, but this is not simulated in our model runs. Finally, we updated the modelled glacier surface by considering ice dynamics, whereas we were unable to homogenize the glaciological mass balance to consider changes in hypsometry over time and its effect on mass balance.

5.5 Effect of glacier ice submergence and emergence on mass balance

Distributed glacier surface mass-balance models can be effectively used to monitor glacier changes with climate. Such models often assume that the glacier geometry, i.e. its surface and extent, typically remains unchanged over the period for which the simulation is performed. A glacier, however, is not a static body. Factors such as glacier bed topography, ice thickness and response time (Jóhannesson and others, 1989) may lead to considerable differences in dynamic adjustments given the same climate forcing. These dynamic adjustments, in turn, affect the mass balance. Reference surface mass-balance data, which omit the geometric adjustments, are better suited for climate change studies, whereas the conventional balance that considers ice dynamics is more relevant for hydrological studies (Elsberg and others, 2001).

To understand the performance of the ice dynamics model, we use yearly surfaces generated by the model to compare mean modelled surface elevations against geodetic elevations (Supplementary Figs S4, S5). For Place Glacier (Supplementary Fig. S4), modelled surface elevations are consistently higher than geodetic surfaces above 2200 m. Below 2200 m, the modelled and geodetic surface elevations at Place Glacier largely agree. At Peyto Glacier (Supplementary Fig. S5), an overestimation of surface elevation can also be seen above 2700 m, with greater agreement below 2700 m. Uncertainty in the geodetic observations overlaps with the modelled surface elevations in nearly all cases, except at higher elevations of 1997 and 2001 for Peyto Glacier due to poor contrast of the air photos in the accumulation regions. Total thickness changes between modelled surfaces in 2020 and geodetic surfaces in 1981 (Place) and 1982 (Peyto) show that the mass-balance and ice dynamics models produce a slight thickening at the highest elevations on Place Glacier, and a thinning at all elevations on Peyto Glacier (Supplementary Fig. S6). Steeper slopes in the upper accumulation area at Peyto Glacier (Supplementary Fig. S7) result in a higher mass flux and a corresponding thinning across the entire elevation range.

While geodetic observations appear to show better correspondence with modelled and glaciological observations at low elevations on both Place and Peyto glaciers (Fig. 9), the comparison is problematic for a number of reasons. First, the time periods considered are not perfectly matched. Second, uncertainties in the geodetic observations have to be considered at higher elevations. Third, the glaciological mass-balance observations do not include dynamic adjustments. In addition, firnification plays an important role in the mass balance at higher elevations. The use of a constant glacier wide density in geodetic mass-balance measurements will contribute to differences in modelled and observed mass changes at high elevations (Li and Zwally, 2011; Pelto and Menounos, 2021). Finally, we considered negligible sliding in our ice dynamics routine which could also lead to some differences with geodetic measurements.

For the full modelling period, overall mass loss is greater (by 21–24%) when the geometry (area and extent) is kept constant and glacier submergence and emergence are ignored. The yearly percentage contribution of mass loss of the area vacated by marginal retreat of the glaciers ranges between 0 and 5% (Supplementary Table S2). This suggests that although the effect is negligible on a yearly basis, the cumulative effect of the area change of the glacier for a longer period on overall volume loss is significant. The surface ice velocity estimates from Gardner and others (2019) show maximum and minimum ice velocities of 19 and 0.05 m a⁻¹, respectively, for Peyto Glacier. For Place Glacier, the data coverage is <30% and the values there are <1.5 m a⁻¹, which is significantly less than the maximum ice velocities for some other glaciers of western Canada (such as 320 m a⁻¹ for Klinaklini Glacier; 190 m a⁻¹ for Bridge Glacier). It can be expected that the effects of ice dynamics on cumulative mass balance will be more pronounced for glaciers having higher rates of ice flux.

5.6 Use of ERA5 climate data to drive glacier mass-balance models

Uncertainties in model forcing data may lead to a poorly calibrated model that can reproduce long-term mass changes, but fail to resolve annual/seasonal mass balances or mass-balance gradients. We evaluated the ERA5 Land dataset using homogenized station data near Place Glacier (Pemberton). Temperature and precipitation records for Pemberton, BC are available for 1912–2020 and 1913–1991 respectively, though there are gaps in the records. Using all available values from 1982, we calculate monthly mean temperature and monthly total precipitation from ERA5 Land data and plot those across the station data. We find excellent correlation for temperature (r = 0.98), and precipitation (r = 0.92) (Supplementary Fig. S8), though ERA5 precipitation is consistently higher and ERA5 temperatures are consistently lower than station values. This is likely due to the differences in ERA5 grid point elevation (1490 m) and station elevation (204 m). The high correlation is perhaps expected as observations are fed through data assimilation (technique by which observations from different sources are processed and adjusted, so that initial conditions are as close as possible to the reality; Asch and others, 2016) to ERA5 Land, which is based on a numerical weather prediction model (Hersbach and others, 2020). Our intention to use ERA5 land to drive a glacier mass-balance model is to also explore how reliable ERA5 might be to estimate changes in glacier mass balance within all of western Canada.

For Peyto Glacier, we compared the ERA5 Land precipitation and temperature with Bow Summit station data. The station data are available at hourly time scale. We converted hourly to daily variables and plotted those across ERA5 daily precipitation and temperature (Supplementary Fig. S9). We observe a weak correlation for precipitation (r = 0.72) but a strong correlation for temperature (r = 0.98). When we convert the data to monthly variables and plot those, we observe a stronger correlation for both precipitation (r = 0.83) and temperature (r = 1) (Supplementary Fig. S10). ERA5 Land slightly overestimates the lower precipitations and underestimates the higher precipitation values (Supplementary Fig. S10). The potential underestimation of precipitation could lead to compensation through other factors (e.g. precipitation and/or temperature gradients in SnowModel). Our assessment of the elevation-binned modelled mass balance versus glaciological and geodetic approaches (Fig. 9) suggests that it may be difficult to identify how model performance is affected, as modelled mass balance at Peyto both overestimates and underestimates glaciological and geodetic approaches, depending on the elevation band being considered.

Based on Mann–Kendall trend test, ERA5 Land winter precipitation does not show any significant trend, but summer temperature shows an increasing trend (Supplementary Fig. S11) for Place Glacier. The ERA5 Land summer temperature and winter precipitation data at Peyto Glacier show no significant trend (Supplementary Fig. S12).

The El Niño Southern Oscillation (ENSO) atmospheric and oceanic circulation system plays an important role in the mass balance of western Canadian glaciers (e.g. Hodge and others, 1998; Moore and others, 2009). We observe strong mass loss in glaciological and modelled mass-balance records at both glaciers during the 1998 El Niño event (Fig. 6). Glacier mass gains or near-balanced conditions are also observed and modelled in strong La Niña years (1996, 1999, 2000). This implies that the ERA5 Land forcing data accurately captures the large-scale changes in temperature and precipitation over western Canada during El Niño (warmer, drier) and La Niña (colder, wetter) events.

6. Conclusions

Our analysis uses high-resolution climate reanalysis data to simulate glacier mass change at two long-term glacier monitoring sites using the physically-based SnowModel coupled with an ice dynamics model. SnowModel is tuned to match geodetic observations of mass change collected from repeat LiDAR surveys, and modelled seasonal and annual mass balances are compared with available in-situ mass-balance data. While cumulative modelled, geodetic and glaciological mass balances at both sites are similar, our analysis suggests that several years of reported glaciological mass balances require careful attention. ERA5 Land captures regional temperature and precipitation variability, and when coupled with SnowModel is able to reconstruct mass-balance series where no other data are available. This also strengthens the reliability of ERA5 data to estimate changes of glacier mass balance within all of western Canada. In the absence of stake records, our LiDAR-based geodetic estimate of 2020 mass balance could be used to fill up the gap in glaciological mass-balance record. Modelled glacier mass loss at Place Glacier has considerably increased after the year 2000 (from − 0.94 m w.e. a⁻¹ in 1981–2000 to − 1.31 m w.e. a⁻¹ in 2000-2020). For Peyto Glacier, modelled mass losses have remained almost similar in the period before (− 0.90 m w.e. a⁻¹ in 1982–2000) and after (− 0.93 m w.e. a⁻¹ in 2000–2020) the year 2000. This is commensurate with the accelerating trend of temperature increase after 2000 based on ERA5 Land data, especially for Place Glacier. Our coupled mass-balance and ice dynamics model shows that ice dynamics play a significant role in the mass balance of glaciers and the effect is perhaps more pronounced for glaciers with higher rates of mass flux.