1. Introduction
Despite accounting for only around 1% of the planet’s total mountain glacier cover (1928 km2 of 706 744 km2 in the Randolph Glacier Inventory (RGI) 7.0), high-altitude tropical glaciers play an important role in sustaining the flow of rivers in heavily populated catchments in countries such as Peru and Bolivia, particularly during the dry season, and most importantly in drought years (Immerzeel and others, Reference Immerzeel2020). The meltwater produced by glaciers in the Outer Tropics is also important for unique but delicate high-altitude ecosystems, which themselves are the host of both arable and pastoral farmland for local populations. However, tropical glaciers have experienced substantial recession (Rabatel and others, Reference Rabatel2013; Seehaus and others, Reference Seehaus2019) and mass loss (Dussaillant and others, Reference Dussaillant2019) in recent decades, and their ability to sustain river flow, particularly during the dry season, is expected to decrease in the coming decades (Huss and Hock, Reference Huss and Hock2018). This will increase water stress experienced by local populations and ecosystems, meaning that appropriate water management strategies must be established.
Glacier modelling is key to the development of robust freshwater management strategies in high-mountain catchments. A variety of ablation processes must be understood, parameterised and calibrated against observations if such models are to produce actionable estimates of future meltwater yield. The ablative processes in operation on tropical glaciers in Peru are heavily influenced by the interplay between short-term (daily to weekly) precipitation dynamics and net shortwave radiation (Fyffe and others, Reference Fyffe2021). High net shortwave radiation fluxes drive substantial ice loss even in the presence of sub-zero air temperatures, particularly when bare-ice conditions present a low albedo glacier surface. Whilst progress has been made recently on understanding the energy and mass balance of tropical glaciers (Perry and others, Reference Perry2017; Fyffe and others, Reference Fyffe2021; Mackay and others, Reference Mackay2025), more detailed examination of specific ablative processes operating at the glacier surface is so far limited, particularly those linked to temporally variable albedo (Mackay and others, Reference Mackay2025).
Supraglacial pond and ice cliff networks on debris-covered glaciers are ablative features that have received considerable recent attention (e.g. Rowan and others, Reference Rowan2021) because they are ablation hotspots over parts of glaciers that would otherwise be protected from substantial melt by the thick debris mantle. On debris-covered glaciers, supraglacial pond and ice cliff networks can account for ∼40% of ice loss in ablation zones despite occupying <10% glacier surface (Thompson and others, Reference Thompson2016). The presence of such features is one of the primary reasons put forward to explain similarities in debris-covered glacier and clean ice mass balance over recent decades, despite the presence of a protective layer of debris on the former (Brun and others, Reference Brun2019). Melt enhancement factors associated with ponds and cliffs on debris-covered glaciers are now increasingly being incorporated into glacier models as a result (Rowan and others, Reference Rowan2021). However, the impact of the development of similar features on clean-ice glaciers has not been studied to the best of our knowledge.
In this study, we examine the impact of the development of supraglacial ice cliffs and ponds on clean-ice glaciers, which are previously unreported in the literature. We focus on selected glaciers in the Cordillera Vilcanota (Fig. 1), a mountain massif in the outer tropics of Peru. We document the extent and geometry of cliffs and ponds and propose a mechanism by which they form and evolve. We also aim to quantify the impact their development has on glacier recession and mass balance. Finally, we relate their development to longer-term changes in glacier surface conditions in the region and discuss their potential role in glacier recession in the future.
The regional setting of the Cordillera Vilcanota in the Outer Tropics of South America (a); the location of the Suyuparina–Quisoquipina catchment within the Cordillera Vilcanota (b); and the three glaciers examined in more detail in the study — Suyuparina, Quisoquipina and Jurcay glaciers (c).
2. Study area
Glaciers cover ∼246 km2 of the Cordillera Vilcanota (Taylor and others, Reference Taylor2022), predominantly above an elevation of around 4500 m, with the region’s highest glacierised peak (Nevado Ausangate) reaching 6384 m a.s.l. Meltwater produced by glaciers in the Cordillera Vilcanota primarily drains into the Vilcanota-Urubamba catchment (Fig. 1a). Substantial reductions in glacier area have occurred in the Cordillera Vilcanota over recent decades. INAIGEM (2023) mapped the extent of 380 glaciers in the Cordillera in 1962 (495.05 km2) and then again in 2020, identifying a reduction in glacier area of 251.76 km2, or 51%, over a 58 year period. Hugonnet and others (Reference Hugonnet2021) and Taylor and others (Reference Taylor2022) derived almost identical estimates of the mean mass balance of glaciers in the Cordillera Vilcanota during 2000–20; −0.47 ± 0.07 and −0.47 ± 0.05 m w.e.a−1, respectively. This estimate is in line with the rate of mass loss (−0.45 m w.e.a−1) experienced by low-latitude glaciers collectively over the same period (Hugonnet and others, Reference Hugonnet2021).
Within the Cordillera Vilcanota, we focus more specifically on clean-ice glaciers within the Suyuparina catchment, an area of the Cordillera which is relied on by local people for freshwater supply and as grazing areas for livestock and native fauna such as Llama (Lama glama), Alpaca (Lama pacos) and Vicuña (Lama vicugna). The hydrology of the catchment has been studied in some detail (e.g. Gribbin and others, Reference Gribbin2024), and several studies have examined the links between glaciers and bofedales, high altitude wetlands, which are prevalent in the area (Davies and others, in review a, in review b). At the head of the catchment (−13.7984°, −70.8878°), the three largest glaciers (Suyuparina, Quisoquipina, Jurcay) now cover 0.97, 2.96 and 1.75 km2, respectively, and cover an altitudinal range of 5155–5657 m a.s.l. These three glaciers have distinct, heavily pitted ablation zones (Fig. 2), which host a substantial number of supraglacial ponds and present highly irregular glacier margins, all of which indicate that they experience spatially focused ablative processes. The occurrence and driving mechanisms of such ablative features form the focus of the rest of the study. Elsewhere in the Cordillera Vilcanota, similar surface features to those in the Suyuparina catchment are found predominantly on glaciers situated on the southern side of the main orographic divide and on the lower reaches of some glaciers draining the Quelccaya Ice Cap (Fig. 1). The features are generally all found on glacier surfaces above an elevation of ∼4500 m a.s.l., are found down-glacier of large crevasse fields, and migrate towards glacier margins.
The ablation zone of Suyuparina Glacier, captured by oblique drone imagery acquired in September 2023 (a), shows a heavily pitted surface and irregular glacier margin. Surface metrics extracted from a 2024 Pléiades DEM covering the wider catchment, including surface slope (b), the hillshade of the DEM (c) and Topographic Position Index (TPI) (d), which were used to aid the delineation of surface features.
3. Methods and data
We examine the evolution of ablative features over the lower reaches of three glaciers in the Suyuparina catchment using a combination of airborne and satellite-derived remotely sensed data spanning the period 1977–2024. We assembled a time series of optical stereo imagery of high enough resolution to allow for the examination of features typically a few tens of metres wide and deep (Table 1). We applied various processing techniques to different imagery types, which are outlined below.
Sources of images used to derive topographic datasets, their resolution and purpose in later analyses.
3.1. Satellite remote sensing data and processing
We generated digital elevation models (DEMs) and orthoimages from pairs of high-resolution stereo images acquired by the American Hexagon KH-9 satellites and by the Pléiades 1B satellite. The Hexagon KH-9 satellite series operated during 1972–84 and carried two cameras, a moderate resolution mapping camera which acquired overlapping images over a footprint of ∼130 × 260 km, and two panoramic cameras arranged in a stereo configuration, which acquired images of contrasting geometry (250 × 20 km) but at much higher resolution (∼0.6 m at the centre of the scene), making this image archive more suitable for this work. Images acquired by cameras aboard the Hexagon KH-9 satellite series were stored as film negatives and jettisoned from the satellites to be retrieved during their return to the Earth’s surface. A copy of this film archive is now maintained by the United States Geological Survey and was made publicly available in stages during 2002–11, with scanning ongoing. We requested the scanning of five Hexagon KH-9 panoramic-camera images (scene IDs D3C1213-200305A001-3 and D3C1213-200305F001/2) at a resolution of 7 microns. We followed the approach of Ghuffar and others (Reference Ghuffar, Bolch, Rupnik and Bhattacharya2022) to generate a DEM and orthomosaic of the panoramic-camera image pair using Sentinel 2 imagery and the Copernicus 30 m global DEM as reference datasets to enable the geolocation of the relative DEM and orthomosaic generated using the Corona Stereo Pipeline (Ghuffar and others, Reference Ghuffar, Bolch, Rupnik and Bhattacharya2022). We generated a 2 m orthoimage and 10 m DEM from the Hexagon KH-9 panoramic-camera image set, which is a suitable choice of resolution for an AOI situated towards the centre of KH-9 panoramic-camera imagery (Ghuffar and others, Reference Ghuffar, King, Guillet, Rupnik and Bolch2023).
The Pléiades satellites (1A and 1B) operate in a synchronous orbit, collecting 0.7 m, 12 bit stereo and tri-stereo panchromatic images of targeted parts of the Earth’s surface every 26 days. The Pléiades image archive has been used to examine a variety of small-scale glacier processes (e.g. Berthier and others, Reference Berthier2014; Miles and others, Reference Miles2018a; Falaschi and others, Reference Falaschi2022; Brun and others, Reference Brun2023) owing to its high resolution and ability to resolve surface details on flat and textureless terrain, such as over glacier accumulation zones. We selected a pair of images from the Pléiades 1B archive over our study site acquired on 7 August 2016, and acquired an additional pair of images over the study area on 11 September 2024. Both of these scenes were acquired towards the end of the dry season (August/September), so they should have captured glacier surface conditions most favourable for DEM generation, i.e. minimal snow and cloud cover and good optical contrast. We generated DEMs and orthoimages from both pairs of Pléiades 1B images without ground control points (GCPs), as accurate, well-distributed point elevation data were unavailable over the study area. Photogrammetric processing was undertaken using CATALYST Earth software (formerly PCI Geomatica), in which we employed the Semi Global Matching algorithm (Hirschmüller, Reference Hirschmüller2005) implemented in the Ortho Engine module to precisely match subsets of the stereo images. Horizontal and vertical shifts of ∼3 m or less can be expected from Pléiades DEMs generated without GCPs (Berthier and others, Reference Berthier2014), although such shifts between pairs of DEMs can be removed using established techniques, described later. We generated 1 m DEMs and orthoimages from Pléiades stereo pairs from August 2016 and September 2024 to allow for further analyses of small-scale (tens of metre size and relief) glacier features. Pléiades DEMs were resampled to 10 m when analysing elevation change between the Hexagon KH-9 panoramic-camera DEM and Pléiades DEMs.
3.2. Drone surveys
The lower reaches of the Suyuparina Glacier (Fig. 1c) were surveyed using a DJI Phantom 4 Pro during the dry season in 2022, 2023 and 2024. Survey areas varied slightly on each occasion but allowed for the detailed study of ice-marginal processes on an annual basis; we therefore use these data to examine the impact of the migration of supraglacial features towards the periphery of the glacier and their impact on ice margin recession. We processed each drone survey using Agisoft Metashape, following a well-established workflow (Smith and others, Reference Smith, Carrivick and Quincey2016) to align sequential images across each full drone survey, produce a sparse point cloud, which we constrained the geolocation of using GCPs (n = 8) extracted from the 2024 Pléiades DEM and orthoimage and then generate a dense point cloud and resulting orthoimage (0.15 m) and DEM (1 m). The root mean squared error of control points following these processing steps in x, y and z directions was 1.28, 1.62 and 1.84 m in the case of the 2022 drone DEM, 1.03, 0.98 and 1.45 m in the case of the 2023 drone DEM and 1.48, 1.42 and 1.89 m for the 2024 DEM, respectively. The flying height of each survey was 230, 182 and 195 m in 2022, 2023 and 2024, respectively.
3.3. Spatial and geometric attributes of supraglacial ponds and cliffs
We manually mapped structures and features present at the surface of glaciers in the Suyuparina–Quisoquipina catchment from orthoimages generated from Hexagon KH-9 Panoramic-Camera scenes acquired in 1977, and Pléiades 1B images acquired in 2016 and 2024. We followed established mapping criteria (Jennings and Hambrey, Reference Jennings and Hambrey2021) to identify and map features including folded primary stratification, crevasses of various arrangements, supraglacial ice cliffs and ponds and crevasse traces (Table 2). We draw inferences about the mechanisms by which different features are formed based on the stress and strain regime they typically represent to be in operation at the glacier surface.
Features and structures present at the surface of Suyuparina–Quisoquipina Glaciers and their appearance on Hexagon KH-9 and Pléiades orthoimages, ordered in their likely sequence of formation.
Spatially variable ablation over glacier surfaces can lead to the development of complex surface topography (e.g. King and others, Reference King, Turner, Quincey and Carrivick2020a), which can be characterised using various topographic metrics to distinguish features attributable to different processes. We primarily use surface slope, aspect and the Topographic Position Index (hereafter TPI) (Weiss, Reference Weiss2001) extracted from satellite and drone derived DEMs to examine the geometry of features which developed over the surface of glaciers in the Suyuparina–Quisoquipina catchment (Fig. 2). The TPI represents differences between the elevation of a central pixel with the average of that of surrounding pixels and is therefore a good expression of local topographic highs and lows. We used window sizes equating to 30 × 30 m for DEMs of different resolutions when calculating the TPI. These topographic metrics also assisted in the identification and mapping of features listed in Table 2.
3.4. Ice loss through DEM differencing
To establish the impact of cliff and pond development on the evolution of glaciers in the Suyuparina–Quisoquipina catchment, we generated a time series of surface elevation change data from the various DEMs generated from satellite and drone imagery. Different combinations of DEMs were best suited to examining the development of these features over different timescales. Elevation change data derived from the Hexagon KH-9 panoramic-camera DEM (1977) and Pléiades DEMs (2016, 2024) covers a period when substantial ice loss has occurred in the Cordillera Vilcanota, thus ice loss captured by these data is likely of greater magnitude than ice loss driven by smaller-scale processes such as cliff and pond development. We therefore use elevation change data derived from Hexagon KH-9 panoramic camera and Pléiades DEMs to illustrate the longer-term behaviour of glaciers in the catchment in response to warming of the regional climate (Supplementary Material, Salzmann and others, Reference Salzmann2013). We focus our analyses of the impact of cliff and pond development on glacier mass balance on the pair of Pléiades DEMs, as the timespan between these two acquisitions (2016–24) is appropriate to capture both ice loss related to climate and ice loss caused by localised processes. The annual timestep of the drone-derived DEMs and orthoimages allowed us to trace the impact of the development of individual features as they migrated towards the margin of the Suyuparina Glacier.
The generation of robust surface elevation change data relies on the removal of a variety of biases that may impact DEMs, particularly on those generated using different methods, as in this study. The characteristics of elevation differences over terrain known to have been stable over the duration of a DEM time series provide the best indication of the presence of any biases. We removed co-registration errors using established methods (Nuth and Kääb, Reference Nuth and Kääb2011) implemented in xDEM (xDEM contributors, 2024, https://zenodo.org/records/11492983), primarily evident between the Hexagon KH-9 PC DEM and Pléiades DEMs owing to the contrasting use of reference data during DEM generation. Elevation differences over stable ground had mean values of −0.16 m (normalised median absolute deviation [NMAD] 1.23 m) in the case of the KH-9 PC and 2024 Pléiades DEM pair, and −0.08 m (NMAD 0.45 m) for the 2016 and 2014 Pléiades DEM pair, indicating good coregistration of these sets of DEMs.
Following DEM co-registration and DEM differencing, the identification and removal of outliers evident in elevation change data is necessary to ensure no subsequent impact on estimates of glacier thinning or mass loss. This was primarily required in the case of elevation change estimates derived from the Hexagon KH-9 PC and 2024 Pléiades DEM. We excluded values outside of the range ± 100 m under the assumption that these are related to DEM errors. We then further refined the dH data by filtering to within ± 3σ within 50 m elevation bands through the full altitudinal range of different glaciers (Gardelle and others, Reference Gardelle, Berthier, Arnaud and Kääb2013). Any resulting gaps were filled with mean values of surface elevation change estimates within the same 50 m elevation band as data gaps. Geodetic glacier mass balance estimates were then derived from the filtered and filled elevation change grids using a conversion factor of 850 kg m−3 (Huss and others, Reference Huss2013). We used the methods developed by Hugonnet and others (Reference Hugonnet2022) to estimate uncertainty associated with dH data and propagated this error in combination with 7% uncertainty associated with the conversion of ice volume change to mass change to calculate uncertainty associated with glacier mass balance estimates (Huss, Reference Huss2013). We do not consider uncertainty associated with glacier area estimates substantial enough to impact the mass balance calculation, due to the high resolution of the datasets used in glacier margin delineation (1–2 m orthoimages) (Table 1).
3.5. Ancillary datasets
To understand how climatic conditions may have aided the development of features present at the surface of the three studied glaciers, we have examined available records of meteorological data gathered by three nearby weather stations. Here we primarily focus on records of radiative fluxes, in particular solar irradiance data from Apogee CS-300 sensors at nearby stations on the Quelccaya Ice Cap (5650 m; Perry, Reference Perry2017) and at Lake Sibinacocha (4890 m) and measured by a Hukseflux NR01 four-component net radiometer in the forefield of the Suyuparina and Quisoquipina Glaciers (5153 m) (Fig. 1). These stations provide observations of incoming shortwave radiation during 2017–24 at Sibinacocha, during 2014–22 close to the summit of Quelccaya and during May 2024–July 2025 in the forefield of Suyuparina and Quisoquipina Glaciers. The available datasets from Sibinacocha and Quelccaya provide hourly average (mean) and maximum incoming and outgoing shortwave flux estimates. The records from Suyuparina and Quisoquipina are available at higher temporal resolution, but we have used hourly average values to allow for direct comparison with Quelccaya and Sibinacocha. We use these to calculate daytime (10:00–21:00) average and maximum shortwave fluxes, which will have provided energy to melt local glacier surfaces over each respective time series. We also illustrate single hourly incident shortwave maxima recorded by each station at each site.
4. Results
4.1. Glacier surface features and structures
We identified 266 individual ice cliffs on glacier surfaces in the Suyuparina catchment in the 1977 Hexagon PC orthoimage (Fig. 3, Table 3): 81 on Suyuparina Glacier, 125 on Quisoquipina Glacier and 60 on Jurcay Glacier. The number of cliffs approximately doubled during 1977–2016, when we identified 525 ice cliffs across the same glaciers: 129 on Suyuparina Glacier, 263 on Quisoquipina Glacier and 112 on Jurcay Glacier. We identified 457 cliffs in the 2024 Pléiades orthoimage, a decrease of 68 overall, 116 of which were on Suyuparina Glacier, 221 on Quisoquipina Glacier, and 120 on Jurcay Glacier. The number and area of supraglacial ponds increased sharply through our study period — in 1977 only two small (576 m2) ponds could be identified on Suyuparina Glacier, whereas in 2016 119 ponds covered 36 984 m2 of glacier surfaces (21 covering 5217 m2 on Suyuparina Glacier, 69 covering 23 181 m2 on Quisoquipina Glacier and 29 covering 8586 m2 on Jurcay Glacier). We identified a smaller total area and number of supraglacial ponds in 2024 (93 covering 29 725 m2) compared to 2016, although primarily only on Quisoquipina Glacier (50 covering 16 396 m2 in 2024), with pond area similar on Suyuparina Glacier (16 covering 5205 m2) and Jurcay Glaciers (27 covering 8124 m2). We note slightly contrasting snow cover conditions during 2016–24 Pléiades orthoimages due to their acquisition dates, with greater snow cover and therefore meltwater availability in the 2016 acquisition, which likely resulted in increased pond area at this point in the study period. Finally, we also saw a substantial increase in the number of mapped crevasses, with 621 spread across glacier surfaces in 1977, 1222 in 2016 and 1409 in 2024. Again, the true number of crevasses present on each glacier at each point in time is obscured by snow cover, which was extensive at the point of the 1977 Hexagon PC images, and slightly more prevalent in the 2016 Pléiades orthoimage compared to that of 2024.
Glacier structures present at the surface of the Suyuparina, Quisoquipina and Jurcay Glaciers (left to right) in 1977 (a) and 2024 (b).
Number, area and % area of glacier surface features mapped from orthoimages in 1977, 2016 and 2024.
The estimation of the total area covered by supraglacial cliffs on a glacier’s surface depends somewhat on the method used to define the spatial limit of such features. Here, we measure cliff area as any ice surface of a slope greater than 30° (derived from the Pléiades DEMs), within a 20 m buffer of manually digitised cliff crests. This method yielded an estimate of cliff area ranging 4.04–5.59% of total glacier area for Suyuparina, Quisoquipina and Jurcay Glaciers in 2024 (Table 3). The area covered by ponds and cliffs increased markedly during 1977–2016 across all three glaciers (Table 3) but decreased on the Suyuparina and Quisoquipina Glaciers during 2016–24, likely due to the melt-out of large cliff/pond features at the glacier margin over this timeframe. Mapped cliff crests had a mean length of 68.3 m in 1977, 79.1 m in 2016 and 75.4 m in 2024. The geometry of supraglacial ice cliffs evolves as they move down-glacier. Over each of Suyuparina, Quisoquipina and Juray Glaciers, cliffs initially predominantly face north in the higher reaches (Fig. 4), but their most common aspect evolves down-glacier, with most cliffs facing broadly south below a certain elevation on each glacier. It is likely that some of the steep, north-facing cliff faces found in the glaciers’ higher reaches have paired south-facing ice faces, given their original formation as crevasse walls. As with our crevasse mapping, it is likely that snow cover has impacted our ability to fully capture south-facing cliff geometry in the highest parts of the three studied glaciers, predominantly in the higher reaches of glaciers in the 2016 Pléiades scenes.
Aspect of supraglacial cliffs on the Suyuparina, Quisoquipina and Jurcay Glaciers, both above and below certain elevations, where their orientation is contrasting. Cliff orientation was derived from the 2024 Pléiades DEM.
4.2. Changes in glacier area and surface elevation
Glaciers draining directly into the Suyuparina catchment (those shown in Fig. 5a) covered an area of 13.94 km2 in 1977. By 2024, the total area covered by the same glaciers had reduced to 8.51 km2, −39% compared to their original 1977 area. The number of individual glaciers in the catchment increased by 2 (9–11) due to a glacier on the western edge of the catchment fragmenting as it has receded. One small glacier, immediately south of Jurcay Glacier, disappeared during 1977–24 (Fig. 5a and b). The differencing of DEMs illustrates the substantial ice loss experienced by all glaciers in the Suyuparina catchment through the study period (Fig. 5). For 1977–2024, glacier surfaces decreased in elevation by a mean of 22.9 m (those shown in Fig. 5a), although thinning over the lower reaches of the larger glaciers in the catchment was typically 80 m or more. The estimated mass balance of Suyuparina, Quisoquipina and Jurcay Glaciers was −0.61 ± 0.04, −0.35 ± 0.03 and −0.53 ± 0.04 m w.e.a−1, respectively, over this period.
Surface elevation change rates (dHdT) over the Suyuparina catchment between (a) 1977 and 2024, (b) 2016 and 2024, and in more detail over the Quisoquipina (c) and Jurcay (d) Glaciers over the 2016–24 periods. Note the different value ranges for panel a.
The differencing of Pléiades DEMs during 2016–24 shows a contrasting pattern of ice loss over the Suyuparina, Quisoquipina and Jurcay glaciers over this shorter contemporary period. Whereas the pattern of thinning is clearly elevation-dependent during 1977–2024 (Fig. 5a), an additional pattern of spatially restricted and enhanced ice loss is evident in the elevation change maps in Fig. 5b. These areas of heightened thinning are coincident with the location of supraglacial cliffs and ponds (Fig. 5c and d). The mean thinning rate within cliff and pond extents (Table 3) was −2.48, −1.83 and −2.31 ma−1 for Suyuparina, Quisoquipina and Jurcay Glaciers, respectively, during 2016–24, compared to −1.13, −0.95 and −1.35 ma−1 over the rest of each glacier’s surface.
On Suyuparina Glacier, during 2022–24 elevation change data from drone surveys (Fig. 6) shows that thinning over the lowermost reaches of the glacier is concentrated preferentially on the northwest-northeast facing slopes of reclined ice cliffs that have developed further up-glacier, as well as the topographic highs between neighbouring cliff and pond complexes (Fig. 5b and c). Thinning rates over ice surfaces oriented between 292.5° (NW) and 67.5° (NE) were 31% higher during 2022–24 (−3.93 ma−1), compared to ice surfaces facing predominantly south (−3.00 ma−1) over the terminus of Suyuparina Glacier. Furthermore, ice surfaces of positive relief relative to their surroundings (TPI < 0) lowered at a mean rate of −12.53 ma−1 over the lower reaches of Suyuparina Glacier, whereas ice surfaces in topographic lows (TPI > 0) lowered at a mean rate of −1.76 ma−1 during 2022–24 (Fig. 6).
Surface elevation change during 2022–24 over the terminus of Suyuparina Glacier (a); TPI (b) and aspect, highlighting predominantly north-facing slopes (c), of the glacier surface over the same area in 2022, and (d) the relationship between elevation change and bins (45° for aspect, 0.5 for TPI) of the metrics shown in (b) and (c). Box plots in (d) show the median, interquartile range (IQR), 1.5* the IQR (whiskers) and values outside of these bounds as crosses.
5. Discussion
5.1. Supraglacial ice cliff formation
Through examining changes in the type and distribution of glacier surface structures (Fig. 3), the geometry of glacier surfaces and spatial patterns of ice loss over time, we propose a mechanism of formation (Fig. 7) of supraglacial cliffs and ponds over the Suyuparina, Quisoquipina and Jurcay glaciers.
Conceptual figure illustrating the proposed mechanism of formation of supraglacial ice cliffs and ponds on clean-ice glaciers. Schematic (left) indicates three apparent stages of formation, illustrated using subsets of the 2024 Pléiades orthoimagery (right).
Structural mapping suggests that supraglacial cliffs on these glaciers originate from fields of large crevasses in each glacier’s higher reaches (Fig. 3) (cf. Kneib and others, Reference Kneib2023). These large crevasses are likely formed when each glacier flows over bedrock steps, which are characteristic of the basalt-dominated landscape of the Cordillera Vilcanota (Davies and others, in review a). Both in each glacier’s lower reaches and immediately down-glacier of large crevasses fields, arcuate 10–125 m long surface depressions form, which possess a steep (<45°), north-facing backwall (Fig. 4). We interpret these features as the earliest instance of supraglacial ice cliffs. The predominance of initially north-facing ice cliffs would suggest that their formation is associated with incident solar radiation. Data from the proximal weather station and elsewhere in the Cordillera Vilcanota capture the extreme nature of incoming radiative fluxes experienced at high altitude (>4800 m) in the study area (Fig. 8). During daylight hours (10:00–21:00), mean incident solar radiation measurements are typically within the range 200–500 W/m2. Hourly average maximum incoming shortwave (hereafter SW) values are regularly sustained at 600–700 W/m2, with single hourly measurements close to or even exceeding the top of atmosphere (ToA) solar constant value of 1370 W/m2 at each of the three local weather stations, which may be a result of the amplification of incident SW radiation by forward scattering of broken clouds (Conway and others, Reference Conway2022; Cordero and others, Reference Cordero2023). At the Suyuparina weather station, hourly maximum SW incident radiation values exceeded the ToA constant on 43 separate occasions during May 2024–July 2025. Such high fluxes of SW can readily drive melt at the glacier’s surface (Fyffe and others, Reference Fyffe2021; Conway and others, Reference Conway2022), particularly where impurities or debris reduce the albedo of the ice surface. We suggest that north-facing crevasse walls are exposed to high radiative fluxes over the glaciers in the Suyuparina catchment, which drives their expansion into the initially arcuate cliffs found mid-glacier.
Incident shortwave radiation measured by weather stations at Sibinacocha (4895 m), in the forefield of the Suyuparina and Quisoquipina Glaciers (5153 m) and at the summit of the Quelccaya ice cap (5650 m). Daily SW maxima represent maximum single-incident SW measurements over hourly intervals from each day in the station records. Daytime max/mean SW values represent mean and maximum hourly values, averaged 10:00–21:00 each day. Top of atmosphere (ToA) values represent incident SW radiation at solar noon each day.
As ice cliffs migrate down-glacier, they become more arcuate and bound by steep ice faces which have variable aspects (Fig. 4), rather than exclusively facing north. At the same time, the proportion of cliffs with supraglacial ponds at their base increases (Fig. 3), suggesting supraglacial meltwater begins to modify their geometry as they reach mid-glacier. The combination of deformation through ice flow and enhanced melt around increasingly prevalent ponds down-glacier likely explains the more arcuate form of cliffs at lower elevations. Supraglacial ponds act as a significant recipient of energy to drive melt within the glacier system on debris-covered glaciers, through thermoerosional ablation (Miles and others, Reference Miles2018b), so it is logical to suggest that meltwater ponding at the base of ice cliffs is responsible for their expansion and evolution towards a circular planform down-glacier (Fig. 2d).
Below ∼5300 m a.s.l. on the Suyuparina and Quisoquipina Glaciers, and ∼5450 m on Jurcay Glacier, it is notable how the predominant aspect of the steepest portions of individual cliffs transitions towards a southerly aspect (Fig. 4). This is likely due to the continued irradiance-driven ablation of north-facing portion of cliffs, which become more gently sloping down-glacier due to melt (e.g. Evatt and others, Reference Evatt2024), with the steepest, south-facing features preserved by topographic shading and steepened by undercutting from resident ponds (Fig. 9). The development of a preferential aspect of ice cliffs (Fig. 4) is akin to those extensively documented on debris-covered glaciers, which predominantly face NNW in Northern Hemisphere mountain regions such as High Mountain Asia (HMA) (Kneib and others, Reference Kneib2023), regardless of glacier aspect and their style of formation.
Individual cliff and pond features continue to expand as they are advected towards glacier margins and may coalesce with neighbouring features, although few survive transit through areas of crevassing lower down on each glacier, with surface fracturing allowing for meltwater drainage and reorganisation of the ice surface. Repeat drone imagery of the margin of Suyuparina Glacier during 2022–24 clearly shows how cliffs and ponds that survive to reach the glacier margin either melt through to bedrock and drain subglacially within ∼100 m of the margin (Fig. 9) or drain once the glacier margin retreats to intersect with each cliff. Subglacial routing of meltwater seems to occur for larger ponds with a larger potential hydraulic gradient, whereas smaller ponds drain from the glacier margin.
Evolution of supraglacial cliffs and ponds at the margin of Suyuparina Glacier during 2022–24. Cliffs hosting large supraglacial ponds ∼100 m from the margin of Suyuparina in 2022 (a, b), the same cliffs and ponds in 2023, partially drained and now within ∼50 m of the margin (c, d) and fully drained ponds in 2024 (e, f), within a few metres of the ice margin, captured by drone surveys.
The mechanisms involved in the subglacial drainage of meltwater ponded close to the ice margin may be similar to those associated with the Marmolada glacier collapse, where the filling and expansion of crevasses up-glacier of the ice avalanche release zone has been identified as a contributing factor to the 2022 event, with meltwater routing subglacially in the absence of efficient drainage routes laterally or englacially (Bondesan and Francese, Reference Bondesan and Francese2023). At Marmolada and the three glaciers studied here, ice margin proximal ice thicknesses are a few tens of metres, suggesting polythermal or fully frozen conditions at the bed of each glacier, and supraglacial features (crevasses at Marmolada, ponds in the Suyuparina catchment) have progressively expanded and filled with meltwater over successive ablation seasons. In the absence of an efficient and linked englacial-subglacial drainage network, the overpressure caused by these water-filled surface features may result in the forced routing of water through subglacial till or pre-existing englacial structures such as foliation surfaces. Hydraulic jacking of ice down-glacier has been inferred at Marmolada (Bondesan and Francese, Reference Bondesan and Francese2023), leading to its collapse, but there is little evidence for the occurrence of such collapse events in the Suyuparina catchment, indicating subglacial routing through a deformable till layer as the most likely pathway for draining supraglacial ponds. Finally, the process of cliff and pond growth and advection has allowed for the development of a highly irregular pattern to the margin of Suyuparina Glacier, which is replicated on the Quisoquipina and Jurcay Glaciers. When cliffs and ponds reach the glacier margin, it instantly retreats by tens of metres (Fig. 9e).
5.2. Impact on glacier state
The impact of supraglacial cliff and pond development on glaciers in the Suyuparina catchment seems twofold. As previously mentioned, contemporary (2016–24) rates of surface lowering are heightened within the extent of cliffs and ponds compared to elsewhere over glacier surfaces (Section 4.2). To quantify the impact of the presence of cliffs and ponds on glacier wide ice loss rates, we masked cliff and pond extent during 2016–24 dH grid and filled resulting gaps using mean elevation change estimates within 50 m elevation bands through the full elevation range of each glacier (local hypsometric filling — McNabb and others. [Reference McNabb, Nuth, Kääb and Girod2019]), retaining and propagating the signal of elevation change caused primarily by climate rather than localised ice loss processes. We then calculated the geodetic mass balance of each glacier using both the unaltered and cliff-masked dH grids. For 2016–24, the mass balance of Suyuparina, Quisoquipina and Jurcay Glaciers was −1.20 ± 0.10, −0.85 ± 0.08 and −1.04 ± 0.09 m w.e.a−1, respectively, including cliffs and ponds. For context, the mean mass balance of glaciers across the low latitudes was −0.38 ± 0.06 m w.e.a−1 during 2000–23 (GlaMBIE, Reference Zemp, Jakob and Dussaillant2025), which is comparable to the ice loss rate of mountain glaciers globally (−0.39 m w.e.a−1 during 2000–19 [Berthier and others, Reference Berthier2024; Hugonnet and others, Reference Hugonnet2021]). With cliffs and ponds masked, we estimate the mass balance of the same three glaciers to have been −1.10 ± 0.10, −0.78 ± 0.08 and −0.94 ± 0.09 m w.e.a−1 over the same period; differences of 10.63, 8.97 and 9.09%. The increased net mass loss caused by cliffs and ponds occurred despite these features accounting for relatively minor portions of each glacier’s total area — 4.45, 4.04 and 5.59% on the Suyuparina, Quisoquipina and Jurcay Glaciers in 2024 (Table 3). The proportion of glacier area occupied by cliffs and ponds increased by two to threetimes during 1977–2016/24 on each glacier, so their impact on overall glacier mass balance has likely become more significant with time in the catchment.
A secondary impact of supraglacial cliff and pond development on the glacier’s state is their impact on glacier surface topography. The spatially focused and heightened melt occurring within pond and cliff extents has created a surface of substantial relief, particularly over the ablation zones of Suyuparina, Quisoquipina and Jurcay Glaciers (Figs. 2 and 6). As previously described (Section 4.2, Fig. 6), steep, north-facing ice surfaces and localised high points in the glacier surface topography of Suyuparina Glacier are recipients of intense incident shortwave radiation and ablate much more quickly than adjacent ice surfaces, so cliff and pond development is facilitating the topographically enhanced melt occurring over the lower portion of these glaciers. An ice surface of such relief has not developed on other glaciers in the catchment where supraglacial cliffs and ponds are absent, and so such topographically focused ice loss has not occurred (e.g. Fig. 5b, glacier immediately south of Jurcay Glacier). The hypsometry of these glaciers, which have a greater portion of their total area above ∼5400 m, is more favourable to snowpack retention, which protects the ice beneath from surface melt.
Numerous studies focusing on debris-covered glaciers have highlighted supraglacial ice cliffs and ponds as ablation hotspots and estimate their contribution to total ice loss budgets to be up to ∼40% (Kneib and others, Reference Kneib2023) despite accounting for a fraction (<10%) of glacier surface area. Our estimates of increased net mass loss attributable to cliffs and ponds are clearly lower, but are made at the glacier scale, rather than considering only the ablation zone, the latter of which is common in studies focusing on debris-covered glaciers. Similarly, thick debris found over glacier ablation zones in regions such as HMA dampens melt rates away from ponds and cliffs, enhancing the local contrast between melt over the two features (Miles and others, Reference Miles2018b). As previously stated, thinning rates away from cliffs and ponds on the Suyuparina, Quisoquipina and Jurcay Glaciers were still substantial during 2016–24 (−1.35, −0.95 and −1.13 ma−1, respectively), and so the melt enhancement factor of cliffs and ponds is lower in our observations on clean-ice glaciers in comparison to studies focused on debris-covered glaciers. Regardless, our work suggests that these are important accelerators of mass loss, enhancing net mass loss by ∼10% on vulnerable and important tropical glaciers (c.f. Immerzeel and others, Reference Immerzeel2020). These have previously been neglected in analyses of glacier mass balance.
5.3. The prevalence of clean-ice cliffs and ponds
Ice cliff and pond networks are prevalent across other glaciers in the Cordillera Vilcanota, as well as other glacierised regions of the Andean outer tropics. We note cliff and pond networks on glaciers within the Kaka, Boopi (Bolivia) and Santa (Peru) River catchments, within the latitudinal range of 8–15°S. Iturizzaga (Reference Iturizzaga2014) noted the development of these features above 4850 m a.s.l on the Kogan Glacier and above 4950 m a.s.l on the Pisco-Norte Glacier in the northern part of the Cordillera Blanca (Santa catchment) in Peru, which represents the lowest latitude example of the features we can establish.
Outside of South America, we also note examples of cliff and pond development on low-latitude, high altitude clean-ice glaciers in the northern Hemisphere, in various parts of HMA, namely Bhutan, north of the main orographic divide close to Mt Everest (Tibet) (Fig. 10) and in the Tien Shan (Ak-Shirak; 41.815, 78.252°). Again, these glaciers are found at high elevation (>5000 m), aside from Ak-Shirak, where features are found as low as 4200 m. The layout of features here again suggests incident isolation being a control, with well-developed features being composed of a north-facing backwall, large ponds and an ‘open’ side which likely reflects a melted-out crevasse wall. Matthews (Reference Matthews2020) also demonstrated exceptionally high incident SW radiative fluxes throughout the higher reaches of glaciers in the Everest region, along with short events where the incoming SW flux may exceed ToA levels.
Examples of supraglacial cliffs and pond networks on other clean-ice glaciers outside of the Suyuparina catchment. (a) Northeast of Mt. Everest on the Tibetan side of the main orographic divide. (b) Bhutan, immediately south of Lugge Tsho glacial lake. (c) Bolivia, due east of Ancohuma. (d) Also in the Cordillera Vilcanota in Peru, west of the Suyuparina catchment. Imagery sources: (a) King (Reference King2020b); (b) USGS Earth Explorer declassified spy satellite image archive image ID DZB1209-500106L002001 in collection declass 2, and Google Earth high-res imagery; (c) declass 3 archive D3C1215-300994A003 and Google Earth; and (d) declass 3 archive D3C1213-200305A002 and Google Earth.
As with glaciers in the Suyuparina catchment, there is some indication that the extent of cliff and pond networks may have increased on glaciers outside our study area over the last few decades. Comparison of historical imagery (1974–84) with contemporary images over those glaciers where we note cliff and pond presence elsewhere in the outer tropics and in HMA (Fig. 10) suggests up-glacier expansion of these features, primarily into areas of glaciers displaying a snow-covered surface in historical imagery. Given the sparse nature of these observations, a much more extensive examination of the prevalence of these features is required before their impact on longer-term glacier behaviour can be established.
5.4. Outlook
Our observations of the formation of supraglacial cliffs and ponds on glaciers in the Suyuparina catchment reiterate the importance of net shortwave radiation and precipitation dynamics (Fyffe and others, Reference Fyffe2021) in the melt regime of high-altitude tropical glaciers. The proliferation of supraglacial cliffs and ponds across the surfaces of Suyuparina, Quisoquipina and Jurcay Glaciers has occurred through a period of consistent temperature increases (Supplementary Material, Salzmann and others, Reference Salzmann2013) and snow cover reductions in the Cordillera Vilcanota (Lamantia and others, Reference Lamantia, Larocca, Thompson and Mark2024), without which the expansion of such features would have been unlikely. Whilst supraglacial cliffs and ponds were already present on the Suyuparina, Quisoquipina and Jurcay Glaciers in the 1970s (Fig. 3), they were restricted to the lower portion of each glacier, with a more extensive snowpack likely protecting ice surfaces at higher elevations from incident SW radiation and lower temperatures limiting the amount of meltwater being produced by the snowpack present over glacier surfaces. As temperatures have risen (Supplementary Material, Salzmann and others, Reference Salzmann2013), the regional snowline altitude has increased (Lamantia and others, Reference Lamantia, Larocca, Thompson and Mark2024), increasing meltwater availability over glacier surfaces, and the precipitation phase has begun to transition from snow to rain at lower elevations (Fyffe and others, Reference Fyffe2021). Snowpack loss has also resulted in a greater portion of each glaciers ice surface being exposed to shortwave radiation driven melt more regularly, particularly over north-facing surfaces and ice-surface topographic highs (Fig. 6). Both processes are now amplified by a lower albedo ice surface and the receipt of energy by supraglacial ponds driving additional supraglacial and englacial melt (Fig. 9) and generating an ice surface of greater relief.
Whilst the impact of shrinking snowcover extent can and has been incorporated into glacier modelling studies (Fyffe and others, Reference Fyffe2021), the impact of supraglacial cliff and pond development on glaciers devoid of debris-cover has not been attempted. The energy balance over such features has been examined in detail in studies focused on debris-covered glaciers (Buri and others, Reference Buri2021; Miles and others, Reference Miles2018b; Reference Miles, Steiner, Buri, Immerzeel and Pellicciotti2022) and even quantified at the catchment level (Miles and others, Reference Miles2018b), suggesting an assessment of the wider impact of cliff and ponds development on clean-ice glaciers is achievable. Similarly, the work of Rowan (Reference Rowan2021) also shows that melt enhancement factors associated with ponds and cliffs on debris-covered glaciers can be effectively represented in glacier models, indicating it is also plausible for those without debris cover.
Before any adaptation of model structure is undertaken to incorporate the impact of cliffs and ponds on clean-ice glaciers, further work should be carried out to address the shortcomings of this study to ensure that the processes driving cliff and pond expansion are understood with greater clarity. Firstly, the contribution of supraglacial cliffs and ponds to the pattern and magnitude of ablation can only be quantified with certainty once the vertical ice fluxes over cliff and pond extents have been resolved, which we have not done in this work. This has been achieved recently in studies of debris-covered glaciers (Brun and others, Reference Brun2018; Bhushan and others, Reference Bhushan2024) where suitably high spatial and temporal resolution data are available. These studies show that the amplifying effect of ponds and cliffs on ablation rates persists after corrections for emergence velocities, although failure to account for ice fluxes over faster-flowing glaciers can lead to considerable uncertainties in the volume change associated with individual cliffs. Future work could therefore prioritise understanding the impact of emergence velocities on ablation of clean-ice glaciers hosting cliffs and ponds in high-altitude settings. Similarly, the temporal resolution of the time series of DEMs and orthoimages we have assembled is low in comparison to those being examined in studies examining cliff and pond behaviour on debris-covered glaciers (e.g. Kneib and others, Reference Kneib2023). Our ability to describe shorter-term variability in cliff and pond extent on clean-ice glaciers, or indeed their persistence over longer timeframes (Supplementary Material), is therefore limited, as is the scope for linking their behaviour to external factors varying over similar timescales, such as regional climatic variability.
Considering these points, an obstacle to improved observations of supraglacial cliff and ponds behaviour on high-altitude, clean-ice glaciers is the need for data of fine spatial and temporal resolution to constrain their development. The largest examples of these features over Suyuparina, Quisoquipina and Jurcay Glaciers reached ∼125 m in length, with a surface relief of up to 40 m. Whilst satellites such as Pléiades and the WorldView sensors are now capable of resolving these features, moderate resolution (>10 m ground sampling distance) optical datasets captured by Sentinel 2A, Landsat 8/9 would only capture the largest features in limited detail. In addition, the reduced footprint of higher resolution sensors such as Pléiades (swath width 20 km) precludes their use in regional scale studies, with most studies utilising these data focusing on sites typically 100–500 km2 (Berthier and others, Reference Berthier2024). Similar limitations apply to the declassified image archive, with the Hexagon KH-9 mapping camera more suited to capturing larger-scale glacier processes or collections of smaller-scale features distributed across glacier surfaces, rather than the evolution of individual features at a more localised scale (e.g. Fig. 9b). As we have shown, Hexagon KH-9 panoramic-camera imagery is more suited to studying fine-scale glacier processes, but success is highly dependent on the situation of an area of interest within the bounds of acquired images, with higher image resolution (and accuracy) found towards the centre of individual scenes. Continuation and expansion of monitoring programmes such as the Pléiades Glacier Observatory (Berthier and others, Reference Berthier2024) likely offer the best opportunity to capture the development of these localised glacier surface features in the future.
Overall, our initial observations suggest that, with rising temperatures and decreased snow cover, supraglacial cliffs and ponds are becoming more prevalent (Fig. 10). These features appear to be hotspots of ice loss, enhancing overall glacier mass loss rates by ∼10%, and accelerating glacier frontal retreat. We tentatively suggest that with continued global temperature increases and less predictable snowfall, these features may become more prevalent across high-altitude and low-latitude mountain regions. As a consequence, in the near future, their role in determining the mass balance of glaciers may become more significant. However, over longer time scales, and at an individual glacier scale, this process may also have a limited timespan. As glaciers retreat, the geometry of ablation zones may become less favourable to their development as glacier thicknesses and flow rates decrease (Dehecq and others, Reference Dehecq2018). Thinner, less dynamic glaciers will be less prone to significant crevassing, which we suggest is the point of formation for supraglacial ice cliffs and ponds on debris-free glaciers, particularly where glaciers flow over bedrock obstacles. Similarly, should a glacier retreat up-valley of former crevasse formation zones, supraglacial cliffs and ponds may cease to develop. Further work refining the conditions of their formation and their longevity is required.
5.5. Conclusions
We have examined the prevalence and influence of supraglacial cliffs and ponds over three high-altitude, clean-ice tropical glaciers in the Suyuparina catchment of the Cordillera Vilcanota. Using high-resolution satellite and drone-derived DEMs and orthoimages, we have combined observations of glacier surface structures, glacier surface geometry and glacier surface elevation change patterns to propose a model of formation of supraglacial cliffs and ponds, which we suggest is primarily a consequence of aspect-dependent surface melt driven by high shortwave radiative fluxes. We identified approximately twice as many supraglacial ice cliffs over the period 2016–24 when compared with 1977, along with a significant increase in the number and area of supraglacial ponds on the same glaciers. We derived surface elevation change data from various combinations of DEMs to quantify rates of thinning and total ice mass loss at the glacier scale. Our results show substantial glacier mass loss in the region since the 1970s (mean rate −0.49 ± 0.04 m w.e.a−1 during 1977–2024), a significant acceleration in the rate of ice mass loss recently (mean rate −1.03 ± 0.07 m w.e.a−1 during 2016–24) and the operation of localised, amplified ice loss processes, coincident with the location of supraglacial cliffs and ponds. Ice melt associated with cliff and pond development accounts for ∼10% of the total ice mass loss budget, despite these features occupying a minor portion (4–5%) of total glacier area. The formation and advection of these features towards glacier margins has resulted in the development of glacier surface topography of substantial relief, which itself is subject to locally amplified melt over north-facing and protruding ice surfaces.
Whilst supraglacial cliffs and ponds have been identified across many debris-covered glaciers and studied extensively in these environments, their occurrence on clean-ice glaciers has received little, if any, concerted attention. Our observations show how the development of these features on clean-ice glaciers can amplify their rates of recession and ice loss. More focused work should be conducted on the prevalence of these features, which, as we show, are not confined to the Cordillera Vilcanota, and their incorporation into glacier models, which will be needed to reliably constrain future freshwater yield from high mountain environments at low latitudes. Such work will be challenged by the requirement for high-resolution observational data and modelling approaches, with moderate resolution satellite imagery (e.g. Sentinel 2a or Landsat) and common model structures (e.g. temperature-index models or energy balance models based on global DEMs) unlikely to resolve these small but abundant ice surface features.
Supplementary material
The supplementary material for this article can be found at 10.1017/jog.2025.10109.
Data availability statement
Datasets documenting changes in glacier extent and surface elevation associated with this study are available on a public repository (10.5281/zenodo.17699108). xDEM is available at https://github.com/GlacioHack/xdem (last accessed: 31 January 2025; DOI: 10.5281/zenodo.8220229).
Acknowledgements
This work was undertaken as part of the Natural Environment Research Council (NERC) Highlights Project ‘Deplete and Retreat: the Future of Andean Water Towers’ (NE/X004031/1). O.King acknowledges an ESA data grant (PP0097129) and the Newcastle University Faculty Research Fund. Tom Gribbin was supported by the CENTA doctoral training programme (NERC grant no. NE/S007350/1). The Quelccaya Ice Cap and Lake Sibinacocha stations were supported by grant AGS-1347179 (CAREER: Multiscale Investigations of Tropical Andean Precipitation) from the U.S. National Science Foundation to L.B. Perry in collaboration with the Universidad Nacional de San Antonio de Abád del Cusco (UNSAAC) and the Peruvian Servicio Nacional de Meteorología e Hidrología (SENAMHI). Maxwell Rado and Sandro Arias coordinated the maintenance and operation of the stations. We thank Sutapa Bhattacharjee for her work towards deriving ERA5 Land temperature and precipitation anomalies.
Open access
