1. Introduction
The Cordillera Blanca (CB) is the most extensively glaciated mountain range in Peru [1] where the present-day climate is characterized by a distinct seasonality, with a dry season from May to September and a wet season from October to April. The seasonality is controlled by the oscillation of the inner-tropic convergence zone [2]. Seasonal temperature variations are small but vary significantly on a daily basis [2]. Wet season snow accumulation at high elevation and meltwater from glaciers are the main source of water [3] for the arid west coast of Peru during the dry seasons. Over the course of the 20th century, glacier mass has significantly contracted (e.g., [1,4,5,6]). Glacier volume loss may influence the stability state and dynamics of a glacier, as a result of a change of the thermal regime, stress state at the bed-ice interface, and associated changes in geometry. Numerous studies suggest that variations in glacier velocities are related to subglacial water pressures, which are controlled by seasonal meltwaters or infiltration following heavy rainstorms [7,8,9,10,11,12,13,14]. The supply of subglacial water is regulated by the extent of Firn coverage [15], as well as the characteristics of the subglacial drainage system. Pressurized subglacial water reduces the effective normal stress at the bed-ice interface, leading to reduced frictional resistance, and thus sliding rates are enhanced [11,15,16,17]. Detailed observational changes in a glacier’s dynamics as well as its geometry, provide insights into whether a potentially hazardous glacier may develop a critical stability state, or conversely remain harmless [8].
2. Materials and Methods
2.1. SAR Offset Tracking
The calculation of displacement fields using Synthetic Aperture Radar (SAR) data is possible using offset tracking methods. Offset tracking is normally used for analyzing rapidly moving objects, such as glaciers, because it overcomes the loss of signal, due to decorrelation, when using InSAR techniques over long time intervals [18,19,20,21]. As the glaciers in the Cordillera Blanca of Peru are relatively small, offset tracking is only feasible using high-resolution SAR images with a resolution of ~2–3 m. Using a unique, continuous time series of archived TerraSAR-X (TSX) radar imagery for the period 15 January 2017–6 April 2020 (3.2 years) we reconstructed ice surface velocities for the Palcaraju glacier. The dataset consists of 51 TSX radar images that were acquired on a descending orbit with a ground resolution of 3 m (Table 1). The radar images were acquired at 10.53 GMT and the distribution of re-visit intervals for the 50 image pairs was 11 days (2 pairs), 22 days (43 pairs), 33 days (1 pair), and 44 days (4 pairs).
Utilizing a normalized cross-correlation of the amplitude component of the SAR images, offsets were measured with rectangular windows at a set of positions uniformly distributed over the scene. To obtain an accurate estimate of subpixel precision of the correlation peak, correlation function values were fitted using a biquadratic polynomial surface. The time interval of the image pairs can be adjusted according to the expected maximum displacement over the glaciers from 11 days to several months. Mismatches or errors were filtered by applying a threshold to the correlation coefficient [22], by (1) iteratively discarding spurious matches based on the angle and size of displacement vectors in the surrounding areas, (2) applying a low-pass filter on the resultant displacement fields, and (3) applying a 2–98% cumulative cutoff to remove potentially uncertain velocity values and outliers. Such values cannot be realistically validated and are likely to be artifacts of the radar data processing method. Hence they have no physical meaning.
Slant range and azimuth offset displacement fields were geocoded and transformed to 3D displacements along the terrain surface using the Advanced Land Observing Satellite (ALOS) World 3D (AW3D30), Digital Elevation Model (DEM) [23]. TSX images were processed in series with offset-tracking procedures [18,19,20,21] to 3D ice surface displacement maps. This involved combining slant-range and azimuth offsets by assuming that flow occurs parallel to the ice surface, as estimated from the DEM. Matching window sizes of 128 × 96 pixels (e.g., 202 m × 192 m) were applied with steps of 16 × 12 pixels (e.g., 25 m × 24 m). The displacement maps in m/a were geocoded to a posting of ca. 60 m. For an estimation of the uncertainties in the ice displacement maps, a precision of 1/10th of a pixel in the offset estimation can be assumed [24,25]. The displacement error of TSX data with pixel sizes in ground-range and azimuth direction of 0.9 m × 2.0 m respectively, and a time interval of 11 days, are thus on the order of 10 m/year. For longer time intervals, the noise level increases, hence, a similar displacement error of about 10 m/year can be assumed for other image pairs. The displacement values were extracted from individual pixels from the spatial displacement maps.
2.2. Satellite Image Analysis
Surface features and glacier extent were mapped in a geographic information software system [26] using high-resolution satellite Pléiades imagery (Table 2). In addition, we visually inspected Google Earth and ESRI World Atlas imagery (Table 2) captured between 1999 and 2020. Aspect and slope were calculated using the AW3D30 DEM, and, together with velocity and elevation data, exported with one data point per pixel for further statistical analyses.
2.3. Statistical Analysis of Time Series
The residuals from a linear fit of the glacier surface velocities over the 3.2 year observation period were used to conduct various correlation analyses. The analyses were based on about 1300 measurement points. Data time series with less than 50% completeness were removed from the statistical analysis. For the remaining data, where gaps were evident, a linear interpolation was applied to ensure consistent time intervals. The filtered and interpolated data allowed us to perform spatial and temporal correlation analysis between the residuals, precipitation, and SST. The Pearson correlation coefficient was determined for all correlation analyses.
As an independent test, the processed time series were also subjected to a time series analysis, without a prior assumption of correlations with climate or weather data. The entire data processing workflow was performed in Python using standard packages for transformation and optimization. The dominant wavelength was extracted for each time series using a Fourier transformation. From this dominant signal, a fit to the time series was performed with a Levenberg-Marquardt or dogleg algorithm, resulting in a single fitted trigonometric function with an estimate of best-fitting values for amplitude, phase, offset, and period, as well as the covariance matrix. Physically unreasonable outliers were removed based on the following criteria:
Period > 2000 days (beyond the observation time; no periods below Nyquist frequency)
Amplitude <–20 and >20 (beyond the data range)
Phase: no obvious outliers
Offset <–4 and >4 (beyond the data range)
In total, 28 data points were excluded based on the above criteria.
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 open-source Python package bayseg, available on
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 Figure 1a and Figure 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 (Figure 1a and Figure 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.
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 (Figure 1 and Figure 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).
3.2. Velocity Variations
Linear regression analysis (see methods section) reveals a coefficient of determination R2 > 0.90–0.99 for all cumulative time series (Figure 4), indicating steady surface movements with no obvious acceleration phases within the observation period. Approximately 2.6% of the analyzed time series showed an R2 ranging between 0.90–0.96, however, this is primarily due to under-sampling of the time series where data points were missing (Figure 5).
Extracted residuals of the time series from linear fitting revealed that the residuals of the glacier surface velocities fluctuate around the long-term (3.2 years) linear trend. Spatial correlation and cluster analysis of the residuals shows that these variations in surface velocity have a spatial and temporal variability with respect to zones in the central, eastern, and western parts of the glacier, as well as zones above 5900 m a.s.l. (Figure 6).
Surface energy and glacier mass balance models (SEMB) across the glaciated Andes [1,36,37] collectively show that the glacier response has a strong spatial variability due to its sensitivity to regional and local meteorological and topographic factors [2,37,38]. In this connection, the observed spatio-temporal variation in our glacier velocity residuals may be interpreted as being closely linked to meteorological and topographical variables, which in turn influence the variability of the glacier’s surface energy, mass balance, meltwater production, and effective stress states at the glacier bed.
3.3. Spatial Correlation and Cluster Analysis
A spatial correlation and cluster analysis following the procedure described in Section 2.3 and Section 2.4 were applied to elucidate details of the factors influencing the spatio-temporal velocity patterns in different parts of the glacier (Figure 7). It is important to note that this analysis is performed independently of the previous correlation evaluations and is performed directly on the fitted time series as an additional investigation using an unsupervised machine learning algorithm. From these analyses, three clusters, each with a distinct behavior, were identified from the residuals (Figure 6b). These clusters may be considered as two end-members; an intra-annual and inter-annual cluster, with a transition cluster between the two (Figure 7c), revealing a behavior similar to the results of the correlation analysis (Figure 7a,b).
4. Discussion
4.1. Intra-Annual Dominant Cluster
The velocities belonging to the intra-annual cluster are found below 5900 m.a.s.l., on slopes with flat to medium steepness (<35°). The slopes are typically south-southeast and steeply southwest exposed. On relatively flat southwest and steep southeast exposed glacier slopes, at altitudes above 5900 m.a.s.l., the intra-annual signal weakens and becomes dominated by a pattern characteristic of the inter-annual cluster.
This cluster is characterized by glacier surface velocities that are in phase with the outer tropical seasons. The velocities generally peak in the dry season and subsequently reach a minimum in the following wet season (Figure 6b). Meteorological data from the region provided by INAIGEM (Palcacocha weather station 4607 m.a.s.l.) show relatively constant mean monthly temperatures and a distinct seasonal fluctuation in monthly precipitation (Figure 6d). Statistical analysis of the residuals from this cluster shows an anti-correlation with monthly precipitation. Higher precipitation during the wet season corresponds to lower surface velocities and vice versa. Palcaraju glacier is predominantly located above the Equilibrium Line Altitude (ELA) (i.e., 4850–4950 m.a.s.l.) where precipitation during the wet season falls as snow [38], and a temperature-sensitive increase of the snowline occurs towards the dry season. For a slightly different aspect (e.g., ~20°), but similar altitude at the Shallap glacier, (i.e., a glacier about 10 km south of the Palcaraju glacier), Gurgiser et al. [38] found surface Albedo up to 80% during the wet season and a significantly reduced Albedo during the dry season. Annual differences in the gross glacier mass balance were solely related to ablation below the ELA, which were driven by net shortwave radiation in both seasons (with stronger surface energy fluxes during the dry season), the temperature-affected altitude of the snow line, and snow quantity in the wet season (i.e., reduced albedo).
Above the ELA, the mass balance was primarily sensitive to annual cumulative precipitation. Gurgiser et al. [38] suggested that above 5000 m a.s.l. seasonal fluctuations of surface and subsurface melting may change gradually with time and altitude between the seasons. They associate this with gradual changes in the snow line, precipitation sum and type, snow cover, snow age, albedo conditions, and the net shortwave radiation budget. High ablation rates during the dry season coincide with the shortwave budget and surface/subsurface melting peaks. Thus, our surface velocity variations in the intra-annual cluster correspond with seasonal variations in ice-melt water production, which in turn modifies subglacial pressure conditions and basal motion.
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,42,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,49,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 (R2 > 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 (R2 > 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 (R2 > 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).
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., R2 > 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
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.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figures and Tables
Enlarge this image.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.
Enlarge this image.Figure 4. Spatial distribution of coefficient of determination (R2) from the linear correlation analysis.
Enlarge this image.Figure 6. (a) Cumulative surface velocities from six representative points, (b) Residuals of the linear fit displaying an intra-annual (seasonal) pattern, from three representative points above 5500 m a.s.l. (No. 143), in the East (No. 2442), and West (No. 752) of the glacier, (c) Residuals of the linear fit displaying an inter-annual pattern, from three representative points for the glacier response in the central, south-exposed part of the glacier (Nos. 702, 1522, 1574), (d) meteorological data; Temperature T, Humidity H, Precipitation P (Palcacocha weather station, altitude 4607 m.a.s.l., latitude 09°24′09.3211″S, 77°23′07.0258″, Instituto Nacional de Investigación en Glaciares y Ecosistemas de Montaña INAIGEM), and (e) Sea Surface Temperature Anomalies (SSTA) shifted by 3 months (see main text) for the NOAA Niño 3.4 region.
Enlarge this image.Figure 7. (a) Map of correlation coefficients between residuals of the long-term linear fit through the velocity time series and precipitation, measured from the Palcacocha weather station, altitude 4607 m.a.s.l., latitude 09°24′09.3211″S, 77°23′07.0258″, INAIGEM, (b) Map of correlation coefficients between residuals of the long-term linear fit through the velocity time series and SST, (c) Results of the cluster analysis show the distribution of intra, inter and transitional clusters.
Enlarge this image.Figure 8. Satellite image of Zone S3 on (a) 20 August 2013 and (b) 3 July 2020 (Google Earth).
TSX SAR offset tracking displacement maps time intervals, perpendicular baseline, and resolution.
| Scene | Image Pair | Date | Days | Perpendicular Baseline | Resolution |
|---|---|---|---|---|---|
| 1 | 6 April 2020 | 3 m | |||
| 2 | 1 | 15 March 2020 | 22 | 74.4 | 3 m |
| 3 | 2 | 31 January 2020 | 44 | 1.5 | 3 m |
| 4 | 3 | 9 January 2020 | 22 | −146.8 | 3 m |
| 5 | 4 | 18 December 2019 | 22 | 104.7 | 3 m |
| 6 | 5 | 26 November 2019 | 22 | −112.4 | 3 m |
| 7 | 6 | 24 October 2019 | 33 | −44.6 | 3 m |
| 8 | 7 | 02 October 2019 | 22 | 136.4 | 3 m |
| 9 | 8 | 10 September 2019 | 22 | −137.0 | 3 m |
| 10 | 9 | 19 August 2019 | 22 | 203.5 | 3 m |
| 11 | 10 | 28 July 2019 | 22 | −125.7 | 3 m |
| 12 | 11 | 14 June 2019 | 44 | 46.0 | 3 m |
| 13 | 12 | 23 May 2019 | 22 | −14.1 | 3 m |
| 14 | 13 | 1 May 2019 | 22 | 29.9 | 3 m |
| 15 | 14 | 9 April 2019 | 22 | −4.2 | 3 m |
| 16 | 15 | 29 March 2019 | 11 | −103.0 | 3 m |
| 17 | 16 | 07 March 2019 | 22 | 6.1 | 3 m |
| 18 | 17 | 13 February 2019 | 22 | −78.5 | 3 m |
| 19 | 18 | 22 January 2019 | 22 | 174.7 | 3 m |
| 20 | 19 | 31 December 2018 | 22 | −88.2 | 3 m |
| 21 | 20 | 09 December 2018 | 22 | 84.8 | 3 m |
| 22 | 21 | 17 November 2018 | 22 | −25.8 | 3 m |
| 23 | 22 | 26 October 2018 | 22 | 71.4 | 3 m |
| 24 | 23 | 4 October 2018 | 22 | −222.0 | 3 m |
| 25 | 24 | 12 September 2018 | 22 | 280.1 | 3 m |
| 26 | 25 | 21 August 2018 | 22 | −164.2 | 3 m |
| 27 | 26 | 30 July 2018 | 22 | 148.7617 | 3 m |
| 28 | 27 | 8 July 2018 | 22 | 26.9 | 3 m |
| 29 | 28 | 16 June 2018 | 22 | −17.8 | 3 m |
| 30 | 29 | 25 May 2018 | 22 | 4.0 | 3 m |
| 31 | 30 | 3 May 2018 | 22 | −66.2 | 3 m |
| 32 | 31 | 11 April 2018 | 22 | 106.8 | 3 m |
| 33 | 32 | 20 March 2018 | 22 | −220.7 | 3 m |
| 34 | 33 | 26 February 2018 | 22 | −65.0 | 3 m |
| 35 | 34 | 13 January 2018 | 44 | 133.6 | 3 m |
| 36 | 35 | 2 January 2018 | 11 | −149.4 | 3 m |
| 37 | 36 | 11 December 2017 | 22 | 30.0 | 3 m |
| 38 | 37 | 19 November 2017 | 22 | −74.9 | 3 m |
| 39 | 38 | 28 October 2017 | 22 | 169.1 | 3 m |
| 40 | 39 | 14 September 2017 | 44 | −5.5 | 3 m |
| 41 | 40 | 23 August 2017 | 22 | 246.2 | 3 m |
| 42 | 41 | 1 August 2017 | 22 | 133.1 | 3 m |
| 43 | 42 | 10 July 2017 | 22 | −206.8 | 3 m |
| 44 | 43 | 18 June 2017 | 22 | 182.6 | 3 m |
| 45 | 44 | 27 May 2017 | 22 | −90.0 | 3 m |
| 46 | 45 | 5 May 2017 | 22 | −87.2 | 3 m |
| 47 | 46 | 13 April 2017 | 22 | 107.1 | 3 m |
| 48 | 47 | 22 March 2017 | 22 | 139.4 | 3 m |
| 49 | 48 | 28 February 2017 | 22 | −111.8 | 3 m |
| 50 | 49 | 6 February 2017 | 22 | −297.0 | 3 m |
| 51 | 50 | 15 January 2017 | 22 | 291.0 | 3 m |
Satellite imagery analyzed in this study.
| Image Capture Date | Image Source | Resolution (m) | Accuracy (m) |
|---|---|---|---|
| 3 July 2020 | Google Earth | ||
| 24 February 2020 | Pleiades | ||
| 22 June 2019 | Google earth | ||
| 08 June 2019 | Pleiades | ||
| 24 June 2018 | Maxar/DigitalGlobe (WV02) * | 0.5 | 8.47 |
| 18 June 2017 | Google Earth | ||
| 30 May 2016 | Google Earth | ||
| 30 May 2016 | Pleiades | ||
| 01 June 2016 | DigitalGlobe (WV02) * | 0.5 | 10.2 |
| 28 April 2013 | Google Earth | ||
| 28 July 2013 | Google Earth | ||
| 5 May 2012 | Google Earth | ||
| 15 January 1999 | Earthstar Geographic (Terracolor) * | 15 | 50 |
| 31 December 1985 | Google Earth | ||
| 1 January 1970 | Google Earth |
* as published on
© 2021 by the authors.
Details
; Amann, Florian 2 ; Strozzi, Tazio 3
; Osten, Julian 2 ; Wellmann, Florian 4
; Jalali, Mohammadreza 2 ; Dufresne, Anja 2 1 Terrasense Switzerland Ltd., 9470 Buchs, Switzerland
2 Department of Engineering Geology and Hydrogeology, RWTH University, 52062 Aachen, Germany;
3 Gamma Remote Sensing Ltd., 3073 Gümligen, Switzerland;
4 Unit of Computational Geoscience & Reservoir Engineering, RWTH University, 52062 Aachen, Germany;