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

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⁻¹ were observed during the winter interval (November 2015–May 2016) and 0.74–5.18 cm d⁻¹ 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.

1. INTRODUCTION

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; Kääb 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; Steiner and 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.

2. STUDY SITE

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).

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.

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 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.

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).

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.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 m².

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.

Table 1. Summary statistics for each ice cliff model

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

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⁻¹) 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 (E _T) 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⁻¹ (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 and 10 November 2015), we take the mean daily displacement of the respective cliff between field campaigns (e.g. 0.0023 m d⁻¹) 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 (E _T) for each cliff model comparison within a field campaign, the individual errors were propagated using (1):

(1) $$E_{T} = \sqrt {{GE}_{{C1}}^{2} {+ GE}_{{C2}}^{2} {+ DE}_{{C}1 - {C}2}^2} $$

Where GE _C1 and GE _C2 are the georeferencing RMS errors associated with clouds C1 and C2, and DE _C1–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, E _norm (%), through refinement of a rescaled measure of the normal scale n(i):

(2) $$n(i) = \; \displaystyle{D \over {\sigma _i (D)}}$$

n(i) is the normal scale D divided by the roughness σ measured at the same scale around i. Where n(i) falls in the range 20–25, E _norm < 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:

(3) $${\rm LOD}_{95\%} \; (d) = \pm 1.96\left( {\sqrt {\displaystyle{{\sigma _1 (d)^2} \over {n_1}} + \displaystyle{{\sigma _2 (d)^2} \over {n_2}}} \; + {\rm reg}} \right)\; $$

where σ ₁ and σ ₂ represent the roughness of each point in sub-clouds of diameter d and size n ₁ and n ₂, and reg is the cloud-to-cloud co-registration error, for which E _T 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. RESULTS

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).

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.

Table 2. Ice cliff characteristics in November 2015

^1,2,3 Indicate the presence of a supraglacial pond during (1) November 2015, (2) May 2016 and (3) October 2016 surveys.

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

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 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.

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.

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

¹ 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).

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

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 adjacent supraglacial pond including cliffs A (1.75 cm d⁻¹), E (4.26 cm d⁻¹) and G (2.62 cm d⁻¹) (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 retreat was 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⁻¹ (winter), and 3.84–5.85 cm d⁻¹ (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⁻¹) 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⁻¹ 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.

Fig. 6. Ice cliff retreat rates shown for winter (November 2015–May 2016) and summer (May 2016–October 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.

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⁻¹ (Fig. 6g). However, the presence of runnels was localised and, when considering the whole cliff, the rate of runnel retreat (~0.47–0.58 cm d⁻¹) was otherwise comparable with the mean retreat rate during winter (0.47 cm d⁻¹). 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⁻¹) compared with southerly faces (~2.62 cm d⁻¹, 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⁻¹ 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.

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.

The profile through Cliff B revealed greater retreat on the southerly face through winter, compared with the northerly 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 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.

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).

5. DISCUSSION

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 (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 supraglacial 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 multi-temporal 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 m² compared with 1313 m² 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⁻¹ (winter), and 0.74–5.18 cm d⁻¹ (summer) and were comparable with those of Brun and others (2016), 0.70–1.20 cm d⁻¹ (winter) and 2.2–4.5 cm d⁻¹ (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–October 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⁻¹ 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⁻¹ (summer), which is much greater than sub-debris melt of 0.15 cm d⁻¹ (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⁻¹ 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⁻¹ 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).

6. CONCLUSIONS

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⁻¹ (winter), and 0.74–5.18 cm d⁻¹ (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.