2.2.1. Multi-SAR-System

All acquired SAR data were processed with a pre-processing framework called Multi-SAR-System. The Multi-SAR-System enables processing of SAR data being acquired with different sensors within the same system and output of the processed data in a uniform format. The pre-processing system is implemented at the Earth Observation Center of the German Aerospace Center, Oberpfaffenhofen, Germany. It contains geocoding, radiometric calibration and image enhancement. During geocoding, radar image distortions induced by the side-looking geometry of SAR-systems are corrected using a digital elevation model (DEM). The geometric quality of the geocoding depends on the height accuracy and the resolution of the elevation model. We used an airborne laser scanning DEM with a spatial resolution of 10 m generated between 2006–2010 (data source: Province of Tyrol—data.tirol.gv.at). The next step, after geometric adjustments, was a coarse radiometric calibration. Local topographic variations and sensor position affect not only location, but also the brightness of the radar return [27]. In mountainous terrain, this effect can be reduced with a gamma correction. This approach compares the ratio of individual area parcels of a synthetic image derived from DEM and orbit information with the illuminated SAR images. Due to the integration of backscatter over the area parcels, one S1 scene require about one full day of processing time. For regions with moderate slope angles such as for the glacier areas investigated here, the approximation of the gamma correction based on the local incidence angle using *γ*0 = *β*<sup>0</sup>*tan*Θ, *γ*0 the backscatter coefficient, *β*0 the backscatter in beta naught and Θ the local incidence angle is sufficient. The processing time using the approximated gamma correction reduces to about 20 min per scene. Our interpolation uses the cubic convolution method on a 17 × 17 pixel raster. The last step of the Multi-SAR-System is the image enhancement to reduce the influence of additive and multiplicative noise contribution. Compared to conventional speckle filtering algorithms, the multi-scale, multi-looking approach applied here adapts the local number of looks to the image content. For heterogeneous surfaces (i.e., mountainous areas), a minimal look number and, consequently, the maximum geometric resolution is necessary to adequately describe surface features. The choice of an appropriate look number is made by the help of a novel perturbation-based noise model that combines both additive and multiplicative noise contributions and automatically adapts to sensor and imaging mode characteristics via the delivered metadata. The result is a very smooth, but detail-preserving multi-looked SAR image representing the local optimal trade-off between geometric resolution and radiometric accuracy. Finally, image values are converted to the standard unit decibel (dB) [28,29].

As a next step, two different masks are applied. The first mask eliminates the radar shadow and layover effects and the second mask removes non-glaciated areas (based on the glacier boundaries of 2009 derived from optical remote sensing data). As a result, we receive SAR data of just the respective glaciated areas described in Figure 1.

#### 2.2.2. Backscatter Variability of Wet Snow Scenes

To determine the variability in *γ*0 for homogenous wet snow surfaces, we calculated the coefficient of variation (CV, *CV* = *μ σ* , *μ* the arithmetic mean and *σ* the standard deviation) per June scene. A low variability in backscatter distribution per glacier areas indicate homogeneous conditions. We set the CV in backscatter to below 20% as criterion for data selection. Table 3 displays the determined CVs for June scenes in cross- and co-polarization.

Almost all cross-polarized SAR scenes in Table 3 show a CV for June acquisitions of 15% or less. The single exception is the S1 scene being recorded in ascending orbit in 2015. All co-polarized scenes, however, result in a variation significantly higher than for cross-polarized acquisitions and show a much larger range. Such larger ranges indicate less consistency. Consequently, we use only on cross-polarized SAR scenes.

A secondary criterion for data selection is time of acquisition. During the main melt season in August, glacier ice and firn surfaces are very much affected by surficial melt water streams during the day, while early in the morning at 5:30 UTC (local time 7:30) nocturnal refreezing has peaked. Unfortunately, all ascending scenes were acquired at 17:10 UTC (local time 19:10), when it is more likely that strong surface wetting and meltwater runoff can lead to misinterpretation for the extent of wet snow covered areas. Hence, in the following, we will present methodology and results only for descending cross-polarized orbits, which were acquired at 5:30 UTC in the morning (Table 2). We excluded the first S1 June scene in Figure 3b (recorded on 2 June 2015) despite a CV value of 14% (Table 3). The median value of this scene is distinctly higher than for the remaining scenes (Figure 3). We expect that higher elevation ranges for all glaciers were still covered by dry snow and, hence, increase the median backscatter to higher values.


**Table 3.** Determined coefficients of variation for all acquired Sentinel-1 scenes in June between 2015–2018. Respective orbits are described in Table 2.

**Figure 3.** Box plots of backscatter distributions of all recorded cross-polarized June acquisitions in descending orbits with calculated coefficients of variation below 20% in Table 3. The red horizontal lines within the boxes represent the median in *γ*0, the boxes frame the interquartile range and the whiskers display extreme values not considered as outliers. Outliers are presented through red crosses. The black horizontal line displays threshold *β*1 = −21.45 dB and the purple line displays threshold *β*2 = −22.41 dB—for details, see Section 2.2.4.

#### 2.2.3. Correction of Systematic Backscatter Offsets

However, cross-polarized wet snow scenes are very homogeneous (Figure 3), and topographic effects very often bias interpretation of satellite data. In addition, very complex topographies, such as mountain glacier regions with various slope aspects, are challenging for space borne acquisitions to interpret. In addition, we relied on a DEM from 2006 for data processing. In particular, at glacier tongues, topography has lowered significantly between time of acquisition and DEM generation. Since we had several acquisitions in June with equal orbits during four years of data acquisition, we could check whether variations to average *γ*0 per scene are fully random or very similar for consecutive years and various melt progresses. We calculated for each June SAR scene deviations from the respective median. Deviations in *γ*0 being constant in location and independent from melt progress can be attributed as systematic due to topography and incidence angle. Next, we calculated the variance of determined deviations for all June scenes for each pixel. Calculated variances for the nine wet snow scenes below a value of 0.5 indicate that offsets from median are systematic and are corrected for. Figure 4 displays the resulting constant deviations used as backscatter offset correction file. To prevent misinterpretation of snow-free glacier tongues, we set corrections for the GPF and HEF glacier tongues to zero. This file was subtracted from each single S1 scene as additional topography correction. Instead of manually discarding snow-free glacier termini, one can apply automated detection algorithms (e.g., [30]) to minimize influences of already snow free glacier parts. We consider platform effects such as variations in perpendicular baseline or image frame as negligible for the here performed analysis. However, in case such deviations would be significant, the effect per pixel would vary with different baselines and, hence, the pixel variances of different scenes increase.

**Figure 4.** Applied backscatter coefficient (*γ*0) correction in decibel (dB), which is subtracted from each scene. Coordinates are given in UTM with datum WGS 1984.

## 2.2.4. Threshold Determination

The above selected nine wet snow scenes (Figure 3) are used to calculate wet snow thresholds. We calculated for each scene the upper quartile (75% percentile, upper frame of the boxes (Figure 3) of the backscatter distribution. The arithmetic mean of those upper quartiles is used to discriminate wet snow and wet firn from dry snow and bare ice areas (threshold *β*1, black line in (Figure 3). In a second step, we used solely the parts of the wet snow scenes with *γ*0 < *β*1 and calculated the 95% percentiles per remaining part of the scenes. The arithmetic mean of those percentiles result in threshold *β*2 (purple line in (Figure 3). The resulting values are *β*1 = −21.45 dB (with standard deviation of 0.76 dB) and *β*2 = −22.41 dB (with standard deviation of 0.18 dB). The determined standard deviations are further used for sensitivity analyses.

For the applied two-step approach, we first analyze for glacier parts being classified as wet using *β*1. After correction for systematic backscatter offsets, we subtracted *β*1 from each S1 scene resulting in a classification as dry/bare ice (positive values) or wet (negative values). In case the classified wet area is less than 50% of the glacier, we further discriminated for firn and wet snow utilizing *β*2. Here, we used again the offset corrected S1 scenes, removed all dry parts classified with *β*1 and subtracted the wet areas by *β*2. The resulting pixels were classified as wet snow (negative values) and firn (positive values). The proportion of each glacier AOI being either wet or wet snow divided by the whole AOI area results in the WSCAF. Glacier parts, which were covered by dry snow, cannot be distinguished from bare ice areas. As a result, we interpret WSCAF per glacier AOI as accumulation area fraction.

#### *2.3. Processing of Optical and Near-Infrared Remote Sensing Data*

Similar to, e.g., Nagler et al. [9] and many other authors, we used L7, L8 and S2 acquisitions to assess the accuracy of the presented SAR workflow on capturing transient snowlines. From now on, we use the term "optical" remote sensing data to abbreviate the usage of Landsat or Sentinel-2 data.

While the SAR data presented here are solely sensitive to wet snow, data within and close to the optical frequency range display brightness differences at the surface. However, cloud coverage often prevents data usage from optical instruments. For instance, in summer 2017, not a single Landsat scene for the here described AOIs could be used for comparison with SAR data. Each Landsat scene for this region was strongly influenced by clouds (a maximum of 20% cloud coverage per glacier was set as cut-off criterion for data usage). For the other observed summer seasons, we were able to make use of eight optical scenes. We used the normalized difference snow index (NDSI) [31–33] for automatic classification of snow covered areas for Landsat data and the Level-2A algorithm for Sentinel-2 data [34]. For L8 scenes, we used a threshold of 0.6 and for L7 this threshold had to be increased to 0.88 to obtain reliable results for snow covered areas. Snow classification thresholds for S2 data were increased to 0.16 (2017) and 0.18 (2018) from the given 0.12 [34] to match best the visually inspected snow boundaries. Since coincidences of optical sensor acquisitions and SAR imagery are rare, the quantification of spatial deviations between both data interpretations is difficult. We defined a range of four days prior and four days after each optical acquisition as being reasonable for direct comparison. For the case of new snow events within this time range, data were not used for analysis.

#### *2.4. Annual Minimum Wet Snow Covered Area Extent*

For remote sensing data, fixed date acquisitions such as annual mass balance assessments per glaciological year (GY) are difficult to achieve due to overpass rates and priority rankings by the satellite operators. However, more importantly, conditions on the glacier can be very different from year to year at fixed dates. In this study, we propose searching for local minimum in a wet snow extent within a range of five weeks prior to 30 September each year and two weeks after that day. Within this time range, we search for the local minimum in WSCAF for SAR scenes. However, recent new snow falls shortly before a SAR acquisition influence the detected WSCAF drastically. A minimum by early October could arise due to the fact that almost the entire glacier was covered by dry snow and, hence, volume scattering drastically increased. To prevent misinterpretation of derived snowlines, we applied a correction algorithm before minimum search based on observed firn and snowlines from optical data. We used the fact that perennial firn is more resistant to melt than seasonal snow. For the last analyzed S2 scene (16 September 2018), we decreased the snow threshold to 0.0 (instead of 0.18) and determined the area per AOI for reflectivities larger than zero. These areas include, at least, a fraction of the perennial firn coverage per AOI. Since those firn areas were formed in previous years from snow reserves, it is impossible that the SAR detected firn and WSCAF (result of step 1) is below those optical results. We defined that any SAR scene with lower "wet" areas than 75% of this optical reference scene as being influenced by recent new snow. Hence, those scenes were excluded from annual minimum snowline searches. The resulting minimum snowline is used as annual AAR from SAR scenes.
