Evolution of Surface Characteristics of Three Debris-Covered Glaciers in the Patagonian Andes From 1958 to 2020

A number of glaciological observations on debris-covered glaciers around the globe have shown a delayed length and mass adjustment in relation to climate variability, a behavior normally attributed to the ice insulation effect of thick debris layers. Dynamic interactions between debris cover, geometry and surface topography of debris-covered glaciers can nevertheless govern glacier velocities and mass changes over time, with many glaciers exhibiting high thinning rates in spite of thick debris cover. Such interactions are progressively being incorporated into glacier evolution research. In this paper we reconstruct changes in debris-covered area, surface velocities and surface features of three glaciers in the Patagonian Andes over the 1958–2020 period, based on satellite and aerial imagery and Digital Elevation Models. Our results show that debris cover has increased from 40 ± 0.6 to 50 ± 6.7% of the total glacier area since 1958, whilst glacier slope has slightly decreased. The gently sloping tongues have allowed surface flow velocities to remain relatively low (<60 m a−1) for the last two decades, preventing evacuation of surface debris, and contributing to the formation and rise of the ice cliff zone upper boundary. In addition, mapping of end of summer snowline altitudes for the last two decades suggests an increase in the Equilibrium Line Altitudes, which promotes earlier melt out of englacial debris and further increases debris-covered ice area. The strongly negative mass budget of the three investigated glaciers throughout the study period, together with the increases in debris cover extent and ice cliff zones up-glacier, and the low velocities, shows a strong linkage between debris cover, mass balance evolution, surface velocities and topography. Interestingly, the presence of thicker debris layers on the lowermost portions of the glaciers has not lowered thinning rates in these ice areas, indicating that the mass budget is mainly driven by climate variability and calving processes, to which the influence of enhanced thinning at ice cliff location can be added.

A number of glaciological observations on debris-covered glaciers around the globe have shown a delayed length and mass adjustment in relation to climate variability, a behavior normally attributed to the ice insulation effect of thick debris layers. Dynamic interactions between debris cover, geometry and surface topography of debris-covered glaciers can nevertheless govern glacier velocities and mass changes over time, with many glaciers exhibiting high thinning rates in spite of thick debris cover. Such interactions are progressively being incorporated into glacier evolution research. In this paper we reconstruct changes in debris-covered area, surface velocities and surface features of three glaciers in the Patagonian Andes over the 1958-2020 period, based on satellite and aerial imagery and Digital Elevation Models. Our results show that debris cover has increased from 40 ± 0.6 to 50 ± 6.7% of the total glacier area since 1958, whilst glacier slope has slightly decreased. The gently sloping tongues have allowed surface flow velocities to remain relatively low (<60 m a −1 ) for the last two decades, preventing evacuation of surface debris, and contributing to the formation and rise of the ice cliff zone upper boundary. In addition, mapping of end of summer snowline altitudes for the last two decades suggests an increase in the Equilibrium Line Altitudes, which promotes earlier melt out of englacial debris and further increases debris-covered ice area. The strongly negative mass budget of the three investigated glaciers throughout the study period, together with the increases in debris cover extent and ice cliff zones up-glacier, and the low velocities, shows a strong linkage between debris cover, mass balance evolution, surface velocities and topography. Interestingly, the presence of thicker debris layers on the lowermost portions of the glaciers has not lowered thinning rates in these ice areas, indicating that the mass budget is mainly driven by climate variability and calving processes, to which the influence of enhanced thinning at ice cliff location can be added.

INTRODUCTION
Driven by climate change, glaciers in most high mountain ranges around the globe are receding and losing mass rapidly (Jomelli et al., 2011;Davies and Glasser, 2012;Zekollari et al., 2019;Davaze et al., 2020;Hugonnet et al., 2021). Whilst glaciers steadily thin and contribute to sea level rise through increased melting (Zemp et al., 2019), another aspect of glacier change and adjustment to drier conditions is the progressive accumulation of surface debris stemming from the surrounding valley walls and lateral moraines (Scherler et al., 2011;Kirkbride and Deline, 2013;Sasaki et al., 2016;Tielidze et al., 2020). Debris cover on glacier surfaces accounts for ∼7% of the mountain glacier area worldwide (Herreid and Pellicciotti, 2020), and is relevant from a surface energy balance point of view, as the distribution and thickness variability of debris layers can alter glacier melt rates (Nicholson et al., 2018;Fyffe et al., 2020). Length adjustments and overall response times of debris-covered glaciers to climate variability are often delayed compared to debris-free glaciers. This delayed signal can occur as ablation rates decrease when debris layers are thicker than the effective thickness (i.e. the thickness above which glacier ice is sheltered and insulated from solar radiation, commonly <0.1 m; Mattson, 2000;Mihalcea et al., 2008;Reznichenko et al., 2010). Normally, the thickness of debris layers increases toward the glacier terminus, resulting in lower or even reversed mass balance gradients (e.g. Benn et al., 2012;Rounce et al., 2018;Bisset et al., 2020).
The spatio-temporal relationship between debris cover, surface mass balance glacier dynamics is nevertheless not simple. The mass balance gradient of debris-covered glaciers may not vary in straightforward form with elevation, and debris-covered glaciers are said to react by thinning or thickening in response to climate forcing rather than terminus recession or advance due to their typically gentle slopes (Bolch et al., 2008;Pellicciotti et al., 2015). Reduced ablation on the debris-covered portions of glaciers can result in a decrease in overall slope, which subsequently favors a reduction in ice flow velocities (Rowan et al., 2015;Dehecq et al., 2019). A feedback mechanism is established whereby rising Equilibrium Line Altitudes (ELA) lead to an earlier melt out and emergence of englacial debris, extending the debris layer further up-glacier. In this way, the overall glacier slope is additionally reduced, leading to stagnation of the glacier tongues (Quincey et al., 2009;Benn et al., 2012;Haritashya et al., 2015). Surface and basal ice melt, coupled with surface lowering, culminate in the formation and expansion of supraglacial lakes, causing further enhancement of glacier ablation (Mertes et al., 2017).
In spite of the above, more than a few studies, particularly in the Himalayas, have shown that debris-covered glaciers can have surface lowering rates comparable to those of debris-free glaciers (Kääb et al., 2012;Ragettli et al., 2016;King et al., 2019). A number of causes for this uncertainty have been identified so far (see e.g. Mölg et al., 2019 and references therein), including the influence of ice cliffs. Ice cliffs are steeply inclined areas of glacier ice, usually covered only by a very thin layer of debris, that vary from sub-circular to largely elongated in shape, and are of varying size (Sakai et al., 2002;Steiner et al., 2019). The occurrence of exposed ice in otherwise predominantly debris-covered glacier surfaces in the form of cliffs and supraglacial ponds (often associated) can boost ablation in comparison to "regular" debris-covered surfaces or even debris-free areas (Sakai et al., 1998;Sakai et al., 2002;Juen et al., 2014). The influence of ice cliffs and supraglacial ponds on the evolution of debris-covered glaciers is presently well acknowledged and the mapping and monitoring of these cryokarst forms has increased in recent years Watson et al., 2017;Herreid and Pellicciotti, 2018;Mölg et al., 2019;Steiner et al., 2019;Ferguson and Vieli, 2020). Bolch et al. (2012) pointed out that understanding the interactions between debris cover, mass balance and glacier dynamics would be a major achievement for envisaging the future evolution of debris-covered glaciers in high mountain areas. Consequently, recent studies have attempted to integrate this subject by assimilating diverse glacier variables (mass balance, debris distribution, thickness and transport, flow rates) and their feedback mechanisms into numerical models (Banerjee and Shankar, 2013;Collier et al., 2015;Rowan et al., 2015;Anderson and Anderson, 2018). Although large-scale assessments on debris-covered ice worldwide distribution exist (Scherler et al., 2018), modeling efforts are currently limited by the relative scarcity of more specific debris-covered glacier data, such as debris thickness and the evolution of ice cliffs and supraglacial ponds. These parameters are available mostly for the Himalaya, Karakorum, (e.g. Han et al., 2010;Thompson et al., 2012;Brun et al., 2018;Nicholson et al., 2018;Rounce et al., 2018;Buri et al., 2021) and to a lesser degree the European Alps (Reid and Brock, 2014;Mölg et al., 2019), but are lacking for other mountain regions where debris-covered glaciers are common (Röhl, 2006). In the Andes of Argentina and Chile, specifics about the evolution of debris layers, and their influence in glacier surface morphology, overall dynamics and mass changes, is particularly limited despite the abundance of debris-covered ice in some sectors such as the Central Andes (see e.g. Bodin et al., 2010;Wilson et al., 2016a;Falaschi et al., 2018;Ferri et al., 2020). In comparison to the Central Andes, however, the Patagonian Andes show a much lower relative concentration of debris-covered glaciers (Falaschi et al., 2013;Masiokas et al., 2015), although debris-covered ice can be relevant in some specific sites. Specifically for the Río Mayer catchment, debriscovered ice accounts for ∼15% of the total glacier area (IANIGLA, 2018). In addition, recent studies have detected signs of increasing proportions of debris-covered ice in the large outlet glaciers from the Patagonian Icefields in recent decades (Glasser et al., 2016).
In this study, we aim to determine the main characteristics of the morphological features of the Río Oro, Río Lácteo and San Lorenzo Sur debris-covered glaciers in the eastern (Argentinian) side of the Monte San Lorenzo massif (Southern Patagonian Andes), namely the distribution and area changes of the debris layer, ice cliffs and supraglacial ponds from 1958 to 2020. In addition, we produce the most complete glacier velocity record from freely available (mostly Landsat) satellite imagery. Whilst the mass balance evolution of these glaciers has been strongly linked to climate variability, proglacial lakes influence glacier evolution to a significant extent as well (Falaschi et al., 2019). Additionally, these three glaciers have one of the longest (geodetic) mass balance records in the Southern Andes and have been studied previously in regards to their geometric evolution. Combining the newly produced datasets and mass balance and climate records previously available, we discuss the spatio-temporal interactions between surface features and glacier velocities, and the likely impact on mass balance evolution. Finally, we evaluate the possible linkage between the observed areal evolution of the debris-covered area in view of the latesummer snowline elevation and climate variability in the area.

STUDY AREA AND PREVIOUS STUDIES
Río Oro, Río Lácteo, and San Lorenzo Sur valley glaciers ( Figure 1; Supplementary Table S1) lie on the Argentinian side of Monte San Lorenzo (47°35′S, 72°18′W; 3,706 m asl), the second highest peak in Southern Patagonia, ∼70 km East of the southernmost tip of the Northern Patagonian Icefield. At these latitudes, the Andes act as an effective orographic barrier to the humid westerly flow (Rivera et al., 2007), producing steep meteorological gradients and a strong rain shadow effect to the east of the main water divide (Bravo et al., 2019;Sauter, 2020). In fact, whilst the Northern Patagonian Icefield has a clear maritime climatic regime and receives extreme precipitation (Lenaerts et al., 2014), the San Lorenzo area lies in a more transitional environment and receives much less precipitation (Damseaux et al., 2020).
The Monte San Lorenzo area is presently well-researched in terms of the number and area of glaciers, the hypsometrical distribution of debris-covered and debris-free ice, and the rates of glacier recession (1% a −1 in the last ∼3.5 decades; Falaschi et al., 2013). The three glaciers investigated in this study have thinned at considerable rates since the mid-20th century in spite of the debris cover, particularly in comparison to nearby debris-free glaciers (Falaschi et al., 2017). A strong link between climate variability (added to the influence of calving dynamics on proglacial lakes) and glacier mass loss has been made (Falaschi et al., 2019). Sagredo et al. (2017) estimated (>1,600 m asl) paleoand modern ELAs for the Tranquilo glacier in the northern flank of Monte San Lorenzo. There is currently, however, no detailed information about how climate change and altitudinal shifts of ELAs with time have promoted the development of debris layers and how the inception of feedback mechanisms between glacier dynamics (e.g. glacier flow velocity), morphological features (ice cliffs, supraglacial ponds) and glacier geometry have affected the evolution of these particular glaciers.
Glacier inventories in Southern Patagonia have shown that the surface of both outlet and mountain glaciers is mostly free of debris, and that the amount of debris-covered ice is usually marginal in comparison (e.g. Glasser and Jansson, 2008;Paul and Mölg, 2014;Masiokas et al., 2015;Glasser et al., 2016). In some specific areas of extremely abrupt relief, high valley walls and lateral moraines can be the source of larger quantities of supraglacial debris, and extensive debris layers and debriscovered glaciers can exist, such as in the Monte San Lorenzo massif. Río Lácteo and San Lorenzo Sur are valley glaciers lying at the very bottom of Monte San Lorenzo's >2,000 m nearly vertical East face, have lake-terminating fronts and noticeable debris layers in an otherwise largely debris-free glacier inventory in the area (Falaschi et al., 2013). On the other hand, Río Oro has a minor proportion of debris-covered area compared to Río Lácteo and San Lorenzo Sur, and its only partially calving-type.
Probably because of their relevance as contributors to sea-level rise (Dussaillant et al., 2019), much scientific research upon the impact of climate change on glaciers has been focused in the nearby Patagonian Icefields. Numerical Modeling and mapping of ELAs and Accumulation Area Ratios (AAR) from remotesensed data has been carried out to test the mass balance sensitivity of glaciers (Barcaza et al., 2009;De Angelis, 2014). In the periphery of the Patagonian Icefields, however, there are still a vast number of smaller mountain glaciers that can provide complementary information for refining our understanding of overall glacier dynamics and glacier-climate interactions.

DATA AND METHODS
For this study we used 1) a set of orthophotos derived from archival aerial images, 2) medium resolution satellite imagery (Landsat 5, 7, 8, SPOT 5, Sentinel 2) available from various archives and 3) high spatial resolution satellite imagery (Pléiades) available from Google Earth. Unfortunately, during the commercial phase of Landsat, Thematic Mapper (TM) data for Southern Patagonia between 1984 and 1997 was not acquired beyond very few local requests (Goward et al., 2006). In general, these scenes are not suitable for glacier mapping (due to cloud cover or seasonal snow), and thus this period had to be omitted from the analysis (except for the 1985-1986 interval). The main characteristics of the implemented imagery (including spatial resolution, usage of each scene/photo is described in Supplementary  Table S2, whereas Supplementary Table S3 shows the time span covered on each of the specific analyses caried out in this study. Glacier-specific geodetic mass balance values for the 1958-2018 period, as well as length and area changes for the investigated glaciers are available from Falaschi et al. (2019). These data are essential for shedding further light into the interaction of debris covers and the dynamic evolution of glaciers, and therefore incorporated to the discussion sections.

Mapping of the Debris-Covered Ice Extent
Glacier outlines for the years 1958, 1981, 2000, 2012, and 2018 were retrieved from Falaschi et al. (2019). For all the remaining survey years, we modified the glacier outlines by visual interpretation of Landsat scenes. In the 1958 and 1981 orthophotos produced by Falaschi et al. (2019), we mapped the extent of debris layers on each individual glacier by onscreen manual digitization. We considered only areas with a fairly continuous debris cover extent, discarding this way small debris patches and very thin debris covered areas. Excluding now the years 1958 and 1981, debris-covered ice area for each survey year was obtained by calculating the debris-free area on the corresponding Landsat scene (from which the glacier outline had been previously generated), using a thresholded ratio image of the NIR/SWIR band (Paul et al., 2015), and substracting the debris free mask from the total glacier area.
So far, a number of (semi) automatic methods for identifying and mapping of debris-covered ice have been proposed (e.g. Shukla et al., 2010;Racoviteanu and Williams, 2012;Rastner et al., 2014). These approaches, which combine multi-spectral and terrain analyses, are a good approximation and reduce processing time, but ultimately require additional manual corrections to produce the final glacier outlines. Paul et al. (2004) stated that after manual editing, the accuracy of the debris-covered glacier outlines can be equal to pixel resolution. We have thus calculated the uncertainty in debris-covered area by implementing ±1 pixel buffers around the debris-covered ice area. Conservatively, the uncertainty in the debris-covered area change between 1958 and 2020 (σ dc.tot ) was calculated by quadratically summing the uncertainties of the initial (σ dc.i ) and final (σ dc.f ) debris-covered areas:

Glacier Surface Velocity
Because glaciers are constantly adjusting and redistributing ice mass in an effort to reach an equilibrium state, ice surface velocity is frequently used as an indicator of glacier mass balance conditions and dynamic state (Wilson et al., 2016b). This relationship assumes that during positive mass balance conditions ice surface velocities increase due to increased ice deformation and vice versa. Here we used Landsat 5-7, ASTER, SPOT 5, and Sentinel-2 images to automatically obtain velocity field maps of the three investigated glaciers for the 1985-1986-2020periods. The 1958-1981and 1981-1985 observation periods were omitted from this analysis due to the lack of corresponding surface features in the imagery used (preventing the use of manual feature tracking techniques e.g. Mölg et al., 2019).
Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 671854 The optical satellite imagery was processed in the CIAS feature-tracking algorithm (Kääb and Vollmer, 2000). This algorithm calculates sub-pixel displacement and direction magnitudes by first coregistering a master (t1 or initial epoch) and slave (t2 or final epoch) images and then comparing the coordinates of the highest correlation pixel block in a search window of the slave image with those of the master pixel block. Given the relatively low speed flow expected for the Río Lácteo and San Lorenzo Sur glaciers, a time separation of ∼1 year between the initial and final scenes was required to detect significant glacier flow that would surpass the total uncertainty in surface velocity (see below). For survey intervals where the initial and final images differed in spatial resolution (e.g. Landsat TM and Landsat ETM+), the highest resolution were resampled to the coarser one before coregistration (Wilson et al., 2016b). Once the CIAS output displacement vectors were generated, a correlation threshold between the initial and final pixel block of r > 0.4 was used to remove abnormal velocities. The resultant displacement points were upscaled to annual velocities and converted to raster grids. The final output cell size of the velocity maps was 30 × 30 and 45 × 45 m for Landsat ETM+ and Landsat TM pairs, respectively. Due to the lack of identifiable features on the debris-free ice, only some portions of the debriscovered tongue of Río Oro can be derived and velocity variations along longitudinal profiles cannot be entirely identified. This prevented a full comparison between the velocity map of Río Oro and the neighboring Río Lácteo and San Lorenzo Sur.
The total uncertainty associated with each glacier velocity map σv tot was calculated for a given time interval as the quadratic sum of the systematic (σv sys ) and random (σv rand ) uncertainties, following the laws of error propagation (Berthier et al., 2003): The systematic uncertainty σv sys depends on the performance quality of the image cross-correlation algorithm and typically, a ±2 pixel accuracy is expected for Landsat scenes: where R is the pixel size and Δt the time separation between the initial and final epochs in days. In turn, the random uncertainty σv rand is measured over stable terrain (off-glacier), in a buffer area around each individual glacier outline minus shadow areas and overlay steep terrain (>25°).

Identification of Ice-Cliffs and Supraglacial Ponds
Mapping of ice exposures (cliffs and supraglacial drainage channels, Supplementary Figure S1) has been carried out in different studies by means of object-oriented and principal components analyses in satellite and aerial imagery (Racoviteanu and Williams, 2012;Mölg et al., 2019), slope thresholding of high-resolution elevation data (Herreid and Pellicciotti, 2018), linear spectral unmixing approaches (Kneib et al., 2021), and of course, by on-screen digitizing of high resolution satellite imagery (Watson et al., 2017;Steiner et al., 2019). For the sake of simplicity, here we mapped ice cliffs and ponds by on screen digitization. We also mapped exposed ice in supraglacial drainage channels, but excluded heavily crevassed areas. The terminal portion of the San Lorenzo Sur glacier was not fully covered in the 2011 imagery available on Google Earth. Since the ice cliffs and supraglacial ponds located on the glaciers investigated here are relatively small in size, only the highest resolution (<1 m) imagery available from our dataset was used (see Supplementary Table S2). The Minimum Mapping Unit (MMU) is the smallest sized object that determines whether an object can be identified and mapped in remotely sensed imagery, and needs to be four contiguous pixels in size. In this way, we mapped ice cliffs ≥4 m 2 , which is the MMU of our lowest resolution data (aerial photos, 1 m resolution). On the other hand, intra-annual variations in both number and area of supraglacial ponds between seasons has shown to be significant . We have minimized this potential bias by mapping ponds on images acquired from very late spring (December 10, 2009) to mid-summer (February 22, 2018).
The overall accuracy of manual mapping of ice cliffs is dependent on several sources of uncertainty, including the spatial resolution and quality of the utilized imagery (optical contrast, shadows), and the user's inherent mapping subjectivity when identifying an ice cliff and its characteristics. Some studies aiming to map ice cliff geometry exhaustively (including information regarding slope, aspect, slope area) have carried out in-depth delineation uncertainty assessments , by checking ice cliff outlines and omissions made by multiple users and slope analyses derived from high resolution DEMs. Whilst such a thorough quality evaluation was not considered here, three independent operators inspected visually all individual ice-cliffs, whilst we also carried out a multiple user digitization of 50 randomly selected cliffs on each scene (Watson et al., 2017). Results of the round-robin digitizations are provided in Supplementary Table S4 and Supplementary Figure S2. We observed two main sources of uncertainties in ice-cliff mapping. One source is related to inherent imagery properties, not necessarily related to spatial resolution, but rather image illumination conditions and contrast (aerial photos were the least adequate in this regard). The second source relates to ice-cliff specifics such as complexity or branchiness, preservation/degradation state and debris layer thickness on ice exposures.
In addition to the ice cliff inventory, we also investigated ice cliff evolution in terms of changes in mean size and total area, and the altitudinal distribution of ice cliffs through time. This was done by calculating the relative proportion of ice-cliff area per 50 m elevation band across each glaciers' debris-covered surface.

Mapping of the Transient Snowline Altitude
We mapped the elevation of the transient snow line (SLA) and Snow Cover Ratio (SCR, defined as the snow-covered area with respect to the total glacier area) for Río Oro, Río Lácteo, and San Lorenzo Sur glaciers. To infer possible ELA changes with time for Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 671854 the three surveyed glaciers, we relied on the assertion that the Snow Line Altitude at the end of the hydrological year provides an effective proxy of the equilibrium line in temperate glaciers (Rabatel et al., 2008 and references therein) for that specific year. Because differences between the theoretical annual ELA and the snowline may exist at the end of summer (due to mass additions or losses e.g. superimposed ice or atypical length of the ablation period), the term end of summer snowline has been coined in the Southern Hemisphere (Clare et al., 2002). The main limitations for deriving end of summer snowlines lies in that suitable satellite imagery for snow mapping on glaciers (no cloud cover or fresh snow on glacier surfaces), DEMs and glacier outlines when dealing with a long time series are not always available for each individual year. To minimize this impact, we implemented a prototype version of the Automatic Snow Mapping on Glaciers (ASMAG) tool (Rastner et al., 2019), running now on the Google Earth Engine (GEE) platform. The algorithm currently maps SCR from Landsat and Sentinel 2 scenes, after converting digital numbers to top of atmosphere reflectance. The histogram approach of Bippus (2011) is then used within a single glacier mask to distinguish snow from ice by implementing an automatic threshold based on the Otsu (1979) thresholding algorithm to the Near Infrared (NIR) band (Rastner et al., 2019). After comparing the automatic derived snow cover area with the manual digitized ones, it became obvious that the ASMAG approach leads to systematically lower values compared to the manual derived ones (mapped on RGB true color compositions). We attribute this to the still missing topographic correction procedure not yet implemented in the GEE procedure so that we used this prototype in this study only for a straightforward selection of the satellite scene with the highest SLA for each single survey year. After knowing all the automatic derived acquisition dates with the highest SLAs (termed "SLAs" hereafter) for each period (ranging between February 1 and April 10), we manually digitized the snow cover area to retrieve SLAs based on three DEMs (SRTM-X, ALOS PALSAR and Pléiades) for three different time periods (2000-2009, 2010-2017, and 2018-2020).

End of Summer Snowline Altitude and Climate
Here we investigate the influence of temperature and precipitation on the end of summer snowline elevations (and consequently ELAs) on Río Oro, Río Lácteo, and San Lorenzo Sur glaciers. For the 2000-2020 period, we performed correlation analyses between mapped snowlines and monthly, annual and seasonal (summer, autumn, winter and spring) temperatures of the Balmaceda meteorological station (45.912°S, 71.694°W, 517 m a.s.l.), retrieved from the GISS-NASA portal. A similar analysis was carried out, relating monthly precipitations for the December and January months prior to the February-March snowline mapping. Results are presented as the individual temperature and precipitation influence on both glacierspecific and glacier-averaged SLAs, as well as the joint temperature and precipitation influence, again at both glacierspecific and glacier-averaged scale. For all statistical analyses, we used a confidence interval of 95%.
Linear fits of the percentage of debris-covered area as a function of time revealed a trend of progressively higher percentages of debris-covered area on Río Lácteo (0.2% a −1 ) and San Lorenzo Sur (0.24% a −1 ) glaciers up to the year 2020, whereas there is no clear trend for Río Oro (Figure 3). In principle, expansion of proglacial lakes affecting the lowest part of the glaciers should lead to a smaller proportion of debris-covered glacier area, yet our results show that yearly debris input has compensated the debris-cover area losses through calving. Indeed, if the yearly debris-covered area is divided by the initial 1958 total glacier area, the debriscovered area increase rates rise to 0.38% a −1 (San Lorenzo Sur), 0.48% a −1 (Río Lácteo), and 0.33% a −1 for Río Oro. The regression analysis of debris-covered area with time for these later cases revealed statistically significant (p < 10 -3 for all glaciers), very high r values of 0.86 (Río Oro), 0.91 (Río Lácteo) and 0.96 (San Lorenzo Sur) at the 95% confidence level.
Concomitantly to the increase in debris-cover area, debris layers have progressively occupied higher elevations on glacier surfaces as well. Between 1958 and 2020, sporadic debris has appeared on accumulation areas, shifting the maximum elevation of debris-covered ice from 1,570 to 3,240 m on Río Oro, 1,515-2,840 m on Río Lácteo, and 1,700-2,840 m on San Lorenzo Sur. Considering the glacier thinning that has taken place during this time interval (Falaschi et al., 2019), mean elevation of the debris covered ice has increased ∼240 m (1,130-1,370 m) on Río Oro, ∼130 m (1,320-1,450 m) on Río Lácteo, and ∼160 m (1,280-1,440 m) on San Lorenzo Sur.
Regarding the thickness of the debris layer on the glaciers investigated, we did not undertake systematic in situ surveys, nor did we estimate it from remote sensing data. Field observations on ice exposures (crevasses, ice cliffs) of Río Lácteo (Supplementary Figure S1), however, revealed a rough estimate of debris thickness between 1 and 2 m in proximity of the glacier terminus. In the lateral margins, boulders up to ∼3 m in diameter are frequent, these possibly resulting from the degradation of the lateral moraines. Thickness of the debris layer progressively decreases in the upward direction and in the uppermost part of the debris-covered tongue, close to the debris-free ice, the debris layer is pebbly in grain size, a few centimeters thick and interspersed with ∼10 cm blocks.

Glacier Velocities
Glacier velocity profiles and maps and were produced for Río Lácteo and San Lorenzo Sur glaciers for the 1985-1986 and 2000-2020 periods (Figures 1, 4, 5). In terms of spatial distribution, the highest velocities for both glaciers were found at the steep icefalls located just below the vertical rock wall of Monte San Lorenzo's East Face (∼90 ± 33 and ∼75 ± 55 ma −1 ) for Río Lácteo and San Lorenzo Sur, respectively). In the case of San Lorenzo Sur, velocities steadily decrease from these icefalls toward the terminus, where the extremely low velocities (<18 ± 29 m a −1 ) indicate that the glacier terminus is practically stagnant. Aside from the icefall area, the velocity maps for Río Lácteo show higher velocities around the central flow line and a flow speed increase toward the calving front which is visible between 2012-2014 and 2018-2020 (maximum velocity of ∼51 ± 34 m a −1 in 2020) (Figures 4, 5).
Overall, we found that the spatio-temporal patterns and magnitudes of ice flows at both glaciers were fairly constant during the survey period . For both glaciers, average velocities only fluctuated within a narrow range of ±36.5 m a −1 along the whole length of the glaciers, and no significant trends of progressive velocity decrease were detected for the longer 1985-2020 time span ( Figure 5). Short-term variations were detected during sub-periods (e.g. speed-ups between the years 2001 and 2007, deceleration between 2015 and 2020), however, these changes were small in magnitude and within the displacement error rate calculated for nearby stable terrain.

Inventory and Area Changes of Ice Cliffs and Supraglacial Ponds
We found that the number of ice cliffs at the individual glacier scale increased considerably during the observation period, with (1,395) ( Table 1). This increase, however, amounts only to a relatively small change in the area occupied by ice cliffs (0.65 ± 0.07 and 0.78 ± 0.02 km 2 ). In absolute terms, the maximum in total ice cliff area (1.02 km 2 ) occurred in the year 2009.
For all surveyed years, San Lorenzo Sur was found to have the highest concentration of ice cliffs amongst the investigated glaciers, accounting for 75-85 and 85-90% of the total cliff number and area, respectively ( Table 1). For this glacier, the total ice-cliff area ranged between ∼0.54 ± 0.06 and ∼0.89 ± 0.23 km 2 , whereas the average area of individual ice-cliffs varied between 670 and 1,900 m 2 . Over time, glacier retreat and proglacial lake expansion removed a number of ice cliffs in the glacier lowest sections, though we still observed a number and area increase of ice cliffs. This possibly means that the relative number increase would be much more if the glaciers had land-terminating tongues. Interestingly, we also found that the maximum and mean area of ice cliffs decreased with time except for the year 1981, when the largest identified ice cliff was 83,918 m 2 on San Lorenzo Sur and caused an anomalously high mean size.  647,000 ± 67,000 627,000 ± 31,000 1,020,000 ± 200,000 556,000 ± 22,000 a 776,000 ± 279,000 774,000 ± 24,000 Share of debris-covered area (%) 3.2 ± 0.3 3.1 ± 0.2 5.6 ± 1.1 3.1 ± 0.1 a 3.9 ± 1.4 3.8 ± 0.1 Mean size (m 2 ) 1,012 1,559 873 741 a 726 556 a Indicates the terminal portion of San Lorenzo Sur was not surveyed for the year 2011 due to lack of coverage (see Figure 6).
Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 671854 Another aspect of ice-cliff evolution is the increasing altitudinal range at which ice cliffs occur (Table 1). Overall, the upper limit of this range has shifted toward higher altitudes on all three glaciers. For Río Oro (1,140-1,190 m) and Río Lácteo (1,270-1,320 m) glaciers, the shift was 40-50 m. The altitudinal shift of ice cliffs is better exemplified at San Lorenzo Sur glacier, where the increase is ∼180 m (1,280-1,480 m). When calculating the percentage share of ice-cliffs within debris-covered areas by elevation (considering 50 m elevation bands), we found that ice-cliffs began to form above 1,250 m at least from the year 2009 onwards, and that ice-cliff area below 1,150 m (the lowermost portion of the glacier tongue) decreased after the same year (Figure 6), probably as a consequence of fast glacier area removal by calving processes.
In comparison to ice-cliffs, the total area and number of supraglacial ponds was found to be consistently smaller, representing merely 0.3% of the total debris-covered area at all times among the three sampled glaciers. The total number and area of supraglacial ponds was shown to be highly variable on all three glaciers across time, with San Lorenzo Sur again hosting the highest concentration (∼70-90% in number). The total pond area on this glacier varied between 0.04 and 0.06 km 2 , whereas the mean size ranged between 350 and 780 m 2 . Except for the two initial survey years (1958 and 1981), more than 50% of supraglacial ponds were associated with an ice-cliff ( Table 2). For the most recent 2018 pond inventory, 38% of the ponds were situated on debris-covered ice surfaces with no relation to a cliff, 36% occur at the bottom of longitudinal or lateral cliffs, 13% are contained within a circular cliff, and a further 13% are associated with complex cliffs.

Changes in Snow Cover Ratio and Snow Line Altitude
Snow Cover Ratios and Snow Line Altitudes were calculated for Río Oro, Río Lácteo and San Lorenzo Sur glaciers. In general, Río Oro has always the greatest SCR values (mean 0.44), whereas the (average) minimum SCR of Río Lácteo and San Lorenzo Sur are 0.23 and 0.16 respectively (Figure 7). Data reveal also that when the minimum SCR values for each survey year are considered, SCR appears to have diminished with time on all three glaciers. Linear fits for the SCR scatter data suggest moderate to high correlation (statistically significant with α 0.05) with time for Río Lácteo (r 0.6; p 0.007), Río Oro (r 0.7; p 0.001) and San Lorenzo Sur (0.76; p 2 × 10 −4 ). The maximum SLA varied strongly on a yearly basis on all glaciers, with altitudinal differences between maximum and minimum values of up to ∼800 m on Río Oro and ∼600 m on Río Lácteo and ∼500 m San Lorenzo Sur. For the full study period, the averaged maximum SLA was 1905 ± 100 m (Río Oro), 1720 ± 85 m (Río Lácteo) and 1815 ± 165 m (San Lorenzo Sur).

Relationships Between Late-Summer Snowline and Climate Variables
At the monthly scale, March temperature is significantly related to SLAs (Figures 8, 9; Table 3). This is not surprising, since most of5 the snow line determinations were based on satellite images collected during March. Although the correlation coefficients between March temperature and SLAs are significant for all  three glaciers, the relationship is stronger when considering the glacier-averaged SLA in all three areas (r 0.64, p < 0.01; see Figure 8A). Based on this result, we show that March temperature explains about 37% of the total variance of SLA in the 2000-2020 interval. The annual temperature, from December of the previous year to November preceding the SLA determination, is also significantly associated, although weaker (r 0.58), to the SLA, indicating that SLA variations are also influenced by snowpack accumulation during the previous year. December-January precipitation is the seasonal arrangement most significant related to SLA. Similar to temperature, the relationships improve when the comparison between December and January precipitation is based on the glacieraveraged SLAs (r 0.68; see Figure 8B). The percentage of total variance in SLAs changes explained by December-January precipitation over the interval 2000-2020 is somewhat higher (R 2 around 44%) than that accounted by temperature.
The joint use of March temperature and December-January precipitation from Balmaceda as predictors of the SLAs at Mount San Lorenzo in a regression analysis shows that these climatic variables can explain more than 50% of the total variance of the glacier-averaged SLAs during the period 2000-2020 (adjusted r 2 55%, F(2,16) 12.031 p < 0.00065; Figure 8 and Table 3). However, we noted that each climatic variable has a distinctive impact on each of the glaciers ( Table 4). For the Río Oro glacier, the total variance explained by both climatic variables is high (adjusted r 2 50%), but only December-January precipitation is a significant (p 0.009) predictor for SLA variations in this glacier. For the Río Lácteo glacier, the explained SLA variance remains significant (adjusted r 2 41%) but none of the predictors are statistically significant at 95% confidence level (p 0.054 and p 0.085 for December-January precipitation and March temperature, respectively). Finally, the SLA variance explained for the San Lorenzo Sur glacier is adjusted r 2 40%, with March temperature being statistically significant (p 0.033) as a regression predictor (Table 5).

Increase in Debris-Covered Area
Overall, the total extent of debris-covered ice area for the three glaciers investigated has increased from 40 ± 0.6% to 50 ± 6.7% from 1958 to 2020. Because glacier recession in this study area is significant due to calving in proglacial lakes (particularly in the cases of Río Lácteo and San Lorenzo Sur), each year a large amount of debris-covered ice is lost (Falaschi et al., 2017). The above estimate of debris-covered area increase represents thus a minimum value of the potential debris-covered ice area increase. To give a better idea of the above, it can be said that if the initial 1958 glacier area is kept for all calculations (i.e., omitting all area losses since 1958 due to calving), the increase would be from 40 to 56%. Throughout the study period, the total debris-covered ice area has increased a steady rate, and we were unable to detect any short-term migration of the upper boundary of the debris layer toward lower elevation. Most debris input probably originates from the headwalls of each glacier, especially the East Face of Monte San Lorenzo, where the accumulation areas of the investigated glaciers are located. Although lateral moraine degradation is relevant on the investigated glaciers (Falaschi et al., 2019), we infer from the field visit to Río Lácteo glacier that they contribute merely to the glacier margins in the form of large boulders. Incidentally, van Woerkom et al. (2019) also showed that debris supply from lateral moraines alone cannot explain thick and continuous debris layers on debris-covered glacier surfaces. On the sampled glaciers here, no large rock avalanches were detected, which would rapidly add to the debris-covered area (see e.g. Berthier and Brun, 2019).
In the nearby Northern Patagonian Icefield, Glasser et al. (2016) found an increase in the number of debris-covered glaciers and a relative increase of ∼4% in debris-covered ice area between 1987 and 2015. Given the much larger size of the outlet glaciers investigated by Glasser et al. (2016) compared to the Monte San Lorenzo glaciers, this seemingly minor relative increase of ∼4% is naturally much larger in terms of the total debris-covered area (139 vs. ∼3 km 2 in San Lorenzo). Our findings nevertheless confirm that the trend of increasing debris-covered area is taking place not only in the large ice masses of the Patagonian 3 | Regression summaries for mean snowline altitudes (mean-SLA) using March temperature, December-January precipitation, and both variables as predictors.

Interaction Between Flow Velocity, Mass Balance, Proglacial Lakes and Ice Cliff Evolution
The flow velocity time series produced here for Río Lácteo and San Lorenzo Sur glaciers is, despite its relatively shortness, the longest one available for Patagonia beyond the outlet glaciers of the Patagonian Icefields (Mouginot and Rignot, 2015). The unavailability of Landsat TM imagery acquired during the second half of the 1980's and most of the 1990's, added to the persistent cloud cover in this region, preclude a longer and high temporal resolution assessment. Whilst some studies have used large boulders visible on old aerial photos to derive isolated and sporadic velocities (e.g. Mölg et al., 2019), this was not possible in our case, as the morphology of the glacier surfaces being observed have transformed greatly between the acquisition of our 1958 and 1981 aerial photographs and the availability of the first satellite images.
A possible explanation for the differences in flow speed evolution between Río Lácteo and San Lorenzo Sur described in Glacier Velocities may be related to the dynamic effect originating at the glacier fronts in response to local lake water depths. Greater proglacial depths enhance calving at the snout of lake-terminating glacier tongues due to increased glacier-lake contact area subjected to subaqueous melt, accelerating flow rates and ice flux toward glacier terminus (Kirkbride and Warren, 1997;Benn et al., 2007). Acceleration of flow rates has been in fact acknowledged at the glacier terminus of several glaciers in the Himalaya Pronk et al., 2021). The deceleration of San Lorenzo Sur could hence be linked to shallower lake depths, whereas the constant velocities on Río Lácteo (and even slight acceleration of the glacier front in recent years) might be connected to the presence of a deeper proglacial lake and consequently to higher calving rates.
The dynamic thinning described above is highly relevant for understanding glacier dynamics, since it has shown to exacerbate glacier retreat and thinning rates of freshwater calving glaciers compared to land-terminating ones King et al., 2019). Ice velocities close to the front of Río Oro (small panels in Figure 4) are similar to the maximum values found at Río Lácteo and San Lorenzo Sur (∼90 m a −1 ), though we interpret them as a result of a higher local slope and not a dynamic acceleration due to the presence of a proglacial lake. Falaschi et al. (2019) found a strong relation between higher rates of proglacial lake expansion and increased mass losses. In particular, the 1958-2018 thinning rates of the lake-terminating Río Lacteo and San Lorenzo Sur (1.87 and 1.94 m a −1 , respectively) appear to be much stronger than the land-terminating Río Oro ( Despite the above limitations in establishing a clear relation between flow speed and thinning rates, it can be argued that the overall mass balance of Río Lácteo and San Lorenzo Sur glaciers is primarily driven by climate variability and calving processes. Indeed, thinning rates increase at lower elevations at all of the investigated glaciers (Falaschi et al., 2019) in spite of the increase in debris layer thickness observed during field trips. Naturally, however, some spatial heterogeneity in the thinning patterns exist due to local variations in debris thickness, which is in turn partly linked to the occurrence of ice cliffs. To quantify the possible influence of ice cliffs on mass balance, we generated buffer zones of 35 m around the mapped ice cliffs (see Mölg et al., 2019) for three time intervals of geodetic mass balance studied in Falaschi et al. (2019). Also included within these buffer areas were supraglacial ponds, as they are closely related to ice cliffs (Watson et al., 2017;Steiner et al., 2019) and actually appear to contribute to their formation. The mean elevation change (dh/ dt) in meters per year for the elevation grid cells that correspond to the buffered areas (including actual ice cliffs and ponds) effectively shows higher thinning rates compared to the debris-covered portions that are not affected by ice cliffs/ ponds on nearly all glaciers and evaluated time periods (Table 6). This signal of increased ablation at ice cliff location might be even attenuated here, since a given ice cliff might have disappeared or migrated downglacier within each surveyed time interval. In addition, DEM differencing using 30-40 m resolution DEMs might be too coarse for the relatively small ice cliffs. High contribution to glacier ablation from ice cliffs (though seasonally variable) has also been acknowledged on several other debriscovered glacier sites like the Himalayas (e.g. Brun et al., 2016;Buri et al., 2021).
Another aspect of glacier velocity in relation to the occurrence of glacier surface features (cliffs and ponds) was brought out by Mölg et al. (2019) for the Zmuttgletscher in the Swiss Alps. These authors showed that time intervals of relatively high flow speeds hampered the emergence of ice cliffs, and concomitantly shifted the upper boundary of debris-covered area downglacier. In principle, our results for the investigated glaciers in Monte San Lorenzo agree with this assertion, since we also found a continuous increase in the debris-covered ice area and a steady shift of the upper altitudinal limit of the ice cliff area toward higher elevations during the study period. Average slopes on Río Lácteo and San Lorenzo, as retrieved from a 1958-2018 multi-sourced DEM time series (see Falaschi et al., 2019), have decreased only very slightly on the debris-covered glacier portions between 1958 and 2018 (9°-8°a nd 11°-9°respectively). The relatively high slopes can explain the still continually flowing glacier tongues, though velocities might be still low enough to have enabled the growth of ice cliff and supraglacial ponds areas at increasingly higher elevation with time. Watson et al. (2017) and Steiner et al. (2019), however, found no clear link between ice cliff spatial density and glacier surface velocity on debris-covered glaciers in the Everest and Langtang regions in the Himalaya, nor a preferential formation sites of ice cliffs on stagnating tongues. Also, Río Oro (3.4% max), and especially San Lorenzo Sur glacier (7.6% max), have a much higher ice cliff share of the debris-covered ice area compared to e.g. Zmuttgletscher (1.5% max), which flows at similar speed compared to San Lorenzo Sur, or a series of five valley glaciers in the Langtang Himalaya (max 3.8%) whose lowermost tongues are almost stagnant . Also, the 16-18% annual survival rate of the ice cliffs in the investigated glaciers (measured for the 2009-2011 and 2017-2018 time intervals) happens to be somewhat higher than for Langtang glacier (9-13% Steiner et al., 2019) despite higher flow velocities.

Influence of Snow Line Altitude/Snow Cover Ratio Variability in Glacier Changes
Our satellite-derived SLAs were highly variable in elevation (various hundreds of meters). This is due to limitations posed by the occasionally unavailable actual late-summer imagery to derive the highest SLA for many survey years (Supplementary  Table S2). In fact, there was only a single suitable scene for SLA mapping in many years. Hence, images acquired earlier in the ablation season will provide a lower SLA compared to actual latesummer scenes, as long as no snowstorms hit the mountains later in the ablation season. Facing the same limitations, Barcaza et al. (2009) recorded differences of up to 300 m in satellite-derived ELAs for singe glaciers in the Northern Patagonian Icefield between 1975 and 2003. With the restrictions imposed by the overall limited satellite imagery availability in mind, we found that SLAs have increased during the 2000-2020 time interval (Figures 8, 9). Increasing SLAs are more evident before the year 2008 than for later years, coinciding with the contemporary temperature and precipitation changes during these earlier years. From 2001 to 2008, December-January monthly precipitation decreases from 120 mm to only 20 mm, whereas there is no clear trend afterward. Likewise, March mean temperature shifts from ∼7°C before 2002 to ∼12°C in 2004 and remained fairly stable since then.
Given that the March air temperature and December-January precipitation can explain more than 50% of the total variance of the glacier-averaged SLAs, a first, preliminary conclusion from Changes in Snow Cover Ratio and Snow Line Altitude is that our mapped glacier-averaged SLAs and their changes through time are fairly representative of the climate variability in the Monte San Lorenzo region. Because in the temperate ice masses of Patagonia superimposed ice is not to be expected (Cuffey and Paterson, 2010), end of summer snowlines are a good approximation of ELAs. Therefore, our findings suggest progressively lower Accumulation Area Ratios (AAR) on the investigated glaciers (Figure 7), which is compatible with glaciers in decidedly disequilibrium state (Bakke and Nesje, 2011). Previously, a paleoclimatic reconstruction by Sagredo et al. (2017) in the nearby Tranquilo glacier (in the northern slope of Monte San Lorenzo), estimated a 130-170 m shift in ELA since the Last Millennium Maximum. Further supportive data to this assertion comes from Glaciar de los Tres, ∼200 km south of San Lorenzo (but still the closest glacier where in-situ measurements are available). Here, glaciological surveys have also shown rapidly increasing ELAs since the late 1990s (WGMS, 2020).
Climate data presented in Falaschi et al. (2019) indicate that there have been clear trends of diminishing annual precipitation and increasing summer air temperatures in the Monte San Lorenzo region, at least since the mid-20th century, promoting rising ELAs in the district. The observed increase in debris-covered area in the sampled glaciers may hence stem not only from the low glacier flow speeds that fail to evacuate debris from the glaciers' surface, but also from the enhanced emergence of englacial debris at higher elevations up-glacier due to higher SLAs with time. In consequence, climate forcing (by way of changes in SLAs through time) needs to be incorporated as a causing factor in the feedback interactions governing the dynamics of the investigated glaciers.
In regards to changes in SLAs and glacier mass balance, we split the 2000-2020 observation period into two sub-periods (2000-2012 and 2012-2020) to match the geodetic mass balance intervals sampled in Falaschi et al. (2019), and found a shift of time-averaged SLAs toward higher elevations on all three glaciers (∼300 m on Río Oro, ∼120 m on Río Lácteo and ∼100 m on San Lorenzo Sur). In principle, rising SLAs (ELAs) should promote higher thinning and ablation rates, as a now larger portion of the glacier is under melting conditions and the glacier surface is exposed to higher temperatures. Given that the largest shift in SLA was found ELA for Río Oro, it is surprising that Falaschi et al. (2019) found that thinning rates of both Río Lácteo and San Lorenzo Sur) have increased between 2000-2012 and 2012-2018 (Table 6), yet the thinning rates of Río Oro have decreased in the last decade approximately.

CONCLUSION
In this study we have reconstructed the evolution of surface features and morphology, debris cover, and a superficial flow velocity of Río Oro, Río Lácteo, and San Lorenzo Sur glaciers in Patagonia, which, taken as a whole, covers a time span of over 60 years. Throughout the study period, glaciers showed 1) an overall debris-covered area increase from 40 ± 0.6 to 50 ± 6.7%  2) slow surface flow velocities  3) an increase in ice cliff area and maximum elevation and 4) a shift toward higher end of summer snowlines .
In the last six decades approximately, summer air temperatures have increased significantly in the San Lorenzo region. Coupled with glacier dynamics linked to freshwater calving processes, this has led to strongly negative mass balance conditions throughout this period, in spite of extensive debris layers on the glaciers surface. In addition to the climate forcing and inherent calving dynamics that drive the mass changes of the sampled glaciers, here we add the abundant presence of ice cliffs as an additional factor leading to enhanced mass loss, especially in comparison to other debris free glaciers in this region. Higher summer temperatures have resulted in increased end of summer snowlines (and presumably ELAs), exposing this way bare ice and supraglacial debris at higher elevations. This has resulted in a ∼10% expansion of the debris cover extent among the three glaciers. The considerably large extent and further expansion of debris cover at higher elevations might suggest an increasing importance of ice cliffs in overall glacier mass budget.
Preliminary in situ observations on Río Lácteo glacier indicate increasing debris layer thickness toward glacier termini, which should hamper ablation and cause a gentler or reversed mass balance gradient. Under these conditions, we hypothesize that ablation maxima at the ice cliff locations, in addition to ice speed acceleration and ice thinning toward the lowermost parts of the glacier tongues have prevented this from happening.
The overall surface flow velocity record of Río Lácteo and San Lorenzo Sur for the last two decades did not show a clear shift toward slower velocities (indicating progressive stagnation) in spite of the (slightly) reduced glacier slope. The overall low velocities might stem from the reduced ice thickness with time rather than the gentle slope of the glaciers. Nonetheless, velocities have not been fast enough to evacuate supraglacial debris, facilitating the up-glacier expansion and thickening of debris layers, nor have they prevented the growth of newly formed ice cliffs and the occurrence of ice cliffs at increasingly higher elevation with time.
The results presented in this study, mostly based on remote sensing data and records, represents a first line approach for understanding the complex relations between climate variability, glacier dynamics, debris cover and ice cliffs of glaciers in southern Patagonia. In this regard, more specific and in situ surveys are necessary (such as glaciological mass balance of annual resolution, ablation and backwasting at ice cliff locations, debris cover thickness measurements) for a better understanding of these feedback mechanisms and their ingestion in larger scale assessments and numerical modeling.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article can be made available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
DF designed the study, analyzed data and wrote the manuscript. AR contributed to the study design and wrote the manuscript. ALR, SM, RV, PR, JZ, and AS analyzed data. All authors contributed to the writing of the final version of the manuscript.

ACKNOWLEDGMENTS
SRTM-X and ALOS PALSAR DEM tiles were obtained from the Deutsches Zentrum für Luft und Raumfahrt (DLR) EOWEB, and the Alaska Satellite Facility portals. The Pléiades data used in this study was provided by the Pléiades Glacier Observatory (PGO) initiative of the French Space Agency (CNES). Pléiades data© CNES 2018, Distribution Airbus D&S. We acknowledge Frank Paul and Nico Mölg from the Geography Department at the University of Zurich for the manual mapping of snowlines and Mariano Castro for his photos and observations on the debris layer of Río Lácteo glacier. Thanks also to Pierre Pitte (IANIGLA) for his insight into ELAs at Glaciar de los Tres. Special thanks to Tobias Bolch (University of St. Andrews, Scotland) for his constructive comments and Ryan Wilson (University of Huddersfield, United Kingdom) for proofreading a previous version of the manuscript. Finally, we thank handling editor Aparna Shukla and two reviewers for their valuable comments and suggestions, which led to a much improved presentation of this study's content.