*2.4. Cluster Analysis of the Time Series*

We employed a Bayesian unsupervised machine learning approach [27], to extract clusters of similar behavior from our surface displacement variations (i.e., residuals). The primary aim of the employed method was to extract clusters from the data based on feature similarity, with consideration of the spatial configuration in the physical space. This combination is enabled through a combination of a Gaussian Mixture Model (GMM) in feature space with a Hidden Markov Random Field (HMRF) model in physical space. Both models are integrated into a Bayesian model, and MCMC sampling is performed to obtain the posterior distribution of the GMM parameters.

The GMM was applied to samples of N data points (image pixels) in M-dimensional feature space to determine a set of multivariate Gaussian distributions (L classes). The distributions were parametrized by their mean, μ, and covariance, Σ for each cluster. As features originate from spatial data sets, the HMRF was then used to consider spatial dependencies of sampling points, with a smoothing coefficient, which parameterized the strength of spatial correlation, independently for each class.

The model parameters (μ, Σ, β), as well as the latent field values, were obtained through a Bayesian optimization, iterative sampling process using a Markov Chain Monte Carlo (MCMC) approach, following an initial Expectation-Maximization (EM) step. The final association of sampling points to clusters were obtained from the Maximum a posteriori (MAP) distribution. In addition, cluster assignment probabilities were obtained from the MCMC chain and combined using information entropy, to obtain a spatial estimate of uncertainty [28], for the assigned cluster values.

The segmentation procedure is described in Wang et al. [27] and has previously been applied to geophysical [29] and arctic sea ice [30] data sets. It is implemented in the opensource Python package bayseg, available on https://github.com/cgre-aachen/bayseg, accessed on 10 June 2021. The clustering workflow is presented in an exemplary form in the supplementary information section (Figure S1a–d).

#### **3. Results**

#### *3.1. Spatio-Temporal Distribution of Glacier Velocities*

We calculated the spatial distribution of average daily surface velocities of the Palcaraju glacier over 50 measurement intervals (i.e., velocity maps with 1300 pixels per interval), for the 3.2 year observation period. As the spatial distribution of individual velocity maps were similar for all 50 intervals, the velocity distribution shown in Figure 1 is an appropriate representation for the observation period. Daily surface velocities averaged over 3.2 years ranged between 0.01 and 0.47 m/day, where both their spatial distribution and magnitude are consistent with results obtained from Sentinel-2 data [31]. The displacement time series, representing the residuals of the linear fit and cumulative surface displacements, for selected points shown in Figure 1, are available in the supplementary materials (Figures S2–S10).

The resulting surface velocities are a combination of basal sliding as well as creep and ablation of the glacier, which cannot be individually quantified without independent measurements. Such measurements are outside the scope of this study. Values from the literature [13] and field measurements [7,9,32] suggest that the contribution of basal sliding to the total surface velocities varies between 50–90%, and is spatially variable across the glacier bed [32].

Higher mean annual surface velocities are observed in areas with a high frequency of crevassing and or steep glacier topography (>30–40◦) (see zones S2, S4a, S4b, E1, E2, and E3 Figures 1a and 2). Based on the occurrence of ridgelines and terrain steps in the 30 m DEM, we infer that steep glacier sections directly reflect the underlying topography of the glacier bed. Although the slope gradient in Zone S4a (Figures 1a and 2) is comparatively flat (20–35◦), such higher velocity rates can be related to larger ice thickness at this location [9,33].

West exposed slopes (Figure 3) typically showed the highest mean annual surface velocities, followed by south and east exposed slopes over the 3.2 year observation period. Surface velocities of south exposed slopes tend to decrease with altitude, while surface velocities of East exposed slopes increase with altitude. West exposed slopes also showed a general increase with altitude, however, the trend is reversed at an elevation range of approximately 5400–5700 m.a.s.l.

**Figure 1.** (**a**) Spatial distribution of mean daily surface velocities with respect to elevation, showing the zones and point selections used in this study (for location see insert map), (**b**) Spatial distribution of contoured mean daily surface velocities on the Palcaraju glacier. Elevation range of zone E1 is 5 400–6050 m a.s.l, and 5 900–6150 m a.s.l for zone S1.

**Figure 2.** Glacier slope map.

**Figure 3.** Glacier aspect map.

Zones with higher surface velocities for all slope aspects were only found below 5700 m a.s.l., suggesting that the steep hanging glacier in E1 and S1 (Figures 1 and 3) are cold glaciers [34], likely to be frozen to the glacier bed. This is supported by estimates of the cold-temperate transition line at approximately 5500 m a.s.l., elevation in the region [35]. Hence, in these areas, the surface velocity field mainly represents creep deformation within the ice body. Below 5700 m a.s.l., the glacier is likely to be polythermal to temperate (i.e., subglacial liquid water is temporarily present), and therefore the measured surface velocities are likely controlled by creep processes, basal sliding along the ice-bed interface, and an indeterminable error due to ablation (e.g., a shift in the slant range geometry between image acquisitions and consequent projection into DEM predating the radar image acquisition, resulting in an over-estimation of surface velocities).
