Scholarly article on topic 'Quantifying ice cliff evolution with multi-temporal point clouds on the debris-covered Khumbu Glacier, Nepal'

Quantifying ice cliff evolution with multi-temporal point clouds on the debris-covered Khumbu Glacier, Nepal Academic research paper on "Earth and related environmental sciences"

Share paper
Academic journal
Journal of Glaciology

Academic research paper on topic "Quantifying ice cliff evolution with multi-temporal point clouds on the debris-covered Khumbu Glacier, Nepal"

Quantifying ice cliff evolution with multi-temporal point clouds on the debris-covered Khumbu Glacier, Nepal


1 School of Geography and water@leeds, University of Leeds, Leeds, LS2 9JT, UK 2Department of Geography, University of Sheffield, Sheffield, S10 2TN, UK 3Lancaster Environment Centre, Lancaster University, Lancaster, LA1 4YQ, UK Correspondence: C. Scott Watson <>

ABSTRACT. Measurements of glacier ice cliff evolution are sparse, but where they do exist, they indicate that such areas of exposed ice contribute a disproportionate amount of melt to the glacier ablation budget. We used Structure from Motion photogrammetry with Multi-View Stereo to derive 3-D point clouds for nine ice cliffs on Khumbu Glacier, Nepal (in November 2015, May 2016 and October 2016). By differencing these clouds, we could quantify the magnitude, seasonality and spatial variability of ice cliff retreat. Mean retreat rates of 0.30-1.49 cm d-1 were observed during the winter interval (November 2015-May 2016) and 0.74-5.18 cm d-1 were observed during the summer (May 2016-October 2016). Four ice cliffs, which all featured supraglacial ponds, persisted over the full study period. In contrast, ice cliffs without a pond or with a steep back-slope degraded over the same period. The rate of thermo-erosional undercutting was over double that of subaerial retreat. Overall, 3-D topographic differencing allowed an improved process-based understanding of cliff evolution and cliff-pond coupling, which will become increasingly important for monitoring and modelling the evolution of thinning debris-covered glaciers.

Keywords: debris-covered glaciers, glacial geomorphology, glaciological instruments and methods, remote sensing, supraglacial lakes


In the coming decades, ongoing mass loss from Himalayan glaciers and changing runoff trends will affect the water resources of over a billion people, including those who require it for agricultural, energy production and domestic usage (Immerzeel and others, 2009, 2010; Lutz and others, 2014; Mukherji and others, 2015; Shea and Immerzeel,

2016). A negative mass-balance regime prevails across glaciers in the central and eastern Himalaya (Bolch and others, 2011; Fujita and Nuimura, 2011; Benn and others, 2012; Kaab and others, 2012, 2015; King and others,

2017), which are widely recognised to be out of equilibrium with current climate (Yang and others, 2006; Shrestha and Aryal, 2011; Salerno and others, 2015). Deglaciation is leading to the development of large proglacial lakes, which may expand rapidly through ice cliff calving (Bolch and others, 2008; Benn and others, 2012; Thompson and others, 2012; Thakuri and others, 2016), and pose potential glacial lake outburst flood hazards (e.g. Carrivick and Tweed, 2013, 2016; Rounce and others, 2016, 2017).

Debris-covered glaciers have a hummocky, pitted surface, caused by variable melt rates under different debris thicknesses, and include extensive coverage of ice cliffs and supraglacial ponds (Hambrey and others, 2008; Thompson and others, 2016; Watson and others, 2016, 2017). Studies using DEM differencing to quantify elevation change over debris-covered tongues have revealed an association

between glacier surface lowering and the presence of ice cliffs and supraglacial ponds (Immerzeel and others, 2014; Pellicciotti and others, 2015; Ragettli and others, 2016; Thompson and others, 2016), confirming historical ice cliff observations (e.g. Inoue and Yoshida, 1980; Sakai and others, 1998; Benn and others, 2001; Sakai and others, 2002). However, raster-based DEMs generally give a poor representation of steep slopes or steeply-sloping topography (Kolecka, 2012) and their differencing incorporates a mixed signal containing surface elevation change related to debris cover, ice cliff dynamics, supraglacial ponds and glacier emergence velocity (Vincent and others, 2016).

Models of glacier evolution do not consider pond dynamics or ice cliff dynamics explicitly, because this requires an understanding of their spatio-temporal distribution (e.g. Sakai and others, 2002; Watson and others, 2017), energy-balance modelling of the ice cliff surface (e.g. Reid and Brock, 2014; Steinerand others, 2015; Buri and others, 2016a, b) and cliff-scale observations of retreat rates (e.g. Brun and others, 2016). Several studies have exploited topographic models derived from unmanned aerial vehicle surveys of Lirung Glacier in the Langtang region of Nepal, to make substantial progress towards understanding ice cliff dynamics (Immerzeel and others, 2014; Steiner and others, 2015; Brun and others, 2016; Buri and others, 2016a, b; Miles and others, 2016a). However, techniques to perform direct comparisons of multi-temporal point clouds without simplification have yet to be exploited.

In this study, we explore ice cliff evolution using multi-temporal point clouds obtained on Khumbu Glacier, Nepal. Specifically we: (1) quantify the retreat of ice cliffs for pre-monsoon and monsoon time periods; (2) compare the spatial variation in retreat across ice cliff faces;and (3) assess the change in ice cliff morphology through time in relation to local topography and the presence of supraglacial ponds.


Field data were obtained on Khumbu Glacier in the Everest region of Nepal during three field campaigns (post-monsoon November 2015, pre-monsoon May 2016 and late-monsoon October 2016). The November 2015-May 2016 and May 2016-October 2016 survey intervals are referred to as 'winter' and 'summer', respectively. The Indian summer monsoon spans the months of June to mid-October (Bollasina and others, 2002; Bonasoni and others, 2008) and is when ~80% of annual precipitation falls (Wagnon and others, 2013).

Khumbu Glacier is ~17 km long, of which the lower 10 km is debris covered (Fig. 1) and the lower ~4 km is

essentially stagnant (Quincey and others, 2009). Supraglacial debris thickness is >2 m in this stagnating region and decreases up-glacier (Nakawo and others, 1986; Rowan and others, 2015). However, the thickness of the debris layer is locally heterogeneous owing to the pitted surface and the presence of ice cliffs and supraglacial ponds. We studied nine ice cliffs on the lower debris-covered glacier (Fig. 1), which is a region of particular interest since supraglacial ponds have begun to coalesce here over the past 5 years (Watson and others, 2016), and a large glacial lake is expected to form (Naito and others, 2000; Bolch and others, 2011; Haritashya and others, 2015).

3. DATA AND METHODS 3.1. Data collection

Terrestrial photographic surveys of nine ice cliffs were carried out during the three field campaigns. Our study cliffs represented ~2% of the total ice cliff extent on Khumbu Glacier, based on the top-edge cliff delineation of Watson and others (2017). We sought to survey cliffs that were broadly representative of the range of cliffs found on

Elevation (m)

| £4,900 _ 4,901 -5,100 fc 5,101 -5,300 | 5,301 -5,500 5,501 -5,700 5,701 - 5,900 ^■ 5,901 -6,000 | £6,000

Ice cliffs | Supraglacial ponds)

- 50 m contour

] Khumbu Glacier

Fig. 1. Ice cliffs and supraglacial ponds on Khumbu Glacier (a), located in eastern Nepal (b). Inset boxes show the location and ID of the ice cliffs surveyed (c). Cliff sites B and D included both northerly- and southerly-facing ice cliffs. The panchromatic background image is from the Pleiades satellite (7 October 2015), and corresponding ice cliffs and ponds are shown. Khumbu Glacier flows in a southerly direction.

Khumbu Glacier, with and without supraglacial ponds, of variable aspect and of variable size, noting the terrestrial survey constraints that precluded surveys of very large cliffs. Four of our nine study cliffs had a supraglacial pond present during the initial survey and the mean length of ice cliffs was 57 m. This compares to the observation that on Khumbu Glacier 47% of ice cliffs were associated with a pond in 2015, and cliffs had a mean length of 54 m (Watson and others, 2017). We note from Watson and others (2017) that cliffs 20-40 m in length were most common, but that some cliffs exceeded 200 m in length.

Two out of the seven individual study sites (Fig. 1) included both northerly- and southerly-facing cliff faces. These southerly-facing cliffs are labelled '-SF' hereafter. Within the first two field campaigns, surveys were conducted at intervals of 7-11 days at cliffs C, D, E and F, which are referred to as 'weekly' surveys. 'Seasonal' surveys refer to those between field campaigns. Each survey typically took <1 h and 122-564 photos were taken of each ice cliff with a highly convergent geometry (Fig. 2a) using a Panasonic DMC-TZ60 18.1 megapixel digital camera. In order to capture the surrounding topography, each photo was taken from a different position but was not necessarily orientated towards the ice cliff. High-contrast temporary ground control points (GCPs) (number of GCPs (n) = 6-15) were distributed around each ice cliff to encompass the survey area extents and surveyed using a Leica GS10 global navigation satellite system (GNSS). Each GCP was occupied in static mode for ~5 minutes. A base station was located on the lateral moraine of the glacier <2 km from our survey sites for the duration of each field campaign and was set to record each day.

3.2. Post-processing

Our GNSS base station data were post-processed against the Syangboche permanent station (27.8142N, 86.7125E) located ~20 km from our field site using GPS and GLObal NAvigation Satellite System (GLONASS) satellites. Our field GCPs were then adjusted with reference to the field base station data following a relative carrier phase positioning strategy. The mean 3-D positional uncertainty was 3.9 mm across all our GCPs (n = 281).

Photographs were input into Agisoft PhotoScan 1.2.3 to derive 3-D point clouds of the ice cliff topography following a Structure-from-Motion with Multiview Stereo (SfM-MVS) workflow (e.g. James and Robson, 2012; Westoby and others, 2012; Smith and others, 2015). First, photographs were aligned to produce a sparse point cloud by matching coincident features. This stage also estimated internal camera lens distortion parameters and scene geometry using a bundle adjustment with high redundancy, owing to large overlapping photographic datasets (Westoby and others, 2012). Only points with a reprojection error of <0.6 were retained and clear outliers (e.g. areas of shadow under overhanging cliffs) were removed manually. Second, GCPs were identified in each photograph to georeference the sparse cloud. GCP placement accuracy was <10 mm (e.g. Fig. 2b). Uncertainties from GCP placement and the post-processed coordinates were used as weights to optimise the point cloud georeferencing to minimise RMSE (Javernick and others, 2014; Stumpf and others, 2015; Smith and others, 2016; Westoby and others, 2016). Third, a dense point cloud was produced using PhotoScan's multiview stereo (MVS) algorithm (Fig. 3). The dense cloud was subsequently edited to remove points that were not on solid surfaces (e.

Fig. 2. The generation and georeferencing of ice cliff point clouds. Photographs of each ice cliff (a) were aligned to produce a sparse point cloud, which was georeferenced using high-contrast pink and yellow markers (b). Dense point clouds were produced and manually edited to remove points not on solid surfaces (e.g. supraglacial ponds) (c).

Fig. 3. Oblique views of the 3-D ice cliff point clouds in November 2015. Cliff IDs correspond to Table 1 and Figure 1. The profiles (red lines) correspond to Figure 7.

g. on supraglacial ponds) and clear outliers. All PhotoScan's processes were run on high quality settings. Georeferencing uncertainty in the final point clouds was <0.035 m (RMSE; Table 1). The final point clouds were sub-sampled using an octree filter in CloudCompare to unify point density across the surveys for each individual ice cliff. Subsequent point cloud densities ranged between 2185 and 14 581 points

per m2.

3.3. Ice cliff displacement

3.3.1. Ice cliff surveys between field campaigns The study cliffs were located down-glacier of the expected active/inactive ice transition on Khumbu Glacier (Quincey and others, 2009). However, small magnitude displacements (<3 m a-1) were observed in this region from dGPS surveys of tagged boulders (Supplementary Fig. 1). Therefore, to correct for ice cliff displacement between field campaigns (November 2015-May 2016 and May 2016-October 2016) we co-registered point clouds using image features that

could be identified in multiple surveys (e.g. identifiable marks on boulders). For each survey, the coordinates were derived for these features (as 'Markers' in PhotoScan), enabling a transform (2-D translation-rotation) to be calculated to co-register the later model with the earlier one. The RMSE between these co-registered features are reported in Table 1, and were used as the total error (ET) in subsequent point cloud differencing. Co-registration errors were subject to sub-debris melt over each differencing period; however, thick debris cover (1->2 m) (Nakawo and others, 1986) and low sub-debris melt rates of 0.0015 m d-1 (Inoue and Yoshida, 1980) over our study area minimised these errors. Glacier emergence velocity could not be calculated for this study;however, this is expected to be low for the slow-moving and gently-sloping debris-covered study area (Nuimura and others, 2012).

3.3.2. Ice cliff surveys within field campaigns To account for cliff displacement between repeat surveys within each field campaign (e.g. Cliff C 3 November 2015

Cliff ID Survey date Georeferencing RMSE m/number of GCPs Tie point RMSE pixels Co-registration RMSE m/number of GCPs

May 2016 to October 2016

November 2015 to May 2016

A 08/11/2015 0.017 (10) 0.96 *

13/05/2016 0.015 (7) 0.75 0.050 (7)* *

05/10/2016 0.014 (8) 0.73 0.262 (8)*

B 03/11/2015 0.025 (14) 0.97 *

16/05/2016 0.013 (14) 0.98 0.070 (12)* *

05/10/2016 0.012 (14) 0.84 0.334 (11)*

C 03/11/2015 0.012 (6) 1.05 *

10/11/2015 0.016 (7) 1.03

16/05/2016 0.011 (7) 0.92 0.025 (7)* *

25/05/2016 0.013 (7) 0.87

03/10/2016 0.016 (6) 0.86 0.183 (6)*

D 01/11/2015 0.016 (8) 1.05 *

10/11/2015 0.012 (8) 0.96

15/05/2016 0.011 (8) 0.80 0.100 (6)* *

25/05/2016 0.011 (8) 0.71

03/10/2016 0.011 (7) 0.68 0.176 (6)*

E 01/11/2015 0.034 (15) 0.87

12/11/2015 0.011 (15) 0.85 *

15/05/2016 0.018 (15) 0.84 0.122 (9)* *

26/05/2016 0.014 (15) 0.81

04/10/2016 0.013 (14) 0.89 0.140 (9)*

F 04/11/2015 0.012 (8) 1.14 *

14/11/2015 0.017 (9) 1.02

14/05/2016 0.008 (9) 1.05 0.041 (8)* *

25/05/2016 0.010 (9) 0.96

04/10/2016 0.007 (8) 1.02 0.186 (8)*

G 07/11/2015 0.017 (7) 0.98 *

15/05/2016 0.004 (6) 0.91 0.031 (7)* *

03/10/2016 0.011 (8) 1.01 0.224 (7)*

* For point cloud differencing, May 2016 data were co-registered with November 2015, and October 2016 were co-registered with May 2016.

and 10 November 2015), we take the mean daily displacement of the respective cliff between field campaigns (e.g. 0.0023 m d-1) and multiply this by the time-separation of the repeat models (e.g. 7 days). The resulting shift (e.g. 0.0161 m) was treated as an additional uncertainty in addition to the respective georeferencing errors. We did not shift these models as described in section 3.3.1, since the expected displacement was <0.04 m for all ice cliffs, which was similar to the uncertainty in identifying coincident features in each model.

To calculate the total error (ET) for each cliff model comparison within a field campaign, the individual errors were propagated using (1):

ET = ^ GE2C1+GE2C2+DE2C1_C2 (1)

Where GEC1 and GEC2 are the georeferencing RMS errors associated with clouds C1 and C2, and DEC1-C2 is the displacement error between clouds C1 and C2.

3.4. Point cloud characteristics and differencing

The mean slope and aspect of ice cliffs were calculated in CloudCompare using the dip direction and angle tool. The aspect of overhanging cliff sections required correction through 180°. The area of cliffs was calculated by fitting a mesh to each point cloud using the Poisson Surface Reconstruction tool in CloudCompare (Kazhdan and Hoppe, 2013).

Cloud-to-cloud differencing was carried out in the open source CloudCompare software using the Multiscale Model to Model Cloud Comparison (M3C2) method (e.g. Barnhart and Crosby, 2013; Lague and others, 2013; Gómez-Gutiérrez and others, 2015; Stumpf and others, 2015; Westoby and others, 2016). M3C2 was created by Lague and others (2013) to quantify the 3-D distance between two point clouds along the normal surface direction and provide a 95% confidence interval based on the point cloud roughness and co-registration uncertainty. The method is therefore ideally suited to quantifying statistically significant ice cliff evolution where the geometry changes in 3-D, and is robust to changes in point density and point cloud noise (Barnhart and Crosby, 2013; Lague and others, 2013).

For point clouds derived from photogrammetric techniques, uncertainty is spatially variable but highly correlated locally, since points in close proximity to one another are derived from the same images (James and others, 2017). Thus, although point-cloud roughness could represent a component of the measurement precision required to derive confidence intervals, it would not include broader photogrammetric contributions. In order to visualise spatially variable photogrammetric and georeferencing precision, we derived 3-D precision maps for May 2016 study cliffs using the method of James and others (2017) (Supplementary Fig. 2). Repeated bundle adjustments implemented in PhotoScan (4000 Monte Carlo iterations) were applied to the sparse point cloud of each ice cliff using GCP and tie

point uncertainties. Point precision estimates were then interpolated onto a 1 m raster grid in sfm_georef v3.0 (James and Robson, 2012; James and others, 2017). Large uncertainties were apparent at the survey edges (e.g. Cliff A, C, D, E), and around supraglacial ponds (e.g. Cliff E) due to poor photograph coverage. In contrast, the mean precision estimates ranged from 7 to 38 mm and were uniform across individual cliff faces. Hence, given the large magnitudes of the measured ice cliff changes, the M3C2-PM (Precision Maps) variant based on photogrammetry-derived precision maps was not required and the native M3C2 algorithm, implemented in CloudCompare, was used throughout.

M3C2 requires two user-defined parameters: (1) the normal scale D, which is used to calculate surface normals for each point and is dependent upon surface roughness and point cloud geometry, and (2) the projection scale d over which the cloud-to-cloud distance calculation is averaged, which should be large enough to average a minimum of 30 points (Lague and others, 2013). We estimated the normal scale D for each point cloud following a trial-and-error approach similar to that of Westoby and others (2016), to reduce the estimated normal error, Enorm (%), through refinement of a rescaled measure of the normal scale n(i):

n(i) =

n(i) is the normal scale D divided by the roughness a measured at the same scale around i. Where n(i) falls in the range 20-25, Enorm < 2% (Lague and others, 2013). In this study, normal scales D ranged from 1.5 to 8 m and the projection scale d was fixed at 0.3 m. The Level of Detection (LoD) threshold for a 95% confidence level is given by:

(d) = ±1.96

a1 (d)2 + a2(d)2 n1 n2

where a1 and a2 represent the roughness of each point in sub-clouds of diameter d and size n1 and n2, and reg is the cloud-to-cloud co-registration error, for which ET is substituted. The error is assumed to be isotropic and spatially uniform across the dataset (Lague and others, 2013).

Distance calculations were masked to exclude points where the change was lower than the Level of Detection (LoD) threshold and were clipped to individual ice cliff faces. Ice cliff retreat rates were divided by respective survey intervals to derive daily retreat rates. Since cliffs often exhibited large changes in geometry between surveys, some cliff normals at time one intersected debris cover at time two. The total retreat of the cliff face was therefore reported, in addition to cliff-to-cliff retreat (i.e. where cliff normals at time one intersected a cliff at time two).

3.5. Other data

Volumetric loss due to ice cliff retreat was estimated from DEM differencing using point clouds gridded at 0.5 m. Cliff retreat rates were also calculated using these volumetric changes for comparison with M3C2 retreat rates, by dividing the volume loss over the respective time period by the cliff area. The zone of ice cliff retreat was defined as the area connected by cliff outlines at respective time periods. Where a

cliff partially or completely degraded and hence was not represented by an outline at time two, the zone of cliff retreat was estimated using M3C2 distance measurements (representing spatially variable ice cliff retreat) from the cliff at time one. These distance measurements were used to define a variable-width buffer in ArcGIS to delineate the maximum extent of ice cliff retreat.

The drainage of the supraglacial pond adjacent to Cliff E provided an opportunity to reconstruct the bathymetry and maximum pond depth using the historic water level from May 2016, with the assumption that subaqueous basal melt and debris inputs to the basin were minimal. Additionally, air temperature at 1 m above the surface was recorded at 20 minute intervals on Khumbu Glacier using a Solinst Barologger Edge, which was sited behind Cliff G. Measurements were recorded from October 2015 to October 2016;however, the station collapsed in August 2016 due to the retreat of Cliff G. The logger was found in an air pocket buried by debris, hence data shown after this collapse revealed a subdued diurnal temperature cycle.


4.1. Summary ice cliff characteristics

Ice cliffs: ranged from 4 to 23 m in height; were all within a 43 m elevation range; had mean slopes of 50° to 73°; and were of variable aspect, with both northerly- and southerly-facing cliffs represented (Fig. 4, Table 2). Of the four cliffs with a supraglacial pond present in November 2015, only Cliff B-SF had a pond remaining in October 2016. Overall, four study cliffs persisted throughout the study period and the other five were buried under debris between May 2016 and October 2016 (Fig. 4c).

The mean slope, maximum cliff height, cliff area and mean cliff aspect were evaluated across our study period (Fig. 4). Southerly-facing cliffs generally featured higher mean slopes, and the greatest changes in cliff slope were observed on southerly-facing cliffs B-SF and E (Fig. 4a). Maximum cliff height reduced for all cliffs over the study period, although this change was generally small for those cliffs that persisted through the study. Cliff E, which lost ~5 m in height over summer (Fig. 4b), was an exception. Notably, persistent cliffs had a starting height >10 m; however, we note that Cliff F decayed, despite a starting height >10 m. With the exception of Cliff B-SF which increased in area, all other cliffs declined in area over the study; however, the rate of this area loss varied from cliff to cliff. Of the four cliffs persisting over the study, two were southerly-facing and two were northerly-facing and all had a supraglacial pond present for part of the study period (Table 2). However, pond dynamics between the observation dates were unknown. The largest changes in cliff aspect were for cliffs C and D-SF, which became increasingly northerly and westerly orientated by 25° and 23°, respectively (Fig. 4d).

4.2. Ice cliff retreat

With the exception of Cliff D-SF, cliff retreat rates were higher over summer than the preceding winter, which corresponds with consistently higher air temperatures during summer (Fig. 5, Table 3). Mean winter temperatures were generally below 0°C, whereas summer temperatures were generally

cliff A —a—dig B Cliff B SF

*-Cliff C —M—Cliff D -«- Cliff D SF

o- Cliff E -x-Cliff F -i-Cliff G

Nov 2015 May 2016 Oct 2016 Nov 2015 Way 2016 Oct 2016

Survey date Survey date

Fig. 4. The evolution of ice cliff mean slope (a), maximum height (b), area (c), and aspect (d) over the study period. Absolute cliff area change is shown in Supplementary Fig. 3.

above 0°C, although several discrete periods of positive air temperature also occurred in winter. Similarly, volumetric losses due to cliff retreat were generally higher during summer, although they were small where cliffs degraded (e.g. D-SF, Table 3). Notably, the M3C2- and DEM-based retreat rates were comparable for most cliffs.

The highest mean retreat rate occurred at Cliff B during summer, although this (along with the retreat at Cliff B-SF) was a combination of subaerial retreat and a large-scale cliff collapse involving a section of the cliff ~30 m in length. Excluding these cliff faces, the highest mean retreat rates observed and the largest seasonal differences in retreat rates were from ice cliffs with an adjacentsupraglacial pond including cliffs A (1.75 cm d-1), E (4.26 cm d-1) and G (2.62 cm

Table 2. Ice cliff characteristics in November 2015

Cliff ID Elevation Surface Maximum Mean Mean

m 2 area m height m slope ° aspect °

A1,2 4959 85 5 56 5

B2 4939 1278 17 54 37

B-SF1,2,3 4933 1313 23 73 235

C 4941 34 4 58 266

D 4933 56 4 57 15

D-SF 4935 37 5 63 210

E1,2 4926 782 15 65 161

F4 4916 757 13 58 132

G1,2 4923 357 10 50 1

1,2,3 Indicate the presence of a supraglacial pond during (1) November

2015, (2) May 2016 and (3) October 2016 surveys.

4 A supraglacial pond was present at Cliff F prior to field surveys based on Pleiades imagery (Fig. 1).

d-1) (Fig. 5b). The retreat rates for weekly surveys were generally higher than seasonal retreat rates, with the exception of Cliff E (Fig. 5). Retreat rates for cliffs that degraded over the study period (A, C, D, D-SF and F) included a transition to sub-debris melt during the summer, and hence were expected to have lower retreat rates and volume losses compared with persistent cliffs. Similarly, where cliffs partially degraded between survey intervals, the cliff-to-cliff retreatwas generally greater than the total retreat rates, which included areas of cliff-to-debris transition (Table 3). Here, retreat rates attributed only to persistent areas of cliff ranged from 0.36-1.68 cm d-1 (winter), and 3.84-5.85 cm d-1 (summer).

Across individual cliff faces, the observed retreat was related to the presence of supraglacial ponds, englacial conduits (expressed as an opening within or below ice cliff faces), cliff aspect, cliff slope and the formation of runnels (Fig. 6). Mean winter ice cliff retreat showed a clear relationship with aspect and peaked at a south easterly aspect of 155°, although this was not the case during summer (Fig. 6h). Maximum retreat rates (10.65 cm d-1) were observed at cliffs B and B-SF (Fig. 6b) where a notable calving event occurred in the summer. This was followed by northerly-facing Cliff G in association with a supraglacial pond, with maximum retreat rates of 6.18 cm d-1 in the zone of thermal undercutting (Fig. 6g). The surface of this pond was frozen between the November 2015 and May 2016 surveys and the pond had drained by October 2016.

Runnels were locations of locally differential retreat on the north-facing cliffs B and G during winter (Figs. 6b, g) and these cliffs also had the lowest mean initial slopes of 54° and 50°, respectively (Table 2). During winter, the relative increase in retreat at Cliff G at locations of runnels was ~0.12-0.24 cm d-1 (Fig. 6g). However, the presence of

3 Air temperature (20 minute intervals} Seven day moving average

Fig. 5. (a) Air temperature at 1 m above the surface recorded at 20 minute intervals with a seven day moving average. Survey intervals are indicated by vertical black lines. The logger mounting collapsed due to ice cliff retreat in August 2016 (shaded area represents data when the logger was partially buried by debris). Mean ice cliff retreat rates for the seasonal (b), and weekly surveys (c). Error bars show one standard deviation.

runnels was localised and, when considering the whole cliff, the rate of runnel retreat (~0.47-0.58 cm d-1) was otherwise comparable with the mean retreat rate during winter (0.47 cm d-1). Evidence of a vertical retreat gradient was apparent on several cliffs during winter (e.g. Figs. 6c, d, f), and was similar during summer, other than where thermal undercutting was apparent (e.g. Cliff G). Cliff B-SF featured the most apparent aspect-related control on retreat during winter with westerly facing ice melting at the slowest rate (~0.84 cm d-1) compared with southerly faces (~2.62 cm d-1, Fig. 6b). Cliff A featured an englacial conduit large enough to enable crouched access into a void behind the cliff face, which was the area of greatest retreat for this cliff (~4 cm d-1 during summer) as the void likely became exposed and degraded (Fig. 6a).

4.3. Ice cliff evolutionary traits

2-D profiles through selected ice cliffs revealed different characteristics of retreat, including ice cliff burial under debris, ice cliff collapse and undercutting by adjacent supraglacial ponds (Fig. 7). Cliff A maintained a similar slope during its retreat over winter, although during summer the slope angle decreased to ~35°, which led to burial under debris (Fig. 7a). There was a small supraglacial pond shallower than 1 m adjacent to Cliff A in the November 2015 and May 2016 surveys, and an associated undercut notch. The pond drained prior to the October 2016 survey, and the steep cliff back-slope led to a relaxation of the cliff slope and hence degradation of Cliff A by October 2016.

The profile through Cliff B revealed greater retreat on the southerly face through winter, compared with the northerly

Table 3. Mean ice cliff retreat rates and volume loss for winter and summer

Cliff November 2015 to May 2016 May 2016 to October 2016

ID ---

M3C2 M3C2 retreat DEM-based DEM-based M3C2 M3C2 retreat DEM-based DEM-based

retreat (cliff only)1 volume loss m3 retreat cm d-1 retreat (cliff only)1 volume loss m3 retreat cm d-1

cm d-1 cm d-1 cm d-1 cm d-1

A 0.46 0.46 68 0.30 1.75 - 420 1.20

B 0.65 0.70 797 0.66 5.18 5.85 12 4262 5.53

B-SF 1.49 1.55 3286 3.05 4.05 3.84

C 0.39 0.58 22 0.36 0.95 - 38 0.97

D 0.30 0.36 28 0.34 1.27 - 63 1.18

D-SF 0.74 0.92 69 0.80 0.74 - 5 0.33

E 1.44 1.68 1779 2.58 4.26 4.70 2950 4.55

F 1.26 1.57 1238 1.32 1.39 - 643 1.93

G 0.47 0.49 315 0.63 2.62 4.10 1786 2.44

1 M3C2 retreat rates (cliff only) represent cliff-to-cliff retreat i.e. excluding areas where ice cliff normals at time one intersected with debris cover at time two (see Supplementary Table 1).

2 Volume loss could not be separated at B and B-SF due to a calving event.

Fig. 6. Ice cliff retreat rates shown for winter (November 2015-May 2016) and summer (May 2016-0ctober 2016). Note the different scale ranges. Distance measurements are clipped to the study cliffs and indicative values are shown for key features. The mean and standard deviation of non-cliff surface elevation changes are reported for winter (w) and summer (s). Ice cliff retreat rate and initial aspect for winter and summer differencing periods are shown in (h), with a sinusoidal regression line (winter). Circled points indicate ice cliffs that disappeared during summer.

face (Fig. 7b). A section of Cliff B collapsed prior to the final survey in October 2016, suggesting that the undercut notch caused the cliff to collapse in a northerly direction. The

supraglacial pond in contact with Cliff B-SF expanded throughout the study and was in contact with the southerly face (although not at the 2-D profile shown), which

Fig. 7. 2-D ice cliff profiles for selected cliffs revealing topographic change over the study period. Ice cliff faces are shown as lines without a transparency, whereas debris-covered areas and water levels are shown with transparency.

exposed a new ice cliff face and caused an increase in cliff area (Fig. 4c). A supraglacial pond was present at the northerly-facing Cliff B in May 2016;however, the water level was well below historic water cut notches.

The two opposing faces of Cliff D and D-SF both became buried over the study (Fig. 7c). The southerly-facing cliff retreated faster than the northerly face during winter; however, the steep back-slope and small area of this southerly face limited the retreat during summer. Debris infill was apparent in the May 2016 profile, caused by the retreat of the southerly face, which had an inwardly sloping cliff top.

Cliff E featured a supraglacial pond over 9.95 m deep, which drained over the summer. There was evidence of deepening towards the cliff faces at profiles 1 and 2 and thermal undercutting of both cliff faces (Figs. 8c, d). The pond at Cliff G also drained over the summer and was also associated with an undercut notch. Cliff G had a gentle back-slope and maintained a similar slope (50°-54°) during retreat (Fig. 7d). The gentle back-slope of Cliff G allowed continued retreat, in contrast to cliffs A and D where a steep back-slope lead to cliff degradation.


5.1. Multi-temporal ice cliff surveys

In this study we have presented the first application of 3-D point cloud differencing to multi-temporal topographic surveys of ice cliffs, revealing evolutionary traits for a variety of ice cliffs present on the lower ablation area of Khumbu Glacier. This method has specific advantages for quantifying the importance of ice cliff retreat and the cliff/ pond interaction when compared with previous approaches. First, the retreat attributed to ice cliffs is calculated along the normal direction of the cliff face, thereby minimising the conflation of topographic change from debris cover, ice cliffs and supraglacial ponds that exists in vertical DEM differencing (e.g. Thompson and others, 2016). However, where a cliff decays between two survey dates, retreat calculations include topographic change related to processes in addition to cliff retreat such as sub-debris melt, hence short survey intervals are preferable. M3C2 allows quantification of the variability of retreat across a cliff face, for example in relation to slope and aspect (Buri and others, 2016a), the presence of runnels (Watson and others, 2017) and supraglacial ponds

Fig. 8. The drainage of a supraglacial pond at Cliff E (a). The drained supraglacial pond provided an opportunity to reconstruct the historic bathymetry (b and c). The data gap at the deepest part of the pond (intersecting with Profile 1) was caused by the remnant presence of water, which had not drained, estimated to be <1 m in depth. Point cloud profiles revealed subaerial ice cliff retreat and thermo-erosional undercutting (d). The yellow star denotes an area of the cliff that was present in November 2015 (a), but had degraded by May 2016 (d).

(Miles and others, 2016a). Second, the mechanism controlling topographic change can be evaluated in three dimensions, revealing the role of ice cliff back-slope in ice cliff persistence and thermo-erosional undercutting by supragla-cial ponds.

5.2. Ice cliff retreat

Observations of ice cliff retreat have previously been obtained from point ablation stake measurements or using static markers on the back-slopes of ice cliffs (e.g. Inoue and Yoshida, 1980; Sakai and others, 1998; Purdie and Fitzharris, 1999; Benn and others, 2001; Sakai and others, 2002; Han and others, 2010; Reid and Brock, 2014; Steiner and others, 2015). Use of ablation stakes restricts assessment of the spatial variation in retreat across an ice cliff face, since stake placements are likely to be aligned vertically down the cliff face and biased towards areas of comparatively safer access. In comparison, Brun and others (2016) used multitemporal fine-resolution topographic surveys to estimate the volumetric mass loss and mean retreat rates from five ice cliffs on the debris-covered Lirung Glacier. Their study cliffs were generally larger (maximum face size of 6441 m2 compared with 1313 m2 in this study), and at ~800 m lower elevation.

5.2.1. Ice cliff retreat through time Four of nine ice cliffs, which all had a maximum height >10 m (Fig. 4b) and featured an adjacent supraglacial pond (Table 2), persisted over 1 year in this study. In contrast, all study cliffs of Brun and others (2016) persisted over 1 year and were all ^9 m high. Mean retreat rates obtained in this study ranged from 0.30-1.49 cm d-1 (winter), and 0.745.18 cm d-1 (summer) and were comparable with those of Brun and others (2016), 0.70-1.20 cm d-1 (winter) and 2.2-4.5 cm d-1 (summer), despite the higher elevation and smaller size of our study cliffs. In our study, ice cliff retreat

rates were generally higher during summer and for southerly-facing ice cliffs, and summer featured the largest variability in retreat rates amongst cliffs (Figs. 5b, c). Lower retreat rates during winter correspond with cooler air temperatures (Fig. 5a). Higher retreat rates during the November 2015 weekly surveys compared with respective seasonal surveys, reflected warmer temperatures before winter. Similarly, the May 2016 surveys generally had higher retreat rates than the May 2016 to October 2016 surveys (Figs. 5b, c). For cliffs C, D and F, this likely reflects cliff degradation before October 2016 and hence a transition from cliff retreat to sub-debris melt for part of the survey interval. The unknown date of cliff degradation is a limitation to assessing mass loss relating only to ice cliff retreat. However, the retreat rates for cliffs that partially degraded between surveys can be quantified by considering the total cliff retreat, which includes areas where the cliff degraded, and the cliff-to-cliff retreat, where cliff normals at time one intersect a cliff at time two (Table 3, Supplementary Table 1). The former is representative of the total cliff evolution, which includes areas of ice cliff degrading and becoming buried by debris, whereas the latter is representative of the retreat rates for persisting areas of cliff.

Cliff B suffered a partial collapse of a ~30 m segment during summer, causing high mass loss due to the effect of this calving event (Fig. 5b, Table 3). Similarly high retreat rates May 2016-0ctober 2016 were observed at cliffs E and G, which both featured supraglacial ponds becoming active (i.e. thawed) during summer. The high variation in retreat rates over summer suggests that more frequent monitoring would be beneficial to assess cliff dynamics such as burial under debris and the interaction with seasonally expanding supraglacial ponds.

5.2.2. Cliff face variation in retreat

Visualising the spatial distribution of retreat across individual

cliff faces revealed vertical and lateral gradients, and

increased retreat attributed to supraglacial ponds during summer (Fig. 6). The influence of cliff aspect is apparent on Cliff B-SF, where a southerly face retreated 1.78 cm d-1 faster than an adjacent westerly face (Fig. 6b). Both the north-facing cliffs B and G displayed evidence of locally enhanced retreat attributed to the presence of vertical runnels observed during the winter surveys (Figs. 6b, g). Runnels were also observed by Watson and others (2017) on other ice cliffs on Khumbu Glacier. The low solar radiation receipt on northerly-facing cliffs during winter may mean that meltwater generated at the less shaded cliff top from sub-debris melt and from melt on the cliff face, is able to incise runnels faster than the background rate of cliff retreat. In contrast, during summer the higher magnitude of retreat likely masks this influence of micro-scale cliff topography (Fig. 6), although the runnels may persist. Runnels also act as preferential pathways for debris slumping from the cliff top and differential retreat may also occur in response to albedo variations across the cliff face. The morphology of the cliff face, including runnel development and self-shading, is therefore likely to locally influence retreat rates as evidenced in this study, but should be explored further with additional 3-D surveys in order to assess their importance seasonally, and at a glacier scale.

Cliff tops exhibited the highest retreat rates in several cases (e.g. Figs. 6c, d, f), which was also observed in the modelled retreat rates at two ice cliffs by Buri and others (2016b). Cliff G also displayed a vertical gradient during the summer; however, this was locally reversed in the area undercut by a supraglacial pond (Fig. 6g).

5.2.3. The influence of aspect

Several studies have observed a prevalence of northerly-facing ice cliffs on debris-covered glaciers in the Northern Hemisphere, suggesting that solar radiation receipt plays a key role in controlling ice cliff development (Sakai and others, 2002; Thompson and others, 2016; Kraaijenbrink and others, 2016b; Watson and others, 2017). Southerly-facing ice cliffs are expected to decay quickly after formation due to high solar radiation receipt, whereas northerly-facing cliffs are more persistent (Sakai and others, 2002). Slope relaxation was apparent on southerly-facing cliffs B-SF and E, which decreased by 14° and 13°, respectively;however, both of these cliffs persisted throughout the study period.

We observed highest ice cliff retreat rates on ice cliffs with a southeasterly aspect (155°) (Fig. 6h), although this trend was only apparent during winter. Also on Khumbu Glacier, Inoue and Yoshida (1980) revealed maximum ice cliff retreat at an aspect of ~190°. Cliff aspect likely has a stronger influence over cliff retreat in the winter due to the low solar angle and cliff self-shading (e.g. Steiner and others, 2015). Additionally, direct solar radiation receipt is reduced during the summer monsoon due to the prevalence of cloud cover (Supplementary Fig. 4). Therefore at this time, diffuse radiation, air temperature and local ice cliff characteristics such as the presence of a supraglacial pond were likely stronger controls on ice cliff retreat than the cliff aspect.

5.3. Local controls on ice cliff evolution

The back-slope of individual ice cliffs influences their longevity, since there is a finite volume of ice for the cliff to retreat into unless accompanied by simultaneous

downwasting of a supraglacial pond. In our study, slope relaxation and cliff degradation (Figs. 4a, 7a, c) were observed on both northerly- and southerly-facing ice cliffs. This contrasts with the observations of Sakai and others (2002) where slope relaxation was a trait of southerly-facing ice cliffs, highlighting the importance of local topography and cliff characteristics, which determine the longevity of individual ice cliffs over an ablation season.

Several studies have observed strong spatial coincidence of ice cliffs and supraglacial ponds (Thompson and others, 2016; Watson and others, 2017) and notable subaqueous melt rates (Sakai and others, 2009; Miles and others, 2016a). The potential importance of ponds for enhancing ice cliff retreat on Himalayan debris-covered glaciers is analogous to the 'wandering lakes' on Antarctic ice-cored moraines (e.g. Pickard, 1983). Thompson and others (2016) observed that 75% of ice cliffs were associated with a supraglacial pond on the Ngozumpa Glacier, and an average of 49% was observed by Watson and others (2017) across 14 glaciers in the Everest region of Nepal;however, these associations are likely to be seasonally variable. Our study revealed greater retreat for ice cliffs associated with a supraglacial pond, and mean retreat rates of pond-contact ice were estimated to be double that of subaerial retreat at Cliff G. However, the pond at Cliff G drained prior to the final survey such that the role of the pond could not be fully isolated from subaerial retreat. All persisting ice cliffs featured a supraglacial pond during their lifespan. We suggest that undercut notches allowed the cliff angle to be maintained during retreat, which promoted cliff persistence (e.g. Figs. 7d, 8d). Therefore our observations, in addition to strong association of cliffs and ponds (e.g. Watson and others, 2017), suggest that supraglacial ponds are likely to play a key role in ice cliff retreat and persistence at a glacier scale. However, quantifying subaqueous retreat using point cloud differencing is hindered by submerged topography, and manual field measurements (e.g. Rohl, 2006) are restricted by falling debris. Additionally, we cannot comment on the spatial variation in the importance of ice cliff retreat, which likely decreases with distance up-glacier using due to thinning debris cover (Thompson and others, 2016; Watson and others, 2017).

5.4. Implications for mass loss at a glacier scale

The ice cliff retreat rates observed in this study support previously observed associations between glacier surface lowering and the presence of ice cliffs and supraglacial ponds (Immerzeel and others, 2014; Pellicciotti and others, 2015; Thompson and others, 2016; Watson and others, 2017). Observed mean ice cliff retreat rates ranged from 0.30 (winter) to 5.18 cm d-1 (summer), which is much greater than sub-debris melt of 0.15 cm d-1 (Aug-1978) observed in a similar region of Khumbu Glacier by Inoue and Yoshida (1980). However, we note that these rates are not directly comparable since our observations represent surface-normal retreat, whereas sub-debris melt represents vertical lowering. The rate of surface lowering related to debris cover ranged from 0.03 to 0.31 cm d-1 on the nearby Ngozumpa Glacier based on the DEM differencing of Thompson and others (2016). However, surface lowering observed from DEM differencing is a function of sub-debris melt and emergence velocity. The latter was not quantified by Thompson and others (2016), or in this study, but was

shown to be +0.37 m w.e a-1 on the debris-covered Changri Nup Glacier (Vincent and others, 2016).

The volumetric loss at ice cliffs is variable and highlights the requirement to up-scale our methodology to the glacier scale in order to capture the full size distribution of ice cliffs present (Table 3). Additionally, knowledge of fine spatio-temporal dynamics of supraglacial ponds is still limited, but reveals potentially large seasonal expansion and contraction of ponds (e.g. Watson and others, 2016; Miles and others, 2016b). This restricts efforts to model the ice cliff/pond interaction (Buri and others, 2016a), or to quantify subaqueous retreat with multi-temporal point clouds. Nonetheless, our results suggest that undercut notches can promote ice cliff persistence by maintaining the slope angle during retreat. However, this requires further investigation at a glacier scale and over longer time periods, with particular attention to the role of undercutting for promoting calving events. A SfM-MVS methodology using time-lapse imagery is one such approach that could provide high temporal resolution.

5.5. Future work

M3C2 offers new opportunities to quantify 3-D topographic change on debris-covered glaciers and this could be used to explore debris redistribution and the formation of ice cliffs, which are currently limiting factors in modelling efforts (Buri and others, 2016a). Similarly, using point cloud data and M3C2 can address several problems related to fine spatio-temporal resolution DEM differencing, including the conflation of several processes contributing to the topographic change signal such as ice cliff retreat, sub-debris melt and supraglacial pond basal melt. Debris thickness estimated along the top edge of the cliff could be accounted for when slumping into supraglacial ponds, which may otherwise be counted as mass loss in DEM differencing. Comparisons of coincident measurements of 3-D and 2-D topographic change would therefore be highly beneficial to fully quantify their limitations. Additionally, the topographic change could be explored further with a greater dataset of ice cliff observations, to quantify specific relationships between cliff retreat and variables such as local slope, aspect and pond presence. However, our dataset demonstrates that ice cliff evolution is highly heterogeneous and that, when considering the dataset as a whole, the relationship between cliff retreat and slope, aspect and pond presence would be highly complex. Moving forward, conceptualising ice cliff evolution requires both local observations as presented in this study, and glacier scale multi-temporal ice cliff datasets (e.g. Watson and others, 2017).

The M3C2 method is not without its own limitations since it is difficult to calculate volumetric mass loss due to the variable alignment of surface normals;however, such methods are likely to become available or can be developed independently for similar applications (e.g. Brun and others, 2016). Non-uniform glacier surface displacement also presents issues when co-registering multi-temporal point clouds; however, this is arguably easier to achieve than using a DEM due to the availability of true-colour point data. However, DEMs and corresponding orthophotos can also be used for this correction (e.g. Kraaijenbrink and others, 2016a). Non-uniform glacier surface displacement is an important consideration if investigating lower magnitude processes such as sub-debris melt, which also requires quantification of glacier emergence velocity (Vincent and

others, 2016). Emergence velocity is expected to be low on slow-moving low gradient debris-covered glacier tongues (Nuimura and others, 2011 );however, positive surface elevation change was observed in this study (e.g. Fig. 6f), which was confirmed by independent dGPS boulder surveys (used in Supplementary Fig. 1). Additionally, point cloud precision estimates based on rigorous photogrammetric processing rather than surface roughness allow improved topographic change detection (James and others, 2017).


We have presented the first multi-temporal 3-D analysis of ice cliff evolution using 3-D point cloud differencing, which was necessary to quantify the spatial heterogeneity of retreat across individual cliff faces and their interaction with supraglacial ponds. Our results revealed the importance of a gentle cliff back-slope to allow continued retreat, and the role of supraglacial ponds in thermo-erosional undercutting, which maintains the cliff angle and delays burial under debris. Mean ice cliff retreat rates observed in this study ranged from 0.30-1.49 cm d-1 (winter), and 0.74-5.18 cm d-1 (summer). Additionally, the four ice cliffs persisting over our 1 year study period were all influenced by supraglacial ponds, and pond-contact ice was associated with a twofold increase in retreat at Cliff G. Our findings add further evidence to the role of ice cliffs as 'hot-spots' of mass loss on heavily debris-covered glaciers and contribute to a previously sparse dataset of ice cliff observations, revealing local controls on cliff retreat, which can be used to validate emerging models of ice cliff evolution (Brun and others, 2016; Buri and others, 2016a).

We observed an aspect-related control on ice cliff retreat during winter; however, local ice cliff characteristics masked any cliff-scale aspect related control on retreat during summer. We observed examples of northerly- and southerly-facing cliffs persisting, but also examples of cliff burial under debris. The controlling factors for ice cliff persistence appeared to be cliffs with a maximum height >10 m and with supraglacial pond influence. Nonetheless, the prevalence of northerly-facing cliffs on debris-covered glaciers in the northern hemisphere (Sakai and others, 2002; Brun and others, 2016; Watson and others, 2017) suggests that over longer timescales (e.g. decadal) the persistence of northerly-facing cliffs is greater in response to self-shading and supraglacial pond association.

M3C2 point cloud differencing was shown to be an effective tool to quantify the spatio-temporal magnitude of retreat across ice cliff faces, and to offer a new opportunity to validate models of ice cliff evolution. It is also more practical than point-based ablation stake measurements. M3C2 could be applied to glacier scale point clouds to enable surface elevation change to be partitioned into surface-normal (ice cliff retreat) and vertical (sub-debris and subaqueous melt) components, and should be compared with the prevailing practice of DEM differencing. These 3-D point cloud data provide a much more realistic representation of surface area compared with a planimetric DEM, and minimise the conflation of different topographic change signals that are common to DEM differencing.


The supplementary material for this article can be found at


C.S.W acknowledges support from the School of Geography at the University of Leeds, the Mount Everest Foundation, the British Society for Geomorphology, the Royal Geographical Society (with IBG), the Petzl Foundation, and water@leeds. The Natural Environment Research Council Geophysical Equipment Facility is thanked for loaning Global Navigation Satellite Systems receivers and technical assistance under loan numbers 1050, 1058 and 1065. Dhananjay Regmi and Himalayan Research Expeditions are thanked for fieldwork support including research permit acquisition, and Mahesh Magar is thanked for invaluable support during data collection. Patrick Wagnon is thanked for providing access to AWS data. This AWS has been funded by the French Service d'Observation GLACIOCLIM, the French National Research Agency (ANR) through ANR-13-SENV-0005-04/05-PRESHINE, and has been supported by a grant from Labex 0SUG@2020 (Investissements d'avenir - ANR10 LABX56). CloudCompare (version 2.8, 2016) is GPL software retrieved from Pleiades images were supplied by Airbus Defence and Space through a Category-1 agreement with the European Space Agency (ID Nr. 32600). We thank two anonymous reviewers for thorough and constructive reviews.


Barnhart T and Crosby B (2013) Comparing two methods of surface change detection on an evolving thermokarst using high-temporal-frequency terrestrial laser scanning, selawik river, Alaska. Remote Sens., 5(6), 2813-2837 Benn DI, Wiseman S and Hands KA (2001) Growth and drainage of supraglacial lakes on debris-mantled Ngozumpa Glacier, Khumbu Himal, Nepal. J. Glaciol., 47(159), 626-638 Benn DI and 9 others (2012) Response of debris-covered glaciers in the Mount Everest region to recent warming, and implications for outburst flood hazards. Earth-Sci. Rev., 114(1-2), 156-174 Bolch T, Buchroithner MF, Peters J, Baessler M and Bajracharya S (2008) Identification of glacier motion and potentially dangerous glacial lakes in the Mt. Everest region/Nepal using spaceborne imagery. Nat. Hazards Earth Syst. Sci, 8(6), 1329-1340 Bolch T, Pieczonka T and Benn DI (2011) Multi-decadal mass loss of glaciers in the Everest area (Nepal Himalaya) derived from stereo imagery. Cryosphere, 5(2), 349-358 Bollasina M, Bertolani L and Tartari G (2002) Meteorological observations in the Khumbu Valley, Nepal Himalayas. Bull. Glaciol. Res., 19, 1-11

Bonasoni P and 21 others (2008) The ABC-pyramid atmospheric research observatory in Himalaya for aerosol, ozone and halo-carbon measurements. Sci. Total Environ., 391(2-3), 252-261 Brun F and 9 others (2016) Quantifying volume loss from ice cliffs on debris-covered glaciers using high-resolution terrestrial and aerial photogrammetry. J. Glaciol., 62(234), 684-695 Buri P and 5 others (2016a) A physically-based 3D-model of ice cliff evolution over debris-covered glaciers. J. Geophys. Res.: Earth Surf., 121(12), 2471-2493 Buri P, Pellicciotti F, Steiner JF, Miles ES and Immerzeel WW (2016b) A grid-based model of backwasting of supraglacial ice cliffs on debris-covered glaciers. Ann. Glaciol., 57(71), 199-211 Carrivick JL and Tweed FS (2013) Proglacial lakes: character, behaviour and geological importance. Quat. Sci. Rev., 78, 34-52 Carrivick JL and Tweed FS (2016) A global assessment of the societal impacts of glacier outburstfloods. Glob. Planet. Change, 144,1-16 Fujita K and Nuimura T (2011) Spatially heterogeneous wastage of Himalayan glaciers. Proc. Natl. Acad. Sci. USA, 108(34), 14011-14014

Gómez-Gutiérrez Á, de Sanjosé-Blasco J, Lozano-Parra J, Berenguer-Sempere F and de Matías-Bejarano J (2015) Does HDR pre-processing improve the accuracy of 3D models obtained by means of two conventional SfM-MVS software packages? The case of the Corral del Veleta rock glacier. Remote Sens., 7(8), p10269 Hambrey MJ and 5 others (2008) Sedimentological, geomorpho-logical and dynamic context of debris-mantled glaciers, Mount Everest (Sagarmatha) region, Nepal. Quat. Sci. Rev., 27(25-26), 2361-2389

Han H, Wang J, Wei J and Liu S (2010) Backwasting rate on debris-covered Koxkar glacier, Tuomuer mountain, China. J. Glaciol., 56(196), 287-296 Haritashya UK, Pleasants MS and Copland L (2015) Assessment of the evolution in velocity of two debris-covered glaciers in Nepal and New Zealand. Geogr. Ann.: Ser. A, Phys. Geogr., 97(4), 737-751

Immerzeel WW, Droogers P, de Jong SM and Bierkens MFP (2009) Large-scale monitoring of snow cover and runoff simulation in Himalayan river basins using remote sensing. Remote Sens. Environ., 113(1), 40-49 Immerzeel WW, van Beek LPH and Bierkens MFP (2010) Climate change will affect the Asian water towers. Science, 328(5984), 1382-1385

Immerzeel WW and 6 others (2014) High-resolution monitoring of Himalayan glacier dynamics using unmanned aerial vehicles. Remote Sens. Environ., 150, 93-103 Inoue J and Yoshida M (1980) Ablation and heat exchange over the

Khumbu glacier. J. Japanese Soc. Snow Ice, 39, 7-14 James MR and Robson S (2012) Straightforward reconstruction of 3D surfaces and topography with a camera: accuracy and geo-science application. J. Geophys. Res.: Earth Surf., 117 F03017, DOI: 10.1029/2011JF002289. James MR, Robson S and Smith MW (2017) 3-D uncertainty-based topographic change detection with structure-from-motion photo-grammetry: precision maps for ground control and directly geor-eferenced surveys. Earth Surf. Process. Landf. DOI: 10.1002/ esp.4125

Javernick L, Brasington J and Caruso B (2014) Modeling the topography of shallow braided rivers using structure-from-motion photogrammetry. Geomorphology, 213, 166-182 Kaab A, Berthier E, Nuth C, Gardelle J and Arnaud Y (2012) Contrasting patterns of early twenty-first-century glacier mass change in the Himalayas. Nature, 488(7412), 495-498 Kaab A, Treichler D, Nuth C and Berthier E (2015) Brief communication: contending estimates of 2003-2008 glacier mass balance over the Pamir-Karakoram-Himalaya. Cryosphere, 9(2), 557-564 Kazhdan M and Hoppe H (2013) Screened poisson surface reconstruction. ACM Trans. Graphics (TOG), 32(3), p29 King O, Quincey DJ, Carrivick JL and Rowan AV (2017) Spatial variability in mass loss of glaciers in the Everest region, central Himalayas, between 2000 and 2015. Cryosphere, 11(1), 407426

Kolecka N (2012) Vector algebra for Steep Slope Model analysis.

Landform Analysis, 21, 17-25 Kraaijenbrink P and 5 others (2016a) Seasonal surface velocities of a Himalayan glacier derived by automated correlation of unmanned aerial vehicle imagery. Ann. Glaciol., 57(71), 103113

Kraaijenbrink PDA, Shea JM, Pellicciotti F, Jong SMd and Immerzeel WW (2016b) Object-based analysis of unmanned aerial vehicle imagery to map and characterise surface features on a debris-covered glacier. Remote Sens. Environ., 186, 581595

Lague D, Brodu N and Leroux J (2013) Accurate 3D comparison of complex topography with terrestrial laser scanner: application to the Rangitikei canyon (N-Z). ISPRS J. Photogramm. Remote Sens., 82, 10-26

Lutz AF, Immerzeel WW, Shrestha AB and Bierkens MFP (2014) Consistent increase in High Asia's runoff due to increasing

glacier melt and precipitation. Nature Clim. Change, 4(7), 587592

Miles ES and 5 others (2016a) Refined energy-balance modelling of a supraglacial pond, Langtang Khola, Nepal. Ann. Glaciol., 57 (71), 29-40

Miles ES, Willis IC, Arnold NS, Steiner J and Pellicciotti F (2016b) Spatial, seasonal and interannual variability of supraglacial ponds in the Langtang Valley of Nepal, 1999-2013. J. Glaciol., 63(237), 1-18

Mukherji A, Molden D, Nepal S, Rasul G and Wagnon P (2015) Himalayan waters at the crossroads: issues and challenges. Int. J. Water Res. Dev., 31(2), 151-160 Naito N, Nakawo M, Kadota T and Raymond CF (2000) Numerical simulation of recent shrinkage of Khumbu Glacier, Nepal Himalayas. In Nakawo M, Raymond CF and Fountain A, eds. IAHS publ. 264 (Symposium at Seattle 2000 - debris-covered glaciers), Seattle, Washington, USA. IAHS Publication, 245-254 Nakawo M, Iwata S, Watanabe O and Yoshida M (1986) Processes which distribute supraglacial debris on the Khumbu Glacier, Nepal Himalaya. Ann. Glaciol., 8, 129-131 Nuimura T and 5 others (2011) Temporal changes in elevation of the debris-covered ablation area of Khumbu Glacier in the Nepal Himalaya since 1978. Arct. Antarct. Alp. Res., 43(2), 246-255 Nuimura T, Fujita K, Yamaguchi S and Sharma RR (2012) Elevation changes of glaciers revealed by multitemporal digital elevation models calibrated by GPS survey in the Khumbu region, Nepal Himalaya, 1992-2008. J. Glaciol., 58(210), 648-656 Pellicciotti F and 5 others (2015) Mass-balance changes of the debris-covered glaciers in the Langtang Himal, Nepal, from 1974 to 1999. J. Glaciol., 61(226), 373-386 Pickard J (1983) Surface lowering of ice-cored moraine by wandering lakes. J. Glaciol., 29(102), 338-342 Purdie J and Fitzharris B (1999) Processes and rates of ice loss at the terminus of Tasman Glacier, New Zealand. Glob. Planet. Change, 22(1-4), 79-91 Quincey DJ, Luckman A and Benn D (2009) Quantification of Everest region glacier velocities between 1992 and 2002, using satellite radar interferometry and feature tracking. J. Glaciol., 55 (192), 596-606

Ragettli S, Bolch T and Pellicciotti F (2016) Heterogeneous glacier thinning patterns over the last 40 years in Langtang Himal, Nepal. Cryosphere, 10(5), 2075-2097 Reid TD and Brock BW (2014) Assessing ice-cliff backwasting and its contribution to total ablation of debris-covered Miage glacier, Mont Blanc massif, Italy. J. Glaciol., 60(219), 3-13 Rohl K (2006) Thermo-erosional notch development at fresh-water-calving Tasman Glacier, New Zealand. J. Glaciol., 52(177), 203-213 Rounce D, Watson C and McKinney D (2017) Identification of hazard and risk for glacial lakes in the Nepal Himalaya using satellite imagery from 2000-2015. Remote Sens., 9(7) 654; DOI: 10.3390/rs9070654 Rounce DR, McKinney DC, Lala JM, Byers AC and Watson CS (2016) A new remote hazard and risk assessment framework for glacial lakes in the Nepal Himalaya. Hydrol. Earth Syst. Sci, 20(9), 3455-3475 Rowan AV, Egholm DL, Quincey DJ and Glasser NF (2015) Modelling the feedbacks between mass balance, ice flow and debris transport to predict the response to climate change of debris-covered glaciers in the Himalaya. Earth Planet. Sci. Lett., 430, 427-438

Sakai A, Nakawo M and Fujita K (1998) Melt rate of ice cliffs on the Lirung Glacier, Nepal Himalayas. Bull. Glaciol. Res., 16, 57-66 Sakai A, Nakawo M and Fujita K (2002) Distribution characteristics and energy balance of ice cliffs on debris-covered glaciers, Nepal Himalaya. Arct. Antarct. Alp. Res., 34(1), 12-19 Sakai A, Nishimura K, Kadota T and Takeuchi N (2009) Onset of calving at supraglacial lakes on debris-covered glaciers of the Nepal Himalaya. J. Glaciol., 55(193), 909-917 Salerno F and 10 others (2015) Weak precipitation, warm winters and springs impact glaciers of south slopes of Mt. Everest (central Himalaya) in the last 2 decades (1994-2013). Cryosphere, 9(3), 1229-1247 Shea JM and Immerzeel WW (2016) An assessment of basin-scale glaciological and hydrological sensitivities in the Hindu Kush-Himalaya. Ann. Glaciol., 57(71), 308-318 Shrestha AB and Aryal R (2011) Climate change in Nepal and its impact on Himalayan glaciers. Reg. Environ. Change, 11, 65-77 Smith MW, Carrivick JL and Quincey DJ (2015) Structure from motion photogrammetry in physical geography. Prog. Phys. Geogr., 40(2), 1-29 Smith MW and 6 others (2016) Aerodynamic roughness of glacial ice surfaces derived from high-resolution topographic data. J. Geophys. Res.: Earth Surf., 121(4), p2015JF003759 Steiner JF and 5 others (2015) Modelling ice-cliff backwasting on a debris-covered glacier in the Nepalese Himalaya. J. Glaciol., 61(229), 889-907 Stumpf A, Malet JP, Allemand P, Pierrot-Deseilligny M and Skupinski G (2015) Ground-based multi-view photogrammetry for the monitoring of landslide deformation and erosion. Geomorphology, 231, 130-145 Thakuri S, Salerno F, Bolch T, Guyennon N and Tartari G (2016) Factors controlling the accelerated expansion of Imja Lake, Mount Everest region, Nepal. Ann. Glaciol., 57(71), 245-257 Thompson S, Benn D, Mertes J and Luckman A (2016) Stagnation and mass loss on a Himalayan debris-covered glacier: processes, patterns and rates. J. Glaciol., 62(233), 467-485 Thompson SS, Benn DI, Dennis K and Luckman A (2012) A rapidly growing moraine-dammed glacial lake on Ngozumpa Glacier, Nepal. Geomorphology, 145, 1-11 Vincent C and 10 others (2016) Reduced melton debris-covered glaciers: investigations from Changri Nup Glacier, Nepal. Cryosphere, 10(4), 1845-1858 Wagnon P and 11 others (2013) Seasonal and annual mass balances of Mera and Pokalde glaciers (Nepal Himalaya) since 2007. Cryosphere, 7(6), 1769-1786 Watson CS, Quincey DJ, Carrivick JL and Smith MW (2016) The dynamics of supraglacial ponds in the Everest region, central Himalaya. Glob. Planet. Change, 142, 14-27 Watson CS, Quincey DJ, Carrivick JL and Smith MW (2017) Ice cliff dynamics in the Everest region of the Central Himalaya. Geomorphology, 278, 238-251 Westoby MJ, Brasington J, Glasser NF, Hambrey MJ and Reynolds JM (2012) 'Structure-from-Motion' photogrammetry: a low-cost, effective tool for geoscience applications. Geomorphology, 179, 300-314 Westoby MJ and 6 others (2016) Interannual surface evolution of an Antarctic blue-ice moraine using multi-temporal DEMs. Earth Surf. Dynam, 4(2), 515-529 Yang X and 6 others (2006) Climate change in Mt. Qomolangma region since 1971. J. Geogr. Sci., 16(3), 326-336

MS received 16 February 2017 and accepted in revised form 25 July 2017