*4.2. Inter-Annual Dominant Cluster*

The inter-annual cluster is characterized by glacier surface velocity variations that are most likely influenced by factors controlling subglacial water pressure conditions over a longer temporal scale. We compared the Sea Surface Temperature (SST) curve for the Niño 3.4 region (National Oceanic and Atmospheric Administration, NOAA) to the residuals of our surface velocity time series. Similar to Maussion et al. [39] and Francou et al. [40], we shifted the SST time series by three months to account for the time lag in the glacier mass balance response, and only consider 0.4 K ≤ SST ≤ –0.4 K, with duration > 6 months [36]. Figure 6e shows the SST curve for the observation period. From January to December 2017 the SST intensity and duration do not qualify as either a La Niña or El Niño event. Between December 2017 and July 2018, the SST intensity and duration qualify this period as a La Niña event, which is followed by an El Niño event commencing in December 2019, extending through to September 2020. The duration of SST > 0.4 K from 20 January, to 15 June 2020, does not qualify as an El Niño event. We note that prior to our observation period, SST was 1.5 K for more than 6 months indicating a strong El Niño event in 2016 (not shown in Figure 6e).

We find that the observed Sea Surface Temperature Anomalies (SSTA) show a strong correlation (r > 0.5, *p* < 0.05) with the residuals obtained from surface velocity time series in the inter-annual cluster (Figure 6c,e). The general decrease in surface velocities at the beginning of the time series corresponds to insignificant SST variations verging into a La Niña event starting in December 2017. From March 2017, SST increases gradually. This is associated with a gradual increase in surface velocities and corresponds to the subsequent El Niño (December 2019–September 2020). Maussion et al. [39] used an SEB/SMB model calibrated against a 4-year temperature time series obtained on the Shallap glacier. Reanalysis of downscaled atmospheric variables were used to estimate the monthly SEB/SMB for the period 1980–2013. Correlation analysis revealed a strong anti-correlation (r < −0.5, *p* < 0.05) between SSTA and individual fluxes of the SEB/SMB. During an El Niño event, the correlation was attributed to increased air temperatures, leading to an elevated snowfall line, an increase in short-wave radiation, and reduced precipitation [39]. At higher altitudes, the ENSO correlation was found to be weaker and mostly related to changes in total precipitation. Therefore, during an El Niño event, meltwater production is strongly enhanced, where subglacial water pressure conditions influence glacier surface velocities, concomitant to the approximate duration of the El Niño event.

Our results also show that steep slopes with both East and West orientations, as well as relatively flat areas (regardless of altitude), are prone to inter-annual glacier surface velocity changes as a result of ENSO. This is supported by the findings of Kaser and Georges [2], where regional-scale differences in the ELA between E- and W-facing glaciers were observed in the CB. They found a zonal asymmetry caused by convective cloud development in the afternoon, which influences the radiation balance and ablation. In the morning, east exposed slopes received a high amount of shortwave solar radiation, whilst during the day, clouds develop, and west exposed glaciers receive a significantly lower amount of shortwave solar radiation. Further, Gurgiser et al. [38] found a topographic signature in surface energy and mass balance above 5000 m a.s.l at the Shallap glacier, where the combination of aspect and increased slope steepness results in a smaller incidence for solar radiation and a slower decrease in Albedo.

#### *4.3. Implications for Avalanching Glacier Instabilities*

Based on our continuous time series of surface velocities, and high-resolution satellite imagery we evaluated whether avalanching glacier instabilities were evident on the Palcaraju glacier for our period of observation. Various authors have suggested a high glacier lake outburst flood (GLOF) susceptibility at Laguna Palcacocha due to the potential of avalanching glaciers, that could induce an impulse wave and result in moraine overtopping [41–43]. A GLOF susceptibility approach, such as those of Wang et al. [44] and Bolch et al. [45] were applied by Emmer and Vilimek [46] for Laguna Palcacocha. Those first-order assessments may only be used as a basis for prioritizing more comprehensive hazard assessments at the local scale [44,45], and require region-specific adjustments to the criteria used [46]. Somos-Valenzuela et al. [42], and Frey et al. [41] defined hazard scenarios for modelling the expected hazard cascade at Laguna Palcacocha. The characteristics of the ice avalanche scenarios, for example, volume and area of initiation, were arbitrarily selected, and lack the detailed analysis concerning the glacier state and dynamics, which are required to understand ice avalanche potential. The scenario definition, in terms of volume for the nearby Laguna 513 can be traced back to those reported by Schneider et al. [47]. The ice avalanche scenarios, selected for the Palcacocha GLOF modeling are therefore poorly constrained, and may ultimately lead to a misrepresentation of the GLOF risk (i.e., particularly when considering pessimistic or worst case hazard scenarios).

We focus our evaluation of potential areas of ice avalanche initiation by identifying deviations from the linear long-term velocity trends that are not related to intra or inter-annual velocities (e.g., indications for an imminent glacier instability), as well as qualitative features of the glacier from high-resolution satellite imagery. The evaluation of ice avalanching potential is primarily concerned with the characteristics of the (thermal) contact between a glacier and the bedrock [48–50] as well as the type of starting zones recognized as being important for their initiation [49]. Type I starting zones have a relatively large area of bedrock with an almost constant slope, where glacier stability depends upon altitude [48,49], and the proportion of the glacier frozen to the bedrock [50]. Failure of cold glaciers (Type IA), which are those that are frozen to the glacier bed, are primarily controlled by intra-glacial rupture processes, and typically occur on bedrock inclinations > 45◦ [8,50,51]. Detachment in steep, cold glaciers is associated with crevasse formation and accelerating surface velocities. For polythermal or temperate glaciers (Type IB), the bedrock inclination for starting zones are >25◦, where several factors, such as adhesion to the glacier bed, proportion of glacier frozen to the bed, lateral support, and subglacial water pressure influence glacier stability [8,49,50]. Type II starting zones are associated with abrupt changes in bedrock inclination [49], where glaciers may develop a steep cliff, and unstable ice slabs form parallel to the cliff.

At the Palcaraju glacier, Type IA starting zones are expected to occur at altitudes above ~5500 m a.s.l. with a slope angle > 45◦. On the Palcaraju glacier, those areas include: (1) the South face of the Palcaraju summit, which is a hanging glacier suggested as a potential ice avalanche scenario by Vilimek et al. [43] (Figure 1, zone S1), (2) the upper part of the steep glaciated ridge in zone S2 (Figure 1), and (3) the hanging glaciers on the SW-slope of the Pucaranra West face (Figure 1, zone E1). At these locations, satellite imagery indicates evidence for relatively small snow avalanches and/or ice falls, some of which reach the relatively flat glacier area above ~5200 m a.s.l. The time series of surface velocities from these locations all show a linear trend (R<sup>2</sup> > 0.97), and the scale of acceleration/deceleration phases of the glacier correspond to seasonal dominated oscillations in S2 and inter-annual dominated oscillations in zone S1 and E2.

The steep (up to 50◦) glaciated ridge between 5200 and 5750 m a.s.l. in zones S2, show mean daily surface velocities of >0.3 m/d. Mean surface velocities of 0.2–0.3 m/d are solely found above 5500 m (Figure 1). Below 5 500 m a.s.l., where the glacier is likely to be polythermal to temperate, daily surface velocities decrease to >0.2 m/d. All velocity time series in this region exhibit a linear trend (R<sup>2</sup> > 0.98), and acceleration/deceleration phases are associated with seasonal-scale dominated oscillations.

The potential ice avalanche release area suggested by Frey et al. [41], is located 200–300 m east of the steep ridge, on a 25–30◦ inclined terrace, at approximately 5430 to 5525 m a.s.l. (Figure 1, zone S3). Satellite imagery shows that the ice thickness and area located on this terrace reduced significantly between 2013 and 2020 (Figure 8). The remaining glacier ice at this location shows mean surface velocities < 0.1 m/d with seasonal dominated oscillation in the West and inter-annual dominated oscillation in the East (Figure S1).

Below ~5400 m a.s.l., the slope of the glacier surface ranges between 1◦–35◦ with some steeper sections in the Northwest and Southeast. For zone S4a (Figure 1), where higher surface velocities are observed (i.e., >0.3 m/d), the slope parallel to the glacier flow lines range between 20–35◦ (<30◦ on average). As mentioned earlier, this zone is characterized by an elevated ice thickness and hence elevated velocities [9,33]. The surface velocity trends over the observation period are linear (R<sup>2</sup> > 0.99), and acceleration/deceleration phases correspond to inter-annual scale-dominated oscillations (Figure S2). The two zones showing higher velocities in the eastern part of the glacier (E2 and lower part of zone E1 in Figure 1) are likely associated with terrain steps in the glacier bed and multiple crevasses perpendicular to the glacier flow line. The maximum average velocity is 0.2–0.4 m/d. Below the higher velocity zones, the average daily velocity decreases to < 0.2 m/d. Linear

regression analysis shows a linear velocity trend (E1 R2 > 0.97, E2 R2 > 0.98), again, characterized by a strong seasonal-scale oscillation (Figures S3 and S4).

**Figure 8.** Satellite image of Zone S3 on (**a**) 20 August 2013 and (**b**) 3 July 2020 (Google Earth).

Somos-Valenzuela et al. [42] suggested a potential ice avalanche starting zone at 5 202 m.a.s.l. in zone S5 (Figure 1), located at the toe of a steep glaciated flank. Within the area delineated by Somos-Valenzuela et al. [42], the slope of the glacier varies between 10–30◦. Satellite imagery shows that the lower end of this zone is in close contact with the surrounding glacier without any visible terrain step, which renders any locally enhanced process kinematically unlikely. Further, the surface velocities are between 0.1–0.2 m/d, exhibiting a linear trend (R2 > 0.97) with acceleration/deceleration phases corresponding to distinct seasonal dominated oscillations (Figure S5).

Despite all accumulated surface velocity time series exhibiting linear trends (e.g., R<sup>2</sup> > 0.95), potential instabilities may form at the glacier terminus, which are characterized by glacier flow directions perpendicular to steep terrain steps (Type IB or Type II ice avalanche starting zones). The velocity-time series for the glacier terminus did not suggest the development of significant ice instabilities within the observation period. However, as our surface velocity field has a 60 m spatial resolution, instabilities smaller than the resolution may not be clearly detectable. Satellite images obtained from 2012 and 2020 show that two distinct zones at the terminus exhibit frequent ice falls or snow avalanches. One zone is located at the terminus of zone S4a and on the terminus of E3 (Figure 1). Due to the dynamics of the glacier, small icefall events at these locations are expected. Those types of events may occasionally reach the Laguna, causing small-scale disturbances, an example of which was captured by cameras operated by INAIGEM (e.g., small icefall on 5 February 2019 and a snow avalanche on 17 January 2021). However, based on available data it is not possible to ascertain whether they are single or multiple icefall events. Evaluation of ice avalanche potential for the Palcaraju glacier, based on our velocity time series indicates that an imminent failure of a glacier instability was not evident during the observation period. However, continuous observations, for which our study provides an important baseline, would be necessary to assess whether the behavior of the glacier changes over time, and to detect indications for an imminent ice avalanche.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/rs13142694/s1, Figures S1–S10 Conceptual representation of the clustering approach and displacement time series residuals for selected points.

**Author Contributions:** Conceptualization, A.K. and F.A.; methodology, A.K., F.A., F.W. and M.J.; software, T.S. and J.O.; investigation, A.K., F.A., A.D., F.W., M.J. and T.S.; data curation, J.O.; writing original draft preparation, A.K. and F.A.; writing—review and editing, A.K., F.A. T.S., J.O., F.W. M.J. and A.D.; visualization, F.W., M.J., and J.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The TerraSAR-X data are copyrighted by the German Aerospace Center (DLR) and were provided through a commercial data purchase. Permission to use the results for this publication was granted by the firm RWE AG. We acknowledge the detailed review and suggestions provided by Martin Funk, VAW ETH Zurich. We thank two anonymous reviewers for their improvements to the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.
