the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Assimilating nearrealtime mass balance stake readings into a model ensemble using a particle filter
Johannes Marian Landmann
Hans Rudolf Künsch
Matthias Huss
Christophe Ogier
Markus Kalisch
Daniel Farinotti
Shortterm glacier variations can be important for water supplies or hydropower production, and glaciers are important indicators of climate change. This is why the interest in nearrealtime mass balance nowcasting is considerable. Here, we address this interest and provide an evaluation of continuous observations of point mass balance based on online cameras transmitting images every 20 min. The cameras were installed on three Swiss glaciers during summer 2019, provided 352 nearrealtime point mass balances in total, and revealed melt rates of up to 0.12 m water equivalent per day ($\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) and of more than 5 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ in 81 d. By means of a particle filter, these observations are assimilated into an ensemble of three TI (temperature index) and one simplified energybalance mass balance models. State augmentation with model parameters is used to assign temporally varying weights to individual models. We analyze model performance over the observation period and find that the probability for a given model to be preferred by our procedure is 39 % for an enhanced TI model, 24 % for a simple TI model, 23 %, for a simplified energy balance model, and 14 % for a model employing both air temperature and potential solar irradiation. When compared to reference forecasts produced with both mean model parameters and parameters tuned on single mass balance observations, the particle filter performs about equally well on the daily scale but outperforms predictions of cumulative mass balance by 95 %–96 %. A leaveoneout crossvalidation on the individual glaciers shows that the particle filter is also able to reproduce point observations at locations not used for model calibration. Indeed, the predicted mass balances is always within 9 % of the observations. A comparison with glacierwide annual mass balances involving additional measurements distributed over the entire glacier mostly shows very good agreement, with deviations of 0.02, 0.07, and 0.24 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$
Glaciers around the world are shrinking. For example, Switzerland has lost already more than a third of its glacier volume since the 1970s (Fischer et al., 2015), glaciers are currently melting at about −0.6 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{\mathrm{1}}$ (water equivalent per year) on average (Sommer et al., 2020), and it is expected that glaciers will continue to lose mass (Jouvet et al., 2011; Salzmann et al., 2012; Beniston et al., 2018; Zekollari et al., 2019). Since glaciers are important for the supply of drinking water or for irrigation and electricity production, there is high interest in nearrealtime glacier mass balance information. Such information has also become important in the context of public outreach, e.g., for demonstrating the consequences of climate change (e.g., Euronews, 2019; Science Magazine, 2019).
A glacier mass balance nowcasting framework assimilating relevant observations could deliver these nearrealtime mass balances whenever required. While nowcasting frameworks exist, for example, for the mass balance of the Greenland Ice Sheet (Fettweis et al., 2013; NSIDC, 2020a), for snow (NSIDC, 2020b; SLF, 2020) or for hydrological purposes (Zappa et al., 2008; Pappenberger et al., 2016; Zappa et al., 2018; WSL, 2020; Hydrique, 2020; Wu et al., 2020), there is no specific framework providing such analyses for mountain glaciers yet. In general, data assimilation is widespread in oceanography, meteorology, hydrology, and snow sciences, “but its introduction in glaciology is fairly recent” (Bonan et al., 2014). Especially regarding glacier mass balance studies, data assimilation and Bayesian approaches have appeared only slowly in published work (Dumont et al., 2012; Leclercq et al., 2017; Rounce et al., 2020; Werder et al., 2020).
In many cases, mass balance analyses are available twice a year and are based on seasonal in situ observations (Cogley et al., 2011). This relatively low frequency is related to the fact that in situ observations are expensive in terms of both time and manpower. Only recently have lowcost and highfrequency monitoring approaches emerged (Hulth, 2010; Fausto et al., 2012; Keeler and Brugger, 2012; Biron and Rabatel, 2019; Carturan et al., 2019; Gugerli et al., 2019; Netto and ArigonyNeto, 2019). However, even with these observations, it is not straightforward to provide analyses at higher frequencies. This is because nearrealtime estimates are often based on ensemble modeling in order to enable a correct quantification of uncertainties. Ensemble modeling is used in glaciology in the context of model intercomparison projects (Hock et al., 2019), future projections for ice sheets and mountain glaciers (e.g., Ritz et al., 2015; Shannon et al., 2019; Golledge, 2020; Marzeion et al., 2020; Seroussi et al., 2020), and also determining the initial conditions for modeling (Eis et al., 2019). However, ensembles are currently not prominent in the calculation of seasonal or daily glacier mass balances. Another reason why calculating higherfrequency glacier mass balance analyses is not straightforward is the lack of knowledge about the shortterm variability in the parameters of the necessary models. Temperature index (TI) models, for example, are parametrizations of the full energy balance equation, and they offset some of the changes occurring in the driving processes through parameter fluctuations (Ohmura, 2001; Lang and Braun, 1990; Hock, 2003; Hock et al., 2005). In a comparison of four TI models and a full energy balance model, Gabbi et al. (2014) showed that all models perform very similarly on a multiyear scale.
In this study, we address the issue of lowfrequency observations, ensemble modeling, and the lack of knowledge about shortterm parameter variability as part of the project CRAMPON (Cryospheric Monitoring and Prediction Online). The latter aims at delivering nearrealtime glacier mass balance estimates for mountain glaciers using data assimilation. To obtain highfrequency data at a relatively low cost, we equipped three Swiss glaciers – Glacier de la Plaine Morte, Findelgletscher, and Rhonegletscher – with seven cameras in total. The cameras were operated in summer 2019 and took images of a mass balance stake marked every 2 cm at 20 min intervals, thus providing estimates of surface point mass balance aggregated to the daily scale. By using a particle filter (e.g., Arulampalam et al., 2002; Beven, 2009; Magnusson et al., 2017), we assimilate these observations into an ensemble of three TI models and one simplified energy balance model.
Ensemble stability and suitability for operational use is ensured by designing the particle filter such that, at any instance, each model has a minimum contribution to the mass balance model ensemble. In particular, models with temporarily bad performance are not excluded from the ensemble and can thus regain weight later. To address parameter uncertainty, we drive the mass balance model ensemble with both Monte Carlo samples of uncertain meteorological input and prior parameter distributions obtained from calibration to past, longerterm seasonal mass balance series. By using an augmented state formulation of the particle filter, we constrain model parameters as well (e.g., Ruiz et al., 2013). We are not aware of glacier mass balance studies that have applied a multimodel ensemble based on a particle filter with the resampling methods we propose, although multimodel particle filters have been used for other applications (e.g., Kreucher et al., 2004; Ristic et al., 2004; Saucan et al., 2013; Wang et al., 2016).
As a result, we demonstrate (1) how a workflow including daily melt observations, ensemble modeling, and data assimilation works in practice, (2) to what extent the assimilated mass balances are able to reproduce cumulative observations, and (3) how the ensemble performs with respect to both reference forecasts and seasonal analyses from in situ measurements.
We use Glacier de la Plaine Morte, Rhonegletscher, and Findelgletscher in summer 2019 as test sites (Fig. 1). The basic morphological characteristics and instrumentations of these glaciers are given in Table 1.
2.1 Continuous in situ mass balance observations
2.1.1 Technical camera station setup
For acquiring daily point mass balances in the field, we use offtheshelf cameras and logger boxes from the company Holfuy Ltd. We mount the cameras to a construction of aluminum stakes that we designed for glacier applications. Figure 2 provides an overview of the camera installation.
The camera observes an ablation stake, which is marked with colored tape at 2 cm intervals. When the surface melts, the aluminum construction holding the camera slides along the mass balance stake. Pictures are taken every 20 min and are sent in realtime to our servers via the Swiss mobile phone network. Differences between subsequent pictures are used to infer daily glacier surface height changes relative to the stake top, which are the basis for ablation measurements (Cogley et al., 2011). All pieces of the construction are lightweight (4 kg for the station plus 4 kg for 8 m of mass balance stakes) and can be mounted by one person.
2.1.2 Camera acquisitions in summer 2019
Figure 3 shows an overview of camera acquisitions and data gaps over the summer 2019.
In total, we obtained 352 daily point mass balance observations between 20 June 2019 and 3 October 2019. The camera longest in the field was on Glacier de la Plaine Morte (91 d between 20 June 2019 and 18 September 2019), while those shortest in the field were two cameras at the tongue of Rhonegletscher (52 d between 13 August 2019 and 2 October 2019). Very few data gaps occurred due to failure of the mobile network over which the data were transmitted.
Once camera images are on our servers, they are read manually to obtain daily cumulative surface height change h(t,z) since camera setup. We assume that the observational error ϵ_{t} of a reading is Gaussian distributed and uncorrelated in time and space. To estimate the standard deviation of the Gaussian error distribution, we performed a roundrobin experiment with seven participants. In this kind of experiment all participants were given the task to read h(t,z) from the same camera images independently, and statistics were made about the degree of agreement between the readings. We found a standard deviation of 1.5 cm, with a range from 0.2 to 1.7 cm. The estimate accounts for reading errors, errors in stake marker positions, and unknown thickness of the melt crust on the ice surface, but it does not account for systematic errors.
The relationship between observations of surface height change since an initial point in time (in our case the time at which a camera is set up) and the cumulative glacier mass balance is given through the simple linear observation operator ℋ:
where b_{sfc}(t,z) ($\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$) is the accumulated surface mass change at elevation z and time t since the day of the first camera observation, ρ_{w}=1000 kg m^{−3} is the water density, and ρ_{bulk} is the temporally weighted bulk density of snow and ice at the camera location (kg m^{−3}). We expect the observation errors to be uncorrelated in time, since every reading is independent from the previous one. To avoid systematic errors in the readings, we exclude the initial, snowcovered phase after camera setup at stations FIN 2 and PLM 1. This is because the camera construction can sink into the snow cover, potentially biasing the snowmelt signal. This “sinking bias” is virtually impossible to distinguish from the actual melt signal. Moreover, the temporally varying density of the melting snow is unknown. Short snow events during the melt season are assigned a density of 150 kg m^{−3}. The calculated snow water equivalent is assigned an uncertainty of 2–3 cm w.e. If a stake reading was impossible, we have resumed with a zero balance after the snowfall melted again. For days without snow, there are three cases that require special attention when reading mass balances: (1) maintenance operations like setup, redrilling, and unmounting of a station, (2) melt that happens during night and that is thus only visible on the next day, and (3) data gaps. Regarding maintenance operations (point “(1)”), we do not consider the observations from days when maintenance has taken place. This is either because those days are not fully covered or because the mass balance stake and the entire station might melt into the ice after redrilling. For nighttime melting (i.e., “(2)”), we equally distribute the overnight melt between the 2 d concerned. This is a tradeoff between warmer temperatures before midnight and colder temperatures but longer time span after midnight. For data gaps (point “(3)”), we experienced only short image transmission outages which were mainly due to a 6 d failure in the mobile network connection on Plaine Morte during September 2019. We excluded the daily readings on these days but were able to reconstruct estimates of cumulative mass balance over the time gaps when acquisitions resumed.
2.2 Meteorological input data
To model glacier mass balance, we employ verified products of daily mean and maximum 2 m temperature T and T_{max}, precipitation sum P, and mean incoming shortwave radiation G from MeteoSwiss as model input (MeteoSwiss, 2018, 2021a, b). These are delivered on grids with approx. 0.2^{∘} spatial resolution, which for Switzerland corresponds to a horizontal resolution of about 2 km.
Temperature uncertainty, given as a rootmeansquare error, varies per season from 0.94 K (May to March) to 1.67 K (December to February) in the Alpine region (Frei, 2014; Christoph Frei, personal communication, 2020). We assume a Gaussiandistributed additive error perfectly correlated in space for a single glacier but independent on different days. The assumption of a perfect error correlation can be justified by the fact that the station network from which the gridded temperature values are interpolated is much sparser than the scale of individual glaciers. The air temperature lapse rate is derived from the 25 closest cells to a glacier outline centroid using a Bayesian estimation based on a linear regression model:
where T_{t,i} is the temperature of the ith grid cell out of the 25 considered cells at time t, e_{t} is the regression line intercept, q_{t} is the regression slope (i.e., $\frac{\partial T}{\partial z}$), h_{i} is the height of the ith grid cell, and ${\mathit{\nu}}_{t,i}\sim \mathcal{N}(\mathrm{0},{\mathit{\sigma}}_{\mathit{\nu},t}^{\mathrm{2}})$ are the residuals independent in space and time. Using a gprior of Zellner (Zellner, 1986), being noninformative in the intercept e_{t} and model noise variance ${\mathit{\sigma}}_{\mathit{\nu},t}^{\mathrm{2}}$ of the regression, we draw samples of the lapse rate q_{t} from the following posterior distribution:
with
Above, p(⋅) means “probability of”, g determines a weighting factor composing the posterior mean (we set g=1), $\widehat{{q}_{t}}$ is the least squares estimator of the slope, q_{0} is the prior mean, which we choose to be an annually varying climatological mean gradient at the respective grid location, $\overline{h}$ is the average height of the 25 grid cells, and ${s}_{t}^{\mathrm{2}}$ is the residual sum of squares. This is, up to a constant, the probability density of a tdistributed random variable with 24 degrees of freedom, shifted by $\frac{(g\widehat{{q}_{t}}+{q}_{\mathrm{0}})}{(\mathrm{1}+g)}$ and multiplied by c. The samples drawn from this distribution are then propagated into the particle filter.
For operational reasons, the precipitation grids contain the 06:00–06:00 local time precipitation sums and thus do not conform with the 00:00–00:00 temperature values (Isotta et al., 2019; MeteoSwiss, 2021b). This might introduce an error, which we cannot quantify. As for temperature, we thus focus again on random errors and pretend for simplicity that the precipitation sum was also from 00:00–00:00. Precipitation uncertainty is generally harder to assess than temperature uncertainty since it involves undercatch errors and skew error distributions. Here, we follow the error assessment by quantiles of precipitation intensity proposed by Isotta et al. (2014) and Christoph Frei, personal communication (2020), who calculated that, for the Alpine region, the standard error at moderate precipitation intensities roughly corresponds to an over or underestimation by a factor of 1.25. The error increases (decreases) towards low (high) precipitation intensities, and it is generally slightly higher in the summertime. We draw samples from a multiplicative Gaussian error distribution and for the same reason as for temperature, we assume perfect precipitation error correlation at the glacier scale. We also derive Bayesian precipitation lapse rates from the surrounding 25 grid cells in the same fashion as we do for the temperature lapse rate. However, to circumvent high errors in the slope calculation due to the boundedness of precipitation towards zero, we (1) calculate the slope on the square root of the precipitation and (2) assign a probability that the reference has actually received precipitation when the reference cell value is zero but other cells have nonzero precipitation.
Shortwave radiation is derived using data from the geostationary satellite series Meteosat. As an uncertainty, Stöckli (2013) gives a mean absolute bias between 9 and 29 W m^{−2}. We assume the errors to be Gaussian and assign a standard deviation of 15 W m^{−2}, perfectly correlated on the glacier scale and independent in time. Shortwave radiation is downscaled from the grid to the glacier with potential radiation (see Sect. 3.1).
2.3 Glacier outlines and measured mass balances
Glacier outlines for the year 2019 are obtained from GLAMOS (Glacier Monitoring Switzerland), and mass balances are calculated over a fixed glacier surface area (Elsberg et al., 2001; Huss et al., 2012).
For calibration and verification, we use mass balance data acquired in the frame of GLAMOS (Glacier Monitoring Switzerland, 2018). First, intermediate stake readings independent from our nearrealtime stations but close to our installations have been made explicitly for this study. Stake locations are depicted in Fig. 1. The reading error for these measurements is usually estimated to be around 5 cm (e.g., Müller and Kappenberger, 1991). Second, we use seasonal, glacierwide mass balances based on in situ observations. These observations are acquired during two field campaigns in April and September, respectively. Glacierwide mass balances are obtained by extrapolating the in situ observations and making the extrapolated values consistent with longterm mass changes. The latter procedure is sometimes referred to as “homogenization” (e.g., Bauder et al., 2007; Huss et al., 2015). For recent years, this homogenization has not been performed since no geodetic mass balances are available yet. The extrapolation method used to infer glacierwide mass balance from point measurements involves an adjustment of the model parameters of an accumulation and TI melt model (Hock, 1999) at locations where observations are available, while mass balances at grid cells without observations are produced using the calibrated model (Huss et al., 2009, 2015). Uncertainties of the glacierwide annual mass balance for the measurement period have been estimated to be 0.09–0.2 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ in six experiments in which GLAMOS (1) model parameters (temperature lapse rate, ratio between melt coefficients, summer precipitation correction) and (2) snow extrapolation parameters have been varied within prescribed ranges and for which (3) mass balance stake reading uncertainty, (4) DEM (digital elevation model) and outline uncertainty, (5) climate forcing uncertainty, and (6) point data availability have been accounted.
3.1 Mass balance modeling
Glacier surface mass balance consists of two components: accumulation and ablation. We model accumulation and ablation on elevation bins whose vertical extent is determined by a ≈20 m horizontal spacing of nodes along the central flow line of the glacier (Maussion et al., 2019). To obtain glacierwide mass balance, node mass balances are weighted with the area per elevation bin. To compute accumulation at different elevations, we employ a simple but widely used accumulation model (e.g., Huss et al., 2008):
where c_{sfc}(t,z) ($\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$) is the snow accumulation at time step t and elevation z, prcp_{scale}(t) is the unitless multiplicative precipitation correction factor, P_{s}(t) is the sum of solid precipitation at the elevation of the precipitation reference cell z_{ref} and time step t, and $\frac{\partial {P}_{\text{s}}}{\partial z}$ is the solid precipitation lapse rate. Following Sevruk (1985), we choose prcp_{scale} to vary sinusoidally by ±8 % around its mean during 1 year, being highest in winter and lowest in summer. This is to account for systematic variations in gauge undercatch depending on the precipitation phase. The water phase change in the temperature range around 0 ^{∘}C is modeled using a linear function between 0 and 2 ^{∘}C; i.e., at 1 ^{∘}C there is 50 % snow and 50 % rain (e.g., Maussion et al., 2019).
Since all three glaciers we investigate are in the GLAMOS measurement program and winter mass balance observations are available, the effect of spatial variations in snow accumulation differing from a linear gradient can be incorporated. This is done by choosing a factor D(z) such that the model mass balance in the elevation bins matches the measured and interpolated distribution of the winter mass balance (Farinotti et al., 2010):
Here, ${C}_{\mathrm{sfc}}^{\mathrm{w}}\left(z\right)$ ($\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$) is the modeled winter surface accumulation, i.e., the sum of individual c_{sfc}(t,z) over the winter period, and ${C}_{\mathrm{sfc},\mathrm{glamos}}^{\mathrm{w}}\left(z\right)$ ($\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$) is the interpolated winter surface accumulation at elevation z.
To model surface ablation, we set up an ensemble of three TI melt models and one simplified energybalance melt model. We choose this ensemble since the individual models differ in the degree of complexity by which they describe the surface energy balance (Hock, 2003). The models range from using only temperature as input for determining melt via employing additionally the potential irradiation to using temperature and the actual shortwave radiation. The ensemble contains the following:

The “BraithwaiteModel” uses only air temperature as input to calculate melt (Braithwaite and Olesen, 1989; Braithwaite, 1995):
$$\begin{array}{}\text{(7)}& {a}_{\mathrm{sfc}}(t,z)={\text{DDF}}_{\mathrm{snow}/\mathrm{ice}}\cdot max\left(T\right(t,z){T}_{\mathrm{melt}},\phantom{\rule{0.33em}{0ex}}\mathrm{0}),\end{array}$$where a_{sfc}(t,z) ($\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) and T(t,z) (^{∘}C) are surface ablation and air temperature at time step t and elevation z, respectively, ${\text{DDF}}_{\mathrm{snow}/\mathrm{ice}}$ ($\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) are the temperature sensitivities (“degreeday factors”) of the surface types (snow and ice), max() is the maximum operator, and T_{melt} (^{∘}C) is the threshold temperature for melt. For this application, we set T_{melt} to 0 ^{∘}C and keep the ratio of ${\text{DDF}}_{\mathrm{snow}}/{\text{DDF}}_{\mathrm{ice}}$ constant at 0.5 (Hock, 2003).

The “HockModel” uses potential incoming solar radiation as an additional predictor for melt (Hock, 1999):
$$\begin{array}{}\text{(8)}& \begin{array}{rl}{a}_{\mathrm{sfc}}(t,z)& =(\text{MF}+{a}_{\mathrm{snow}/\mathrm{ice}}\cdot {I}_{\mathrm{pot}}(t,z\left)\right)\\ & \cdot max\left(T\right(t,z){T}_{\mathrm{melt}},\mathrm{0}),\end{array}\end{array}$$where MF ($\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) is the temperature melt factor, ${a}_{\mathrm{snow}/\mathrm{ice}}$ ($\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{W}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}$) are the radiation coefficients for snow and ice, respectively, I_{pot}(t,z) (W m^{−2}) is the potential clearsky direct solar radiation at time t and elevation z, T_{melt} is set again to 0 ^{∘}C, and the ratio of ${a}_{\mathrm{snow}}/{a}_{\mathrm{ice}}$ is 0.8 (Hock, 1999; Farinotti et al., 2012). I_{pot}(t,z) is computed at 10 min intervals following the methods described in Iqbal (1983), Hock (1999), and Corripio (2003) and by using swissALTI3D (swisstopo, 2020) as a background elevation model. Daily values are then obtained by averaging these values, and values for the different glacier elevations are aggregated. We assume equal uncertainties for both actual and potential incoming shortwave radiation G and I_{pot}.

The “PellicciottiModel” employs surface albedo and actual incoming shortwave solar radiation (Pellicciotti et al., 2005):
$$\begin{array}{}\text{(9)}& {a}_{\mathrm{sfc}}(t,z)=\left\{\begin{array}{l}\text{TF}\cdot T(t,z)+\text{SRF}\cdot (\mathrm{1}\mathit{\alpha}(t,z\left)\right)\\ \cdot G({I}_{\mathrm{pot}},t,z),\phantom{\rule{0.25em}{0ex}}\text{for}T(t,z){T}_{\mathrm{melt}}\\ \mathrm{0},\phantom{\rule{0.25em}{0ex}}\text{for}T(t,z)\le {T}_{\mathrm{melt}}\end{array}\right),\end{array}$$where TF ($\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) is the temperature factor, SRF (${\mathrm{m}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{W}}^{\mathrm{1}}$) is the shortwave radiation factor, and α(t,z) and $G({I}_{\mathrm{pot}},t,z)$ (W m^{−2}) are the albedo and incoming shortwave radiation at time t and elevation z, respectively. Note that in this case, T_{melt}=1 ^{∘}C (Pellicciotti et al., 2005).
Albedo is approximated according to the combined decay equation for deep and shallow snow in Brock et al. (2000):
$$\begin{array}{}\text{(10)}& \begin{array}{rl}\mathit{\alpha}(t,z)& =(\mathrm{1}{e}^{(\text{swe}(t,z)/{\text{swe}}^{*})})\\ & \cdot ({p}_{\mathrm{1}}{p}_{\mathrm{2}}\cdot {\mathrm{log}}_{\mathrm{10}}({T}_{\mathrm{acc}}(t,z)\left)\right)\\ & +\phantom{\rule{0.33em}{0ex}}{e}^{(\text{swe}(t,z)/{\text{swe}}^{*})}\cdot \left({\mathit{\alpha}}_{u}\right(t,z)\\ & +{p}_{\mathrm{3}}\cdot {e}^{{p}_{\mathrm{4}}\cdot {T}_{\mathrm{acc}}(t,z)}),\end{array}\end{array}$$where swe(t,z) is the snow water equivalent at time t and elevation z, ${\text{swe}}^{*}=\mathrm{0.024}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ is a scaling length for swe, p_{1}=0.713, p_{2}=0.155, p_{3}=0.442, and p_{4}=0.058 are empirical coefficients as given in Brock et al. (2000), α_{u} is the albedo of the underlying firn/ice below the snow, and T_{acc}(t,z) is the accumulated daily maximum temperature >0 ^{∘}C since a snowfall event at elevation z. To avoid infeasible albedo values, α(t,z) is clipped as suggested in Brock et al. (2000).

The “OerlemansModel” calculates melt energy as the residual term of a simplified surface energy balance equation (Oerlemans, 2001):
$$\begin{array}{}\text{(11)}& {a}_{\mathrm{sfc}}(t,z)={\displaystyle \frac{{Q}_{\text{m}}(t,z)\phantom{\rule{0.125em}{0ex}}\mathit{\delta}t}{{L}_{f}\phantom{\rule{0.125em}{0ex}}{\mathit{\rho}}_{\mathrm{w}}}},\end{array}$$where
$$\begin{array}{}\text{(12)}& \begin{array}{rl}{Q}_{\text{m}}(t,z)& =(\mathrm{1}\mathit{\alpha}(t,z\left)\right)\cdot G({I}_{\mathrm{pot}},t,z)\\ & +{c}_{\mathrm{0}}+{c}_{\mathrm{1}}\cdot T(t,z).\end{array}\end{array}$$Here, Q_{m}(t,z) (W m^{−2}) is the melt energy at time t and elevation z, δt=1 d is a time step, ${L}_{f}=\mathrm{3.34}\times {\mathrm{10}}^{\mathrm{5}}$ (J kg^{−1}) is the latent heat of fusion, and c_{0} (W m^{−2}) and c_{1} ($\mathrm{W}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}$) are empirical factors. The albedo α is calculated as well according to Eq. (10).
3.2 Mass balance model calibration
For the data assimilation procedure described in Sect. 3.3, we need a prior estimate for the model parameter values of Eqs. (5), (7), (8), (9), and (12). To obtain this, we calibrate the three investigated glaciers against the GLAMOS glacierwide mass balances (Sect. 2.3) and use an iterative procedure similar to Huss et al. (2009), illustrated in Fig. 4. Additionally, we calibrate the snow redistribution factor D(z) annually.
Following Huss et al. (2009), all model parameters are initially set to values reported in the literature (Hock, 1999; Oerlemans, 2001; Pellicciotti et al., 2005; Farinotti et al., 2012; Gabbi et al., 2014), and a twostep calibration procedure is then applied: first, the precipitation correction factor is tuned so that the winter mass balance of a given year is reproduced. In this step, the melt factors are held constant at their initial values. In a second step, the calibrated precipitation factor is kept constant, and the melt factors are optimized to reproduce the annual mass balance. The two steps are repeated alternately, and both precipitation correction and melt factors converge with every iteration. We terminate the iteration after the absolute difference to the winter and annual mass balance drops below 1 $\mathrm{mm}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$
Once the model parameters have been optimized, we determine the value of D(z) that matches the interpolated winter mass balance. Since this may result in changes in the required model parameters, the iterative procedure is applied one more time as a final step.
3.3 Particle filtering
To ensure that all mass balance model predictions stay within the observational uncertainty, we perform data assimilation. In particular, we employ a particle filter since it does not restrict the class of state transition models and error distributions. Especially when temperatures are around the melting point, the system becomes nonlinear since melt occurs above but not below this point. As a consequence, the distributions we deal with are not necessarily Gaussian. The facts that (a) the temperature chosen to parametrize the melting point is not the same for all four models, (b) the individual model prior distributions are combined to obtain the ensemble prediction, and (c) there can also be accumulation contributing to the overall mass balance add further complexity. We do not use other data assimilation approaches, such as variational methods or ensemble Kalman filtering, because variational methods encounter difficulties when dealing with nonGaussian priors (van Leeuwen et al., 2019), whilst the ensemble Kalman filter in its original form is not designed for multimodel applications as we use in our case. Overall, particle filtering is a very flexible, generalizable, and readily implementable data assimilation method.
Some extensions of the common particle filter framework allow model parameters and model performance to be estimated over time. With this, we aim at providing optimal, daily mass balance estimates at the glacier scale.
3.3.1 General framework
The general framework for data assimilation consists of a system whose state x_{t} evolves according to a model, but only partial and uncertain observations y_{t} of the state are available.
Here, x_{t−1} is the state at the previous time step t−1, g(⋅) is the state transition function, ℋ is the observation operator as introduced in Eq. (1), ϵ_{t} is the observation error vector at time t, and β_{t} is a random variable that describes model uncertainties. The term for β_{t} does not need to be strictly additive, and it should include uncertainties stemming from the model input variables. The goal of data assimilation is to compute conditional distributions of the system state x_{t} based on observations ${\mathit{y}}_{\mathrm{1}:t}=({\mathit{y}}_{\mathrm{1}},{\mathit{y}}_{\mathrm{2}},\mathrm{\dots},{\mathit{y}}_{t})$ sequentially for $t={t}_{\mathrm{0}},{t}_{\mathrm{0}}+\mathrm{1},\mathrm{\dots}{t}_{\mathrm{end}}$, where t_{0} and t_{end} are the time steps with the first and last observations, respectively. In our case, these conditional distributions describe the cumulative mass balance of a glacier, given all available camera observations.
To put this general framework into practice, we use the particle filter, which is a sequential Monte Carlo data assimilation method. Instead of handling conditional distributions of x_{t} analytically, the particle filter approximates the conditional distribution of a state x_{t} at time t given the observations y_{1:t} by a weighted sample of size N_{tot} (e.g., van Leeuwen et al., 2019).
Here δ(⋅) is the Dirac delta function, the elements x_{t,k} of the sample are called “particles”, and the weights w_{t,k} associated with the particles x_{t,k} sum to unity.
Usually, particle filtering comprises three repeated steps: the predict step, the update step, and the resampling step. In our case, these steps mean the following: during the predict step, particles holding possible mass balance states are propagated forward in time using the state transition in Eq. (13), where g(⋅) represents the ensemble prediction of mass balance Eqs. (5)–(10). This acts as a prior estimate of the mass balance distribution. In the update step, the weights of the propagated particles are recalculated based on Bayes' theorem. This accounts for the information of the next camera observation. In the last step, particles are resampled according to the updated weights. This step is necessary to restore particle diversity that is reduced during the update step. Resampling avoids socalled particle degeneracy, in which all weights collapse on only a few particles. Beyond the common threestep scheme, we additionally estimate model parameters with the particle filter by augmenting the state vector with model parameter values. In this way, we add an additional fourth step to the particle filter scheme, in which we evolve model parameters temporally according to a defined memory parameter. This prevents a collapse of the ensemble due to overconfidence, meaning that model parameter variability would become too low over time.
3.3.2 Application of the framework
The flowchart in Fig. 5 visualizes how the particle filter is implemented in our modeling framework. Figure 6 sheds light on how we perform the individual particle filter steps.
The temporal dynamics of the glacier mass balance state can be described by the accumulation model in Eq. (5) combined with the four different melt models in Eqs. (7), (8), (9), and (11). A priori, it is not known which model performs best, and each model has a set of unknown parameters. To take these two uncertainties into account, we augment the state vector by the model index ${m}_{t}\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}\mathit{\}}$ and the model parameters θ_{t}. In this way, a model weight and the model parameter values are also estimated based on the observations. Although the unknown parameters are different for each model, we do not use an additional model index for θ_{t}. Instead, we ensure that for all particles, θ_{t,k} is always the parameter vector associated with model m_{t,k}. This means that, when following a given particle backwards in time, its entire dynamics is governed by one single model only. In the forward direction, a particle can “die off” and be replaced by another particle with different model index. In this case, both the model index m_{t,k} and the entire past trajectory are changed to the new model.
As the state has to provide all information that is needed to predict the next observation, we also include the surface albedo and the snow water equivalent on the ice in our state vector. Hence, the state vector is defined as
Here, ξ_{t} is called the physical state.
3.3.3 Predict step
During the predict step, the explicit temporal evolution of the physical state ξ_{t} involves the randomized error accounting for uncertainties in the meteorological input (Sect. 2.2). Here, we call these errors η_{t} and set an additional scalar subscript to indicate that the errors are different for each meteorological variable. We first predict c_{sfc}(t,z), T_{acc}(t,z), and swe(t,z).
Based on Eqs. (7), (8), (9), and (11), the predicted mass balance is then
where the errors η_{t} are independent in time but partly perfectly correlated in space for reasons described in Sect. 2.2. By introducing both the meteorological uncertainty η and the parameter uncertainties, we shift the majority of the uncertainty contained in β_{t} to these variables. Since the remaining uncertainty for β_{t} is small and hard to quantify, we set β_{t}=0 for simplicity. With this assumption, we neglect some additional uncertainty contained in β_{t}, being aware that this might lead to “jumps” in the temporal evolution of the model performance. Finally, the observations y_{t} depend only on the cumulative mass balance at the elevation z of the camera as specified in Eq. (1).
We use a total of N_{tot}=10 000 particles and set the weights at time t_{0} (i.e., the time when the first camera observation is available) to $\mathrm{1}/{N}_{\mathrm{tot}}$. Since at t_{0} all models have equal probabilities, ${N}_{\mathrm{tot}}/\mathrm{4}$ particles are assigned to each of the four models. The initial value of b_{sfc}(t_{0},z) is set to zero for all particles, whereas α(t_{0},z) is determined by the maximum air temperature since the last snowfall before t_{0}, and swe(t_{0},z) depends on the cumulative mass balance before t_{0}. Finally, the initial calibration parameter values ${\mathit{\theta}}_{{t}_{\mathrm{0}},k}$ of the particles with model index j are obtained by drawing Monte Carlo samples from a normal distribution fitted to the logarithmized parameters of model j, as they were calibrated in the past (see Sect. 3.2). Table 2 shows the means and standard deviations for the input parameters of the three glaciers.
3.3.4 Update step
In the update step, all particles are then reweighted by multiplying the density of the observations y_{t} given the state of individual particles x_{t,k} with their respective weights at t−1 and normalizing the weights to sum to unity (van Leeuwen et al., 2019).
In our case, ${\mathit{y}}_{t}=h(t,z)$, and $p({\mathit{y}}_{t}\mid {\mathit{x}}_{t,k})$ is the normal density with mean ${b}_{sfc}(t,z{)}_{k}/{\mathit{\rho}}_{\mathrm{bulk}}$ and standard deviation σ_{ϵ}, evaluated at h(t,z). After updating the model predictions with the observations, we are interested in (a) the posterior model probabilities π_{t,j}, (b) the posterior distribution of model parameters θ, and of course (c) the posterior distribution of the physical state given all observations y_{1:t}. These quantities can be decomposed from the approximation with weighted particles in Eq. (14). The posterior model probability is given by
where π_{t,j} is the approximation of the posterior model probability at time t and model j. The posterior distribution of the parameters of model j is approximated by
The posterior distribution of the physical state takes the model uncertainty into account. It combines the posterior distributions under the different models j according to the law of total probability, in which we can insert Eqs. (21) and (22):
As the observations only measure the mass change since the installation of a camera, a difficulty occurs if several cameras are installed on different days at different elevations of the same glacier. We elaborate on the technical details for these cases in Appendix A.
3.3.5 Resampling
During the resampling step, the updated weights are used to choose a new set of N_{tot} particles with equal weights. To achieve equal weights, particles with low weights are removed, whereas those with high weights are duplicated. Because there is no stochasticity in the evolution of m_{t} though, when the particle index k is fixed, for some models only a few particles with the according model index survive after a couple of iterations. If this occurs, the respective model has little chance to become better represented at later time steps, which is unfavorable, since the model might give better predictions on average.
To overcome this problem, we assign a minimum contribution to each model of the ensemble, regardless of the model's performance at a certain time step. To compensate for the potentially too high resampling rate of a poor prediction, we lower the weights of all particles of a model whose contribution has been deliberately increased to match the chosen minimum contribution. In turn, we increase the weights of all other particles to compensate for their underrepresentation. This means that on average, the original weights per model remain unchanged. For technical details of the resampling procedure, see Appendix B.
3.3.6 Parameter evolution
The dynamics of the augmented state is defined such that the model index does not change over time. However, parameters are evolved temporally such that after a long period without observations, θ is distributed according to the prior parameter distribution.
Here, μ_{0} and Σ_{0} are the prior mean and the prior covariance of θ at the starting time t_{0}, which we determine from the calibration procedure in Sect. 3.2, and $\mathit{\rho}\in [\mathrm{0};\mathrm{1}]$ is a memory parameter that we choose to be 0.9. This step accounts for the fact that parameters are not necessarily constant in time, and it also ensures the reintroduction of parameter diversity which is lost during the resampling step.
3.4 Validation scores
To validate the daily mass balance predictions made with the particle filter, we use the CRPS (continuous ranked probability score). The CRPS is designed to estimate the deviation of a probabilistic forecast from an observation. It takes into account both the deviation of the median forecast from the actual observation (forecast reliability) and the spread of the forecast distribution (forecast resolution). This means that a forecast close to the observation median can still receive a poor CRPS if the forecast distribution spread is high, and the other way around. Lower values of the CRPS correspond to better forecasts. The minimum value is zero, corresponding to a perfect, deterministic forecast of the observation.
The CRPS is defined as (Hersbach, 2000)
where P_{f}(⋅) is the forecast mass balance cumulative probability distribution, and H(⋅) is the Heaviside function. The usual choice for P_{f} is the weighted ensemble distribution of the particles from the predict step, i.e., a discrete step function with jumps of height ${w}_{t\mathrm{1},k}$ at the positions ℋ(b_{sfc}(t,z)_{k}), where b_{sfc}(t,z)_{k} are the prediction particles. Note that this setting does not account for the observation error of h(t,z), implying that the score is not “proper”; i.e., it does not always return the best value when the prediction distribution is the true distribution (Ferro, 2017; Brehmer and Gneiting, 2019). To obtain a proper score, one can use the forecast of the camera reading h(t,z), which is the Gaussian mixture with weights ${w}_{t\mathrm{1},k}$, mean values ℋ(b_{sfc}(t,z)_{k}), and common variance ${\mathit{\sigma}}_{\mathit{\u03f5}}^{\mathrm{2}}$. Despite being proper, it still has some theoretical shortcomings (Ferro, 2017). Since for our data the values of the two scores do not differ much, we use only the proper score in all result figures but give also the value of the common CRPS in square brackets in the text.
4.1 Mass balance observations
Figure 7 shows the observed cumulative mass balance at the individual cameras, an example of meteorological conditions at station FIN 1 (providing the longest time series), daily mass balance rates at FIN 1, and four example camera images.
Considering all stations, we have observed ice melt rates of up to 0.12 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$ and a cumulative mass balance of about −5.5 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ in 81 d close to the terminus of Findelgletscher (FIN 1). Different camera stations reveal different melt rates and total ablation, which generally depend on the station's elevation. However, stations at different elevations can have similar melt rates as well. For example, station RHO 4 at 2589 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ experienced an average melt rate of −0.047 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$, while the average melt rate at FIN 2 (3015 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) was −0.043 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$ during the same period (we count only days with net ablation). Further, station PLM 1 had the lowest average melt rate despite not being at the highest elevation. This might be due to the meteorological conditions, such as the formation of local cold air pools, and the socalled “Massenerhebung effect” (Barry, 1992). The latter describes the tendency of higher temperatures to occur at the same elevation in the inner Alps as on their outer margins. For all stations, average melt rates during July and August are similar (0.073 ± 0.012 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$ in July vs. 0.062 ± 0.011 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$ in August) and 0.02–0.03 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$ lower (i.e., 0.044 ± 0.014 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) in September. On Glacier de la Plaine Morte, the difference is most pronounced with a drop of 0.06 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$ between August and September. Again, this is probably caused by local effects. On average, the difference between minimum and maximum melt measured at different stations on a particular day was 0.035 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$. Over the observational period, this difference ranged from 0.005 to 0.081 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$. The highest difference (0.081 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) occurred on 1 September 2019 in connection with the passage of a convergence line and cold front (German Meteorological Service, 2019): while Glacier de la Plaine Morte was already under the influence of cooler weather, Findelgletscher and Rhonegletscher experienced another meltintensive day. The variability at individual stations, measured as standard deviation of a 14 d running mean, was generally low during July and August (0.016 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) and increased at the beginning of September (0.026 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$). We attribute this increase to the onset of intermittent snowfalls at individual sites.
As shown by the pictures from station FIN 1 (Fig. 7d), summer 2019 is characterized by a variety of events, ranging from very hot, meltintensive days to snowfalls at high elevations. The time series of normalized mean daily temperature and shortwave radiation at station FIN 1 (Fig. 7b) illustrate that two heat waves occurred at the end of June and end of July 2019. The total amount of water released by snowmelt and ice melt on Swiss glaciers during these heat waves was 0.8 km^{3} (Swiss Academy of Sciences, 2019), which approximately equals the annual amount of drinking water consumed in the country. These extreme phases are also mirrored in the melt observed at our stations (Fig. 7): for FIN 1, daily melt rates peaked between 0.09 and 0.12 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$. For days with a rangenormalized temperature exceeding 0.8 (9 d in total, Fig. 7b), the average melt rate is 0.1 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$ at that station. During these nine days, modeled, glacierwide melt indicates the release of 6×10^{6} m^{3} of water. Another phase with very high melt occurred at the end of August 2019. Here, normalized temperature and radiation are average (mean values of 0.6 and 0.5, respectively). The exact causes for this strong melt event are unclear, and we speculate that it might be related (at least in part) to rain events that were not captured by the meteorological input despite being visible in our camera images between 28 and 31 August 2019. Summer melt phases were also interrupted by two snowfalls of different strengths: small amounts from 15 July 2019 (image 2 of Fig. 7d) and larger amounts, summing up to 0.25 m snow height in total, in early September (image 4 of Fig. 7).
4.2 Particle filter mass balance validation
Besides the direct observations presented above (Sect. 4.1), our framework enables predictions of daily mass balance. In this section, these predictions are (i) validated against reference forecasts (Sect. 4.2.1), (ii) crossvalidated against test subsets of observations (Sect. 4.2.2), and (iii) compared against glacierwide mass balances reported by GLAMOS (Sect. 4.2.3).
4.2.1 Validation against reference forecasts
We consider two types of reference forecasts: first, a forecast with (i) mean glacierwide melt parameters as obtained from past calibration (Sect. 3.2) and (ii) the precipitation correction factor prcp_{scale} constrained by the 2019 GLAMOS winter mass balance; and second, a forecast with a partially informed model including the same constraint for prcp_{scale} but also a tuning of the melt parameters to reproduce one further intermediate point measurement. The latter measurement is the cumulative mass balance between September 2018 and 2019 at the mass balance stake closest to each camera (locations in Fig. 1). Since there are up to four stake readings per glacier, we calculate single parameter sets tuned to reproduce all possible combinations of stake readings per glacier. This results in 19 CRPS values in total, for which we calculate the median. We also distinguish between the case in which the uncertainties in the meteorological inputs are taken into account and the case in which they are not.
Finally, we calculate the CRPS for the two reference forecasts by inserting two different values into the CRPS equation: (a) the mass balance of each day separately and (b) the cumulative mass balance. Note that for the particle filter, there is no need to make this distinction. Indeed, the daily deviation from a mass balance observation also equals the deviation from the cumulative observation. Figure 8 shows the results of the validation.
For the particle filter, daily and cumulative melt observations are generally reproduced well, with an average CRPS of 0.012 [0.012] m (proper CRPS outside and standard CRPS inside the square brackets). At the end of the assimilation period, Rhonegletscher has an average CRPS of 0.017 m, which is almost double the CRPS for Findelgletscher (CRPS = 0.01 m) and Glacier de la Plaine Morte (CRPS = 0.008 m). The high value of Rhonegletscher is related to the switchon of cameras RHO 1 and RHO 3. Indeed, the glacier also has CRPS ≈ 0.01 m before that. Poor predictive performances also occur after snowfalls, probably related to the uncertainties by which the mass balance stake can be read during these times. We run experiments in which the particle filter is limited to using mean parameters and/or single models instead of parameter distributions and the model ensemble. In the experiments, the resulting average CRPS values are higher than the average CRPS obtained with the ensemble and timevariant parameters. The lowest single values occur for specific combinations when running the particle filter with the BraithwaiteModel and OerlemansModel and flexible parameters on Glacier de la Plaine Morte. Note that if no probabilistic temperature and precipitation lapse rate is used, the resulting CRPS values from the experiments with mean parameters and/or only one model are even higher than the highest CRPS obtained using the ensemble and timevariant parameters. The experiments thus show that it is beneficial to include all four models and parameter uncertainty into the particle filter.
Comparing the CRPS of the particle filter with the reference forecasts, the performance closest to the particle filter is delivered by the forecast produced with mean melt parameters and no uncertainty in the meteorological input (mean CRPS = 0.013 [0.015] m). When the CRPS is calculated from the cumulative mass balance produced with mean melt parameters, the CRPS increases to 0.333 [0.243] m on average. This is because the mean parameters do not adapt to the meteorological conditions over time, and in this case, the cumulative mass balance can temporarily be under or overestimated or even diverge completely over time. Somewhat counterintuitively but for the same reason, the CRPS is similar when parameters have been tuned to match nearby stake readings. For the cumulative deviation, we find CRPS = 0.25 [0.251] $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ when considering meteorological uncertainty, and CRPS = 0.294 [0.28] $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ when not doing so. Compared to both the particle filter prediction and the prediction with mean melt parameters, the CRPS of daily mass balances produced without considering meteorological uncertainty is slightly higher (median CRPS: 0.023 [0.025] $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$).
In general and for the individual glaciers, the particle filter improves the CRPS of the reference forecasts by 95 % to 96 %. For the daily forecasts, the performance of the particle filter is only partly better, with improvements in CRPS between 8 % and 48 %. Along with the performance, a further important advantage of the particle filter is that it provides daily estimates for the results' uncertainties without need for further calculations. Indeed, this information can be essential, especially for the operational application of our framework.
4.2.2 Crossvalidation
A different approach for validating the particle filter is to only use subsets of the available camera observations as input and to evaluate the predicted mass balances at the remaining locations (crossvalidation). We do so by splitting the available observations into training and test subsets of cameras, i.e., by keeping the time series of a given station together (as opposed to splitting individual time series). Our test sets always contains one time series; i.e., we perform a leaveoneout crossvalidation. Figure 9 shows the temporal evolution of the CRPS at the test locations, i.e., at the stations not used by the particle filter.
We find that, in general, the cumulative mass balance at the test locations follows the cumulative observation curve well but not as closely as when the test location's data are assimilated with the particle filter. This shows the benefit of having several cameras per glacier mounted at different elevations. For Findelgletscher, we find 8.8 % average deviation (median CRPS of 0.071 and 0.108 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ for the two stations) when comparing the cumulative mass balance curve with the particle filter prediction. For Rhonegletscher, the average deviation at the test locations is 9.0 % average deviation (median CRPS 0.14, 0.148, 0.067, and 0.178 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ depending on the station). The highest CRPS values for Rhone stem from the period after the middle of August when two additional cameras were set up, but the values still outperform the reference forecasts of Sect. 4.2.1.
The temporal pattern evident in Fig. 9 includes an increasing CRPS through time but at different rates depending on the crossvalidation subset. The individual pattern originates from (1) a stations' representativity for the given elevation band it is located in, (2) the combination of stations in the crossvalidation subsets, and (3) cumulative error characteristics, since we observe cumulative mass balance over time. Station RHO 3, for example, can generally be modeled with lower errors compared to other stations. We speculate this being related to its location, which is in a relatively flat area with little crevasses. The other stations are instead either in the vicinity of crevasses (RHO 4) or influenced by shadows from the surrounding terrain, dark glacier surface, or steep ice (RHO 1 and RHO 2). RHO 1 and RHO 2 also show that even neighboring stations can exhibit different melt. This affects the results of the crossvalidation whenever one of these two stations is excluded from the training dataset.
The above results show the ability of the particle filter to also predict melt at locations without observations, albeit with a lower performance when compared to the situation in which all observations are assimilated. The results also show that even with an augmented particle filter, it is demanding to find a unique, glacierwide parameter set that correctly reproduces the mass balance at all locations.
4.2.3 Comparison to GLAMOS glacierwide mass balances
We compare our assimilated model ensemble predictions to the glacierwide annual mass balance reported by GLAMOS at the autumn field date of the mass budget year 2019. We do so by running the model from the field campaign date in autumn 2018. Figure 10 illustrates the different model and parameter settings used during the simulation.
During the period preceding the installation of our cameras, we calculate mass balance with the parameters calibrated in Sect. 3.2. This results in about 45 distinct model runs, which we call “free model runs”. We use this first period to provide initial conditions for the particle filter period, which lasts from the first camera setup on a respective glacier either until cameras are retrieved or until the autumn field date is reached (whatever comes first). To achieve a connection between the free model run and the period during which the particle filter is used, we sample 10 000 times from the initial conditions at the first camera setup date. We refer to this procedure as “particle filtering without preselection (of initial conditions)”. Not all free model runs have to be used, though: they can also be preselected based on the cumulative mass balance observed at the stakes closest to the camera stations. For this case, we select model runs that reproduce these observations within an estimated reading uncertainty of ± 0.05 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ (“particle filtering with preselection”). The cumulative mass balances calculated with these two procedures are compared to the GLAMOS analyses in Table 3.
For particle filtering without preselection of initial conditions, the difference from the GLAMOS analyses is 0.67 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ for Rhonegletscher, 0.2 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ for Findelgletscher, and 0.05 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ for Plaine Morte. With preselection, in contrast, the absolute difference changes by −0.07, −0.24, and +0.02 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$, respectively. Consequently, including the stake mass balance readings improves the match to the GLAMOS analyses for Rhonegletscher and Plaine Morte, while it has only little effect for Findelgletscher. A reason for this can be either that the mass balance stakes are not at the observation locations or that the mass balance gradients of the preselected runs are unfavorable. Overall, the differences from the GLAMOS analyses can be explained by (1) the difference in the approaches used to calculate glacierwide mass balances from point observations, (2) the use of only 1–4 point observations located in the ablation zone and covering <30 % of the glacier elevation range, compared to the complete network of 5–14 stakes over the entire elevation range used in the GLAMOS analyses, (3) the lack of representativeness of the camera observations for the accumulation zone of the glaciers, i.e., biased vertical mass balance gradients, (4) the lack of representation of individual winter accumulation measurements in our glacier model, or (5) a problem with representing the mass balance of the glacier with only one parameter set. Also note that 91 %–99 % of the total uncertainty for the model runs with data assimilation stem from the period before the particle filter can be initialized, i.e., before the installation of the first camera station. Figure C1 in the Appendix shows the evolution of the mass balance state over the assimilation period by the example of Findelgletscher.
4.3 Individual model performance
We analyze model performance by considering the temporal evolution of the model probabilities π_{t,j} and model particle numbers N_{t,j} for the four melt models over time at individual glaciers. High model performance is indicated by high probabilities and large particle numbers over long time periods.
Figure 11 shows the model performance of all four melt models and at all three glaciers.
In general, we find that the model probabilities are sensitive to the ensemble input, such as the parameter priors, and the prescribed meteorological uncertainty. This is an indication of the ensemble choosing the model combination that best reproduces the observations at any time. Note that none of the models is removed from the ensemble in the resampling step even when the model performs poorly. During the assimilation period, indeed, models can recover and can show good performances at a later stage (see, for example, the HockModel for Rhonegletscher or the PellicciottiModel for Findelgletscher). This shows the utility of the resampling procedure introduced in Sect. 3.3.5.
During the assimilation period of an individual glacier, often one model dominates the ensemble for a given amount of time (“model dominance” being the case in which the model probability is >0.5). Model dominance, and especially fast switches between dominant models, can be indicative of a mode collapse, resulting from either an overconfident likelihood or prior, or both, operating in an Mopen framework (Bernardo and Smith, 2009), i.e., the case in which the “true” model is not a choice amongst the available models. In our case, we believe that the ensemble prior might be overconfident on average since we have chosen the observational error conservatively; i.e., we have chosen the largest errors emerging in the roundrobin experiment (Sect. 2.1.2). This would lead to a model preferably obtaining high weights, which had already dominated on the previous days. However, when the likelihood is overconfident or there is strong evidence that a previously wellperforming model now performs worse, the filter might switch back and forth between individual models that best describe the observations. We accept this model dominance and the fast switching as a sign that the overall ensemble performance is improved. Averaged over all glacier and time steps, the PellicciottiModel has the highest model probability (0.39), while the BraithwaiteModel has an average model probability of 0.24, the OerlemansModel of 0.23, and the HockModel of 0.14. The relatively high probabilities assigned to the PellicciottiModel can have various reasons, and we suspect that two are of particular importance in our case: first, the calibration might have led to a broad prior parameter distribution, allowing for the model to adapt to various combinations of meteorological input and observed melt. Second, using the actual solar radiation G instead of the potential irradiation I_{pot} might provide a further advantage since this accounts for partly cloudy conditions and diffuse radiation, which the potential irradiation is not able to cover. The fact that the second highest probability is assigned to the OerlemansModel (which uses G as well), supports this possible explanation.
In terms of the temporal evolution, the model dominance for Rhonegletscher and Glacier de la Plaine Morte is determined already within the first few days and changes only a little after that. Changes in model dominance can be observed for Rhonegletscher and Findelgletscher, instead. In the case of Rhonegletscher, for example, the model dominance switches from the PellicciottiModel to the OerlemansModel and later to the HockModel. For Findelgletscher, however, there is a transition from the OerlemansModel to the PellicciottiModel. This transition is particularly noticeable between 7 and 8 August 2019 (Fig. 12). The causes for it are not entirely clear, and we speculate that it might be related to the precipitation event starting on 6 August. Perhaps surprisingly, the model dominance seems to be little influenced by snowfall events (e.g., from 9 to 17 September on Findelgletscher or from 5 September to 11 September on Rhonegletscher), even if surface albedo is taken into account very differently by the individual models.
Figure 13 shows the evolution of the distribution of individual model parameters during the assimilation period. The example refers to Findelgletscher.
Three phases of quick parameter changes can be observed. First, the parameters change rapidly on the first days of the assimilation period. This means that the prior parameter distributions do not match the exact parameter distributions needed to model the mass balance at the camera locations. This is due to both the calibration time span (seasonal calibration vs. daily application) and the low sample size of the calibrated parameters. A second rapid change can be observed after the second camera has been switched on, i.e., on 24 July 2019. Here, an adjustment in the parameters is needed in order to accommodate the mass balance at both stations equally well. The third rapid change starts when ablation at station FIN 1 is highest but when radiation and temperature are not at their maximum. Here, the change might be due to the model being forced to yield high ablation rates despite only moderate meteorological forcing. This shows the advantage of employing the model ensemble as opposed to, for example, a single model with deterministic parameters: the ensemble also reproduces system states which cannot be explained by the uncertain meteorological input.
In this study, we mounted seven cameras on three Swiss glaciers, delivering 352 point mass balance observations throughout summer 2019. At the camera locations, we observed daily melt rates of up to 0.12 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$ and cumulative melt of up to ∼5 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ in 81 d. To calculate nearrealtime mass balances, we used an ensemble of three TI models, a simplified energy balance model, and meteorological model input. The camera observations were assimilated into the model ensemble by using a specifically developed particle filtering scheme. The particular focus was put on delivering a stable ensemble capable of reproducing glacier mass balance throughout the summer. Variability in the model parameters, as well as particle filter stability, was considered. For the former, a prior parameter distribution obtained from calibration against past seasonal glaciological mass balances was used as input to an augmented particle filter capable of estimating parameters while assimilating observations. For the latter, the particle filter was designed such that models with temporarily poor performance can recover at a later stage. For the mass budget year 2019, we calculate cumulative mass balances of −1.79 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$, −0.48 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$, and −0.84 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$ for Glacier de la Plaine Morte, Findelgletscher, and Rhonegletscher, respectively.
The mass balances given by the particle filter were closer to the cumulative observations (CRPS = 0.012 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{w}.\mathrm{e}.$) than two reference forecasts which either assumed no measurements to be available or only used one intermediate set of stake readings. Measured with the CRPS for cumulative mass balances, the particle filter improves the performance of reference forecasts by 95 % to 96 %. As a further advantage, the particle filter delivers direct uncertainty estimates. A leaveoneout crossvalidation procedure showed that the cumulative mass balance predicted with the particle filter is within 9 % of the observations at any location. In an analysis of the individual model performance, we found that our technique to prevent models from being removed from the ensemble is useful since models can recover at a later stage. In terms of model ensemble, the TI model by Pellicciotti et al. (2005) obtained the highest average model probability (0.39). None of the four models has an average probability <10 %, and even if individual models can temporarily perform poorly, our technique preventing models from being removed from the ensemble completely allows them to recover at a later stage. Fast temporal switches between model probabilities are attributed to overconfident likelihood or prior distributions, or both. As a future venture, we envision an extension of the particle filter in which glacier mass balances and model parameters are further constrained by remotely sensed observations of albedo and snow lines. These measurements are indirect but have the potential to (1) complement the camera observations extensively and to (2) overcome the limited knowledge about the spatial and temporal extrapolation of glacier mass balances and model parameters.
Assume that camera i is installed at elevation z_{i} on day t_{i−1}, where ${t}_{\mathrm{0}}<{t}_{\mathrm{1}}<{t}_{\mathrm{2}}\mathrm{\dots}$ (to be coherent with earlier notation that the first camera is installed at time t_{0}). From time t_{i−1} onwards, we include ${b}_{\mathrm{sfc}}({t}_{i\mathrm{1}},{z}_{i})$ in the state vector as a component which remains constant. Then the observations at time $t>{t}_{i\mathrm{1}}$ are functions of the state at time t:
The true value of ${b}_{\mathrm{sfc}}({t}_{i\mathrm{1}},{z}_{i})$ is unknown, and the uncertainty is represented by the values ${b}_{\mathrm{sfc},k}({t}_{i\mathrm{1}},{z}_{i})$ of the particles. Thus at time t, the contribution from the observation h(t,z_{i}) to the weight of particle k is proportional to
Although ${b}_{\mathrm{sfc},k}({t}_{i\mathrm{1}},{z}_{i})$ never changes during the propagation step, it will change in the resampling steps. Thus the uncertainty about ${b}_{\mathrm{sfc}}({t}_{i\mathrm{1}},{z}_{i})$ will decrease as time proceeds. This is presumably not realistic, but the effect of small errors in the baseline also diminishes as time proceeds.
The technical details of the resampling procedure in Sect. 3.3.5 are the following: if, after prediction and update, N_{t,j} denotes the number of particles with model index j, we prevent models from not being resampled by choosing a minimum model contribution $\mathit{\varphi}<\frac{\mathrm{1}}{\mathrm{4}}$ to the ensemble. This ensures that the resampling step preserves a minimum particle number ${N}_{t,j}\ge \mathit{\varphi}{N}_{\mathrm{tot}}$ representing model j. For our application, we choose ϕ=0.1. If the posterior probability of model j (Eq. 21) is smaller than the minimum contribution ϕ, an unweighted sample that represents π_{t,j} correctly must have less than ϕN_{tot} particles with model index j. To ensure our minimum contribution condition though, we generate a weighted sample $({\stackrel{\mathrm{\u0303}}{\mathit{x}}}_{t,k},{\stackrel{\mathrm{\u0303}}{w}}_{t,k})$ such that each model index j appears at least ϕN_{tot} times, and the weights are as close to uniform as possible. We select the particles ${\stackrel{\mathrm{\u0303}}{\mathit{x}}}_{t,k}$ in a twostep resampling procedure: first, the number N_{t,j} of particles with model index j is chosen to be ${N}_{t,j}=\mathit{\varphi}{N}_{\mathrm{tot}}+{L}_{t,j}$, where L_{t,j} are excess frequencies. We obtain these frequencies by sampling a total of N_{tot}(1−4ϕ) model indices from $\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}\mathit{\}}$ with weights proportional to how much a model probability exceeds the chosen minimum contribution, i.e., $max(\mathrm{0},{\mathit{\pi}}_{t,j}\mathit{\varphi})$. In a second step, we draw for each model a resample of size N_{t,j} with weights ${w}_{t,k}/{\mathit{\pi}}_{t,j}$ from the particles with model index j. The combined set of the N_{tot} resampled particles gives the new filter particles ${\stackrel{\mathrm{\u0303}}{\mathit{x}}}_{t,k}$.
However, introducing a restriction on the minimum number of particles per model can lead to biased estimates as poor models with probability ${\mathit{\pi}}_{t,j}\le \mathit{\varphi}$ are overrepresented in the ensemble. To compensate for poor models occurring too often among the resampled particles (and the other models not often enough), the following weight has to be given to ${\stackrel{\mathrm{\u0303}}{\mathit{x}}}_{t,k}$:
These weights sum to unity and preserve the original weights w_{t,k} on average. Since they can become very small though, we work with the logarithm of the weights to avoid numerical underflow. It should be noted that we insert ${\stackrel{\mathrm{\u0303}}{w}}_{t\mathrm{1},k}$ for ${w}_{t\mathrm{1},k}$ in Eqs. (14) and (20). In order to see that the weights we choose for ${\stackrel{\mathrm{\u0303}}{\mathit{x}}}_{t,k}$ are correct, we denote the number of times the particle x_{t,k} is selected in the resampling procedure by ${\stackrel{\mathrm{\u0303}}{M}}_{t,k}$. This means that the resampling gives x_{t,k} the random weight $\frac{{\stackrel{\mathrm{\u0303}}{M}}_{t,k}}{{N}_{\mathrm{tot}}}$, which is then multiplied by the additional weight ${\stackrel{\mathrm{\u0303}}{w}}_{t,k}$. Hence, x_{t,k} receives the total weight
If ${m}_{t,k}=j$, it holds that
i.e., on average the new weights ${w}_{t,k}^{\prime}$ are equal to the original weights.
The camera observations are available under the following DOI: https://doi.org/10.3929/ethzb000508515 (Landmann, 2021). The meteorological data can be obtained as a paid service from https://www.meteoschweiz.admin.ch/home/klima/schweizerklimaimdetail/raeumlicheklimaanalysen.html (last access: 28 October 2021), and the glacier outlines and mass balances are available free of charge from the GLAMOS web site at https://doi.glamos.ch/data/inventory/inventory_sgi2010_r2010.zip (last access: 28 October 2021) and https://doi.glamos.ch/data/massbalance/massbalance_observation_elevationbins.csv (last access: 28 October 2021). The code used to produce results and figures can be obtained from the authors upon request.
Time lapse videos of all camera observations used in this study are available as videos under the following DOIs: PLM 1: https://doi.org/10.5446/48826 (Landmann, 2020a); FIN 1: https://doi.org/10.5446/48824 (Landmann, 2020b); FIN 2: https://doi.org/10.5446/48825 (Landmann, 2020c); RHO 1: https://doi.org/10.5446/48820 (Landmann, 2020d); RHO 2: https://doi.org/10.5446/48821 (Landmann, 2020e); RHO 3: https://doi.org/10.5446/48822 (Landmann, 2020f); RHO 4: https://doi.org/10.5446/48823 (Landmann, 2020g).
JML had the particle filter idea, implemented all models, did all figures, and wrote the paper. HRK supervised the particle filter methodology, brought in the method to prevent models from disappearing from the ensemble, and reviewed the paper. MH commented on the method, reviewed the paper, and mounted some of the stations. CO prepared and mounted most of the stations. MK commented on the particle filter and reviewed the paper. DF did the overall supervision, proposed to use data assimilation in JML's doctorate, commented on the method, reviewed the paper, and acquired the funding.
The authors declare that they have no conflict of interest.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We would like to acknowledge the funding that we got from GCOS (Global Climate Observing System) Switzerland and the extensive support that we got from the manufacturer of the cameras and transmitter boxes, Holfuy Ltd (in particular Gergely Mátyus). We would like to thank the teachers of the Joint ECMWF and University of Reading Data Assimilation Training Course that helped to significantly improve Johannes Marian Landmann's knowledge on data assimilation, in particular Javier Amezcua. Further, we would like to thank Anastasia Sycheva and Emmy Stigter for testreading the methods chapter for reader friendliness. We appreciate the help from all people that conducted the field work apart from the authors, namely Małgorzata Chmiel, Amaury Dehecq, Lea Geibel, Katja Henz, Serafine Kattus, Johanna Klahold, Claudia Kurzböck, Amandine Sergeant, and Michaela Wenner, and we thank all people that took part in the roundrobin experiment apart from the authors, namely Amaury Dehecq, Eef van Dongen, Elias Hodel, Jane Walden, and Michaela Wenner. Kudos to Dominik Gräff for helping to “TCify” Sect. 3.4 after the review comment by Douglas Brinkerhoff. We would also like to thank Bertrand Cluzet and colleagues for having changed the acronym of their project which coincidentally was the same as ours (CRAMPON).
This paper was edited by Marie Dumont and reviewed by Douglas Brinkerhoff and one anonymous referee.
Arulampalam, M. S., Maskell, S., Gordon, N., and Clapp, T.: A tutorial on particle filters for online nonlinear/nonGaussian Bayesian tracking, IEEE T. Signal Proces., 50, 174–188, https://doi.org/10.1109/78.978374, 2002. a
Barry, R.: Mountain Weather and Climate, Physical Environment Series, Routledge, 1992. a
Bauder, A., Funk, M., and Huss, M.: Icevolume changes of selected glaciers in the Swiss Alps since the end of the 19th century, Ann. Glaciol., 46, 145–149, https://doi.org/10.3189/172756407782871701, 2007. a
Beniston, M., Farinotti, D., Stoffel, M., Andreassen, L. M., Coppola, E., Eckert, N., Fantini, A., Giacona, F., Hauck, C., Huss, M., Huwald, H., Lehning, M., LópezMoreno, J.I., Magnusson, J., Marty, C., MoránTejéda, E., Morin, S., Naaim, M., Provenzale, A., Rabatel, A., Six, D., Stötter, J., Strasser, U., Terzago, S., and Vincent, C.: The European mountain cryosphere: a review of its current state, trends, and future challenges, The Cryosphere, 12, 759–794, https://doi.org/10.5194/tc127592018, 2018. a
Bernardo, J. M. and Smith, A. F.: Bayesian theory, vol. 405, John Wiley & Sons, Hoboken, New Jersey, USA, 2009. a
Beven, K.: Environmental modelling: An uncertain future?, Routledge, New York, 2009. a
Biron, R. and Rabatel, A.: SmartStake: an autonomous measurement station for high resolution glacier ablation monitoring, Data sheet, available at: https://bubblyflow.com/wpcontent/uploads/2021/02/SmartStakev062020.pdf (last access: 28 October 2021), 2019. a
Bonan, B., Nodet, M., Ritz, C., and Peyaud, V.: An ETKF approach for initial state and parameter estimation in ice sheet modelling, Nonlin. Processes Geophys., 21, 569–582, https://doi.org/10.5194/npg215692014, 2014. a
Braithwaite, R. J.: Positive degreeday factors for ablation on the Greenland ice sheet studied by energybalance modelling, J. Glaciol., 41, 153–160, https://doi.org/10.3189/S0022143000017846, 1995. a
Braithwaite, R. J. and Olesen, O. B.: Calculation of glacier ablation from air temperature, West Greenland, in: Glacier fluctuations and climatic change, edited by: Oerlemans, J., Springer, Amsterdam, 219–233, https://doi.org/10.1007/9789401578233, 1989. a
Brehmer, J. R. and Gneiting, T.: Properization: constructing proper scoring rules via Bayes acts, Annals of the Institute of Statistical Mathematics, pp. 1–15, 2019. a
Brock, B. W., Willis, I. C., and Sharp, M. J.: Measurement and parameterization of albedo variations at Haut Glacier d'Arolla, Switzerland, J. Glaciol., 46, 675–688, https://doi.org/10.3189/172756500781832675, 2000. a, b, c
Carturan, L., Cazorzi, F., dalla Fontana, G., and Zanoner, T.: Automatic measurement of glacier ice ablation using thermistor strings, J. Glaciol., 65, 188–194, https://doi.org/10.1017/jog.2018.103, 2019. a
Cogley, J., Hock, R., Rasmussen, L., Arendt, A., Bauder, A., Braithwaite, R., Jansson, P., Kaser, G., Möller, M., Nicholson, L., and Zemp, M.: Glossary of glacier mass balance and related terms, IHPVII technical documents in hydrology No. 86, IACS Contribution No. 2, Tech. rep., International Association of Cryospheric Sciences (IACS), Paris, 2011. a, b
Corripio, J. G.: Vectorial algebra algorithms for calculating terrain parameters from DEMs and solar radiation modelling in mountainous terrain, Int. J. Geogr. Inf. Sci., 17, 1–23, https://doi.org/10.1080/713811744, 2003. a
Dumont, M., Durand, Y., Arnaud, Y., and Six, D.: Variational assimilation of albedo in a snowpack model and reconstruction of the spatial massbalance distribution of an alpine glacier, J. Glaciol., 58, 151–164, https://doi.org/10.3189/2012JoG11J163, 2012. a
Eis, J., Maussion, F., and Marzeion, B.: Initialization of a global glacier model based on presentday glacier geometry and past climate information: an ensemble approach, The Cryosphere, 13, 3317–3335, https://doi.org/10.5194/tc1333172019, 2019. a
Elsberg, D. H., Harrison, W. D., Echelmeyer, K. A., and Krimmel, R. M.: Quantifying the effects of climate and surface change on glacier mass balance, J. Glaciol., 47, 649–658, https://doi.org/10.3189/172756501781831783, 2001. a
Euronews: From Siberia to Switzerland, scorching August leads to more fires, less ice, available at: https://bit.ly/2BLLrfU (last access: 28 October 2021), 2019. a
Farinotti, D., Magnusson, J., Huss, M., and Bauder, A.: Snow accumulation distribution inferred from timelapse photography and simple modelling, Hydrol. Process., 24, 2087–2097, https://doi.org/10.1002/hyp.7629, 2010. a
Farinotti, D., Usselmann, S., Huss, M., Bauder, A., and Funk, M.: Runoff evolution in the Swiss Alps: projections for selected highalpine catchments based on ENSEMBLES scenarios, Hydrol. Process., 26, 1909–1924, https://doi.org/10.1002/hyp.8276, 2012. a, b
Fausto, R. S., Van As, D., Ahlstrøm, A. P., and Citterio, M.: Assessing the accuracy of Greenland ice sheet ice ablation measurements by pressure transducer, J. Glaciol., 58, 1144–1150, 2012. a
Ferro, C. A.: Measuring forecast performance in the presence of observation error, Q. J. Roy. Meteor. Soc., 143, 2665–2676, 2017. a, b
Fettweis, X., Franco, B., Tedesco, M., van Angelen, J. H., Lenaerts, J. T. M., van den Broeke, M. R., and Gallée, H.: Estimating the Greenland ice sheet surface mass balance contribution to future sea level rise using the regional atmospheric climate model MAR, The Cryosphere, 7, 469–489, https://doi.org/10.5194/tc74692013, 2013. a
Fischer, M., Huss, M., and Hoelzle, M.: Surface elevation and mass changes of all Swiss glaciers 1980–2010, The Cryosphere, 9, 525–540, https://doi.org/10.5194/tc95252015, 2015. a
Frei, C.: Interpolation of temperature in a mountainous region using nonlinear profiles and nonEuclidean distances, Int. J. Climatol., 34, 1585–1605, https://doi.org/10.1002/joc.3786, 2014. a
Gabbi, J., Carenzo, M., Pellicciotti, F., Bauder, A., and Funk, M.: A comparison of empirical and physically based glacier surface melt models for longterm simulations of glacier response, J. Glaciol., 60, 1140–1154, https://doi.org/10.3189/2014JoG14J011, 2014. a, b
German Meteorological Service: Frontal Analysis Europe 20190901, available at: http://www1.wetter3.de/archiv_dwd_dt.html (last access: 28 October 2021), 2019. a
Glacier Monitoring Switzerland: Swiss Glacier Mass Balance, release 2018, GLAMOS Data [data set], https://doi.org/10.18750/massbalance.2018.r2018, 2018. a
GLAMOS: GLAMOS web page, Web site, available at: https://www.glamos.ch/en/, last access: 6 August 2020. a
Golledge, N. R.: Longterm projections of sealevel rise from ice sheets, WIREs Clim. Change, 11, e634, https://doi.org/10.1002/wcc.634, 2020. a
Gugerli, R., Salzmann, N., Huss, M., and Desilets, D.: Continuous and autonomous snow water equivalent measurements by a cosmic ray sensor on an alpine glacier, The Cryosphere, 13, 3413–3434, https://doi.org/10.5194/tc1334132019, 2019. a
Hersbach, H.: Decomposition of the Continuous Ranked Probability Score for Ensemble Prediction Systems, Weather Forecast., 15, 559–570, https://doi.org/10.1175/15200434(2000)015<0559:DOTCRP>2.0.CO;2, 2000. a
Hock, R.: A distributed temperatureindex iceand snowmelt model including potential direct solar radiation, J. Glaciol., 45, 101–111, https://doi.org/10.3189/S0022143000003087, 1999. a, b, c, d, e
Hock, R.: Temperature index melt modelling in mountain areas, J. Hydrol., 282, 104–115, 2003. a, b, c
Hock, R., Jansson, P., and Braun, L. N.: Modelling the Response of Mountain Glacier Discharge to Climate Warming, 243–252, Springer Netherlands, Dordrecht, 2005. a
Hock, R., Bliss, A., marzeion, b., Giesen, R. H., Hirabayashi, Y., Huss, M., Radic, V., and Slangen, A. B. A.: GlacierMIP – A model intercomparison of globalscale glacier massbalance models and projections, J. Glaciol., 65, 453–467, https://doi.org/10.1017/jog.2019.22, 2019. a
Hulth, J.: Using a drawwire sensor to continuously monitor glacier melt, J. Glaciol., 56, 922–924, https://doi.org/10.3189/002214310794457290, 2010. a
Huss, M., Farinotti, D., Bauder, A., and Funk, M.: Modelling runoff from highly glacierized alpine drainage basins in a changing climate, Hydrol. Process., 22, 3888–3902, https://doi.org/10.1002/hyp.7055, 2008. a
Huss, M., Bauder, A., and Funk, M.: Homogenization of longterm massbalance time series, Ann. Glaciol., 50, 198–206, https://doi.org/10.3189/172756409787769627, 2009. a, b, c, d
Huss, M., Hock, R., Bauder, A., and Funk, M.: Conventional versus referencesurface mass balance, J. Glaciol., 58, 278–286, https://doi.org/10.3189/2012JoG11J216, 2012. a
Huss, M., Dhulst, L., and Bauder, A.: New longterm massbalance series for the Swiss Alps, J. Glaciol., 61, 551–562, https://doi.org/10.3189/2015JoG15J015, 2015. a, b
Hydrique: Example hydrological nowcast, available at: https://fribourg.swissrivers.ch/appSite/index/site/fribourg, last access: 6 October 2020. a
Iqbal, M.: An introduction to solar radiation, Academic Press, https://doi.org/10.1016/B9780123737502.X50010, 1983. a
Isotta, F. A., Frei, C., Weilguni, V., Perčec Tadić, M., Lassègues, P., Rudolf, B., Pavan, V., Cacciamani, C., Antolini, G., Ratto, S. M., Munari, M., Micheletti, S., Bonati, V., Lussana, C., Ronchi, C., Panettieri, E., Marigo, G., and Vertačnik, G.: The climate of daily precipitation in the Alps: development and analysis of a highresolution grid dataset from panAlpine raingauge data, Int. J. Climatol., 34, 1657–1675, https://doi.org/10.1002/joc.3794, 2014. a
Isotta, F. A., Begert, M., and Frei, C.: Longterm consistent monthly temperature and precipitation grid datasets for Switzerland over the past 150 years, J. Geophys. Res.Atmos., 124, 3783–3799, 2019. a
Jouvet, G., Huss, M., Funk, M., and Blatter, H.: Modelling the retreat of Grosser Aletschgletscher, Switzerland, in a changing climate, J. Glaciol., 57, 1033–1045, https://doi.org/10.3189/002214311798843359, 2011. a
Keeler, M. L. and Brugger, K. A.: A method for recording ice ablation using a lowcost ultrasonic rangefinder, J. Glaciol., 58, 565–568, 2012. a
Kreucher, C., Hero, A., and Kastella, K.: Multiple model particle filtering for multitarget tracking, in: Proceedings of the Twelfth Annual Workshop on Adaptive Sensor Array Processing, 16–18 March 2004, Massachusetts, USA, 2004. a
Landmann, J. M.: Time lapse video melt at Holfuy station PLM1, TIB AVPortal [video supplement], https://doi.org/10.5446/48826, 2020a. a
Landmann, J. M.: Time lapse video melt at Holfuy station FIN1, TIB AVPortal [video supplement], https://doi.org/10.5446/48824, 2020b. a
Landmann, J. M.: Time lapse video melt at Holfuy station FIN2, TIB AVPortal [video supplement], https://doi.org/10.5446/48825, 2020c. a
Landmann, J. M.: Time lapse video melt at Holfuy station RHO1, TIB AVPortal [video supplement], https://doi.org/10.5446/48820, 2020d. a
Landmann, J. M.: Time lapse video melt at Holfuy station RHO2, TIB AVPortal [video supplement], https://doi.org/10.5446/48821, 2020e. a
Landmann, J. M.: Time lapse video melt at Holfuy station RHO3, TIB AVPortal [video supplement], https://doi.org/10.5446/48822, 2020f. a
Landmann, J. M.: Time lapse video melt at Holfuy station RHO4, TIB AVPortal [video supplement], https://doi.org/10.5446/48823, 2020g. a
Landmann, J. M.: Glacier mass balance stake readings and videos from automated realtime cameras in summer 2019, ETH Zürich [data set], https://doi.org/10.3929/ethzb000508515, 2021. a
Lang, H. and Braun, L.: On the information content of air temperature in the context of snow melt estimation, IAHS Publ., 190, 347–354, 1990. a
Leclercq, P., Aalstad, K., Elvehøy, H., and Altena, B.: Modelling of glacier surface mass balance with assimilation of glacier mass balance and snow cover observations from remote sensing, in: EGU General Assembly Conference Abstracts, vol. 19 of EGU General Assembly Conference Abstracts, p. 17591, 2017. a
Magnusson, J., Winstral, A., Stordal, A. S., Essery, R., and Jonas, T.: Improving physically based snow simulations by assimilating snow depths using the particle filter, Water Resour. Res., 53, 1125–1143, https://doi.org/10.1002/2016WR019092, 2017. a
Marzeion, B., Hock, R., Anderson, B., Bliss, A., Champollion, N., Fujita, K., Huss, M., Immerzeel, W., Kraaijenbrink, P., Malles, J.H., Maussion, F., Radić, V., Rounce, D. R., Sakai, A., Shannon, S., van de Wal, R., and Zekollari, H.: Partitioning the Uncertainty of Ensemble Projections of Global Glacier Mass Change, Earth's Future, https://doi.org/10.1029/2019EF001470, 2020. a
Maussion, F., Butenko, A., Champollion, N., Dusch, M., Eis, J., Fourteau, K., Gregor, P., Jarosch, A. H., Landmann, J., Oesterle, F., Recinos, B., Rothenpieler, T., Vlug, A., Wild, C. T., and Marzeion, B.: The Open Global Glacier Model (OGGM) v1.1, Geosci. Model Dev., 12, 909–931, https://doi.org/10.5194/gmd129092019, 2019. a, b
MeteoSwiss: Daily, monthly and yearly satellitebased global radiation, Tech. rep., MeteoSwiss, available at: https://www.meteoswiss.admin.ch/content/dam/meteoswiss/en/climate/swissclimateindetail/doc/ProdDoc_SIS.pdf (last access: 28 October 2021), 2018. a
MeteoSwiss: Documentation of MeteoSwiss GridData Products: Daily Mean, Minimum and Maximum Temperature: TabsD, TminD, TmaxD, available at: https://www.meteoswiss.admin.ch/content/dam/meteoswiss/de/serviceundpublikationen/produkt/raeumlichedatentemperatur/doc/ProdDoc_TabsD.pdf last access: September 2021a. a
MeteoSwiss: Daily Precipitation (final analysis): RhiresD, Tech. rep., MeteoSwiss, available at: https://www.meteoswiss.admin.ch/content/dam/meteoswiss/de/serviceundpublikationen/produkt/raeumlichedatenniederschlag/doc/ProdDoc_RhiresD.pdf (last access: 28 October 2021), 2021b. a, b
Müller, H. and Kappenberger, G.: ClaridenfirnMessungen 19141984: Daten und Ergebnisse eines gemeinschaftlichen Forschungsprojektes, Verlag d. Fachvereine, Zürich, 1991. a
Netto, G. T. and ArigonyNeto, J.: Opensource Automatic Weather Station and Electronic Ablation Station for measuring the impacts of climate change on glaciers, HardwareX, 5, e00053, https://doi.org/10.1016/j.ohx.2019.e00053, 2019. a
NSIDC: Greenland Ice Sheet Today, available at: https://nsidc.org/greenlandtoday/, last access: 4 September 2020a. a
NSIDC: Snow Today, available at: https://nsidc.org/snowtoday, last access: 22 May 2020b. a
Oerlemans, J.: Glaciers and climate change, Balkema Publishers, Rotterdam, 2001. a, b
Ohmura, A.: Physical Basis for the TemperatureBased MeltIndex Method, J. Appl. Meteorol., 40, 753–761, https://doi.org/10.1175/15200450(2001)040<0753:PBFTTB>2.0.CO;2, 2001. a
Pappenberger, F., Pagano, T. C., Brown, J. D., Alfieri, L., Lavers, D. A., Berthet, L., Bressand, F., Cloke, H. L., Cranston, M., Danhelka, J., Demargne, J., Demuth, N., de SaintAubin, C., Feikema, P. M., Fresch, M. A., Garçon, R., Gelfan, A., He, Y., Hu, Y. Z., Janet, B., Jurdy, N., Javelle, P., Kuchment, L., Laborda, Y., Langsholt, E., Le Lay, M., Li, Z. J., Mannessiez, F., Marchandise, A., Marty, R., Meißner, D., Manful, D., Organde, D., Pourret, V., Rademacher, S., Ramos, M. H., Reinbold, D., Tibaldi, S., Silvano, P., Salamon, P., Shin, D., Sorbet, C., Sprokkereef, E., Thiemig, V., Tuteja, N. K., van Andel, S. J., Verkade, J. S., Vehviläinen, B., Vogelbacher, A., Wetterhall, F., Zappa, M., Van der Zwan, R. E., and Thielendel Pozo, J.: Hydrological Ensemble Prediction Systems Around the Globe, Springer Berlin Heidelberg, Berlin, Heidelberg, 1–35, https://doi.org/10.1007/9783642404573_471, 2016. a
Pellicciotti, F., Brock, B., Strasser, U., Burlando, P., Funk, M., and Corripio, J.: An enhanced temperatureindex glacier melt model including the shortwave radiation balance: development and testing for Haut Glacier d'Arolla, Switzerland, J. Glaciol., 51, 573–587, https://doi.org/10.3189/172756505781829124, 2005. a, b, c, d
Ristic, B., Arulampalam, S., and Gordon, N.: Beyond the Kalman filter: Particle filters for tracking applications, vol. 685, Artech house, Boston, 2004. a
Ritz, C., Edwards, T. L., Durand, G., Payne, A. J., Peyaud, V., and Hindmarsh, R. C.: Potential sealevel rise from Antarctic icesheet instability constrained by observations, Nature, 528, 115–118, 2015. a
Rounce, D. R., Khurana, T., Short, M. B., Hock, R., Shean, D. E., and Brinkerhoff, D. J.: Quantifying parameter uncertainty in a largescale glacier evolution model using Bayesian inference: application to High Mountain Asia, J. Glaciol., 66, 1–13, https://doi.org/10.1017/jog.2019.91, 2020. a
Ruiz, J. J., Pulido, M., and Miyoshi, T.: Estimating Model Parameters with EnsembleBased Data Assimilation: A Review, Journal of the Meteorological Society of Japan. Ser. II, 91, 79–99, https://doi.org/10.2151/jmsj.2013201, 2013. a
Salzmann, N., Machguth, H., and Linsbauer, A.: The Swiss Alpine glaciers' response to the global “2 ^{∘}C air temperature target”, Environ. Res. Lett., 7, 044001, https://doi.org/10.1088/17489326/7/4/044001, 2012. a
Saucan, A.A., Chonavel, T., Sintes, C., and Le Caillec, J.M.: Interacting multiple model particle filters for side scan bathymetry, in: 2013 MTS/IEEE OCEANS – Bergen, 1–5, https://doi.org/10.1109/OCEANSBergen.2013.6608125, 2013. a
Science Magazine: Europe's record heat melted Swiss glaciers, available at: https://bit.ly/2VpvAL3 (last access: 28 October 2021), 2019. a
Seroussi, H., Nowicki, S., Payne, A. J., Goelzer, H., Lipscomb, W. H., AbeOuchi, A., Agosta, C., Albrecht, T., AsayDavis, X., Barthel, A., Calov, R., Cullather, R., Dumas, C., GaltonFenzi, B. K., Gladstone, R., Golledge, N. R., Gregory, J. M., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huybrechts, P., Jourdain, N. C., Kleiner, T., Larour, E., Leguy, G. R., Lowry, D. P., Little, C. M., Morlighem, M., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Reese, R., Schlegel, N.J., Shepherd, A., Simon, E., Smith, R. S., Straneo, F., Sun, S., Trusel, L. D., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., Zhao, C., Zhang, T., and Zwinger, T.: ISMIP6 Antarctica: a multimodel ensemble of the Antarctic ice sheet evolution over the 21st century, The Cryosphere, 14, 3033–3070, https://doi.org/10.5194/tc1430332020, 2020. a
Sevruk, B.: Systematischer Niederschlagsmessfehler in der Schweiz, Der Niederschlag in der Schweiz, Kommissionsverlag Geographischer Verlag Kümmerly + Frey, Bern, 1985. a
Shannon, S., Smith, R., Wiltshire, A., Payne, T., Huss, M., Betts, R., Caesar, J., Koutroulis, A., Jones, D., and Harrison, S.: Global glacier volume projections under highend climate change scenarios, The Cryosphere, 13, 325–350, https://doi.org/10.5194/tc133252019, 2019. a
SLF: WSL Institute for Snow and Avalanche Research (SLF) Operational SnowHydrological Service, available at: https://www.slf.ch/en/snow/snowasawaterresource/snowhydrologicalforecasting.html, last access: 10 June 2020. a
Sommer, C., Malz, P., Seehaus, T. C., Lippl, S., Zemp, M., and Braun, M. H.: Rapid glacier retreat and downwasting throughout the European Alps in the early 21 st century, Nat. Commun., 11, 1–10, 2020. a
Stöckli, R.: The HelioMont Surface Solar Radiation Processing, Tech. Rep. 93, MeteoSwiss, Zürich, 2013. a
Swiss Academy of Sciences: Press Release on Glacier Melt 2019, available at: https://bit.ly/2UK6YfD (last access: 18 October 2021), 2019. a
swisstopo: Swisstopo Swissalti3D, available at: https://shop.swisstopo.admin.ch/de/products/height_models/alti3D, last access: 8 June 2020. a, b
van Leeuwen, P. J., Künsch, H. R., Nerger, L., Potthast, R., and Reich, S.: Particle filters for highdimensional geoscience applications: A review, Q. J. Roy. Meteor. Soc., 145, 2335–2365, https://doi.org/10.1002/qj.3551, 2019. a, b, c
Wang, R., Work, D. B., and Sowers, R.: Multiple model particle filter for traffic estimation and incident detection, IEEE T. Intell. Transp., 17, 3461–3470, 2016. a
Werder, M. A., Huss, M., Paul, F., Dehecq, A., and Farinotti, D.: A Bayesian ice thickness estimation model for largescale applications, J. Glaciol., 66, 137–152, https://doi.org/10.1017/jog.2019.93, 2020. a
WSL: Swiss Federal Institute for Forest, Snow and Landscape Research (WSL) platform for drought monitoring drought.ch, available at: http://www.drought.ch/Messungen/index_DE#, last access: 8 June 2020. a
Wu, W., Emerton, R., Duan, Q., Wood, A. W., Wetterhall, F., and Robertson, D. E.: Ensemble flood forecasting: Current status and future opportunities, WIREs Water, 7, e1432, https://doi.org/10.1002/wat2.1432, 2020. a
Zappa, M., Rotach, M. W., Arpagaus, M., Dorninger, M., Hegg, C., Montani, A., Ranzi, R., Ament, F., Germann, U., Grossi, G., Jaun, S., Rossa, A., Vogt, S., Walser, A., Wehrhan, J., and Wunram, C.: MAP DPHASE: realtime demonstration of hydrological ensemble prediction systems, Atmos. Sci. Lett., 9, 80–87, https://doi.org/10.1002/asl.183, 2008. a
Zappa, M., van Andel, S. J., and Cloke, H. L.: Introduction to Ensemble Forecast Applications and Showcases, pp. 1–5, Springer Berlin Heidelberg, Berlin, Heidelberg, https://doi.org/10.1007/9783642404573_451, 2018. a
Zekollari, H., Huss, M., and Farinotti, D.: Modelling the future evolution of glaciers in the European Alps under the EUROCORDEX RCM ensemble, The Cryosphere, 13, 1125–1146, https://doi.org/10.5194/tc1311252019, 2019. a
Zellner, A.: On assessing prior distributions and Bayesian regression analysis with gprior distributions, Bayesian inference and decision techniques, in: Bayesian Inference and Decision Techniques: Essays in Honor of Bruno de Finetti, edited by: Goel, P. and Zellner, A., Elsevier Science Publishers, Inc., New York, 233–243, 1986. a
 Abstract
 Introduction
 Study sites, data, and field instrumentation
 Methods
 Results and discussion
 Conclusions
 Appendix A: Handling of multiple cameras
 Appendix B: Resampling procedure
 Appendix C: Temporal evolution of the mass balance state with the example of Findelgletscher
 Code and data availability
 Video supplement
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Review statement
 References
 Abstract
 Introduction
 Study sites, data, and field instrumentation
 Methods
 Results and discussion
 Conclusions
 Appendix A: Handling of multiple cameras
 Appendix B: Resampling procedure
 Appendix C: Temporal evolution of the mass balance state with the example of Findelgletscher
 Code and data availability
 Video supplement
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Review statement
 References