## 3 Materials and Methods
### Accuracy and Uncertainties  
The uncertainty in elevation change was estimated on stable ground outside the glacier area. We assumed that these areas did not change in elevation during our study period. In order to remove outliers, a mask which contains only elevation points ranging from 2500 to 3900 m a.s.l. was generated. For all LiDAR post-processing we used the 2015 LiDAR data as reference. Both LiDAR data surveys of 2009 and 2015 were acquired with the same nominal vertical accuracy (\(\pm\)0.3 m) (Table 1). The comparison on stable areas on slopes <40\({}^{\circ}\) after the horizontal and vertical co-registration procedure resulted in a relative vertical error of \(\pm\)3.0 m between 1955 and 2000, \(\pm\)4.2 m between 2000 and 2009, and \(\pm\)0.3 m between 2009 and 2015.  
Since previous studies using SRTM and TanDEM-X data for the Southern Hemisphere acquired in the austral summer for the Southern Hemisphere reported minimal snow and ice penetration, and due to elevation bias due to SAR signal penetration in snow and ice [57], we compared SRTM-C and X band elevation over 20 glaciers located further south of the Echaurren Norte Glacier after respective co-registration of the DEMs [55]. We did not observe elevation difference that we could attribute to radar signal penetration (Figure S2). This indicated that melt conditions prevailed in the Central Andes during X- and C-band summer SRTM acquisitions on glacier surfaces up to 4000 m, preventing significant signal penetration (e.g., [57; 58; 59]) (Figure S2). Landsat imagery from the year 2000 shows a minimal snow cover in the areas surrounding the Echaurren Norte Glacier, confirming that all the winter snow had ablated (e.g., [21]). We assumed similar acquisition conditions for the TanDEM-X dataset with negligible radar signal penetration. This assumption was also supported by the increase of debris cover until 2015.  
Hence, the mass balance inaccuracy was calculated using the equation of error propagation (Equation (1)), including the following uncertainties (e.g., [53]):  
\[d\mathrm{M}=|\mathrm{M}|^{*}\sqrt{\left[\frac{\delta_{th}}{\frac{dh}{dt}\ total}\right]^{2}+\left[\frac{\delta_{A}}{A}\right]^{2}+\left[ \frac{\delta_{\rho}}{\rho}\right]^{2}} \tag{1}\]  
Here, (\(d\mathrm{M}\)) is the uncertainty in mass balance (M), (\(\delta_{th/dt}\) total) the mean error of elevation change over stable ground. This mean uncertainty corresponds to the residual after co-registration on stable areas. (A) denotes the glacier area and (\(\delta_{A}\)) the uncertainty for the glacier extent mapping (Table 1). We considered a density uncertainty (\(\delta_{\rho}\)) of \(\pm\)7%, which corresponds to \(\pm\)60 kg m\({}^{-3}\)[56].  
For TanDEM-X, error assessment was done using the error approach given by [PERSON] et al. [53]. This used the Median Absolute Deviation (MAD) calculated for each 5\({}^{\circ}\) slope bin:  
\[MAD_{Aw}=\frac{\sum MAD_{slpi}*A_{slpi}}{A} \tag{2}\]  
Equation (2) shows the performance for the area weighted (\(MAD_{Aw}\)), where (\(MAD_{slpi}\)) corresponds to each \(MAD\), (A) corresponds to the total analysed area and (\(A_{slpi}\)) is the area covered by each slope bin. The \(MAD_{Aw}\), is equivalent to \(\delta dh/dt\ total\)[53].