Our analysis spans the periods 1970–2013 and 2013–2019, corresponding to two spatial scales (regional/local). Satellite data used for the various time steps are summarized in Table 1. For the 1970s decade, we obtained declassified imagery from the Corona KH-4B mission (Dashora et al., 2007). Corona KH-4B satellites acquired data from September 1967 and May 1972 using a pair of counter-rotating panoramic cameras with a 30° separation angle (one camera tilted 15° forward and one 15° backwards), yielding stereo imagery (Dashora et al., 2007). We obtained three stripes of stereo images from the mission 1,112 acquired on Nov 19, 1970. This date corresponds to the end of the ablation season in this part of the Himalaya (Oct/Nov). Image stripes were scanned at 7 microns directly from the original film stripes and USGS; the digital scans were obtained at no cost from the US Geological Survey (https://earthexplorer.usgs.gov). The scenes were cloud-free over the glaciers.

For the recent time period, we obtained high-resolution images stripes at 5 m from Planet’s RapidEye constellation (Planet_Labs, 2021). RapidEye consists of multispectral data (five spectral bands in the visible and near infrared) with a spatial resolution of 6.5 m on the ground and a positional accuracy of <10 m (Table 1) (Planet Labs, 2016). We used Level 3A data, which consist in stripes of multispectral, radiometrically corrected orthotiles processed at 5 m spatial resolution. These represent surface reflectance computed using at-sensor radiance based on the 6S radiative transfer model (Vermote et al., 1997) and MODIS data, which accounts for atmospheric effects (Planet Labs, 2021). We mosaicked these using nearest neighbour to output a single image covering our study area (Figure 1A). The Nov 20, 2019 image was cloud-free and was acquired at the same time of the year as the Corona data (Nov 19, 1970) and a few days after our field campaign (Nov 12–16th, 2019). The 2019 image had some seasonal snow outside the glaciers in some areas, so we used a second scene from Oct 20, 2015 to refine the glacier extents around rock outcrops and debris cover, and to eliminate seasonal snow.

The High Mountain Asia (HMA) DEM at 8 m spatial resolution (Shean, 2017) was chosen as reference to represent the modern glacier surface for the 2010s decade. The HMA DEM tiles are generated from very high-resolution satellite imagery from Digital Globe Inc., from 2002 to 2016, processed using the NASA AMES Stereo Pipeline as described in Shean et al. (2016). We obtained six DEM tiles derived from along-track Worldview two imagery from the NSIDC DAAC High Mountain Asia (HMA) collection (https://nsidc.org/data/highmountainasia) with dates ranging from November 16, 2013 to November 25, 2013. Although DEM tiles were available from subsequent years, the coverage of the study area was incomplete and in particular large voids were present over some of the glacier accumulation areas. The geolocation accuracy of the HMA DEM is estimated at < 5 m, with relative vertical errors of 1–2 m. Larger uncertainties occur due to clouds, steep slopes image saturation and/or correlation failures (Shean et al., 2016). We also obtained the ALOS Global Digital Surface Model (AW3D30) version 2.2 (30 m) (JAXA, 2019), constructed from data acquired from 2006 to 2011. The vertical accuracy of this dataset is estimated at ∼10 m in Eastern Nepal (Tadono et al., 2014). The AW3D30 DEM was used as regional elevation reference for ground control points in the Corona DEM extraction process (see section Corona DEM Generation and Co-Registration).

In November 2019, we conducted a UAV survey of the lower ablation area of Ponkar Glacier (∼4 km², 3,646–4,238 m a.s.l.) (Figure 1B) using two DJI Phantom four Pro UAVs, with an 80% forward and 70% side overlap, covering an area of 3.36 km². A total of 1,110 images were orthorectified using nine ground control points (GCPs) and were processed using a dense stereo matching algorithm in Pix4D Mapper software to produce a 0.50 m orthoimage and 1 m DEM. The accuracy of the orthorectification process was estimated as 1 m horizontally and <1.5 m vertically (2 and 3 times the ground sampling distance, respectively). Both the orthoimages and the resulting DEM were used to refine the geomorphological mapping (see section Mapping of Supraglacial and Geomorphological Features).

The lower ablation area of Ponkar Glacier has been surveyed in the field for mass balance and hydro-geochemical studies since 2016 (Shrestha et al., 2021). Ice ablation measurements on the lower area of Ponkar Glacier have been carried out yearly from 2016 to 2019. A Kovacs hand driller was used to drill the debris-covered ice up to 2.5 m to install the ablation stake at four locations under different debris thicknesses ranging from 11 to 20 cm. Four stakes (Stakes 1–4) were installed in March 2016 and resurveyed in different dates during field trips from 2016 to 2019. Their locations are shown on Figure 1B and their are characteristics listed in Table 2.

Glacier outlines for the year 1970, including debris cover, were delineated using on-screen manual digitization complemented with image thresholding from the Corona panchromatic image. Shaded relief terrain derived from the 1970 Corona DEM was used to identify rock outcrops. The low contrast of the black and white Corona image made it difficult to distinguish ice from snow in a few areas. Therefore, Corona outlines were compared with the ICIMOD inventory (Mool et al., 2002), which was constructed by manual mapping on topographic maps at 1:50,000 scale published in the 1950s–1970s by the Survey of India (Survey of India, 2005). The latest glacier inventory in this area is the GAMDAM glacier inventory (Sakai, 2019), constructed based on semi-automated mapping from Landsat imagery and manual editing. However, most of the GAMDAM glacier outlines in our study area dated from 1994, except for those in the northern part (China) which dated from 2010, and were thus not suitable for analysis. For the year 2019, glacier outlines were derived by manual digitization on the RapidEye image; areas of seasonal snow were removed manually. When glaciers disconnected or disintegrated, we counted them as special glacier entities as per guidelines established by the Global Land Ice Measurements from Space community (Raup and Khalsa, 2007; Racoviteanu et al., 2009). Debris-covered glacier sections were cross-checked visually against high-resolution SPOT imagery from Oct 24, 2019 visible on the Google Earth platform. For both the 1970 and 2019 time steps, ice divides were constructed on the basis of the AW3D30 DEM complemented with the HMA DEM using the buffer method (Granshaw and Fountain, 2006; Bolch et al., 2010). These were kept constant for the two time steps to facilitate the glacier change detection analysis and to minimize uncertainties related to differences in the upper part of the accumulation area. Rock outcrops were delineated for 1970 and 2019 separately.

We calculated glacier surface slope, mean altitude, and altitudinal range for all the glaciers in the study area based on the AW3D30 DEM (for the regional trends) and the HMA DEM (for the in-depth study on Ponkar Glacier). We also calculated glacier surface velocities in the Manaslu region using a combination of Sentinel-2 and Landsat-8 imagery, following Millan et al. (2019). The ice velocity maps were calculated using a feature tracking algorithm on all possible image pairs, with repeat cycles ranging from the sensor nominal repeat cycle (i.e., 5 days for Sentinel-2, 16 days for Landsat-8) up to more than a year (Millan et al., 2019). Ice velocity mosaics are finally assembled at a sampling resolution of 50 m, by computing a pixel-by-pixel weighted average (Millan et al., 2019). We also obtained yearly velocity data from 2000 to 2017, generated using auto-RIFT (Gardner et al., 2018) and provided by the NASA MEaSUREs ITS_LIVE project (Gardner et al., 2019).

The Corona image stripes were set up as a single block file in ERDAS IMAGINE Photogrammetry Suite (ERDAS, 2010) using the UTM projection system (zone 45N), with WGS84 vertical and horizontal datums. Since Corona parameters are not available, we used a non-metric camera model based on the actual flying height, focal length and the nominal scale of the photos (Table 1). The interior and exterior orientation of the stripes were solved in the aerial triangulation step using Ground Control Points (GCPs) and Tie Points (TPs). For the northern part of the study area (four stereo image stripes), TPs were generated automatically in ERDAS IMAGINE Photogrammetry and were complemented with manually generated points. This resulted in a total of 3,047 TPs with a mean uncertainty of 6.6 m and a standard deviation of 9.1 m. A total of 255 GCPs were manually digitized on stable terrain from the RapidEye scene and identified on each of the stereo image stripes. Elevations of the ground control points were extracted from the AW3D30 DEM. The total root mean square error for the triangulation (RMSE) was 1.9 pixels, with ground x, y and z accuracy of 14.9, 11.2 and 19.3 m, respectively. Point clouds were generated using the Dense Point Matching with the Enhanced Automatic Terrain Extraction (eATE) in ERDAS Imagine Photogrammetry module using the Normalized Cross-Correlation routine with a 9 × 9 window size, spike interpolation setting to remove extreme values and a point threshold of 2.5 (number of points for interpolation). No filtering was applied to the point clouds. For the southern part of the study area where Thulagi Glacier is located, two 1970 image stripes had been processed similarly in a previous study (Robson et al., 2018) using 75 GCPs and 1837 tie points with a RMSEz <5 m. To extract the DEMs, we utilised both the eATE routine in ERDAS and PCI Geomatica. Each DEM was extracted as a point cloud which was then gridded using a finite difference interpolation method with 10 iterations within Geomatica Banff (PIC Geomatics Inc., 2020).

Two DEMs were produced: one for the northern stereo pairs, and one for the southern stereo pairs. Each DEM was co-registered independently before being mosaiced together. The DEM co-registration was performed in a two-step process using the Numpy and ArcPy python packages. Linear co-registration biases were removed using the method outlined by Nuth and Kääb (2011) which minimizes the root mean square residuals of the elevation bias over stable (i.e non-glacial) terrain. Terrain with slope values less than 10^o and greater than 40^o was excluded. The co-registration was repeated iteratively until the mean elevation bias over stable terrain decreased by <2%. In both cases, this led to the process being run 8 times. Corona imagery contains significant panoramic lens distortions which can result in non-linear elevation biases of up to 20 m (Bolch et al., 2011; Robson et al., 2018). These elevation biases were removed by fitting a sixth order polynomial for along track and third order for across track and elevation to the stable terrain elevation biases (Figure 2). A summary of the co-registration shifts is given in Table 3. The process was iterated a total of six times until no noticeable trend in stable ground polynomials was visible. The co-registered surface elevation rasters were then mosaicked together to one file. The same co-registration routine was applied to co-register the 2019 UAV data to the 2013 DEM with the linear-co-registration being repeated four times and the non-linear co-registration repeated three times.

Surface elevation changes and the resulting geodetic mass balance were computed at two time steps for the two spatial domains using the geodetic (DEM differencing) method (Berthier et al., 2004). For the large spatial domain (Manaslu region), we used the 1970 Corona DEM and the 2013 HMA DEM. For the smaller domain (the lower area of Ponkar Glacier), we used the 2013 HMA DEM and the 2019 UAV data. DEMs were subtracted on a cell-by-cell basis for each of the time periods. For both periods, the elevation difference rasters resulting from the co-registration and mosaicking were cleaned using the method set out by Gardelle et al. (2013) where the mean surface elevation change per 50 m altitudinal band was calculated. Any pixels that were not within three standard deviations of the mean value were excluded. Voids were filled by fitting third-order polynomials to the mean elevation change per 50 m altitudinal band, after McNabb et al. (2019). A uniform glacier ice density of 850 ± 60 kg m⁻³ (cf. Huss, 2013) was used to convert volume changes into mass changes for the geodetic mass balance calculation.

We mapped the glacier surface features and associated geomorphology in 1970 and 2019 at the local scale (Ponkar Glacier). Ice-surface features mapped included glacier flow units, crevasses and ogives; mapped landforms included terminal and lateral moraines and large areas of outwash which complement the glacier extents. These were mapped manually using on-screen digitization on the Corona and RapidEye images on the basis of standard criteria for identifying and mapping ice-structural and glacial landforms (adapted from Glasser et al., 2005) and Goodsell et al. (2005)) (Table 4). Surface features (ice cliffs, supraglacial ponds and ice-surface vegetation) were mapped from Corona and RapidEye scenes using an image segmentation workflow implemented in the ENVI Feature Extraction Module (Harris Geospatial, 2017). This consisted in applying an edge algorithm to identify the pond, ice cliffs and vegetation segments using a texture kernel size of 3 pixels. Segments were merged manually to prevent over-segmenting and to combine different segments into one pond, ice cliff or vegetation patch.

The uncertainty of the elevation change was calculated based on the method set out by Falaschi et al. (2019). The uncertainty of the volume change ( E Δ v i ) (in m³) was determined by summing up the standard error ( E Δ h i ) per 50 m altitudinal band, multiplied by the area of each altitudinal band ( A i ) in order to account for the hypsometry: E Δ v i =   ∑ i n E Δ h i   ×   A i (1) Where the standard error ( E Δ h i ) is derived from the standard deviation over stable ground ( σ s t a b l e ), divided by the effective number of observations (N): E Δ h i =   σ s t a b l e N (2)

N is calculated using the number of pixels (N_tot) in the DEM differencing, the pixel size (PS), and the distance of spatial autocorrection, which following Bolch et al. (2011) we took to be equal to 20 pixels: N =   N t o t × P S 2 d (3)

In addition to E Δ v i we also considered the uncertainties relating to the volume to mass conversion ( E ρ ), which following Huss (2013) was taken to be ±60 kg m⁻³ and the error in delineating the glacier outline ( E a ) which was assumed to be 5%. E Δ v i was combined with E ρ and E a to derived the total uncertainty, E Δ v t o t : E Δ v t o t = E 2 Δ v i + E ρ 2 + Ε a 2 + Ε r 2 (4)

Surface elevation change raster data tends to contain data voids over steep or shadowed terrain, as well in areas of low image contrast. The polynomial fit used to interpolate each glacier was inspected, and any glacier where the fit was poor was not included in the analysis of surface change. In the absence of comparable in-situ data, we estimated DEM uncertainties in the based on the literature combined with biases measured over the stable terrain using a conservative approach. Magnusson et al. (2016) demonstrated that such techniques typically overestimate uncertainties as they neglect the spatial dependence of DEM errors.

Uncertainties in the 1970 and 2019 glacier areas and geomorphology units were determined using the buffer method (Granshaw and Fountain, 2006; Pekel et al., 2016; Paul et al., 2017). This consists of applying a buffer inside and outside the glacier polygons and calculating the standard deviation as an uncertainty estimate. Following the method of Granshaw and Fountain (2006), we used ∼½ of the RMSE of the geolocation errors of both images for the size of the buffer (5 m buffer for RapidEye and 6.5 m buffer for Corona). Geolocation errors were calculated as the RMSE of each satellite image orthorectification. The uncertainty of the area change analysis was estimated as the RMSE of the two error uncertainties (1970 and 2019). The error of the supraglacial vegetation was estimated using a ½ pixel buffer (1 and 2.5 m for Corona and RapidEye images, respectively). For ice cliffs and ponds mapping, we used uncertainty values from Steiner et al. (2019) study also based on manual digitization of high resolution imagery. The uncertainty in the glacier terminus recession was estimated using two RapidEye 2 pixels and 5 Corona pixels.

Remote sensing estimates of glacier area, area changes, and geomorphological features are subject to the quality of the images as well as the accuracy of the glacier delineation procedure (Paul et al., 2013). Using the buffer method, we obtained a glacier area uncertainty of 2.7% for the RapidEye imagery and 2.1% for the Corona imagery. While Corona imagery is panchromatic and is likely to have more uncertainty, the total error is slightly less than RapidEye because the buffer method takes into account the higher spatial resolution of Corona compared to RapidEye. The uncertainty of the area change analysis, estimated as the RMSE uncertainties of the two dates was 3.4%. The uncertainty in the glacier terminus recession was estimated at 10 m (i.e., two RapidEye pixels and 5 Corona pixels). The error of the supraglacial vegetation amounted to 9.7 and 8.2% for 1970 and 2019, respectively. Supraglacial pond and ice cliffs were difficult to identify on the Corona image as they sometimes mixed with shadow and lithology. We used uncertainty values of 26 and 21% respectively (cf. Steiner et al., 2019). The feature detection routine helped mitigate some of this challenge and helped automate the delineation of supraglacial ponds and ice cliffs, but we remain cautious about the area estimates for these two features.

With regards to surface elevation changes, all of our regional mass balance estimates, along with 89% of the individual glacier mass balance estimates, and 94% of the surface lowering estimates were significantly different from zero. The non-linear co-registration biases (Figure 2) show that there were deviations from the polynomial of up to 5 m. These biases are also visible on Figure 4, and predominantly occur on steep slopes as well as areas of low-image contrast, such as shadows. In most cases, the surfaces of the glaciers in the study area had good image contrast, especially over the debris-covered surfaces. As such, our error estimates of mass balance are conservative. At the glacier scale, the removal of the non-linear co-registration biases on the UAV data was complicated by the lack of stable terrain in the dataset. The only terrain that could be used to identify the biases was on the steeply sloped eastern side of the terminus of Ponkar Glacier. As such we cannot be confident that the polynomials used to remove the biases are equally applicable across the entire study area; this most likely explains the 0.2–0.5 ma⁻¹ elevation bias that can be observed on the lateral moraine of Ponkar Glacier in Figure 4. Nevertheless, the changes visible on the glacier surface are greater than the biases remaining on the stable ground so we are confident that the elevation change patterns observed are reliable.

In this study, we find that glaciers in the Manaslu region of Nepal have been receding and losing mass over the last five decades. The region-wide mean rate of area loss (−0.23% a⁻¹, 1970–2013) is slightly lower than those reported for the Everest area over the last four decades (Thakuri et al., 2014: −0.27 ± 0.06% a⁻¹, 1962–2011; King et al., 2020: −0.39 ± 0.13 m a⁻¹, 1962–2018). The area changes over the entire period 1970 to 2019 is much lower than the area loss over recent periods reported in the same area in a previous study (Robson et al., 2018: −1.2% a⁻¹, 2001–2005; −0.4 a⁻¹, 2005–2013). This suggests an accelerated area loss in the last decade, also noted in other areas of the Himalaya (Ren et al., 2006; Thakuri et al., 2014; King et al., 2020). However, we interpret the difference in rates of area change between our study and Robson et al. (2018) with caution, because of the different glacier mapping methods used. Regional area change estimates going back to the 1960s and 1970s remain limited in the Himalaya, making it difficult to directly compare the behaviour of glaciers in the Manaslu region to other regions. Here we compare our findings with existing remote sensing studies for similar timescales.

For the period 1970 to 2013, our estimated mean geodetic mass balance (−0.17 ± 0.03 m w. e.a⁻¹) for the Manaslu region is lower than the larger central-west Himalaya domain where are study area is located (King et al., 2019: an average of ∼ −0.25 ± 0.11 w. e.a⁻¹ from ∼1974–2015), as well as than Himalaya-wide estimates (Maurer et al., 2019: −0.31 ± 0.13 m w. e. a⁻¹, 1975–2016). Both cited studies used similar methodology to ours (declassified Hexagon KH-9 imagery combined with SRTM and ASTER DEMs, respectively). However, the slight differences in the size of the domains, the time span of ±3 years and the type of imagery used may account for some of the differences between our study and the cited studies. With respect to temporal evolution of the mass balance, the rate of mass loss over the period 1970–2013 in the Manaslu area is lower than rates reported for the last decade (2000–2015/2016) in the central Himalaya (Brun et al., 2017: −0.28 ± 0.08 m w. e. a⁻¹; King et al., 2019: −0.26 ± 0.11 m w. e a⁻¹). This suggests an acceleration of mass loss in this area in the last decade compared to the previous decades, which is consistent with findings from other studies (Maurer et al., 2019; Hugonnet et al., 2021).

In addition to mass loss, our study shows that glaciers are slowing down as they experience negative mass balance, a general pattern noted in other studies across the Himalaya. Dehecq et al. (2019) reported a slowdown of glaciers in High Mountain Asia for the period 2000–2016. For instance, in the Manaslu region, glaciers experience moderate thinning along with moderate glacier velocities compared to other regions such as Nyainqntanglha and Lahaul Spiti regions (Dehecq et al., 2019). Coupled with the low supraglacial pond and ice cliff coverage in this region (Racoviteanu et al., 2021), this indicates that glaciers in this area might be in the initial phase of retreat, characterized by melting and slow downwasting associated with slowly retreating glacier termini (Kirkbride, 1993). The presence of only three lake-terminating glaciers out of 213/274 glaciers in 1970/2019 indicates that the third phase in the evolution of a glacier, marked by the development of deepening proglacial lakes (Kirkbride, 1993; Benn et al., 2012) is not under way in this area. This is different than the eastern Himalaya, for example Sikkim, which is characterized by accelerated growth of supraglacial and proglacial lakes (Basnett et al., 2013; Racoviteanu et al., 2015; Shukla et al., 2018).

While the area and elevation change in the Manaslu region are moderate compared to other regions, and rather homogenous (clean vs. debris-covered glaciers), differences exist on a glacier-by-glacier basis. The influence of supraglacial debris on glacier area and surface changes in the Himalayas has been addressed in recent studies (Brun et al., 2017; Dehecq et al., 2019; Shean et al., 2020). In this discussion, we focus on the contrasted behaviour of the lake-terminating Thulagi Glacier and the land-terminating Ponkar Glacier and we infer links between surface morphology (supraglacial pond and ice cliff coverage, slope, and debris cover), surface velocity and mass balance of these two glaciers.

These glaciers are both covered with debris but behave differently with respect to area and elevation losses. The lake-terminating Thulagi Glacier had higher area loss, terminus retreat and thinning rates compared to the non-lake terminating Ponkar Glacier from 1970 to 2013/2019. The proglacial Thulagi Lake, previously identified as a hazardous glacial lake, has doubled in size since the 1970s based on our estimates. The growth of Thulagi Lake occurred concomitantly with the thinning of Thulagi Glacier noted in this study as well as in previous studies (Pant and Reynolds, 2000; Robson et al., 2018; Maskey et al., 2020; Watson et al., 2020). Ponkar Glacier, in contrast, did not develop a proglacial lake at its terminus. In terms of surface elevation changes, Ponkar Glacier stands out as an exception to the regional pattern noted earlier (i.e. thinning towards the terminus with stable or thickening trends in the upper areas, vs. consistent thinning across the debris-covered part for most of the glaciers in the region). The thickened upper ablation zone shown in Figures 3, 4 is probably best explained by the velocity data in Figures 7, 8, which indicates that this part of the glacier is still dynamically active. The simplest explanation is therefore that there is still some ice actively flowing from the accumulation area into the upper ablation area but that the glacier is then more or less dynamically inactive below this (i.e. in the middle and lower ablation zone). We infer that the recent prolonged period of negative mass balance led to reduced ice velocities and enhanced ice-surface lowering, especially near the glacier terminus.

Using the example of Ponkar Glacier, we note that changes in the ice-surface morphology based on satellite imagery are commensurate with the recession and thinning of the glacier between 1970 and 2013/2019. For example, we note an increase in the area of ice-surface vegetation over the 1970–2019 period. This supports the idea that Ponkar Glacier is currently in a stagnation phase, in comparison with Thulagi Glacier which is at the stage of lake development and calving. Furthermore, the ice-flow contribution from Flow Unit 1 to the main trunk (Flow Unit 2) decreased by 2019. Ogives, which form beneath icefalls on glacier surfaces indicate active glacier flow; crevasses generally open and form beneath actively-flowing ice. The decrease in number and extent of both of these ice-surface features, coupled with an increase in the size of lateral moraines and expanded outwash plains further indicate that Ponkar Glacier is currently less dynamically active than in 1970. Crevasses in Flow Unit 1 are only visible on that section of the glacier in the 1970 image and some crevasses appear to have closed up by 2019. While this can be a function of glacier recession (i.e., glacier recession means there is more available space for outwash plains to expand on the valley floor), the expansion of outwash plains could also indicate an increase in ablation and associated fluvial activity between 1970 and 2019.

The large lateral moraines present along the glacier margins in 1970 indicate that glacier thinning was already underway five decades ago. Moraines expanded in size between 1970 and 2019 as the glacier continued to thin and recede during this period. The most conspicuous change observed is near the snout in the “embayment” area, on the true left-hand side near the terminus, where moraines can be seen actively forming at the 1970 ice margin. By 2019, these moraines separated from the glacier itself and a new lateral lake had formed within the moraines in the embayment area (Figure 9B).

On Ponkar Glacier, supraglacial ponds covered a small percent of the debris-covered glacier area in both 1970 and 2019 (0.25% and 0.45% of the debris-covered area, respectively). Ice cliffs covered 1–1.35% of the debris-covered area. The low supraglacial pond coverage is a pattern over the entire Manaslu area as noted in a previous study (Racoviteanu et al., 2021). With respect to regional ice cliff coverage, we cannot compare directly the local estimates with the regional ones due to the difference in the spatial resolution of the data used in the cited study vs. this study (30 vs 5 m). However, we note that the pond and ice cliff coverage that we obtained for the post-monsoon (dry) season in 1970 and 2019 is about half that of the area reported from other glaciers in the Nepal Himalaya for the same season. For example, Steiner et al. (2019) reported 2.5% ice cliff coverage and 0.9% pond coverage relative to the debris-cover area of five glaciers in the Langtang valley, about 115 km east of our study area in October 2015. Brun et al. (2018) found a 7–8% ice cliff coverage on the Changri Nup Glacier in the Khumbu Himalaya during the dry season, east of our study area, based on Pléiades imagery (2 m). This is ∼7 times larger than the ice cliff coverage we estimated on Ponkar Glacier for the dry season (November 1970 and 2019). The increased number of ice cliffs from 1970 to 2019 on Ponkar Glacier appears to be consistent with the thinning patterns we noted earlier (see section Changes at Glacier Scale: Ponkar and Thulagi Glaciers). The doubling of the number and surface area of supraglacial ponds on Ponkar Glacier from 1970 to 2019 also supports the idea that the glacier was less dynamically active in 2019 than it was in 1970. The slight increase in the number of ice cliffs on this glacier suggest that ice cliffs are in the “development” phase as defined by Sakai et al. (2002). Ponkar Glacier represents a good example of the relationship between glacier mass balance, glacier dynamics and changes in supraglacial debris cover for a glacier which has not reached the transition to the calving stage. In terms of local ablation, we found that surface elevation changes between 2013 and 2019 were noticeably higher within the vicinity of supraglacial ponds and exposed ice (up to −4.5 ma⁻¹) (Figure 4B), with implications for the future behaviour of this glacier. Other studies in the Himalaya clearly showed that while these surface features cover only several percent of the total ablation area of glaciers (Steiner et al., 2019), they can enhance local ablation rates by up to three times (Brun et al., 2016; Irvine-Fynn et al., 2017; Miles et al., 2018b; Buri et al., 2021).

In this study, we have shown that glaciers in the Manaslu region have undergone thinning and frontal recession in recent decades in response to regional climate change. During this period, the glacier ablation zones have accumulated debris at an increased rate. Debris naturally accumulates in a supraglacial position during surface lowering of glaciers in mountain areas such as the Himalaya for two main reasons. First, material falls onto the ice surface from the collapsing inner faces of the surrounding lateral moraines (Hambrey et al., 2008). These moraines are progressively de-stabilised during ice-surface lowering as the supporting ice is removed. This is especially the case near the terminus of glaciers, where ice velocity is low or where the ice is stagnant and therefore supraglacial material cannot be evacuated. Second, there is a tendency for an increase in supraglacial debris cover simply because debris becomes concentrated on the ice surface during ablation (i.e., as the proportion of ice in an ice-rock mixture decreases during ablation, the proportion of debris increases). The ∼12% increase in supraglacial debris cover from 1970 to 2019 noted in this study is consistent with the increase in debris cover reported in other studies in the Himalaya and the Alps (Mölg et al., 2019; Xie et al., 2020). Again, this is especially the case near glacier termini where velocities tend to be low during periods of negative mass balance (Benn and Evans, 1998).

The regional pattern of surface lowering leads to an increase in debris cover for the two main reasons noted above, which in turn reduces ablation. This accounts for the large and increasing number of debris-covered glaciers in the Himalaya (Scherler et al., 2018; Herreid and Pellicciotti, 2020). This poses the question: how will this, and other glaciers in the region, evolve as climate warming continues? There are perhaps two obvious answers to this question: first, all the glaciers could melt and lose mass through the rest of this century and beyond, leaving behind small niche glaciers at high altitudes and plateau icefields on broad mountain tops; second, debris-covered glaciers such as Ponkar Glacier might undergo a transition to a rock glacier as debris cover accumulates and glacier snouts reorganize and undergo increased downslope creep (Harrison et al., 2021). This transition process has been hypothesized to occur in the arid Andes (e.g. Monnier and Kinnard, 2015), in the Himalayas (Jones et al., 2019) and elsewhere globally (Knight et al., 2019).

If Ponkar Glacier does undergo a future transition to a rock glacier then we would envisage this to occur in two stages. In Stage 1, climate change leads to increased air temperature and associated reduction in solid precipitation, which influences ice flow and drives glacier recession. While few long-term climate datasets exist in the high mountain region of western Nepal, there are data from valley locations in the region that show a recent rise in temperature in the area allied to a reduction in precipitation (Basnet et al., 2020; Paudel, 2020). If these trends exist at high altitudes, we would expect these to induce negative mass balance of the glacier. As a result, the glacier thins and accumulates surface debris via rockfalls from valley sides, mass-wasting of lateral moraines and transport of debris from subglacial and englacial sources (Benn and Owen, 2002; Knight and Harrison, 2009; Kirkbride and Deline, 2013; van Woerkom et al., 2019). In Stage 2, accumulation of surface debris occurs to the extent that mass balance of the terminus becomes more positive as ice melt is reduced. In conjunction, basal shear stress increases as does internal deformation of the terminus ice mass as debris load increases. This induces flow and deformation of the terminus which transitions to an ice-cored rock glacier. However, the processes outlined in Stage 2 have clearly not occurred on Ponkar Glacier, where the terminus area shows low velocities rather than the increased velocity hypothesised in this transition model (see Figure 8B), despite the clear increase in supraglacial debris in this zone.