## 1 Introduction  
Glacier mass balance is a key measure to determine the contribution of glaciers to regional hydrology [1, 2] and global sea level [3, 4, 5]. As the number of glaciers with direct measurements of mass balance (using the glaciological method) is limited (e.g., [6]) and their representativeness for the larger mountain ranges often unknown [7], it is appealing to obtain glacier mass changes from remote sensing data. Over larger regions the differencing of digital elevation models (DEMs) from two points in time (geodetic method) has been applied widely at approximately decadal resolution to determine total volume and mass changes of glaciers [8; 9; 10; 11]. However, for several geoscientific and socio-economic applications, the availability of region-wide glacier mass changes with annual resolution would be beneficial. As direct measurements reveal a high correlation of glacier mass balance with (a) the snow cover on a glacier and (b) the elevation of the snow line altitude (SLA), mapping snow cover (SC) on glaciers from satellite images offers a proxy for glacier mass balance (e.g., [12; 13; 14]), whereby the remotely sensed snow cover ratio (SCR) is taken as a proxy for the accumulation area ratio (AAR) [12; 13; 15] and the elevation of the snow line at the end of the ablation period [12; 16] as a proxy for the equilibrium line altitude (ELA). To compensate for local variability in elevation, the ELA is defined mathematically as the elevation where the vertical mass balance profile obtained by direct measurements crosses zero [17]. Mapping SC extent has also been applied to reconstruct missing mass balance measurements [18] on the basis of manually selected satellite scenes coinciding with the end of the ablation season, revealing a high correlation between satellite-derived snowlines and ELAs measured in the field. Snow and (clean) glacier ice together can be easily mapped from multispectral satellite images by application of various band ratios like near infrared (NIR) / shortwave infrared (SWIR) [19], red / SWIR [20], or the normalized difference snow index (NDSI), which is given as (green - SWIR) / (green + SWIR) [21], and a manually selected threshold. In contrast, the classification of snow on glaciers is difficult with these band ratios, as the shape of the spectral curves of ice and snow is very similar and ratios thus result in about the same values for both facies. A robust and widely used approach for mapping SC on glaciers is based on the use of a threshold on either single-band reflectance values, preferably in the NIR to avoid saturated values over snow with 8-bit Landsat data [12; 19; 22; 23; 24; 25; 26; 27; 28; 29], or albedo products that are readily available [25]. Mapping snow extent on glaciers has thus come a long way, starting in the 1970s with early tests using contrast-enhanced Landsat multispectral scanner (MSS) satellite images [22; 23], to the recent use of near-daily MODIS data to derive time series of surface albedo and SC maps for even small glaciers (e.g., [24; 25; 30]). For repeat application over large regions, it would be beneficial to derive the snow-covered area on glaciers automatically, including determination of the snow line altitude (SLA) near the end of the ablation season (as a proxy for the ELA). This is made challenging because of the following: (i) Due to the strong impact of terrain orientation (defined by its normal vector) on surface reflectance in mountain topography, deriving any measure of SC from reflectance requires implementation of a topographic normalisation using a DEM; (ii) SC on glaciers is often patchy and snow lines are in general not parallel to elevation contour lines; and (iii) highly variable atmospheric conditions (e.g., due to cloud shadows) imply that the same threshold might not work equally well on all glaciers.