1. Introduction
Glaciers in the European Alps are melting about twice as fast as the worldwide average (Zemp and others, Reference Zemp2019; Hugonnet and others, Reference Hugonnet2021). In the two extreme years 2022 and 2023, Swiss glaciers lost about 10% of their total ice volume (GLAMOS, 2023). Such drastic loss increases the interest in monitoring short-term glacier changes too, especially when considering the significance that glacial runoff has in mitigating the impacts of drought during such hot and dry summers (Zappa and Kan, Reference Zappa and Kan2007; Van Tiel and others, Reference Van Tiel, Van Loon, Seibert and Stahl2021; Pelto and others, Reference Pelto, Dryak, Pelto, Matthews and Perry2022, Ultee and others, Reference Ultee, Coats and Mackay2022, Cremona and others, Reference Cremona, Huss, Landmann, Borner and Farinotti2023). Numerous studies investigated long-term glacier evolution at the regional scale (Marzeion and others, Reference Marzeion, Hofer, Jarosch, Kaser and Mölg2012, Gabbi and others, Reference Gabbi, Carenzo, Pellicciotti, Bauder and Funk2014, Compagno and others, Reference Compagno, Eggs, Huss, Zekollari and Farinotti2021, Schuster and others, Reference Schuster, Rounce and Maussion2023), and operational monitoring set-ups to detect glacier changes at the daily to weekly scale are advancing (Fausto and others, Reference Fausto2021, Cremona and others, Reference Cremona, Huss, Landmann, Borner and Farinotti2023, Voordendag and others, Reference Voordendag, Prinz, Schuster and Kaser2023). The glacio-meteorological processes underlying these events are well understood too (Greuell and Böhm, Reference Greuell and Böhm1998, Greuell and Smeets, Reference Greuell and Smeets2001, Sauter and Galos, Reference Sauter and Galos2016, Mott and others, Reference Mott, Stiperski and Nicholson2020, Goger and others, Reference Goger, Stiperski, Nicholson and Sauter2022, Voordendag and others, Reference Voordendag2024), but comprehensive quantitative analyses remain limited. In this context, Landmann and others (Reference Landmann, Künsch, Huss, Ogier, Kalisch and Farinotti2021) developed CRAMPON (for ‘Cryospheric Monitoring and Prediction Online’), an operational glacier mass-balance model for calculating regional mass balance at the daily scale by relying on the infrastructure provided by the Open Global Glacier Model OGGM (Maussion and others, Reference Maussion2019). The application of CRAMPON has the potential to systematically quantify glacier melt, i.e. one of the major runoff contributors during hot and dry periods, and thus to enable the implementation of optimized water-allocation strategies (Anghileri and others, Reference Anghileri, Botter, Castelletti, Weigt and Burlando2018, Landmann and others, Reference Landmann, Künsch, Huss, Ogier, Kalisch and Farinotti2021).
One of the main challenges when aiming at such type of modelling is the lack of observations for model calibration (Hock and others, Reference Hock, Hutchings and Lehning2017). The Glacier Monitoring Switzerland (GLAMOS) program, for example, currently monitors the seasonal mass balance of about 20 out of the 1400 glaciers in Switzerland (Huss and others, Reference Huss, Dhulst and Bauder2015). This means that for most glaciers, the only available information to constrain model parameters is the long-term geodetic mass balance (Bauder and others, Reference Bauder, Funk and Huss2007, Fischer and others, Reference Fischer, Huss and Hoelzle2015, Huss and others, Reference Huss2021). Several studiesused such glacier-specific geodetic mass balances to calibrate glaciological models (Huss and others, Reference Huss, Bauder, Funk and Hock2008a, Zekollari and others, Reference Zekollari, Huss and Farinotti2019, Rounce and others, Reference Rounce2023, Schuster and others, Reference Schuster, Rounce and Maussion2023). However, this type of calibration is subject to parameter equifinality, meaning that different parameter combinations can result in the same ice volume change (Finger and others, Reference Finger, Pellicciotti, Konz, Rimkus and Burlando2011, Gabbi and others, Reference Gabbi, Carenzo, Pellicciotti, Bauder and Funk2014). Indeed, solely relying on geodetic mass balances is insufficient for separating the components of mass change—i.e. accumulation and melt. The equifinality problem becomes particularly important when aiming to model glacier changes at the sub-seasonal to daily time scale, and therefore, additional observations are required.
In recent decades, remote sensing has achieved significant advances in this respect, providing a wealth of new observations of the mountain cryosphere (e.g. Taylor and others, Reference Taylor, Quincey, Smith, Baumhoer, McMillan and Mansell2021). In the context of glacier monitoring, observations of the glacier surface are particularly valuable as direct correlations have been observed between glacier mass balance and surface characteristics such as surface albedo, snowline altitude or snow-covered area fraction (SCAF) (Brun and others, Reference Brun2015, Paul and others, Reference Paul, Winsvold, Kääb, Nagler and Schwaizer2016, Davaze and others, Reference Davaze2018, Naegeli and others, Reference Naegeli, Huss and Hoelzle2019, Rastner and others, Reference Rastner2019). Such remotely sensed observations of the glacier surface thus have considerable potential to constrain glacier mass balance at the sub-seasonal time scale (Dyurgerov and others, Reference Dyurgerov, Mikhalenko, Kunakhovitch, Ushnurtsev, Liu and Xie1994, Pelto and others, Reference Pelto, Kavanaugh and McNeil2013) and indeed, a number of studies combined SCAF observations and modelling to improve mass balance monitoring (Hock and others, Reference Hock, Kootstra and Reijmer2007, Hulth and others, Reference Hulth, Rolstad Denby and Hock2013, Huss and others, Reference Huss2013). More recently, Barandun and others (Reference Barandun2018) used transient snowline observations from satellite imagery and terrestrial automated cameras to constrain a glacier mass-balance model. In an application to Central Asia, the approach also demonstrated its efficacy in evaluating glacier mass balance at a regional scale (Barandun and others, Reference Barandun2021), even though the focus was on annual and not on sub-seasonal mass balance.
In this study, we investigate the potential of using SCAF observations from satellite imagery to constrain the temporal variability of modelled glacier mass balance at the sub-seasonal to daily time scale for ten selected glaciers in the Swiss Alps. The SCAFs are derived from Sentinel-2 observations over the period 2015–23 and are combined with long-term observations of glacier mass balance during model calibration. More specifically, we minimize the difference between modelled and observed SCAFs throughout the melt season while reproducing the measured long-term mean annual mass balance. The study then evaluates the benefit of using such a calibration strategy by benchmarking it against observed seasonal mass balances. Finally, we highlight and discuss the application potential of the approach on glaciers without in situ observations and emphasize its potential for providing more accurate estimates of glacier mass balance at the regional scale.
2. Study sites and data
The study sites comprise ten glaciers in the Swiss Alps (Fig. 1, Table 1). These glaciers are located in different climatological regions and range from 2 to 80 km2 in size. The glaciers are covered by the Sentinel-2 tiles R065-T32TSM, R065-T32TNS and R108-T32TLS (at 10 m resolution) from which up to nine SCAF observations are derived per melt season (June–September), yielding a total of 20–47 observations during the period 2015–23 (see Section 3.2 for derivation). The selected glaciers are also included in the long-term mass-balance monitoring programme of GLAMOS, providing seasonal glacier-wide mass balance inferred from a network of in situ point observations (Huss and others, Reference Huss, Dhulst and Bauder2015, GLAMOS, 2023). The glacier-wide mass balances serve to (i) derive mean annual mass balance over the period 2015–23, used during model calibration, and (ii) evaluate the present approach at the seasonal scale.
Figure 1. Study sites. (a) Overview of the study region. Glaciers are shown in violet, the ten selected study sites being highlighted in dark violet and outlines refer to the Swiss Glacier Inventory 2016 (Linsbauer and others, Reference Linsbauer2021) whereas Swiss boundaries are provided by Swisstopo (2024). The extent of the Sentinel-2 tiles R108-T32TLS, R065-T32TMS and R065-T32TNS are shown in red, yellow and blue, respectively. Panels (b–d) show close-ups of the four regions containing the study sites.
Table 1. Glacier area, elevation range and number of SCAF observations per year used during calibration (range and average), for the considered glaciers. Glacier area and elevation range refer to the Swiss Glacier Inventory 2016 (Linsbauer and others, Reference Linsbauer2021)
The meteorological inputs used to model glacier mass balance are the verified grid products (e.g. the analysis product incorporating information from station measurements) at a resolution of ca. 2 km of daily mean 2 m temperature (T), precipitation sum (P) and daily mean incoming shortwave radiation (G) (MeteoSwiss, 2018, 2021a, 2021b).
The glacier geometry is described by the glacier outlines from the Swiss Glacier Inventory 2016 (SGI2016 Linsbauer and others, Reference Linsbauer2021) and the underlying SwissALTI3D digital elevation model (DEM) (Swisstopo, 2020). The glacier geometry is assumed to remain constant over the study period 2015–23.
3. Methods
3.1. Glacier surface classification
To determine the temporal evolution of the SCAF for the ten target glaciers, we first derive glacier surface classes from selected Sentinel-2 multi-spectral scenes over the melt seasons 2015–23. Only Sentinel-2 scenes affected by less than 15% cloud cover according to the associated metadata are considered for the classification. Clouds from the Sentinel-2 scenes are masked using the s2cloudless python module (Sentinel Hub, 2024). Based on Sentinel-2 multi-spectral L1C data, spectral surface reflectance is calculated by combining the radiative transfer model 6S (Vermote and others, Reference Vermote, Tanré, Deuze, Herman and Morcette1997) driven by atmospheric parameters from the Copernicus Atmosphere Monitoring Service (CAMS) with the approach of Sandmeier and Itten (Reference Sandmeier and Itten1997). A series of thresholds is applied to the resulting spectral surface reflectances to classify snow, clean glacier ice and snow free debris cover on glaciated areas (Schwaizer and others, Reference Schwaizer, Nemec, Nagler, Mölg and Paul2023). On selected dates, based on a thorough visual check of different spectral band combinations, considering only scenes acquired between late July and early September, an additional class including firn, old snow and bright glacier ice is added. Clouds over glaciated areas and non-glaciated areas are masked.
3.2. Snow-covered area fraction
Starting from the glacier surface classification product, the SCAF observations are computed in three steps (Fig. 2). In the first step, the glacier surface classification products are cropped with the outline of every glacier (Fig. 2, step 1). In the second step, a reclassification is performed to make the observations and the modelled surface types comparable. In particular, the classes ‘cloud’, ‘debris’ and ‘firn/old snow/bright ice’ are reclassified into either ‘snow’ or ‘ice’ (Fig. 2, step 2). To perform the reclassification, we divide the DEM of the glacier into elevation bands of 20 m and for each band, the predominant class among ‘ice’ and ‘snow’ is determined. If the predominant class of the band is ‘ice’, the cells that belong to the classes ‘cloud’, ‘debris’ or ‘firn/old snow/bright ice’ in that band are assigned to the class ‘ice’. On the contrary, if the predominant class of the band is ‘snow’, cells that belong to the classes ‘cloud’, ‘debris’ or ‘firn/old snow/bright ice’ in that band are assigned to the class ‘snow’. In elevation bands without any cell of the classes ‘ice’ or ‘snow’, all grid cells are categorized based on the predominant class of the closest neighbouring bands that contain ice or snow. On average over the ten sites, the percentage of cells that were reclassified is 7%. Furthermore, to minimize uncertainties of the reclassification, we exclude any observation with less than 80% of the glacier area belonging to the classes ‘ice’ or ‘snow’. In the third step, we calculate the SCAF by summing the area of all cells of the class ‘snow’ and dividing it by the total glacier area (Fig. 2, step 3).
Figure 2. Deriving the observed snow-covered area fraction (SCAF). Step 1: The glacier facies product is cropped with the glacier outline. Step 2: Cells that do not belong to the class ‘snow’ or ‘ice’ (white areas within the glacier outline in this panel) are reclassified into one of the two classes. Step 3: The SCAF is calculated by dividing the snow-covered area by total glacier area.
3.3. Mass-balance modelling
To model glacier surface mass balance, we use CRAMPON, an operational mass-balance model for calculating glacier-specific snow accumulation and snow/ice melt at the daily scale (Landmann and others, Reference Landmann, Künsch, Huss, Ogier, Kalisch and Farinotti2021). The model relies on the infrastructure provided by the Open Global Glacier Model (Maussion and others, Reference Maussion2019), which discretizes a given glacier into elevation bands with a horizontal distance along the flowline of 20 m, resulting in a vertical resolution of about 2–20 m. At every elevation band, accumulation is modelled according to concepts proposed in Huss and others (Reference Huss, Farinotti, Bauder and Funk2008b):
\begin{equation}
{c_\mathrm{sfc}}(z) = {c_\mathrm{prec}} \cdot P \cdot \left[ 1 + (z - {z_\mathrm{ref}}) \cdot \frac{\partial P}{\partial z} \right],
\end{equation} where 
${c_\mathrm{sfc}}(z)$ is the snow accumulation 
$\left( \mathrm{m\ w.e.\ {d}^{-1}} \right)$ at elevation z, 
${c_\mathrm{prec}}$ is the precipitation correction factor, P is the sum of solid precipitation at the reference elevation 
${z_\mathrm{ref}}$ and 
$\frac{\partial P}{\partial z}$ is the solid precipitation lapse rate.
CRAMPON includes four melt models (Landmann and others, Reference Landmann, Künsch, Huss, Ogier, Kalisch and Farinotti2021); however, the Braithwaite (Braithwaite, Reference Braithwaite1995) and the Oerlemans (Oerlemans, Reference Oerlemans2001) model revealed to be less reliable in reproducing the course of the SCAF over the melt season, which caused convergence issues for the calibration approach. For this study, we thus only use the two melt models by Hock (Reference Hock1999) and Pellicciotti and others (Reference Pellicciotti, Brock, Strasser, Burlando, Funk and Corripio2005). The Hock model computes melt as:
\begin{equation}
{a_\mathrm{sfc}}(z) = (\mathrm{MF} + {a_\mathrm{snow/ice}} \cdot I_{\mathrm{pot}}(z)) \cdot \max(T(z) - T_{\mathrm{melt}}, 0),
\end{equation} where 
$\mathrm{MF}$ is the temperature melt factor 
$\left( \mathrm{m\,w.e.\,K}^{-1}\, \mathrm{d}^{-1} \right)$, 
${a_\mathrm{snow/ice}}$ are two radiation factors for snow and ice 
$\left( \mathrm{m\,w.e.\,m}^{2} \,\mathrm{d}^{-1}\, \mathrm{W}^{-1}\, \mathrm{K}^{-1} \right)$ and 
$I_{\mathrm{pot}}(z)$ is the potential clear-sky direct solar radiation 
$\left( \mathrm{W\,m}^{-2} \right)$. We assume 
$T_{\mathrm{melt}} = 0^{\circ}\mathrm{C}$ and a ratio 
${a}_\mathrm{\mathrm{snow}}/{a}_\mathrm{\mathrm{ice}}$ of 0.75 (Hock, Reference Hock1999, Farinotti and others, Reference Farinotti, Usselmann, Huss, Bauder and Funk2012). The Pellicciotti model computes melt as:
\begin{equation}
{a_\mathrm{sfc}}(z) = \left\{
\begin{array}{ll}
\mathrm{TF} \cdot T(z) + \mathrm{SRF} \cdot (1 - \alpha(z)) \cdot G(z) & \mathrm{if\,} T(z) \gt T_{\mathrm{melt}} \\
0 & \mathrm{if\,} T(z) \leq T_{\mathrm{melt}},
\end{array}
\right.
\end{equation} where 
$\mathrm{TF}$ is the temperature factor 
$\left( \mathrm{m\,w.e.\,K}^{-1}\, \mathrm{d}^{-1} \right)$, 
$\mathrm{SRF}$ is the shortwave radiation factor 
$\left( \mathrm{m}^3\, \mathrm{d}^{-1}\, \mathrm{W}^{-1} \right)$, 
$\alpha(z)$ is the albedo 
$\left( \mathrm{W\,m}^{-2} \right)$ parametrized following Brock and others (Reference Brock, Willis and Sharp2000) and G(z) is the incoming shortwave radiation 
$\left( \mathrm{W\,m}^{-2} \right)$. Note that for this model, 
$T_{\mathrm{melt}} = 1^{\circ}\mathrm{C}$ (Pellicciotti and others, Reference Pellicciotti, Brock, Strasser, Burlando, Funk and Corripio2005). The outcomes from the two melt models are combined by taking the mean between the two. More detailed information about the CRAMPON model set-up is provided by Landmann and others (Reference Landmann, Künsch, Huss, Ogier, Kalisch and Farinotti2021).
To model the evolution of the snow cover, we assume the initial snow depth at the start of the study period to be zero over the entire glacier. Consequently, the snow water equivalent (SWE) is updated for each elevation band depending on the climatic mass balance, i.e. the modelled daily accumulation and melt. To extract the modelled SCAF at each time step, we sum the area of the elevation bands in the model with a SWE larger than zero with respect to the beginning of each hydrological year and divide it by the glacier area.
3.4. Model calibration
To tackle parameter equifinality during calibration, we combine the information provided by the SCAF observations with observed long-term mean mass balance according to GLAMOS. For every glacier, the mean mass balance over the period 2015–23 is available from annual in situ observations at a network of sites extrapolated to the entire surface (Huss and others, Reference Huss, Dhulst and Bauder2015, GLAMOS, 2023). This is used to constrain the melt parameters (cf. Eqns 2 and 3), whereas the SCAF observations over the same period are used to constrain the precipitation correction factor 
$c_{\rm{prec}}$ (Eqn 1). This is achieved in three steps using a brute force approach to select the optimal parameter set within a pre-defined parameter space for 
$c_{\rm{prec}}$ (Fig. 3).
Figure 3. Model calibration on average mass balance and SCAF observations. Step 1): For each of the two melt models j = [1, 2], 
$c_{\rm{prec,i}}$ is varied within the interval 0.5–5.0 in steps of 0.1, and the melt parameters are calibrated to match the long-term mean mass balance. Step 2): The models are run with every value of 
$c_{\rm{prec},i}$, and for every model run, the root-mean-square error (RMSE) between the average SCAF and the observed SCAF is calculated. Step 3): The 
$c_{\rm{prec},i}$ with the lowest RMSE is selected as the optimal precipitation correction factor (
$c_{\rm{prec,opt}}$).
In the first step, the precipitation correction factor 
$c_{\rm{prec}}$ is varied in the interval 0.5 – 5.0 at equally spaced steps of 0.1, and for each of the values, the parameters of the melt models are tuned to agree with the average observed mass balance over 2015–23. This provides a total of 45 parameter sets per model. In the second step, the root-mean-square error (RMSE) between the modelled and observed SCAF is calculated for every parameter set. To avoid spurious results caused by potential outliers (for example, very high SCAF after a summer-snowfall event or unreliable SCAFs due to misclassification of shadows and clouds), the RMSE is only calculated accounting for SCAF observations that are contained within a given interval. This interval is defined by the course of the SCAF throughout the summer season that would be expected to occur after a very snow-rich and a very dry winter. More specifically, we set the upper bound of the interval for realistic SCAF values to be the course of the SCAF that is modelled for a winter balance of 2.1 m w.e., i.e. the average of the maximum winter balance values observed for each of the ten glaciers during 2015–23. Similarly, we set the lower bound of the interval to be the course of the SCAF modelled for a winter balance of 0.8 m w.e., i.e. the average of the observed minimum winter balances. Although these values are arbitrary to some degree, we found them useful to exclude spurious SCAF observations with strong leverage on the computed RMSE. This substantially reduces the biasing effect of individual SCAF observations on the selected 
$c_{\rm{prec}}$ value. In the third step, the 
$c_{\rm{prec},i}$ with the lowest RMSE is selected as the optimal precipitation correction factor 
$c_{\rm{prec,opt}}$. Together with the corresponding melt parameters, this defines the optimal model parameter set.
To evaluate the benefit of our approach, we also calibrate the models with a second strategy that serves as a benchmark. In this strategy, the models are calibrated using observed mean mass balance but no snow-cover information. In this case, the precipitation correction factor is set to 
$c_{\rm{prec}}$=1.6, i.e. the median value obtained for the ten glaciers and the period 2015–23 when calibrating the model with seasonal mass-balance observations as described in Landmann and others (Reference Landmann, Künsch, Huss, Ogier, Kalisch and Farinotti2021). This scenario acts as the reference, pretending that SCAF observations were not available but rough constraints on snow accumulation based on mean observed values are given.
By deriving the mean absolute deviation (MAD) between the modelled and the observed seasonal glacier-specific mass balance from GLAMOS (2023) for each of the two calibration strategies, we evaluate the benefit of constraining glacier mass balance with sub-seasonal SCAF observations.
The outcomes of the two calibration strategies are compared in more detail for Silvrettagletscher due to its comprehensive long-term observational network, allowing thorough control of winter and annual mass-balance distribution. In particular, we compare the daily cumulative mass balance and the daily cumulative glacier storage change, i.e. the glacier-specific mass balance multiplied by the glacier area, resulting from the two calibration strategies over different time periods. The periods are thereby chosen to be (i) the hydrological year, (ii) the melt season defined to last from 18 May (the long-term spring maximum date) to 30 September and (iii) heatwave periods with extreme glacier melt that occurred in summer 2022. This allows us to evaluate the relevance of the approach at the seasonal and sub-seasonal scale, as well as during dry and hot periods of a few days, when extreme glacier melt is particularly important. The heatwave periods are defined according to Cremona and others (Reference Cremona, Huss, Landmann, Borner and Farinotti2023) and are based on temperature time series from five meteorological stations distributed across Switzerland. Following this approach, the first heatwave occurred between 17 and 21 June 2022, the second between 14 and 26 July and the third between 31 July and 6 August.
4. Results and discussion
4.1. Performance evaluation
As an example, Fig. 4 illustrates the application of our approach for Rhonegletscher. The progressive melting of the snow cover during the melt season of 2018 is depicted in Fig. 4a. Until the first half of June, the glacier was completely snow-covered, as evidenced by the scene on 16 June, when the melt-out just started (SCAF = 0.97, Fig. 4a.2). The snow cover gradually depleted after that, causing a gradual decrease of the SCAF, which reached a value of 0.55 on the 20 August 2018 (Fig. 4a.6). Due to cooler temperatures and occasional snowfalls, the SCAF decrease slowed down in September, stabilizing at 0.53 by the 9th and 19th of September (Fig. 4a.7–8). The modelled SCAF and cumulative glacier-specific mass balance for Rhonegletscher are very consistent with the observations over the study period. The MAD between modelled and observed SCAF over 2015–23 is 4%. The MAD between modelled and observed seasonal mass balance from GLAMOS (2023) over 2015–23 is 0.32 m w.e. for winter and 0.24 m w.e. for summer mass balance. This comparison between model output and observations is shown for the SCAF (Fig. 4b) and for the mass balance (Fig. 4c) for three selected years, i.e. 2018 (an average year), 2021 (a year with favourable glacier mass balance conditions) and 2022 (very negative mass balance).
Figure 4. Constraining sub-seasonal mass balance with SCAF observations on Rhonegletscher. (a) Evolution of the observed snow cover over summer 2018. (b) Comparison between modelled SCAF (blue lines) and observed SCAF (black points with uncertainty bars, cf. Section 4.3) over the summer seasons of 2018, 2021 and 2022. (c) Comparison between modelled daily cumulative glacier-specific mass balance (blue lines) and observed seasonal glacier-specific mass balance (red dots with bars representing the observation uncertainty of ±0.25 m w.e. (Huss and others, Reference Huss2021)). The blue-shaded areas in panels (b) and (c) indicate the model uncertainties, derived as described in Section 4.3.
A similarly good agreement between modelled and observed SCAF and glacier-specific mass balance is observed at the other nine study sites too (Table 2). For the ten study sites combined, the MAD between modelled and observed SCAF is 6% (Fig. 5), while the MAD between modelled and observed glacier-specific mass balance is 0.28 m w.e. both for the winter and the annual mass balance (Fig. 6). Such relatively small deviations align with the findings by Barandun and others (Reference Barandun2018) and illustrate the potential that SCAF observations hold for model calibration when used in combination with long-term mean mass balances, from in situ observations or geodetic mass balances from the comparison of repeated DEMs.
Figure 5. Comparison between modelled and observed SCAF for the ten study sites (depicted by different colours). The mean absolute deviation (MAD) is 6%. See Table 1 for glacier names.
Figure 6. Comparison between modelled and observed mass balance. The blue line shows the average of the two melt models while the shaded area shows the uncertainty (cf. Section 4.3). The seasonal mass-balance observations from GLAMOS (2023) are shown by the red points with bars representing the observation uncertainty of ±0.25 m.w.e (Huss and others, Reference Huss2021). The availability of SCAF observations is provided at the top of each panel (vertical lines). The colour of the lines indicates whether the observation was used during calibration (green) or discarded according to the filtering criterion described in Section 3.4 (red).
Table 2. Mean absolute deviation between modelled and observed (i) winter mass balance, (ii) annual mass balance and (iii) SCAF over the period 2015–23
4.2. Relevance for determining seasonal and sub-seasonal mass balances
To evaluate the benefit of our approach when determining glacier mass balances at the seasonal and sub-seasonal scale, we calibrate the models with two strategies. In the first strategy, which serves as benchmark, we pretend that the SCAF observations are not available and only calibrate the models against the 2015–23 mean mass balance. In the second strategy, instead, we follow the approach presented in Section 3.4, i.e. we make use of the SCAF observations.
A comparison of the results obtained with these two calibration strategies is shown in Fig. 7. Over the ten considered glaciers, integrating the SCAF observations during model calibration reduces the MAD of the winter mass balance from 0.36 m.w.e to 0.28 m.w.e, i.e. by 22%. The MAD of the annual mass balance, instead, remains virtually unchanged: 0.29 m.w.e when ignoring and 0.28 m.w.e when including the SCAF observations. This latter result is expected, as both calibration strategies are constrained by the mean annual mass balance over the period 2015–23. When considering glaciers individually, the added value of our approach is even more evident: for Silvrettagletscher, where the benefit is most important, the MAD of the winter mass balance is lowered as much as by 70%, i.e. from 0.92 m w.e. to 0.26 m w.e. This demonstrates the capacity of SCAF observations to constrain glacier-specific accumulation more accurately, especially when site-specific characteristics, such as local topography or meteorological conditions, strongly affect snowfall amounts, or the precipitation dataset used is strongly biased.
Figure 7. Benefit of the calibration strategy using SCAF observations. The mean absolute deviation between modelled and observed seasonal mass balance is compared for both winter (blue) and annual (rose) mass balances. ‘cal. 
$\overline{\dot B}$ no SCAF’ stands for the benchmark calibration strategy in which no SCAF observations are used while ‘cal. 
$\overline{\dot B}$ with SCAF’ stands for the calibration strategy presented in Section 3.4. The boxplots show the mean value (red line) together with the 25%/75% (boxes) and 2.5%/97.5% (bars) quantiles.
To evidence the importance of accurately reproducing the winter mass balance and to show the impact it has on computed summer melt rates, we compare glacier mass balance and storage change for Silvrettagletscher resulting from the two calibration strategies. Figure 8a shows the 2015–23 average evolution of the cumulative glacier mass balance over the course of a year for the two calibration strategies, as well as the difference between them. At the end of winter, the mass balance calibrated with the presented approach is about twice as high (i.e. 0.7 m w.e. more positive) as in the benchmark calibration and agrees far better with the in situ observations from GLAMOS. Figure 8b shows the 2015–23 average evolution of the cumulative glacier storage change during the melt season, again discerning the two calibration strategies. At the end of the melt season, the model calibrated with the SCAF observations indicates 1.85 · 106 m3 additional glacier mass change compared to the model calibrated without SCAF observations, i.e. a surplus of 32% relative to the annual meltwater yield. This again highlights the importance of accurately constraining winter accumulation, also when being interested in mass changes and water runoff dynamics that happen over the ablation season.
Figure 8. Effect of the calibration strategy on the cumulative mass balance and the cumulative glacier storage change for the example of Silvrettagletscher. Violet colours refer to the benchmark calibration while orange colours refer to the presented approach relying on SCAF observations. The shading reflects the variability between the two melt models and the fluctuations in mass balance over 2015–23. (a) Evolution of the 2015–23 average daily cumulative mass balance (top) and difference between the two calibration strategies (bottom). The observed mass balance from GLAMOS (2023) is shown with red markers. (b) Same as panel ‘(a)’ but for the daily cumulative glacier storage change over the melt season. (c) Daily cumulative glacier storage change during the three summer heat waves of 2022 (top) and difference between calibration strategies (bottom).
Finally, we compare the glacier storage change occurring during the 2022 summer heat waves for both calibration approaches (Fig. 8c). In the first, second and third heat wave (i.e. during 17–21 June, 14–26 July and 31 July–6 August), the difference in the cumulative glacier storage change calculated with the two calibration strategies corresponds to 1.5 · 105 m3 (34%), 4.8 · 105 m3 (40%) and 2 · 105 m3 (38%), respectively. Also these differences indicate the importance of accurate constraints on winter accumulation, even when dealing with extreme melt rates over relatively short periods. The finding is of particular relevance when considering that the frequency of such events is expected to increase in the future (Fischer and others, Reference Fischer, Sippel and Knutti2021), and when considering that the glacier melt contributions during such periods help in mitigating the effects of potential droughts (Zappa and Kan, Reference Zappa and Kan2007, Pelto and others, Reference Pelto, Dryak, Pelto, Matthews and Perry2022, Ultee and others, Reference Ultee, Coats and Mackay2022, Cremona and others, Reference Cremona, Huss, Landmann, Borner and Farinotti2023). We thus suggest that our calibration strategy could be valuable when aiming at estimating the water availability in such particularly critical periods.
The advantage of our approach lies in its potential for application at the regional scale as the method can separate the components of the mass balance, i.e. accumulation and ablation, without the need of laborious in situ seasonal mass-balance observations. Thus, the method successfully addresses the challenge of equifinality between accumulation and melt parameters by combining the information content of mass change and SCAF observations (Fig. 9). For the example of Rhonegletscher, Fig. 9a illustrates that several combinations of accumulation and melt parameters may reproduce the observed mean mass balance with the same goodness-of-fit (GOF), taken here as the normalized absolute deviation between the observed and modelled mass balance. All parameter combinations on the diagonal of the explored parameter space (dark green colours in Fig. 9a) are equally valid, i.e. result in the same GOF. However, when additionally considering the agreement between observed and modelled SCAF, measured by the normalized RMSE between observed and modelled quantities, the RMSE reveals a clear minimum (Fig. 9b–d). This demonstrates that the equifinality between accumulation and melt parameters is successfully resolved. Note that a certain degree of equifinality still remains between the two melt parameters, i.e. the parameters separately accounting for the effects of temperature and solar radiation (cf. Eqns 2 and 3). With our approach, this equifinality cannot be resolved, as two parameters are calibrated on one observation (the mean mass balance). We argue that this equifinality is much less relevant for modelling sub-seasonal mass balance than the equifinality between accumulation and melt parameters, especially because the relation between the two melt parameters is better bounded by physical constraints.
Figure 9. Illustration of parameter equifinality for the investigated glaciers. (a) Normalized absolute deviation between modelled and observed mean glacier-wide mass balance for different parameter sets for Rhonegletscher. (b) Normalized root-mean-square error between modelled and observed SCAF for the parameter sets with similar performance regarding mean mass balance. The red dot indicates the parameter set with the lowest normalized RMSE, i.e. the optimal parameter set. In panels (a) and (b), dark green indicates a good fit, whereas dark red indicates a poor fit. GOF stands for goodness-of-fit. (c) Visualization of the RMSE with respect to observed SCAF along the black transect in panel (b) for Rhonegletscher. (d) RMSE with respect to observed SCAF for the ten study sites, with the optimal parameter set represented by the red dot.
The potential for regional-scale assessments is confirmed when comparing the precipitation correction factor obtained with the present approach with the one resulting from the calibration on seasonal mass-balance observations (e.g. Huss and others, Reference Huss, Bauder and Funk2009, Landmann and others, Reference Landmann, Künsch, Huss, Ogier, Kalisch and Farinotti2021) over the period 2015–23 (Fig. 10). The high correlation with an MAD of 0.16 demonstrates that our approach is able to resolve actual accumulation amounts through remotely sensed observations, i.e. without the need of direct in situ observations, therefore suggesting a considerable application potential in areas where direct measurements are difficult or unavailable.
Figure 10. Comparison between the precipitation correction factor obtained by the calibration approach presented in our study and the precipitation correction factor constrained by seasonal mass-balance observations. The MAD is 0.16.
4.3. Sensitivity analysis
Uncertainties in the modelled mass balance were quantified by conducting a sensitivity analysis. This analysis consisted in re-calibrating the model while randomly varying the SCAF observations and the observed average mass balance within given uncertainty ranges. Variations in SCAF observations were imposed by assuming a normally distributed uncertainty with a mean of ‘zero’ and a standard deviation of 0.05. This results in SCAF uncertainties potentially exceeding ±0.1, which we consider to be a conservative estimate when compared to values reported in previous studies (Huss and others (Reference Huss2013) and Kenzhebaev and others (Reference Kenzhebaev, Barandun, Kronenberg, Chen, Usubaliev and Hoelzle2017), for example, estimate SCAF uncertainties of 0.025 and 0.1, respectively). For the mean mass balance, instead, we assumed the uncertainty proposed by Huss and others (Reference Huss2021) for decadal average mass balances, i.e. ±0.1 m.w.e y−1. We performed ten simulations by varying the mean mass balance by +0.1 m.w.e y−1 and superimposing the random uncertainty on the SCAF, and ten simulations by varying the mean mass balance by −0.1 m.w.e y−1, again superimposing the uncertainty on the SCAF. This results in a total of 20 simulations. To evaluate the overall model uncertainty, we use the interval between the minimum and maximum output resulting from the mean of the two melt models for all 20 model realizations, and we account for the spread between the two melt models. When averaged over the 20 simulations, the ten study glaciers, and the period 2015–23, the uncertainty in the modelled SCAF over the melt season reaches a maximum of 0.16 in individual years, with a mean of 0.03. The average uncertainty of the seasonal mass balance, instead, is 0.19 m.w.e for winter and 0.13 m.w.e for summer. From this assessment, we infer that our approach is generally more sensitive to variations in the mean mass balance as compared to variations in the SCAF. This is an important aspect to consider when applying the method at the regional scale relying on geodetic mass balance estimates.
Our approach assumes a constant glacier geometry over the study period, motivated by the relatively short time span considered. To investigate the effect of changing glacier area on model results, we performed two simulations assuming fixed rates of glacier area reduction. Given that Linsbauer and others (Reference Linsbauer2021) reported an average area reduction rate of −0.8% per year for Swiss glaciers and allowing for some variability in this value, we assumed a reduction rate of either −0.5% or −1% per year. To update the glacier area every year, the ablation area of each glacier is reduced by the specified amount. The ablation area is thereby assumed to be the portion of the glacier below the glacier’s median elevation. For the two cases, differences in the modelled SCAF compared to the reference case with constant geometry (averaged over the ten study sites) correspond to 0.013 and 0.018 on average, respectively, but may reach up to 0.066 for individual years. On average over the ten study sites, differences for the seasonal mass balance compared to the reference case are nearly equivalent for summer and winter balance and correspond to 0.04–0.07 m.w.e, for area reduction rates of −0.5% and −1% per year, respectively.
To investigate the effect of the selected melt models on the outputs, we performed two additional simulations in which the approach was first re-run only using the Hock model, and then only using the Pellicciotti model. On average over the study sites, differences in the modelled SCAF over the melt season between the two simulations reach a maximum of up to 0.30 with a mean of 0.07. The average difference of the seasonal mass balance between the two models, instead, is 0.13 m.w.e for winter and 0.14 m.w.e for summer and thus within the range of uncertainty. By averaging the results of the two melt models, we aim to mitigate potential inconsistencies due to uncertainties of the individual models and to enhance the overall reliability of our approach.
4.4. Limitations and challenges
While the proposed approach allows for better constraining seasonal mass-balance dynamics of unmeasured glaciers, thus offering potential for regional mass balance assessments, it also has some limitations. Between about October and May, for example, the glaciers are often completely snow-covered and the SCAF values hold virtually no information on the mass balance. Furthermore, glacier size and elevation range can impact the method’s effectiveness. SCAF observations are most useful when the SCAF slowly decreases throughout the melt season, as happens for large glaciers spanning an extensive elevation range. For small glaciers or for glaciers with limited elevation ranges, instead, the transition between completely snow-covered conditions to snow-free conditions is often fast. This presents a challenge for the approach, as the observed signal is less reliable because few observations (if any) are present during the melt-out period. In addition, local effects, such as peculiar spatial patterns in the snow cover distribution, for example due to wind drifted snow, play a more significant role on smaller glaciers. Such effects are difficult to capture with our model set-up describing glaciers in elevation bands and thus increase the likelihood of inaccurate parameter estimates. Still, our method showed to be effective for relatively small glaciers such as Ghiacciaio del Basodino (1.64 km2) and Adlergletscher (2.10 km2) as well as for glaciers with small elevation ranges such as Glacier de la Plaine Morte, providing results comparable with larger glaciers (Table 2, Fig. 5). This is an important consideration since it indicates that the method can maintain its potential also in the future, with ongoing climate change causing the glaciers to become smaller and to extend over more limited elevation ranges. This notwithstanding, we emphasize the importance of considering the specific characteristics of individual glaciers when applying our method, and to keep in mind the method’s limitations, especially on small glaciers.
5. Conclusions
In this study, we presented a method for constraining glacier mass balance at the sub-seasonal scale relying on Sentinel-2-derived SCAF observations and long-term mean mass balance. We showcased the potential of the method for ten glaciers in the Swiss Alps over the period 2015–23. The combination of SCAF observations with long-term mass balances was proven to be effective in tackling the problem of parameter equifinality, i.e. in separating the total glacier balance in its accumulation and ablation components. The method holds promise for regional-scale applications, as it relies on remote-sensing products while avoiding the need for in situ measurements. On average over the ten considered study sites, the approach provides estimates of winter mass balance that are 22% more accurate than using an average value of the precipitation correction factor, while improvements of up to 70% were found in individual cases. Furthermore, the study evidenced the importance of constraining winter accumulation for accurately capturing the evolution of mass balance over the melt season, leading to important differences in the cumulative glacier storage change at the end of the melt season (in our example at Silvrettagletscher, the difference almost reached 2 million m3 of water, i.e. 32% of the glacier storage change produced over the melt season). Such differences were shown to be significant also for shorter periods (e.g. weeks) and for periods of extreme melting (e.g. during heatwaves), emphasizing the relevance of our approach in the context of short-term glacier monitoring. We thus suggest that the integration of remotely sensed SCAF observations into monitoring approaches has a substantial potential to provide more accurate mass-balance estimates at the sub-seasonal scale, which can, in turn, be important for the optimization of water resources in the Alps, especially during dry and hot summers.
Acknowledgements
The study is supported by the ESA project ‘Glacier Science in the Alps’ (Grant: 4000133436/20/I-NB9), and by the project ‘Machine Learning aided ForecastIng of drought related eXtremes (MaLeFiX)’, conducted in the frame of the programme ‘eXtremes’ of the Swiss Federal Institute for Forest, Snow and Landscape Research (WSL). Swiss Meteorological data were.
Open access
