1. Introduction
Over the past decade, large-scale simulations of glacier mass changes predominantly use positive degree-day (PDD) models (e.g. Radić and Hock, Reference Radić and Hock2014; Clarke and others, Reference Clarke, Jarosch, Anslow, Radić and Menounos2015; Huss and Hock, Reference Huss and Hock2015; Maussion and others, Reference Maussion2019; Rounce and others, Reference Rounce2023). These models, which estimate glacier melt based on the empirical relationship between temperature and melt, are relatively straightforward to apply but have several known limitations: they are often not transferable across different spatial and temporal scales (MacDougall and others, Reference MacDougall, Wheler and Flowers2011; Rounce and others, Reference Rounce, Khurana, Short, Hock, Shean and Brinkerhoff2020), exhibit excessive sensitivity to temperature shifts (Pellicciotti and others, Reference Pellicciotti, Brock, Strasser, Burlando, Funk and Corripio2005; Huss and others, Reference Huss, Farinotti, Bauder and Funk2008) and lack the capacity to provide process-based insights into glacier responses to climate change. As a result, large-scale projections of glacier mass change still contain significant uncertainties, particularly at the scale of individual glaciers and glacierized watersheds, where the impacts of deglaciation on local freshwater supplies are most critical (Anderson and Radić, Reference Anderson and Radić2020). While advances in nonlinear melt models, such as deep learning approaches (Bolibar and others, Reference Bolibar, Rabatel, Gouttevin, Galiez, Condom and Sauquet2020), overcome some of the limitations of PDD models, they also fall short in delivering process-based understanding of glacier-climate interactions. Although progress continues to be made in enhancing ice-dynamics modeling with more physics-based approaches (Rounce and others, Reference Rounce2023), the melt modeling component continues to rely on empirical methods, highlighting the need for physics-based melt models, such as surface energy balance (SEB) models.
Using SEB models in large-scale glacier simulations presents substantial challenges due to the limited availability of crucial climatic inputs, including shortwave and longwave radiation, wind speed, specific humidity and surface characteristics like albedo and roughness length, all of which need to be resolved accurately at the glacier scale. Consequently, large-scale SEB models are rare and often depend on semi-empirical methods that require calibration with mass-balance observations (Giesen and Oerlemans, Reference Giesen and Oerlemans2013; Huss and Hock, Reference Huss and Hock2018; Gunnarsson and others, Reference Gunnarsson, Gardarsson and Pálsson2023) or statistical downscaling techniques, which are challenging to generalize across varying timescales (Noël and others, Reference Noël, Berg, Lhermitte, Wouters, Schaffer and Van den Broeke2018; Shannon and others, Reference Shannon2019). Even regional applications of complex SEB models, such as COSIPY (Sauter and others, Reference Sauter, Arndt and Schneider2020), require calibration of parameters such as albedo with mass-balance observations to achieve comparable performance to PDD models (Temme and others, Reference Temme2023). Until recently, mass-balance observations for individual glaciers covered fewer than 0.1% of glaciers worldwide. The availability of glacier-specific remote sensing estimates from 2000 onward (Hugonnet and others, Reference Hugonnet2021) has facilitated broader model calibration over multi-year timescales (e.g. Compagno and others, Reference Compagno2022; Rounce and others, Reference Rounce2023; Schuster and others, Reference Schuster, Rounce and Maussion2023; Temme and others, Reference Temme2023; Hanus and others, Reference Hanus, Schuster, Burek, Maussion, Wada and Viviroli2024; Zekollari and others, Reference Zekollari2024). While these estimates align well with in situ mass-balance observations for long-term trends, significant discrepancies persist at annual and sub-annual timescales (Hugonnet and others, Reference Hugonnet2021).
Recent advances in climate reanalysis products and the growing availability of remote sensing data have created new opportunities to address some limitations in large-scale SEB modeling. For instance, the ECMWF Reanalysis v5 (ERA5) reanalysis dataset (Hersbach and others, Reference Hersbach2020) reasonably represents key meteorological inputs for SEB models at glaciers in western Canada (Draeger and others, Reference Draeger, Radić, White and Tessema2024). However, while ERA5 provides surface albedo data, it poorly represents albedo at the glacier scale (Draeger and others, Reference Draeger, Radić, White and Tessema2024). In contrast, Moderate Resolution Imaging Spectroradiometer (MODIS) data (Hall and Riggs, Reference Hall and Riggs2021) enables regional mapping of snow and ice albedos and are commonly used in SEB studies on ice caps (Gascoin and others, Reference Gascoin2017; Gunnarsson and others, Reference Gunnarsson, Gardarsson and Pálsson2023) and ice sheets (Stroeve and others, Reference Stroeve2005; Box and others, Reference Box, Fettweis, Stroeve, Tedesco, Hall and Steffen2012). However, for mountain glaciers, MODIS-derived albedo presents larger errors and more frequent data gaps (Davaze and others, Reference Davaze2018; Xiao and others, Reference Xiao, Ke, Fan, Shen and Cai2022), and its direct application in regional-scale SEB models remains unexplored.
This study builds on recent advances by leveraging ERA5 climatic data and MODIS-derived albedo to investigate physics-based melt modeling with minimal downscaling. Specifically, we evaluate an elevation-dependent SEB model, with varying calibration levels, for glacier mass-balance simulations and compare its performance against the widely used elevation-dependent PDD model. A central focus is on representing glacier albedo in the SEB model, contrasting a simple empirical albedo model with a novel machine learning approach trained on MODIS data. With the limitations of mass-balance observations at sub-annual scales such as those provided in the Hugonnet and others Reference Hugonnet(2021) database in mind, we investigate the impact of using seasonal versus annual mass-balance data for model calibration. We evaluate model performance across a diverse set of Canadian glaciers with in situ mass-balance observations, spanning varied climatic regimes—from small temperate valley glaciers in the Canadian Rockies and the maritime climate of the Coast Mountains to large polythermal glaciers in the Arctic.
2. Data and methods
2.1. Study glaciers and their mass-balance observations
We select 23 glaciers across Canada (Fig. 1, Table 1), each with a minimum of 3 years of elevation-based mass-balance records since 1979 (WGMS, 2022). On average, these glaciers have 13 years of mass-balance data post-1979, with individual records ranging from 3 to 39 years. The glaciers are clustered into four geographic subregions: West Coast, Rockies, East Coast and Arctic (Fig. 1). Devon Ice Cap North West in the Arctic is the largest glacier in our sample at 765.44 km2, though it represents an outlier in size, with the next largest, Bridge Glacier, covering 81.98 km2. The smallest glacier, Abraham Glacier, located on Canada’s East Coast, has an area of 0.503 km2. Of the 23 glaciers 17 also have seasonal mass-balance measurements and mass-balance profiles, allowing for more precise calibration of the elevation-dependent mass-balance model and a more comprehensive comparison of model outputs than possible with annual mass-balance data alone. Further details for each glacier are available in Table S3.
Figure 1. Geographic distribution of glaciers included in this study, organized by subregion. Marker colors and shapes represent subregions: red circles for West Coast, green squares for Rockies, orange diamonds for East Coast and blue stars for Arctic glaciers. The number of glaciers in each subregion is indicated by corresponding colored labels. Red dots denote locations of on-glacier meteorological observation sites used for ERA5 data comparison. Map projection: WGS84.
Table 1. Details of glaciers used in the study
2.2. Climate data and downscaling
We use surface-level, hourly ERA5 reanalysis climate data at a 
$30 \times 30\,\mathrm{km}$ grid resolution (Hersbach and others, Reference Hersbach2020) from 1979 to 2019 in this study. For each glacier, we select the reanalysis data from the nearest ERA5 grid point to its central coordinates. Depending on the variable, hourly data are either averaged or summed to produce daily values. The PDD model utilizes mean daily near-surface temperature, while the SEB model relies on mean daily values for near-surface temperature, dewpoint temperature, wind speed, surface pressure, as well as daily totals for incoming shortwave and longwave radiation, and precipitation. Both daily precipitation and mean near-surface temperature are used in the accumulation model. To demonstrate the differences in these variables among the four subregions, we show their time series as the area-weighted mean across all glaciers per subregion (Figure S1).
Among the meteorological inputs used in this study, we apply downscaling only to temperature by lapse-rate correcting ERA5 temperature to account for elevation differences between the ERA5 grid cell and the glacier site. Following Draeger and others Reference Draeger, Radić, White and Tessema(2024), we calculate the lapse rate using multi-level monthly ERA5 temperature data between the surface and 400 hPa, spanning 17 levels. This monthly lapse rate is linearly interpolated to obtain a daily lapse rate, lr (∘C m−1), which we then apply to correct the temperature at each elevation band of each glacier:
\begin{equation}
T_h = T_\mathrm{ERA} + lr(h - h_\mathrm{ERA}) ,
\end{equation} where Th (∘C) is the corrected near-surface temperature for the given elevation, h (m a.s.l.), while 
$T_\mathrm{ERA}$ (∘C) is the ERA5 temperature at the elevation of the grid cell, 
$h_\mathrm{ERA}$ (m a.s.l.).
2.3. Albedo data
We use the MOD10A1 MODIS/Terra Snow Cover Daily L3 Global 500 m SIN Grid, Version 6 dataset (Hall and Riggs, Reference Hall and Riggs2021) to derive daily snow albedo values for our 23 study glaciers during the summer months (June–September) from 2000 to 2019. The MODIS Terra satellite, operating in a near-polar, sun-synchronous orbit, captures up to nine daily observations of high Arctic glaciers (Williamson and others, Reference Williamson, Copland, Thomson and Burgess2020). Data processing utilizes the Normalized Difference Snow Index (Hall and others, Reference Hall, Salomonson and Riggs1995), which calculates albedo based on the visible-to-infrared radiation ratio and identifies snow-covered surfaces. The algorithm dynamically assigns snow albedo values to snow/ice-covered grid cells, leaving non-snow-covered cells with null values. Daily snow albedo data are provided on a 500×500 m grid, with the highest quality observation each day selected based on illumination, satellite angle, cloud cover and fractional snow cover (Hall and Riggs, Reference Hall and Riggs2021).
To ensure data quality, we filter the MODIS albedo values by removing physically unrealistic observations, specifically those >0.99 and <0.05 (Williamson and others, Reference Williamson, Copland, Thomson and Burgess2020). As our model operates at an elevation-band level, we convert gridded daily MODIS data into an average daily albedo for each elevation band of each glacier (detailed in the next section). Due to cloud cover, the dataset is often incomplete; not all days contain observations, and available data may cover only select elevation bands. For some bands, no albedo observations are available. Across the study glaciers, albedo observations are present in at least one elevation band 17%–77% of the time, with an average observation coverage of 44%.
2.3.1. Comparison between MODIS-derived and in situ albedo
MODIS-derived albedo often differs from in situ albedo measurements at glacier sites (Davaze and others, Reference Davaze2018), but such comparisons are limited due to the scarcity of field observations. To estimate these differences for our study glaciers, we compiled in situ albedo observations from four glacier sites in or near our study area (see Fig. 2). Radiometer data, which record incoming and reflected shortwave radiation, are available for two study glaciers in the Rockies, covering multiple summer seasons (Fitzpatrick and others, Reference Fitzpatrick, Radić and Menounos2017; Reference Fitzpatrick, Radić and Menounos2019). For the Arctic subregion, radiometer observations are accessible from nearby glaciers, specifically Belcher Glacier (two summer seasons; Sharp and others Reference Sharp, Duncan, Wolken and Boon2011) and Prince of Wales Glacier (one summer season; S. Marshall, personal communication, May 2020). We derive daily average albedo values from these radiometer data by calculating the ratio of integrated hourly reflected to incoming shortwave radiation, excluding hours with solar zenith angles exceeding 65∘ to minimize error and bias (Fitzpatrick and others, Reference Fitzpatrick, Radić and Menounos2017). Additional details on the radiometer measurements are provided in Table S1. In processing MODIS data to derive an average albedo for each glacier elevation band, we select the MODIS-derived albedo from the band nearest to each glacier’s site with daily in situ albedo values (Fig. 2, Table 2). The comparison shows that MODIS-derived albedo consistently underestimates albedo at all four sites, with a mean bias error (MBE) of −0.14 across sites.
Figure 2. Daily MODIS-derived albedo (α) compared to in situ albedo from radiometer measurements at on-glacier sites over multiple summer seasons for four glaciers (listed in the legend; see text for site details).
Table 2. Performance metrics (root-mean-square error (RMSE), mean bias error (MBE) and correlation coefficient (r)) evaluating daily MODIS-derived albedo against in situ albedo observations at four glaciers across multiple summer seasons. The number of observations per glacier is listed in the final column. The last row summarizes metrics calculated across all sites combined.
2.4. Models
To simulate mass changes of our study glaciers, we use an elevation-dependent mass-balance model, commonly used for large-scale simulations (e.g. Radić and Hock, Reference Radić and Hock2014; Huss and Hock, Reference Huss and Hock2015; Rounce and others, Reference Rounce2023). Our focus is on modeling the climatic or reference mass balance, which assumes constant glacier area and hypsometry over time, set to the values from the initial year of the mass-balance observational record for each glacier. Consequently, the model does not account for ice flow contributions to glacier geometry changes or the feedback between changes in glacier thickness, area and mass balance. For consistency with the mass-balance data, the hypsometry is sourced from the WGMS dataset, where glacier area and mass balance are provided for 100 m elevation bands.
For each elevation band, we calculate daily specific mass balance (in m water equivalent, w.e.), as the difference between modeled daily accumulation and ablation. The glacier-wide mass balance is then computed as the area-weighted mean of specific mass balance across all elevation bands for each glacier. Summer mass balance is derived as the integrated mass balance over the summer season, defined as the period from 01 May to 01 October, which roughly corresponds to the reported observation dates for the study glaciers (WGMS, 2022). Winter mass balance is calculated over the remaining months of the year. We maintain a constant definition for summer and winter periods to allow for fair comparison between glaciers and over the model period, as well as to maintain consistency with other regional/global approaches to mass-balance modeling (e.g. Radić and Hock, Reference Radić and Hock2014; Shannon and others, Reference Shannon2019).
2.4.1. Accumulation model
We adopt a relatively simple accumulation model, commonly used in previous large-scale glacier mass-balance simulations (e.g. Radić and Hock, Reference Radić and Hock2014; Rounce and others, Reference Rounce2023), where daily accumulation, c (in m w.e.), at elevation h (m a.s.l.), the center elevation of each elevation band, is calculated as:
\begin{equation}
c(h) = \begin{cases}
0 & T_{h} \gt 2\\
k_p P_\mathrm{ERA} \times (1+ d_p(h - h_\mathrm{ERA})) & T_{h} \leq 2 ,
\end{cases}
\end{equation} where Th (∘C) is the lapse-rate corrected ERA5 temperature for the given elevation band, 
$P_\mathrm{ERA}$ (m d−1) is ERA5 precipitation and kp and dp are calibration-dependent parameters: a multiplicative precipitation factor (kp) to correct the ERA5 precipitation, which often underestimates precipitation over mountainous regions (Tarek and others, Reference Tarek, Brissette and Arsenault2020), and a precipitation gradient (dp, 100 m−1) to account for orographic precipitation increases not captured by ERA5’s original spatial resolution. Despite not being mass conserving, this correction is often necessary to achieve a better match between modeled and observed profiles of winter glacier mass balance (Radić and Hock, Reference Radić and Hock2014). Both parameters are calibrated against the mass-balance observations for each of our study glaciers; further details are provided later in this section.
2.4.2. Melt model: PDD approach
We use a PDD melt model following Radić and Hock Reference Radić and Hock(2014), which assumes an empirical relationship between temperature and glacier’s daily melt, a (in m w.e.), given by:
\begin{equation}
a(h) = f_\mathrm{ice/snow} \int \mathrm{max}(T_h,0) \mathrm{d}t,
\end{equation} where 
$f_\mathrm{ice/snow}$ is the calibration-dependent melt factor for ice or snow (m w.e. d−1 
$^\circ\mathrm{C}^{-1}$); both are constant values that are individually calibrated for each glacier. The melt factor for snow (
$f_\mathrm{snow}$) is used above the assumed firn-line elevation regardless of snow cover, while below the firn-line elevation we apply the melt factor for ice (
$f_\mathrm{ice}$) when snow depth is zero and 
$f_\mathrm{snow}$ when the snow depth is greater than zero. Following Huss and Hock Reference Huss and Hock(2015) and Rounce and others Reference Rounce, Khurana, Short, Hock, Shean and Brinkerhoff(2020), the firn-line elevation is set to the elevation where the 5 year running mean of each elevation band’s annual mass-balance changes from positive to negative.
2.4.3. Melt model: SEB approach
We use a simple SEB model, keeping only the key contributors to melt energy as is done in Draeger and others Reference Draeger, Radić, White and Tessema(2024). Daily surface melt, a (in m w.e.) is calculated for each elevation band as:
\begin{equation}
a(h) = \frac{Q_\mathrm{melt}}{L_f\rho_w} ,
\end{equation} where 
$Q_\mathrm{melt}$ is the total daily energy available for surface melt (J m−2), Lf is the latent heat of fusion (
$334000 \text{J kg}^{-1}$) and ρw is the density of water (
$1000\,\text{kg m}^{-3}$).
To calculate the energy available for melt, we use:
\begin{equation}
Q_\mathrm{melt} = Q_S ^\downarrow (1 - \alpha) + Q_L ^\downarrow - Q_L ^\uparrow + Q_H + Q_E .
\end{equation} 
$Q_S ^\downarrow$ represents incoming shortwave radiation, α is surface albedo, simulated by albedo models described below, 
$Q_L ^\downarrow$ and 
$Q_L ^\uparrow$ represent incoming and outgoing longwave radiation, respectively, and QH and QE denote the sensible and latent heat fluxes, respectively. All fluxes are integrated over one day, thus are expressed in J m−2. We intentionally neglect (i) empirical correction schemes applied to shortwave radiation fluxes, such as separating direct and diffuse components, (ii) conductive heat flux into the ground and (iii) energy supplied by rain. These terms generally rely on calibration-dependent parameterizations and typically play a minimal role in determining seasonal melt, depending on glacier type and location (Hock, Reference Hock2005; Fitzpatrick and others, Reference Fitzpatrick, Radić and Menounos2017; Reference Fitzpatrick, Radić and Menounos2019). Our rationale for excluding these terms is to simplify the SEB model while retaining the most impactful terms for surface melt.
$Q_L ^\downarrow$ is taken directly from ERA5 without modifications. 
$Q_L ^\uparrow$ is computed using the Stefan–Boltzmann law for blackbody radiation, assuming the ice/snow surface remains at a constant temperature of 0∘C year-round (i.e. the ice/snow surface is consistently at the melting point) (Hock, Reference Hock2005; Cuffey and Paterson, Reference Cuffey and Paterson2010; Fitzpatrick and others, Reference Fitzpatrick, Radić and Menounos2017). We also tested an approach where the ice/snow surface temperature is set to the lapse-rate-corrected 2 m air temperature (Th), when it falls below freezing, but observed negligible effects on the modeled annual 
$Q_\mathrm{melt}$ relative to the former approach.
The turbulent fluxes, QH and QE, are calculated using the bulk aerodynamic method, specifically following the C log approach described by Fitzpatrick and others Reference Fitzpatrick, Radić and Menounos(2017):
\begin{equation}
\begin{aligned}
Q_H(h) &= \frac{c_p \, \rho_a \, C_H \, p \, u \, (T_h - T_0)}{p_0} \\
Q_E(h) &= \frac{0.622 \, l_v \, \rho_a \, C_E \, p \, u \,(e_h - e_0)}{p_0} ,
\end{aligned}
\end{equation} where cp is the specific heat capacity of air (
$1005\,\text{J kg}^{-1} \text{K}^{-1}$), p is the surface air pressure (hPa), p 0 is the air pressure at sea level (
$1013.25\,\text{hPa}$), ρa is the density of air at sea level (
$1.29\,\text{kg m}^{-3}$), lv is the latent heat of vaporization of water (2.514 
$\times 10^6\,\text{J kg}^{-1}$) and T 0 is the temperature of the ice/snow surface (assumed to be at 0∘C over a summer season).
The mean daily wind speed, u (m s−1), is taken directly from ERA5 data (
$u = u_{\text{ERA}}$) and is also tested with a calibration-dependent multiplicative wind correction factor, ku, such that 
$u = k_u \times u_{\text{ERA}}$. The vapor pressures e 0 and eh represent the surface vapor pressure (assumed to be 
$6.1078\,\text{hPa}$ at 0∘C) and the vapor pressure 2 m above the surface, respectively. The latter is calculated using ERA5 data via the August–Roche–Magnus formula (Alduchov and Eskridge, Reference Alduchov and Eskridge1996), with relative humidity derived from ERA5 temperature and dewpoint temperature.
The bulk transfer coefficients for heat (CH) and moisture (CE) are dimensionless and are parameterized following Fitzpatrick and others Reference Fitzpatrick, Radić and Menounos(2017), without stability corrections. These coefficients depend on the roughness lengths for momentum (
$z_{0,v}$), temperature (
$z_{0,T}$) and humidity (
$z_{0,q}$) at the glacier surface. For this study, we use roughness length values of 
$z_{0,v} = 10^{-3}\,\text{m}$, 
$z_{0,T} = 10^{-5}\,\text{m}$ and 
$z_{0,q} = 10^{-5}\,\text{m}$, consistent with typical observations at mid-latitude glaciers (Hock, Reference Hock2005; Fitzpatrick and others, Reference Fitzpatrick, Radić and Menounos2019).
2.4.4. Albedo models
During the winter months, the albedo for each glacier is assigned its specific snow albedo value, calculated as the mean annual maximum MODIS albedo over the summer months during the 2000–2019 period (see Table S2 for individual glacier values). Across the 23 study glaciers, the average MODIS-derived albedo values are 0.85 for snow, 0.52 for firn and 0.23 for ice. The evolution of summer albedo in each glacier’s elevation bands is simulated using two distinct modeling approaches. The first approach, denoted as 
$\alpha_\mathrm{Dec}$, is based on the empirical model of Hirose and Marshall Reference Hirose and Marshall(2013), developed for a glacier in the Canadian Rockies. This model employs a logarithmic expression to represent the decay of albedo following snowfall events. The daily summer albedo, 
$\alpha_\mathrm{Dec}(t)$, is simulated for each elevation band as:
\begin{equation}
\alpha_\mathrm{Dec}(t) = \alpha_\mathrm{snow} - A \ln\sum \text{PDD}(t),
\end{equation} where 
$\alpha_\mathrm{snow}$ is the albedo of fresh snow, A is a decay factor calibrated individually for each glacier and PDD(t) represents the cumulative positive degree-days (in ∘C) up to time t. The sum of PDD(t) is calculated daily, starting from the onset of the melt season and resetting to zero after each fresh snowfall. Hirose and Marshall Reference Hirose and Marshall(2013) calibrate this model against in situ albedo observations and find that A = 0.05668. We optimize A using the Nelder–Mead method (Nelder and Mead, Reference Nelder and Mead1965) to minimize the root-mean-square error (RMSE) between the MODIS-derived and modeled daily albedo over the 2000–2019 period. Optimization of A is performed individually for each glacier, with RMSE calculated using two metrics: the daily glacier-wide albedo and the mean summer albedo for each elevation band, both weighted equally. The optimized A values for each glacier are provided in Table S2. To ensure physical realism, we constrain albedo in the model to avoid values below a minimum threshold: above the firn-line, this threshold is set to firn albedo, and below the firn-line, it is set to ice albedo. For most glaciers, we find A to be below 0.5, except for the Arctic glaciers (see Table S2).
The second albedo model, denoted as 
$\alpha_\mathrm{MLP}$, uses a neural network multilayer perceptron (MLP) to simulate daily summer albedo. The MLP model is trained on MODIS-derived albedo values for each elevation band, aggregated across glaciers within the same subregion rather than on individual glaciers. This approach assumes that a model trained on a broader dataset within a subregion is more robust, benefiting from the increased amount of training data, as machine learning models typically improve with larger datasets (Hsieh, Reference Hsieh2009). The MLP is a feed-forward neural network consisting of an input layer, one or more hidden layers and an output layer, with neurons in each layer connected via nonlinear activation functions. During training, the model’s parameters (weights) are optimized to minimize the mean square error between the predicted and observed MODIS-derived daily albedo values for each elevation band.
For each subregion, the MLP model is trained using data from any study glacier with available MODIS albedo observations in a given elevation band, implicitly incorporating elevation information. When multiple glaciers have observations in the same elevation band, their data are combined for model training. The model inputs include the day of the year, as well as ERA5 temperature and precipitation data for the target glacier, spanning the current day and the preceding 6 days. Each input variable is normalized to a range of 0–1 based on its minimum and maximum values. The paired input data and MODIS albedo observations are randomly shuffled and split into training and testing datasets, ensuring data from each glacier within the subregion are represented in both. We use 80% of the data for model training and validation and reserve 20% for testing. Within the training set, 80% is allocated for training and 20% for validation.
For the model runs, we test various configurations with 1–2 hidden layers, each containing between 1 and 15 neurons. The model with the best performance (i.e. lowest mean square error) is selected. This process is repeated 20 times using different initial weights, creating an ensemble of the 20 best-performing MLP models for each elevation band and subregion. We experimented with several model architectures, including training separate models for each glacier and using different combinations of input parameters. Ultimately, we chose the described model setup as it minimized the need for interpolation and effectively reproduced the MODIS observations. The average RMSE across all glaciers between modeled and MODIS-derived albedo is 0.11.
The model is implemented in Python (Van Rossum and Drake, Reference Van Rossum and Drake2009), utilizing the Scikit-Learn library (Pedregosa and others, Reference Pedregosa2011). We use the rectified linear activation unit as the activation function and optimize the model with the Adaptive Moment Estimation (Adam) method. Adam combines the advantages of both moving average gradients and adaptive learning rates to update the neural network weights, effectively minimizing the objective function (Kingma and Ba, Reference Kingma and Ba2014; Sun and others, Reference Sun, Cao, Zhu and Zhao2019). For further details on these algorithms, refer to Kingma and Ba Reference Kingma and Ba2014.
After calibrating or training both albedo models, the resulting albedo values are incorporated into the SEB model. Given the known biases between MODIS-derived and in situ albedo, as well as discrepancies in ERA5 incoming shortwave radiation compared to observations at glacier surfaces (Draeger and others, Reference Draeger, Radić, White and Tessema2024), we introduce an additional calibration parameter as an additive correction to the modeled albedo. The final albedo is then computed as:
\begin{equation}
\alpha = \alpha_\mathrm{mod} + \Delta\alpha ,
\end{equation} where 
$\alpha_\mathrm{mod}$ represents the albedo value from either the empirical model (
$\alpha_\mathrm{Dec}$) or the machine learning model (
$\alpha_\mathrm{MLP}$), and 
$\Delta\alpha$ is the additive correction factor. The melt model that uses the empirical representation of albedo is denoted SEB
$_\mathrm{Dec}$ and the one that uses the MLP albedo representation is denoted SEB
$_\mathrm{MLP}$.
Since neither of the two albedo models fully capture short-term albedo variations due to fresh snowfall events during the summer, we account for this process separately for each study glacier. On days when daily snowfall, as estimated from the accumulation model (Eqn (2)), exceeds 0.5 mm in a given elevation band, we override the modeled albedo by setting it to the snow albedo value derived from MODIS data. This adjustment is applied regardless of the albedo value predicted by the empirical or machine learning models. Finally, to ensure physical plausibility, the albedo time series for each elevation band is constrained such that the values do not exceed 0.9 or fall below 0.1.
2.4.5. Model calibration
As described in previous sections, several calibration-dependent parameters are optimized during the calibration of the mass-balance model. These include the precipitation correction parameters (kp and dp) in the accumulation model and the melt factors (
$f_\mathrm{snow}$ and 
$f_\mathrm{ice}$) in the PDD model. While the SEB model was originally designed to be calibration-independent, we also test model runs that incorporate a calibration-dependent wind correction (ku) and albedo correction (
$\Delta\alpha$). Since the melt and accumulation models must operate together to generate mass-balance estimates, it is not possible to optimize them independently. All calibration procedures aim to minimize the RMSE between the model output and observed mass-balance data (WGMS) for each study glacier. The RMSE is computed from both the seasonal mass-balance time series and the average seasonal mass-balance profile, with both metrics given equal weight. For comparison, we also perform a calibration where the RMSE is calculated from annual mass-balance observations instead of seasonal ones.
The calibration parameters are optimized using the Nelder–Mead optimization scheme (Nelder and Mead, Reference Nelder and Mead1965). The bounds for the precipitation correction parameter (dp) range from 0 to 3% per 100 m elevation increase, and for the precipitation factor (kp), they range from 1 to 3, in line with the mass-balance model used by Radić and Hock Reference Radić and Hock(2014) for Canadian glaciers. Following Hock Reference Hock(2003), the bounds for the snow melt factor (
$f_\mathrm{snow}$) are set between 1 and 9 mm w.e. d−1 
$^\circ\mathrm{C}^{-1}$, and for the ice melt factor (
$f_\mathrm{ice}$), between 3 and 11 mm w.e. d−1 
$^\circ\mathrm{C}^{-1}$, with the condition that 
$f_\mathrm{ice}$ must be greater than 
$f_\mathrm{snow}$. The bounds for the albedo correction parameter (
$\Delta\alpha$) are between 0 and 0.3, with the upper bound determined by the difference between the maximum physically reasonable albedo for snow (0.9) and the lowest MODIS-derived snow albedo (0.6 for the Nordic glacier). The multiplicative wind correction factor (ku) is bound between 0 and 4, based on the work of Shannon and others Reference Shannon(2019).
2.4.6. Model intercomparison and sensitivity tests
Our primary goal is to evaluate and compare the performance of the models in replicating the observed (WGMS) interannual variability of glacier-wide summer, winter and net mass balance, as well as the average mass-balance profile over the observational period. Given the limited sample size, we do not apply calibration-validation methods (e.g. leave-one-out cross-validation); thus, our assessment focuses on how well each model fits the available data. We begin by analyzing model runs where the SEB model parameters remain uncalibrated. We then proceed to evaluate runs with an optimized albedo correction (
$\Delta \alpha$), and finally, we assess runs that include both optimized albedo and wind corrections.
We also evaluate the impact of temperature perturbations on the simulated mass balance over the 1980–2019 period. To do this, we apply uniform temperature adjustments to the ERA5 input temperatures for each glacier and across all models, including the albedo models. The temperature perturbations are set at 
$-2^\circ$C, 
$-1^\circ$C, 
$+1^\circ$C and 
$+2^\circ$C for the entire 40 year period. The impact is quantified by calculating the mass-balance sensitivity (
$\frac{\partial B}{\partial T}$) for each glacier, defined as the difference in the mean modeled mass balance between the adjusted temperature scenario and the original scenario, divided by the temperature perturbation.
3. Results
3.1. Calibrated parameters
Here we provide a summary of the calibrated parameter values across the study glaciers, with detailed values available in Table S2. The precipitation correction parameters (kp and dp) are calibrated separately for each of the three melt models (SEB
$_\mathrm{MLP}$, SEB
$_\mathrm{Dec}$ and PDD). The multiplicative factor kp had mean values (with relative variance in parentheses) of 1.76 (19%) for the SEB
$_\mathrm{MLP}$ model, 1.94 (16%) for the SEB
$_\mathrm{Dec}$ model and 1.87 (27%) for the PDD model. The precipitation gradient dp averaged 1.11 (48%) 100 m−1 for the SEB
$_\mathrm{MLP}$ model, 0.90 (72%) 100 m−1 for the SEB
$_\mathrm{Dec}$ model and 0.85 (75%) 100 m−1 for the PDD model.
In the PDD model, we observe substantial variability in the melt factors (
$f_\mathrm{snow}$ and 
$f_\mathrm{ice}$) across the study region, though some consistency is noted within certain subregions. On the West Coast, glaciers generally displayed minimal differences between 
$f_\mathrm{snow}$ and 
$f_\mathrm{ice}$, with a mean 
$f_\mathrm{snow}$ of 3.68 mm w.e. d−1 
$^\circ\mathrm{C}^{-1}$ (5.8% variance) and a mean 
$f_\mathrm{ice}$ of 3.84 mm w.e. d−1 
$^\circ\mathrm{C}^{-1}$ (1.4% variance). In contrast, East Coast glaciers exhibited greater variability and larger differences between melt factors, with a mean 
$f_\mathrm{snow}$ of 1.75 mm w.e. d−1 
$^\circ\mathrm{C}^{-1}$ (16.7% variance) and a mean 
$f_\mathrm{ice}$ of 4.84 mm w.e. d−1 
$^\circ\mathrm{C}^{-1}$ (69.0% variance). For Arctic glaciers, where only two glaciers are included in the analysis, the 
$f_\mathrm{ice}$ values spanned a wide range, from 3.77 to 10.79 mm w.e. d−1 
$^\circ\mathrm{C}^{-1}$.
In the SEB models, the albedo correction factor (
$\Delta\alpha$) ranged from 0 to 0.3 across the study glaciers. The mean 
$\Delta\alpha$ is 0.11 for the MLP model and 0.19 for the logistic decay model. Analyzing the subregional differences, the mean 
$\Delta\alpha$ values (with relative variance) for the MLP model are as follows: 0.06 (141.4%) in the Arctic, 0.22 (66.7%) on the East Coast, 0.09 (74.6%) on the West Coast and 0.08 (85.1%) in the Rockies. For the logistic decay model, the mean values are 0.18 (89.1%) in the Arctic, 0.28 (11.1%) on the East Coast, 0.16 (51.1%) on the West Coast and 0.18 (36.1%) in the Rockies.
3.2. Modeled versus observed albedo
To assess the performance of the two albedo models in simulating seasonal albedo evolution, we compare time series of modeled and MODIS-derived daily glacier-wide albedo, averaged over all available MODIS data days and then averaged across glaciers in each subregion (Fig. 3, a–d). Both the MLP and logistic decay models successfully capture the gradual decrease in albedo during the early melt season, showing the best agreement with MODIS-derived albedo for Arctic and East Coast glaciers. However, for glaciers in western Canada, the MLP model more accurately follows the seasonal pattern of MODIS-derived albedo compared to the logistic decay model, which exhibits a faster decline in albedo at the onset of the melt season. In contrast, the MLP model tends to produce higher albedo values in late summer relative to MODIS-derived albedo.
Figure 3. (a–d) Time series of daily average albedo, shown as a 5 day running mean, averaged across glaciers within each subregion. Values are derived from MODIS observations (black), the machine learning (MLP) model (blue) and the logistic decay model (green), calculated as the mean of all days with MODIS observations during the 2000–2019 period. (e–f) Average summer (June–August) albedo versus glacier elevation for each subregion, derived from MODIS and the two models, calculated as the mean of all days with MODIS observations for each elevation band of the study glaciers during the same period.
We further compare modeled and MODIS-derived average albedo across elevation bands in each subregion, averaged over all available MODIS observation days (Fig. 3, e–h). Both the modeled and MODIS-derived albedo generally display the expected increase with elevation. However, relatively low absolute values are observed at the highest elevation bands in both MODIS data and the models. Except for Arctic glaciers, the average summer albedo for the highest elevations rarely exceeds 0.5, whereas a range of 0.6–0.7 would be more physically realistic for the accumulation zones of glaciers without debris cover (Hock, Reference Hock2005; Cuffey and Paterson, Reference Cuffey and Paterson2010). More detailed results of the elevation profiles for individual glaciers are provided in Figure S2, where these subregional patterns are consistently confirmed.
The relatively low MODIS-derived albedo in glacier accumulation zones is further supported by comparisons with in situ albedo observations from multiple glacier sites (Fig. 4). We closely examine the time series of daily albedo values from both in situ measurements and MODIS observations, alongside outputs from the two albedo models, for two glaciers in the Rockies where in situ albedo data are available. The modeled albedo includes a calibrated correction factor (
$\Delta \alpha$) and adjustments for fresh snowfall events. Both albedo models successfully capture the gradual decline in albedo during the early melt season and the sharp increase as the melt season ends. Although neither model fully replicates short-term albedo variations due to fresh snowfall, they generally indicate the presence of these events, albeit with some discrepancies in timing compared to observations. Across all sites, the RMSE between modeled and in situ albedo time series ranges from 0.12 to 0.16, suggesting a reasonable fit. However, MODIS-derived albedo is observed to be consistently lower than in situ measurements, particularly in the accumulation zone (e.g. Station 3 on Conrad Glacier), where MODIS-derived albedo is on average 0.2 units lower than the in situ observations over the observational period.
Figure 4. Time series of daily albedo at multiple sites on two glaciers (Conrad and Nordic) in the Rockies, shown for various summer seasons. The albedo values are derived from MODIS observations (black), automatic weather station (AWS) measurements (gray), the logistic decay model (green) and the machine learning (MLP) model (blue). Both the logistic decay and MLP models include the calibrated albedo correction factor, 
$\Delta \alpha$ (refer to text for calibration details). MODIS albedo values are taken from the grid cell nearest to the AWS location, while modeled albedo values correspond to the elevation band closest to each AWS site. Data from Conrad Glacier are labeled as: Conrad Stn 1 (a and b), Conrad Stn 2 (c) and Conrad Stn 3 (d), while data from Nordic Glacier are labeled as Nordic Stn (e). Further site details are available in Table S1.
Finally, we examine whether the calibrated albedo correction factor (
$\Delta \alpha$) is associated with the tendency of both MODIS data and the models to underestimate albedo values at higher glacier elevations. Our analysis shows a strong, statistically significant negative correlation between the magnitude of the calibrated albedo correction and the maximum MODIS-derived summer albedo across the elevation range—typically in the accumulation zone for most glaciers (Fig. 5). This finding suggests that lower albedo values in the accumulation zone necessitate a larger correction to align the modeled mass balance with observations during the calibration process.
Figure 5. Albedo correction (
$\Delta \alpha$) from the SEB
$_\mathrm{MLP}$ model versus the maximum MODIS-derived albedo along the glacier elevation range for each study glacier. Points are color-coded by subregion, with the linear regression slope and Pearson’s correlation coefficient shown in the top right corner.
3.3. Model intercomparison
We compare the performance of the PDD and SEB modeling approaches using box-and-whisker plots that illustrate the distribution of normalized root-mean-square error (NRMSE) across the study glaciers (Fig. 6). The NRMSE is calculated by dividing the RMSE, which measures the discrepancy between observed and modeled mass-balance time series as well as the averaged mass-balance profiles (using the same RMSE metric applied during model calibration), by the range of observed values (i.e. the difference between the maximum and minimum values). This normalization provides a standardized error metric, enabling direct comparison across glaciers of different sizes. We evaluate NRMSE separately for annual (net), summer and winter mass balances. The annual NRMSE is calculated for all 23 study glaciers, while the seasonal NRMSE is evaluated for a subset of 17 glaciers.
Figure 6. Box plots of normalized root-mean-square error (NRMSE) for annual (a), summer (b) and winter (c) mass balance across 23 glaciers, comparing SEB
$_\mathrm{MLP}$ (blue), SEB
$_\mathrm{Dec}$ (green) and PDD (pink) models under different calibration scenarios. Scenarios include no calibration (SEB
$_\mathrm{MLP}$, SEB
$_\mathrm{Dec}$), precipitation-only (SEB
$_\mathrm{MLP}$ P, SEB
$_\mathrm{Dec}$ P), precipitation and albedo (SEB
$_\mathrm{MLP}$ P,α, SEB
$_\mathrm{Dec}$ P,α) and precipitation, albedo and wind corrections (SEB
$_\mathrm{MLP}$ P,α,u, SEB
$_\mathrm{Dec}$ P,α,u). PDD model results are shown for melt factor calibration (PDD) and calibration of melt factor and precipitation (PDD P). Solid lines represent median NRMSE, dashed lines indicate mean NRMSE and median values are labeled above the whiskers. Boxes span the interquartile range, with whiskers extending to the 5th and 95th percentiles.
First, we examine the results from the uncalibrated SEB models and the calibrated PDD model, where melt factors are individually optimized for each glacier (Fig. 6). In these runs, no precipitation corrections are applied (i.e. 
$d_p = 0$ and 
$k_p = 1$). The performance for winter mass balance is similarly poor across all three models, with median NRMSE values of 51.5%, 54.0% and 42.4% for the SEB
$_\mathrm{MLP}$, SEB
$_\mathrm{Dec}$ and PDD models, respectively. Statistical analysis using an independent two-sample t test reveals no significant differences among these distributions, as shown by the box plots. In contrast, the PDD model shows significantly better performance than the SEB models for summer and net mass balance. For annual (net) mass balance, the median NRMSE of the PDD model is 17.8%, substantially lower than the 47.4% and 72.7% recorded by the SEB
$_\mathrm{MLP}$ and SEB
$_\mathrm{Dec}$ models, respectively.
Next, we consider the model intercomparison using SEB models calibrated with varying sets of calibration-dependent parameters (Fig. 6). In these runs, all models include calibrated precipitation correction parameters (kp and dp). Compared to the previous results without precipitation correction, incorporating only these calibrated parameters leads to lower median NRMSE values for winter mass balance across the study glaciers. For summer mass balance, we observe a median NRMSE reduction of 9.2% for the SEB
$_\mathrm{MLP}$ model, 2.2% for the SEB
$_\mathrm{Dec}$ model and 0.5% for the PDD model. Notably, the NRMSE distributions for summer and annual mass balance in the SEB models, when only precipitation factors are calibrated, show statistically significant differences from those of the PDD model, as indicated by an independent two-sample t test.
Incorporating the albedo correction into the SEB models results in the most significant improvement in summer mass-balance performance, reducing the median NRMSE from 24.6% to 13.0% for the SEB
$_\mathrm{MLP}$ model and from 44.9% to 13.1% for the SEB
$_\mathrm{Dec}$ model. Similarly, for annual mass balance, the albedo correction decreases the median NRMSE from 18.5% to 11.4% in the SEB
$_\mathrm{MLP}$ model and from 38.5% to 10.7% in the SEB
$_\mathrm{Dec}$ model. With this correction, the NRMSE distributions of the SEB models are no longer significantly different from those of the PDD model. Adding a wind correction alongside precipitation and albedo corrections provides only minor, statistically insignificant improvements, lowering the median summer NRMSE by 1.3% for the SEB
$_\mathrm{MLP}$ model and by 1.6% for the SEB
$_\mathrm{Dec}$ model.
Across all subregions, we observe a consistent reduction in NRMSE values with the introduction of calibrated parameters. In each subregion, incorporating the albedo correction significantly reduces NRMSE for both the SEB
$_\mathrm{MLP}$ and SEB
$_\mathrm{Dec}$ models, while adding the wind correction yields more modest improvements. The East Coast glaciers typically have the highest mean annual NRMSE values, starting at 30.3% for the uncalibrated SEB
$_\mathrm{MLP}$ model, which decreases to 25.6% after applying the albedo correction and results in a mean annual NRMSE of 21.8% for the PDD model. In contrast, the Rockies glaciers exhibit the lowest NRMSE values, with the SEB
$_\mathrm{MLP}$ model achieving a mean annual NRMSE of 16.2%, which further decreases to 11.9% with the albedo correction applied, while the PDD model has a comparable mean annual NRMSE of 12.3%.
At the scale of individual glaciers, we find that the PDD and SEB
$_\mathrm{MLP}$ models, with calibrated precipitation and albedo corrections, produce similarly good fits to the mass-balance observations (annual comparison of NRMSE is shown in Figure S3). Specifically, the NRMSE for summer mass balance across the 17 glaciers is strongly and significantly correlated between the two models (r = 0.88; Fig. 7). This correlation in NRMSE is stronger for the SEB
$_\mathrm{MLP}$ model than for the SEB
$_\mathrm{Dec}$ model, which exhibits substantially larger NRMSE values than the PDD model for several glaciers. Overall, the correlation of NRMSE between the SEB and PDD models is stronger for winter (r in the range 0.86–0.89) than for summer mass balance (r in the range 0.70–0.88). Notably, the SEB
$_\mathrm{MLP}$ model consistently achieves substantially lower NRMSE values than the PDD model for winter mass balance across multiple glaciers (Figure S4).
Figure 7. Normalized root-mean-square error (NRMSE) for summer mass balance: SEB
$_\mathrm{MLP}$ versus PDD model (a) and SEB
$_\mathrm{Dec}$ versus PDD model (b). SEB models include calibrated precipitation and albedo corrections, while the PDD model uses calibrated melt factors and precipitation. Results are shown for 17 glaciers with seasonal observations, colored by subregion.
Finally, we examine the effects of calibrating the models using annual versus seasonal mass-balance observations (Fig. 8). For consistency, we limit this analysis to the 17 glaciers with available seasonal mass-balance data. We compare model runs for the SEB models with calibrated precipitation and albedo corrections, and the PDD model that has calibrated precipitation and melt factors. For all three models (SEB
$_\mathrm{MLP}$, SEB
$_\mathrm{Dec}$ and PDD), calibration based on annual mass-balance observations leads to lower median NRMSE values for annual mass balance; however, the differences are not statistically significant. In contrast, for seasonal mass balance, median NRMSE values are either similar or higher, with a notably larger spread in NRMSE across the 17 glaciers. When calibrated solely on annual mass balance, the models for many glaciers tend to overestimate both summer melt and winter accumulation, balancing these errors to achieve a better fit to the net mass balance (see Figure S5). In contrast, calibration using seasonal mass-balance observations effectively reduces biases in winter and summer mass balances, leading to more accurate seasonal simulations.
Figure 8. Box plots of normalized root-mean-square error (NRMSE) for annual, summer and winter mass balance across 17 glaciers with seasonal observations. Results compare SEB
$_\mathrm{MLP}$ (blue), SEB
$_\mathrm{Dec}$ (green) and PDD (pink) models under two calibration scenarios: seasonal (‘S’) and annual (‘A’) mass-balance observations. SEB models include calibrated precipitation factors and albedo corrections, while the PDD model includes calibrated melt and precipitation factors. Solid lines indicate median NRMSE, dashed lines show mean NRMSE and median values are labeled above the whiskers. Boxes represent the interquartile range, with whiskers extending to the 5th and 95th percentiles.
3.4. Modeled versus observed in situ mass balance
In addition to evaluating NRMSE, we assess how well the models replicate the in situ (stake) mass-balance observations at the study glaciers. Although these in situ measurements are not directly used for model calibration, their aggregates, such as glacier-wide mass-balance time series and average mass-balance profiles, are considered. For each glacier, we evaluate the ‘goodness of fit’ between modeled and observed (WGMS) seasonal and annual in situ mass balances, using the modeled value from the elevation band nearest to the stake measurement. The performance metrics include RMSE, MBE, Pearson correlation coefficient (r) and Nash–Sutcliffe efficiency (NSE) coefficient, as detailed in Table 3. For conciseness, we present results only for the SEB model runs incorporating calibrated precipitation and albedo corrections.
Table 3. Performance metrics for simulations of annual (bold), summer and winter mass balance from the PDD model and two SEB models, summarized for each subregion. Results are shown for SEB models with calibrated precipitation and albedo corrections and the PDD model with calibrated precipitation and melt factors. Metrics include root-mean-square error (RMSE), mean bias error (MBE), correlation coefficient (r) and Nash–Sutcliffe efficiency (NSE). Reference data are in situ specific mass-balance measurements from WGMS (WGMS, 2022). The number of observations per subregion is listed in the final column
Across all subregions and models, there is a statistically significant correlation (at the 95% confidence level) between modeled and observed in situ annual mass balance (Table 3). For summer mass balance, the correlation coefficients are consistently higher than for winter mass balance, with all correlations statistically significant (p-value <0.05), except for the winter mass balance in the Arctic subregion. Despite this, the RMSE for winter mass balance is generally smaller than for summer mass balance, likely due to the lower interannual variability in winter mass balance. When comparing RMSE and correlation coefficients, the PDD model gives a slightly better fit to observations than the SEB
$_\mathrm{MLP}$ model in the East Coast subregion. In contrast, the SEB
$_\mathrm{MLP}$ model either performs similarly to or marginally better than the PDD model in the other subregions. A detailed comparison of the modeled versus observed in situ mass balance for each study glacier is provided in Figures S7 and S8.
3.5. Sensitivity to temperature perturbations
To quantify mass-balance sensitivity to temperature changes, the models were first run with unperturbed ERA5 data over the 1980–2019 period, followed by runs in which ERA5 temperatures were systematically perturbed. When forced with unperturbed ERA5 data, the PDD and SEB
$_\mathrm{MLP}$ models produced closely resembling mass-balance reconstructions across all study glaciers. Both models captured a consistent negative long-term trend over the 40 year period, with the SEB
$_\mathrm{MLP}$ model generally exhibiting a slower decline in mass balance compared to the PDD model (Figure S6). To assess whether these differences in long-term trends were statistically significant, we performed an analysis of covariance (ANCOVA). The results indicate that none of the glacier-specific trends differ significantly at the 95% confidence level (noting that ANCOVA considers both trend slopes and intercepts).
In addition to mass-balance trends, both models yielded similar interannual variability in mass balance. On average, the variance of the PDD model’s mass-balance time series was higher than that of the SEB model, with a mean variance ratio (
$\sigma^2_\mathrm{PDD}/\sigma^2_\mathrm{SEB}$) of 1.49 ± 0.36 across all study glaciers. However, when averaged across glaciers, differences in trends and interannual variability among SEB
$_\mathrm{MLP}$, SEB
$_\mathrm{Dec}$ and PDD models were negligible (Fig. 9).
Figure 9. Modeled summer, winter and annual (net) mass balance from 1980 to 2019, averaged across all study glaciers, comparing SEB
$_\mathrm{MLP}$, SEB
$_\mathrm{Dec}$ and PDD models. (a) Results using unperturbed ERA5 temperature data. (b) Results with a +1∘C temperature perturbation. SEB model results include calibrated precipitation and albedo correction factors, while the PDD model incorporates calibrated precipitation and melt factors.
When a +1∘C temperature perturbation was introduced to ERA5 data, the PDD model exhibited a substantially more negative mass-balance response than the SEB model (Fig. 9). Across all models, winter mass balance was less affected by the temperature perturbation compared to summer mass balance. These results are further summarized in box plots for all temperature perturbation scenarios, illustrating the sensitivity of annual, winter and summer mass balance separately (Fig. 10).
Figure 10. Box plots of the temperature sensitivity of mass balance (
$\frac{\partial B}{\partial T}$) for annual (a), summer (b) and winter (c) mass balance across the 23 study glaciers under different temperature perturbations. The sensitivity 
$\frac{\partial B}{\partial T}$ is determined by comparing the mass-balance time series from the original model runs (SEB
$_{MLP}$, SEB
$_{Dec}$ and PDD) to those from model runs with applied temperature perturbations. Results are shown for the SEB models with calibrated precipitation and albedo correction factors, and for the PDD model with calibrated precipitation and melt factors. Solid horizontal lines indicate the median values, while dashed lines represent the mean values of 
$\frac{\partial B}{\partial T}$. The median value is labeled above the top whisker of each distribution. Box edges denote the first and third quartiles, and whiskers represent the 5th and 95th percentiles.
The results for winter mass-balance sensitivity to a 2∘C temperature increase show that the SEB models have a median 
$\frac{\partial B}{\partial T}$ value of 
$-0.03\,\mathrm{m\,w.e.\,yr}^{-1}$ 
$^\circ\mathrm{C}^{-1}$, whereas the PDD model shows a higher sensitivity with a median value of 
$-0.11\,\mathrm{m\,w.e.\,yr}^{-1}$ 
$^\circ\mathrm{C}^{-1}$ (Fig. 10). In contrast, the summer mass balance exhibits greater sensitivity to temperature increases across all models. For the SEB
$_\mathrm{MLP}$ model, the median 
$\frac{\partial B}{\partial T}$ values range from −0.20 to 
$-0.21\,\mathrm{m\,w.e.\,yr}^{-1}$ 
$^\circ\mathrm{C}^{-1}$, while for the SEB
$_{Dec}$ model, they range from −0.28 to 
$-0.29\,\mathrm{m\,w.e.\,yr}^{-1}$ 
$^\circ\mathrm{C}^{-1}$. The PDD model, however, shows significantly higher sensitivity, with median values ranging from −0.50 to 
$-0.53\,\mathrm{m\,w.e.\,yr}^{-1}$ 
$^\circ\mathrm{C}^{-1}$ across the temperature perturbations. In all cases, the PDD model’s sensitivity is significantly different from that of the SEB models (as confirmed by an independent two-sample t test) and exhibits a wider spread of 
$\frac{\partial B}{\partial T}$ values compared to the SEB models.
At the scale of individual glaciers, the PDD model consistently produces higher annual mass-balance sensitivities to a +1∘C temperature change compared to the SEB models (Fig. 11). The annual mass-balance sensitivities across the 23 glaciers are significantly correlated between the PDD and SEB
$_\mathrm{MLP}$ models (r = 0.54) and between the PDD and SEB
$_\mathrm{Dec}$ models (r = 0.61). While the SEB models maintains relatively stable sensitivities across different magnitudes of temperature perturbations, the PDD model exhibits substantial variability. Specifically, the lowest sensitivity for all models, but particularly the PDD model, occurs under a 
$-2^\circ$C perturbation and increases progressively with temperature perturbations up to 
$+2^\circ$C. This effect is most pronounced for glaciers with a large difference between the snow melt factor (
$f_\mathrm{snow}$) and ice melt factor (
$f_\mathrm{ice}$), where the PDD model shows the most significant shifts in mass-balance sensitivity, especially for winter mass balance.
Figure 11. Mass-balance sensitivity to a 1
$^\circ\mathrm{C}^{-1}$ temperature increase (
$\frac{\partial B}{\partial T}$) derived from the PDD model compared to the SEB
$_\mathrm{MLP}$ model (a) and the SEB
$_\mathrm{Dec}$ model (b). Results are shown for all 23 study glaciers, with points colored by subregion. SEB models incorporate calibrated precipitation factors and albedo corrections, while the PDD model includes calibrated melt and precipitation factors.
The mass-balance sensitivity to temperature change varies significantly across the four subregions. The Arctic region shows the lowest mean annual 
$\frac{\partial B}{\partial T}$ values for all three models (PDD, SEB
$_\mathrm{MLP}$ and SEB
$_\mathrm{Dec}$), with values ranging from −0.02 to 
$-0.11\,\mathrm{m\,w.e.\,yr}^{-1}$ 
$^\circ\mathrm{C}^{-1}$ for the SEB
$_\mathrm{MLP}$ model and −0.03 to 
$-0.18\,\mathrm{m\,w.e.\,yr}^{-1}$ 
$^\circ\mathrm{C}^{-1}$ for the PDD model. In contrast, the West Coast glaciers exhibits the highest mean annual mass-balance sensitivity, with 
$\frac{\partial B}{\partial T}$ values ranging from −0.21 to 
$-0.61\,\mathrm{m\,w.e.\,yr}^{-1}$ 
$^\circ\mathrm{C}^{-1}$ for the SEB
$_\mathrm{MLP}$ model and from −0.24 to 
$-0.76\,\mathrm{m\,w.e.\,yr}^{-1}$ 
$^\circ\mathrm{C}^{-1}$ for the PDD model.
4. Discussion
Our results demonstrate that the PDD model, calibrated with glacier-specific melt factors, initially outperforms the SEB models when no additional parameters are adjusted. Introducing calibrated precipitation correction factors (dp and kp) to the models leads to significant improvements in simulating winter and net mass balance. However, even with these corrections, the SEB models underperform relative to the PDD model. It is only after incorporating a calibrated albedo correction in the SEB models that their performance matches or slightly exceeds that of the PDD model. The median NRMSE values for net mass balance are 12.3% for SEB
$_\mathrm{MLP}$ and 11.6% for SEB
$_\mathrm{Dec}$, compared to 12.8% for the PDD model. These findings highlight the pivotal role of accurately simulating glacier albedo, as net shortwave radiation is the dominant contributor to melt energy at our sites, accounting for ∼70% of the total melt energy.
A key focus of this study is the use of MODIS-derived albedo data to train the MLP model and calibrate the logistic decay model. However, for many study glaciers, the MODIS-derived albedo exhibits relatively low summer values, including unrealistic estimates (below 0.5) in glacier accumulation zones. These underestimated albedo values, when used as inputs in the SEB models, lead to an overestimation of melt and consequently higher NRMSE values. The calibrated albedo correction mitigates this bias by increasing the albedo by an average of 0.11, significantly improving the model’s fit to observed mass balance. Introducing a calibrated wind speed correction, in addition to the precipitation and albedo corrections, yields only marginal enhancements in model performance.
Our analysis also reveals a consistent tendency for MODIS to underestimate albedo when compared with in situ albedo observations from several glacier sites. This finding is consistent with other studies that have compared MODIS-derived albedo with in situ snow albedo measurements across various mountainous regions, where results often vary widely (e.g. Stroeve and others, Reference Stroeve, Box and Haran2006; Sorman and others, Reference Sorman, Akyurek, Sensoy, Sorman and Tekeli2007; Williamson and others, Reference Williamson, Copland and Hik2016; Calleja and others, Reference Calleja, Corbea-Pérez, Fernández, Recondo, Peón and de Pablo2019). These discrepancies highlight the site-specific performance of MODIS in assessing snow and ice albedo. According to the MODIS user guide, snow albedo estimates are generally within 10% of in situ observations under optimal conditions. However, the guide also notes that errors can be significantly higher over steep mountainous terrain (Hall and Riggs, Reference Hall and Riggs2021). Our findings underscore the need for applying albedo corrections to MODIS data or to models trained on this data, particularly given the variability in MODIS performance across different regions and local conditions. We observe that the greater the underestimation of summer albedo in the glacier accumulation zone, the larger the required albedo correction.
Additional factors to consider when using MODIS-derived albedo data include atmospheric conditions, shadowing effects and surface type changes. Shadowing from local terrain and variations in surface type primarily impacts grid cells near glacier edges, introducing potential errors in the observed MODIS albedo (Williamson and others, Reference Williamson, Copland, Thomson and Burgess2020). Furthermore, extensive cloud cover often coincides with snowfall events, preventing albedo measurements from being recorded. Consequently, MODIS observations may miss the rapid, short-term increases in albedo over summer following fresh snowfall on exposed glacier surfaces. To address this discrepancy, we adjust the albedo models by increasing albedo to fresh snow values on days with modeled snowfall during summer. This adjustment has a minimal impact on the simulated mass balance, resulting in an average difference of <2% in the NRMSE of summer mass balance across our study glaciers. While the frequency of fresh snowfall events has been shown to significantly influence summer albedo and, consequently, summer mass balance (Oerlemans and Klok, Reference Oerlemans and Klok2004), our findings suggest that accurately resolving the average summer albedo across the glacier elevation range is more critical than capturing short-term albedo variability from individual snowfall events.
While the albedo data simulated by both models—MLP and logistic decay—ultimately result in similar simulations of summer mass balance, there are critical differences between these models that can affect their suitability for regional and global-scale applications. Notably, without the calibrated albedo correction, the SEB
$_\mathrm{MLP}$ model showed a significantly better fit to mass-balance observations compared to the SEB
$_\mathrm{Dec}$ model, with a median summer NRMSE of 24.6%, compared to 44.9% for the SEB
$_\mathrm{Dec}$. One key difference lies in their sensitivity to temperature: the logistic decay albedo model relies solely on temperature as an input, making it more responsive to temperature perturbations, as evidenced by our mass-balance sensitivity tests. In contrast, the MLP albedo model incorporates additional features and is designed to be elevation-specific within subregions rather than glacier-specific, potentially enhancing its spatial transferability. Moreover, while we utilized a neural network (MLP) approach for regional albedo simulation in this study, this is just one of many possible machine learning strategies. More advanced architectures, such as recurrent neural networks or long short-term memory models, have yet to be explored and could offer further improvements in simulating albedo on regional scales.
Since our SEB model is elevation-dependent, it does not explicitly account for the effects of glacier slope and aspect on incoming shortwave radiation. In contrast, a fully distributed SEB model (e.g. Hock and Holmgren, Reference Hock and Holmgren2005; Hill and others, Reference Hill, Dow, Bash and Copland2020) that incorporates these topographic factors would likely offer a more accurate depiction of the net radiative balance but would introduce substantially greater model complexity. Instead, the calibrated albedo correction we apply to help mitigate biases in ERA5 incoming radiation (Draeger and others, Reference Draeger, Radić, White and Tessema2024) and inaccuracies in MODIS-derived albedo data also indirectly compensates for the combined effects of glacier slope and aspect on shortwave radiation.
Recognizing the limitations of sub-annual mass-balance observations, particularly from remote sensing data such as those provided in the Hugonnet and others Reference Hugonnet(2021) database, we explore the effects of using seasonal versus annual mass-balance data for model calibration. Our analysis shows that calibrating both the PDD and SEB models using annual mass-balance observations results in a similarly good fit to annual mass balance as when seasonal observations are used. However, relying solely on annual data for calibration introduces notable biases in the seasonal components; overestimated (or underestimated) winter accumulation is compensated by overestimated (or underestimated) summer melt, effectively masking the errors in net mass-balance simulation. In contrast, calibrating the models with seasonal mass-balance data helps to optimize the accumulation and melt processes separately, leading to more accurate simulations of both winter and summer mass balance. These findings underscore the critical importance of in situ seasonal mass-balance measurements, especially as the field increasingly relies on remote sensing data (Compagno and others, Reference Compagno2022; Rounce and others, Reference Rounce2023; Hanus and others, Reference Hanus, Schuster, Burek, Maussion, Wada and Viviroli2024; Zekollari and others, Reference Zekollari2024), which currently have limited accuracy at annual and sub-annual scales. Accurate seasonal calibration ensures that simple models are better tuned to represent the distinct processes driving glacier mass changes, a crucial consideration for reliable mass-balance projections in future climate scenarios. Our results highlight the importance of seasonal-scale calibration while recognizing the value of glacier-specific calibration based on long-term trends, as implemented in Rounce and others Reference Rounce(2023) and Zekollari and others Reference Zekollari(2024). These findings highlight the need for further improvements in global-scale glacier modeling.
We attribute the strong performance of the SEB models (with calibrated precipitation and albedo corrections) in fitting mass-balance observations to the use of ERA5 reanalysis data, which reasonably captures local meteorological conditions at multiple glacier sites in western Canada (Draeger and others, Reference Draeger, Radić, White and Tessema2024). The only bias correction applied to the ERA5 input data is a lapse rate adjustment for 2 m air temperature. This correction proved to be crucial: omitting it increased the median NRMSE by 24.6% for the SEB
$_\mathrm{MLP}$ model and by 49.9% for the PDD model using the same melt factors. The PDD model shows a greater sensitivity to the lapse rate correction because temperature directly drives its melt calculations. In contrast, the SEB model is less affected by temperature biases; near-surface temperature primarily influences turbulent heat fluxes, which contribute <30% to seasonal melt energy at the study sites. Additionally, temperature had a minimal effect on simulated albedo evolution in the MLP model, resulting in an average difference of 0.08 in summer albedo across the study glaciers.
The higher sensitivity of the PDD model to temperature changes is further evident in our sensitivity tests, where temperature is uniformly adjusted across the entire modelled period. The PDD model consistently exhibits a greater response, with median summer 
$\frac{\partial B}{\partial T}$ values around 
$-0.52\,\mathrm{m\,w.e.\,yr}^{-1}$ 
$^\circ\mathrm{C}^{-1}$, compared to 
$\sim-0.20\,\mathrm{m\,w.e.\,yr}^{-1}$ 
$^\circ\mathrm{C}^{-1}$ for the SEB
$_\mathrm{MLP}$ model. The spatial pattern, with the lowest sensitivities observed in Arctic glaciers and the highest in West Coast glaciers, aligns with previous findings showing greater sensitivities for maritime glaciers and lower sensitivities for Arctic glaciers (Radić and Hock, Reference Radić and Hock2011). Additionally, our sensitivity values fall within the range reported for static mass-balance sensitivities across various glaciers globally (e.g. de Woul and Hock, Reference de Woul and Hock2005). The temperature sensitivity of the PDD model is also influenced by the difference between 
$f_\mathrm{snow}$ and 
$f_\mathrm{ice}$; glaciers with larger differences between these melt factors show a more pronounced, nonuniform response to varying magnitudes of temperature changes. Given that many large-scale glacier models primarily rely on projections of temperature and precipitation (e.g. Marzeion and others, Reference Marzeion2020; Rounce and others, Reference Rounce2023), the differing temperature sensitivities of the PDD and SEB models could lead to substantially different predictions of future glacier mass changes in climate projections where temperature is the only changing variable. Unlike the PDD model, which responds solely to temperature, the SEB model accounts for multiple climatic variables, providing additional degrees of freedom in its response. Consequently, it exhibits lower sensitivity to temperature alone. While these findings do not invalidate the use of PDD models for glacier melt projections, they highlight the limitations of relying exclusively on temperature as the sole meteorological input for melt calculations.
Given our focus on evaluating a relatively simple SEB model, we deliberately excluded components that would require more complex and calibration-intensive frameworks. This choice was guided by the principle that increasing model complexity does not necessarily yield better regional-scale glacier simulations (Temme and others, Reference Temme2023) and can complicate the diagnostic assessment of model performance. However, some of our simplified assumptions may not be universally applicable across glaciers from different climatic regimes. For instance, processes such as meltwater refreezing within the snowpack and heat flux into the glacier surface are more critical for SEB and mass balance in Arctic glaciers compared to those in mid-latitudes (Woodward and others, Reference Woodward, Sharp and Arendt1997; Hock, Reference Hock2005). Similarly, rain heat flux has been shown to contribute significantly (up to 20%) to daily melt energy during extreme rainfall events on some glaciers in the Rockies (Fitzpatrick and others, Reference Fitzpatrick, Radić and Menounos2017). However, the uncertainty in the models used to estimate rain heat flux remains relatively high (Fitzpatrick and others, Reference Fitzpatrick, Radić and Menounos2017).
Despite its simplicity, our SEB model—enhanced with precipitation and albedo corrections—achieves a ‘goodness of fit’ comparable to or better than those reported in studies using more complex, heavily parameterized SEB models, such as Shannon and others Reference Shannon(2019), which are part of global-scale simulations. Specifically, our model consistently produces lower RMSE values (0.20–0.64
$\,\mathrm{m\,w.e.\,yr}^{-1}$ compared to 0.96–1.73
$\,\mathrm{m\,w.e.\,yr}^{-1}$ in Shannon and others Reference Shannon(2019)), higher correlation coefficients, lower MBEs and NSE coefficients closer to one—metrics evaluated using approximately the same in situ mass-balance observations across Canadian glaciers. While differences in model complexity, climate forcing and downscaling approaches exist between our study and Shannon and others Reference Shannon(2019), the superior fit of our model to in situ mass-balance observations underscores its potential for large-scale mass-balance modeling.
5. Conclusions
Large-scale glacier mass-balance models commonly rely on PDD melt models, which, while widely used, have well-recognized limitations. This study investigated the ability of a relatively simple, elevation-dependent SEB model, requiring minimal downscaling of climate input data, to replicate seasonal and annual mass-balance observations. We compared mass-balance models utilizing SEB and PDD melt components against observations from 23 glaciers across Canada, spanning diverse climatic regimes. The models were driven by ERA5 data, with downscaling limited to lapse-rate corrections of 2 m air temperature for better local representation. We initially assessed SEB models without calibrated parameters and then introduced calibrated correction factors for precipitation, albedo and wind to enhance the alignment between modeled and observed mass balance. A central aspect of this study was the improved representation of glacier albedo within the SEB model, contrasting a simple empirical albedo approach (SEB
$_\mathrm{Dec}$) with a machine learning-based method trained on MODIS-derived albedo data (SEB
$_\mathrm{MLP}$). Our principal conclusions are:
• The accuracy of ERA5 input data, especially the lapse rate correction for 2 m temperature, was crucial for the performance of both SEB and PDD models. Initially, the PDD-based mass-balance model outperformed the SEB-based model when no calibrated parameters or only precipitation correction factors were applied, emphasizing the efficiency of glacier-specific melt factor calibration in PDD models.
• Incorporating a calibrated albedo correction significantly improved the performance of SEB models, bringing them on par with or marginally better than the PDD model in simulating glacier mass balance. The calibrated SEB
$_\mathrm{MLP}$ model outperformed the SEB
$_\mathrm{Dec}$ model in most scenarios, demonstrating the importance of advanced albedo parameterization.• MODIS-derived albedo data tended to underestimate summer albedo at the study glaciers, particularly in the accumulation zones, resulting in overestimated melt. Calibrated albedo corrections in the SEB models significantly improved their alignment with observed mass balance, underscoring the importance of calibrating albedo or, more broadly, net shortwave radiative fluxes to enhance model accuracy.
• Calibration using seasonal mass-balance data, as opposed to annual mass balance, led to more accurate simulations of summer and winter mass balance. This finding underscores the importance of using detailed seasonal data for calibrating mass-balance models to avoid compensating biases in accumulation and melt during the calibration process.
• Although the calibrated SEB and PDD models produced similar reconstructions of glacier mass balance over the 1980–2019 period in terms of long-term trends and interannual variability, the PDD model exhibited twice the sensitivity to temperature changes when temperature was the only altered climate variable. Since climate change affects multiple climatic variables, these findings highlight the limitations of the PDD model in relying solely on temperature as the primary input for melt calculations.
While a more complex, fully distributed SEB model incorporating glacier slope and aspect could yield more accurate simulations, the simpler, elevation-dependent SEB model assessed in this study provided a good balance between computational efficiency and model performance. This suggests that for regional applications, simplified models like the SEB model can be effective without unnecessary complexity, as long as key parameters such as albedo are properly calibrated.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0022143025100531.
Data availability statement
All data used in this study were sourced from publicly available data repositories (see Methods for further details).
Acknowledgements
We thank the two anonymous reviewers, as well as the editors David Rounce and Hester Jiskoot, for their helpful comments and feedback throughout the review process. We acknowledge funding from the Natural Sciences and Engineering Research Council of Canada Discovery Grant and the Polar Continental Shelf Program to Valentina Radić and the General Sir John Monash Foundation scholarship to Hannah Phelps.
Competing Interests
The authors declare that they have no conflicts of interest.
Open access
