1. INTRODUCTION
Debris-covered glaciers are notable features in the Pamir-Karakoram-Himalaya (PKH) region (e.g. Scherler and others, Reference Scherler, Bookhagen and Strecker2011; Bolch and others, Reference Bolch2012), covering ~10% of the glacierized area (Bolch and others, Reference Bolch2012), and are of importance for melt and mass balance because they normally lay at low elevations. Despite an insulating effect of debris when thick enough (e.g. Østrem, Reference Østrem1959; Nicholson and Benn, Reference Nicholson and Benn2006; Reid and Brock, Reference Reid and Brock2010; Lejeune and others, Reference Lejeune, Bertrand, Wagnon and Morin2013), recent large-scale studies based on remote sensing have provided evidence that they might be loosing mass at rates comparable with those of debris-free ice (e.g. Kääb and others, Reference Kääb, Berthier, Nuth, Gardelle and Arnaud2012; Nuimura and others, Reference Nuimura, Fujita, Yamaguchi and Sharma2012; Gardelle and others, Reference Gardelle, Berthier, Arnaud and Kääb2013; Holzer and others, Reference Holzer2015), which has been referred to as the debris-covered glacier anomaly (e.g. Pellicciotti and others, Reference Pellicciotti2015).
A number of pioneering (e.g. Sakai and others, Reference Sakai, Takeuchi, Fujita and Nakawo2000, Reference Sakai, Nakawo and Fujita2002) and more recent studies (Steiner and others, Reference Steiner2015; Buri and others, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016; Miles and others, Reference Miles2016) have suggested that supraglacial lakes and cliffs are responsible for larger than expected volume losses. These features have low albedo and are exposed directly to the atmosphere, enhancing the radiative transfer and turbulent energy fluxes at the glacier surface (Buri and others, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016; Miles and others, Reference Miles2016). Due to their steep slopes (commonly >45° and sometimes overhanging) and the complex topography of debris-covered glaciers, the cliffs also receive additional longwave radiation emitted from the surrounding debris-covered terrain, thus further enhancing their melt (Reid and Brock, Reference Reid and Brock2014; Steiner and others, Reference Steiner2015; Buri and others, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016).
However, actual quantitative estimates of the total contribution of ice cliffs backwasting (i.e. the volume loss at cliff scale) to the total mass balance of debris-covered glaciers do not exist. Modelling studies are still in their infancy and cannot be yet used to quantify the total contribution of cliffs to the glacier mass balance. Most numerical models (based on calculation of the energy balance at the cliff surface) have been developed at the point scale (Sakai and others, Reference Sakai, Nakawo and Fujita1998, Reference Sakai, Takeuchi, Fujita and Nakawo2000; Han and others, Reference Han, Wang, Wei and Liu2010; Reid and Brock, Reference Reid and Brock2014; Steiner and others, Reference Steiner2015) and some have quantified the total contribution of cliffs to melt by extrapolation of point estimates of backwasting (Sakai and others, Reference Sakai, Nakawo and Fujita1998, Reference Sakai, Takeuchi, Fujita and Nakawo2000; Han and others, Reference Han, Wang, Wei and Liu2010; Reid and Brock, Reference Reid and Brock2014). Most of these models were validated only against point scale measurements and over short periods of time given the difficulty of maintaining ablation stakes on cliffs (Reid and Brock, Reference Reid and Brock2014; Steiner and others, Reference Steiner2015) or against cliff edge retreat measurements (Sakai and others, Reference Sakai, Nakawo and Fujita1998; Han and others, Reference Han, Wang, Wei and Liu2010). A recent attempt has been made to develop a fully distributed model of cliff backwasting that could provide estimates of total cliff backwasting volumes and patterns (Buri and others, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016), but that model focuses on mass losses due to the energy exchange with the atmosphere and still does not include all relevant processes contributing to cliff volume changes. In addition, numerical models need validation through distributed ground data that, to our knowledge, are not yet available.
Thus, this paper's main aim is to develop a method that allows calculation of total backwasting volumes of supra-glacial cliffs from dense point clouds of glacier and cliff topography obtained via terrestrial and aerial photogrammetry, using Structure from Motion (SfM) techniques (Westoby and others, Reference Westoby, Brasington, Glasser, Hambrey and Reynolds2012) or, alternatively using dGPS measurements of the cliff outline, which are comparable with the type of data that may be obtained from very-high-resolution optical satellite photogrammetry. The volumetric change estimates obtained in this way can be used both to quantify total volume losses and to validate the models currently being developed.
In formulating the new method proposed here, we focus on an approach to provide total volume losses, but not the spatial pattern of backwasting. Several methods have been proposed to estimate temporal changes of complex surfaces. One of the most refined methods is the M3C2 algorithm implemented in the Cloud Compare open-source software (Lague and others, Reference Lague, Brodu and Leroux2013). This method enables computation of distances between two point clouds. Its main strength is that it does not require a priori information on the surface orientation, as it calculates local normals to the surface using neighboring points. However, a consequence of the method's practicality is that the surface change is measured perpendicular to the surface, which presents a challenge for ice cliff backwasting, where sloped features undergo a linear translation. Given the limitation of tools like M3C2, we developed a specific method to estimate volumetric change associated with ice cliff backwasting, which we applied to selected cliffs on the debris-covered tongue of Lirung Glacier in the Nepalese Himalaya. We used an existing dataset of triangulated irregular networks (TINs) obtained from unmanned aerial vehicle (UAV) surveys (Immerzeel and others, Reference Immerzeel2014; Kraaijenbrink and others, Reference Kraaijenbrink2016) and a novel TIN dataset derived from terrestrial photogrammetric survey. We calculate cliff backwasting over two monsoon seasons (May–October 2013 and May–October 2014) and one winter season (October 2013–May 2014). We then compare these estimates with those obtained with a simpler approximation similar to the method of Han and others (Reference Han, Wang, Wei and Liu2010) and discuss whether simplified geometry assumptions can produce a realistic estimate of backwasting. We also document the cliff evolution as inferred from our detailed estimates of cliff geometry and volumetric loss and discuss the possible processes responsible for the observed evolution.
2. STUDY AREA
Lirung Glacier (28.24°N; 85.56°E) is a debris-covered glacier located in the Upper Langtang Valley of the Central Nepalese Himalayas (Fig. 1, inset). It flows from Langtang Lirung peak (7234 m a.s.l.) down to ~4000 m a.s.l. (e.g. Sakai and others, Reference Sakai, Nakawo and Fujita1998; Immerzeel and others, Reference Immerzeel2014; Buri and others, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016). In this study, we focus on the lower ablation zone (ranging between 4000 and 4450 m a.s.l.), which is entirely covered by a heterogeneous debris layer and is detached from the glacier upper area (Ragettli and others, Reference Ragettli2015).
Fig. 1. Map of Lirung Glacier tongue showing the cliffs surveyed in May 2014 (black polygons and dots) and in October 2014 (red polygons and dots). The blue polygons show the areas where terrestrial photographs were taken to measure ice cliff backwasting. The background is an orthophoto taken in October 2013 (Immerzeel and others, Reference Immerzeel2014). The black thick line represents the glacier outline, manually delineated based on the May 2013 UAV-based orthophoto (Immerzeel and others, Reference Immerzeel2014). All the coordinates are in UTM 45N/WGS84 (m). The inset shows the location of the Langtang catchment in Nepal.
Ice cliffs on Lirung Glacier have been monitored during four field campaigns in May 2013, October 2013, May 2014 and October 2014. The monitored cliffs have a height ranging between 10 and 30 m, a width ranging between 40 and 200 m and a mean slope ~45° (Table 1). For the remainder of this study, we focus on the cliffs indicated as cliff 1–5 (Fig. 1; Table 1). Cliff 4’ indicates the western side of cliff 4, which was observed for the first time in October 2013 (this area was previously covered by debris) and has been surveyed since May 2014.
Table 1. Characteristics of the studied cliffs
For all cliffs we present the mean value of multiple years of observations (Fig. 7a for dates). The mean slope and aspect is calculated from 20 cm DEMs derived from the cliff TINs and therefore are only indicative (because for an equal area a steep slope is represented by fewer pixels).
* Cliff 4’ indicates the western side of cliff 4, which was observed for the first time in October 2013 and has been surveyed since May 2014.
3. DATA
3.1. Differential GPS data
The dGPS measurements are based on simultaneous use of two GPS devices: a base station and a rover. The relative accuracy between both devices is usually better than 5 cm, in the horizontal and vertical (Immerzeel and others, Reference Immerzeel2014). During the May 2014 field campaign, 13 cliff outlines were mapped by dGPS on Lirung Glacier (Fig. 1); in October 2014 two cliffs were fully mapped and three were partially mapped (i.e. only the upper edge was mapped; Fig. 1). To map a cliff outline, one dGPS point was recorded every 5 m for all debris/ice margins. The cliff edge was measured as accurately as possible, but it was not always possible to map the exact edge, it being too hazardous due to unstable terrain and rockfall. In cases where a supraglacial pond was present at the foot of the cliff, the cliff bottom was not surveyed.
The dGPS points have two main uses: to reconstruct the cliff surface from its outline and to serve as ground control points (GCPs) for the photogrammetric models.
3.2. High-resolution topographic data
The topographic data were utilized as TINs from terrestrial or airborne (UAV) SfM photogrammetry. All the terrestrial and aerial photographs were processed using the commercial Agisoft software (Agisoft, 2013), which uses an SfM workflow (e.g. Westoby and others, Reference Westoby, Brasington, Glasser, Hambrey and Reynolds2012; Immerzeel and others, Reference Immerzeel2014; Passalacqua and others, Reference Passalacqua2015; Kraaijenbrink and others, Reference Kraaijenbrink2016). The first step of SfM workflow is based on keypoint (i.e. matching features) extraction from the images collected in the field. This is achieved by feature recognition algorithms, such as the Scale Invariant Feature Transform method (e.g. Snavely and others, Reference Snavely, Seitz and Szeliski2008). These keypoints allow for the derivation of a first low-density point cloud and constrain the camera pose (so-called bundle block adjustment). The cloud densification is done by a pair-wise depth map computation algorithm (Furukawa and Ponce, Reference Furukawa and Ponce2009). It is finally triangulated to obtain an irregularly triangulated mesh, also called a TIN.
In May and October 2013, the topographic data were obtained by UAV photogrammetry (Immerzeel and others, Reference Immerzeel2014; Kraaijenbrink and others, Reference Kraaijenbrink2016). Despite a high viewing angle, there was enough overlap in the photographs to resolve the overhanging sections of the cliffs on the final TINs. The flight altitude was chosen in order to obtain images with a ground resolution ~4–7 cm (Kraaijenbrink and others, Reference Kraaijenbrink2016). A step-by-step description of the processing workflow is available in Kraaijenbrink and others (Reference Kraaijenbrink2016). The extraction of the cliff surfaces from the TIN representing the whole glacier is based on visual inspection of the TIN. The resulting TINs have an original point density ranging from 2.3 to 7.8 points m−2. The TINs representing cliff 2 and cliff 3 in May 2013 were linearly re-sampled at 10 points m−2 to facilitate the manual delineation of the cliff.
In May and October 2014, the topography of the cliffs was obtained from terrestrial photographs using SfM. We processed between 10 and 200 pictures to generate each cliff TIN. We filtered the sparse point clouds by removing ~15–20 % of the points with the highest reprojection error and reconstruction uncertainty. The TINs were derived from the dense point clouds using a number of triangle elements with the same order of magnitude as the number of points in the cloud. The final TINs from the terrestrial photogrammetry had a much higher point density than the UAV TINs, ranging from 14 to 243 points m−2. We used 4–8 distinguishable features (such as a high point on the cliff, or a cliff corner) as GCPs for each cliff in order to properly georeference the SfM output. The GCPs were chosen among the dGPS points forming the cliff outline. We then manually extracted the cliff outline from the TIN, which we compared with the full set of dGPS observations for validation of the SfM TIN (Fig. 2; Table 2).
Fig. 2. Distance (3-D) between the field dGPS mapped cliff outlines (dots) and the generated TIN cliff outlines (black). The dots are coloured as a function of the distance to the closest point in the TIN outline. All the coordinates are in UTM 45N/WGS84. Note the difference in horizontal scale between the different cliffs. Cliff 4 is poorly aligned for October 2014, partially due to the lack of dGPS points at the bottom of the cliff. dGPS points are very helpful in aligning the 3-D models and constraining the photogrammetric processing.
Table 2. Differences between dGPS and TINs cliff outlines
All the values correspond to the absolute difference between the cliff TIN outline and dGPS outline. Table based on the same data as Figure 2.
3.3. Displacement data
As the avalanche-fed and nearly-flat tongue of Lirung Glacier is disconnected from its steep accumulation zone, the ice flows very slowly. Immerzeel and others (Reference Immerzeel2014) measured surface velocities between ~0 m a−1 (almost stagnant terminus) to 5 m a−1.
Between May 2014 and October 2014, marked rocks near the cliffs were monitored with a dGPS to track the glacier surface displacement. The surface displacement data used in this study are listed in Table 3. All cliff outlines are translated horizontally and vertically down-glacier to the October 2014 base to remove glacier motion from the volume calculations (Fig. 3). The local displacement for each cliff is assumed to be equal to the median of the measured velocities in a 50 m buffer around the cliff outlines, excluding the cliffs themselves.
Fig. 3. Cliff outlines from May 2013 to October 2014. The light coloured outlines are the original cliff outlines and the darker outlines are after correction for the glacier flow according to the October 2014 base. All the coordinates are in UTM 45N/WGS84 system (m). Note the difference in the horizontal scale.
Table 3. Glacier surface displacement corresponding to each cliff
Dx and Dy are increasing towards the east and north, respectively. Dz is the difference in elevation. The data from May–October 2014 are obtained by tracking marked rocks with dGPS. The May–October 2013 data and the October 2013–May 2014 come from feature tracking on UAV orthophotos and DEMs (Immerzeel and others, Reference Immerzeel2014).
4. METHODS
4.1. Underlying principle
Ice cliff retreat can be idealized as in Figure 4, where cliff backwasting between t i and t i+1 is the volume represented by the shaded area (cross-sectional view). In this case, the cliff backwasting can be seen as the volume sandwiched between the cliff surfaces at different times after correction of the cliff position for glacier flow. We developed two methods to calculate this volume: (1) the ‘dGPS method’ triangulates the dGPS outline of the cliffs to produce two interpolated surfaces, and the outlines of the two surfaces are then triangulated to enclose a volume; (2) the ‘TIN method’ is based on the extraction of the cliff triangulated surface (i.e. TIN) obtained by photogrammetry. Then the outlines of the cliff at t 1 are triangulated with those at t 5 to enclose a volume (Figs 5, 6).
Fig. 4. Idealized cross-section views of a retreating cliff between times t 1 and t 5. The cliff backwasting is represented by the grey shaded areas (projections of the backwasting volumes).
Fig. 5. Chart of the workflow.
Fig. 6. Schematic of the principle of the triangulation (a) of a cliff outline (used in the dGPS method) and (b) between two cliff outlines to calculate the melted ice volume. For the sake of readability, the density of points is reduced. In real cases, the points are closer and more evenly distributed along each outline.
The TIN method is expected to provide the most realistic volume estimates, as it does not require assumptions on the cliff surface geometry. Nevertheless, the dGPS method has the advantage that it can be applied to very high-resolution satellite images or airborne data, and thus offers greater potential for future applications at regional scale.
4.2. Surface reconstruction and volume calculation
4.2.1. dGPS method
The main assumption we make to triangulate the surface from the outlines is that the cliff surface is mostly flat so that it is consistent to link the upper edge of the cliff with the bottom edge. As the dGPS points are not equally spaced (especially when parts of the outlines cannot be reached because they are too steep or bordered by a lake), we interpolate along the surveyed cliff edge to obtain cliff outlines of 500 nodes on the upper and lower edges for 1:1 correspondence between outlines. The number of points interpolated between two adjacent dGPS points is a function of the distance between these points normalized by the median distance between adjacent dGPS points along this cliff edge. The upper and lower cliff edges are separated manually. Once the outline is interpolated, we triangulate it according to the rule illustrated in Figure 6a. The point number (i − 1) of the bottom edge (y i−1) is linked with the point number i of the bottom edge (y i ) and with the point number i of the upper edge (x i ). Then x i is linked with y i and x i+1.
4.2.2. TIN method
For the TIN method the cliff surface is delineated and extracted based on visual inspection of the coloured TIN and of the raw photographs (the cliff has different color and texture than the surrounding terrain).
In both cases we obtain a TIN that contains the cliff surface. These TINs are linearly translated down-glacier to the position where they would have been in October 2014 if the glacier was motionless (Fig. 3). We used for this the velocity data of Table 3.
The triangulation of the ice volume loss is similar to the triangulation of the surface from the dGPS outline (Fig. 6b). For the TIN method, the original cliff outline extracted from the TIN is also interpolated to produce a cliff outline of 1000 points.
The final three meshes (cliff at t 1, cliff at t 2 displaced to October 2014 base and the mesh obtained by joining the two cliff outlines) are imported into the open-source software MeshLab (Cignoni and others, Reference Cignoni, Corsini and Ranzuglia2008). In Meshlab, the face normals are oriented outwards and merged into a single mesh, from which the volume is calculated according to Mirtich (Reference Mirtich1996).
4.3. Uncertainty assessment
The main sources of uncertainties in our methods are: (1) the emergence velocity, which is not taken into account, (2) uncertainty of the cliff outline position as determined from the TIN vertices, (3) uncertainty and spatial variability of the glacier surface displacement field and (4) uncertainty of the georeferencing.
It is difficult to assess the uncertainty associated with the first and second sources of errors. The emergence velocity can be estimated by calculating the ice flow through a cross section of the glacier. This can be performed by measuring the ice surface velocity (and assuming a given ice velocity profile) and the ice thickness along the cross section. The surface velocity can be estimated from our data (either from the UAV or the dGPS data), but the ice thickness is unknown. However, the emergence velocity was previously estimated to be ~0.18 m a−1 (Naito and others, Reference Naito1998), which is low compared with the 1.1 m a−1 average lowering of the glacier surface (Immerzeel and others, Reference Immerzeel2014). The discrimination between cliff and non-cliff surfaces is a potential source of error ((2) above), but it is usually straightforward to distinguish between vertical bare ice and the surrounding rocks because they have different colors and the cliff surface is smooth and flat compared with the surrounding environment.
Uncertainty on cliff displacement due to glacier motion ((3) above) is not easy to assess. Our method assumes that the local displacement field is homogeneous within the cliff. This is a reasonable assumption for small cliffs (e.g. cliff 1) but is more questionable for larger cliffs, even though they are mostly perpendicular to the flow. The larger cliffs (almost as large as the glacier width) are very likely to experience some differential flow. The glacier-wide boulder tracking study of (Immerzeel and others, Reference Immerzeel2014) showed that the 147 boulders have a mean displacement (horizontal and vertical) of 1.1 m with a standard deviation of 0.7 m over the study period. Assuming the uncertainty in the displacement is equal to the ratio of the standard deviation to the average, we obtain an uncertainty of 60%. The displacements between May and October, and October and May are, on average 1.4 and 0.4 m, respectively (Table 3; for May–October we use the average over the two periods May–October 2013 and May–October 2014). This leads to uncertainties in displacement, Δdisplacement, of 0.8 m between May and October and 0.2 m between October and May.
The main uncertainty remains the location of the cliff outlines as derived from the SfM workflow and by dGPS measurement ((4) above). The absolute median distance between each point of the cliff dGPS outline and of the closest point of the cliff TIN outline is higher in z (0.6 m) than in x and y (0.4 and 0.4 m, respectively; Table 2). We also take into account the intrinsic uncertainty on dGPS location, which is assumed to be ±0.3 m. The total uncertainty on the cliff location can be calculated as the quadratic sum of these errors, which leads to an uncertainty in outline location, Δoutline, of 0.9 m.
Assuming that both errors in displacement and location are independent, the overall uncertainty in the location is calculated as:
$$\Delta _{{\rm overall}} = \sqrt {\Delta _{{\rm displacement}}^{\rm 2} + \Delta _{{\rm outline}}^{\rm 2}} = \left\{ {\matrix{ {1.2} \hfill & {m\;{\rm between}\;{\rm May}\;{\rm and}\;{\rm October}} \hfill \cr {0.9} \hfill & {m\;{\rm between}\;{\rm October}\;{\rm and}\;{\rm May}} \hfill \cr}} \right.$$
It is difficult to assess whether this overall error in location has preferential directions, so we use it as the upper bound for the worst-case scenario (i.e. we take the error to be exactly in the backwasting direction). Using typical values 6 m of backwasting (D bw) over the monsoon season, and 2 m over the winter, obtained with our method (Table 4), we obtain the following relative uncertainties:
$$\Delta _{{\rm bw}} = \displaystyle{{\Delta _{{\rm overall}}} \over {D_{{\rm bw}}}} = \left\{ {\matrix{ {20\%} \hfill & {{\rm between}\;{\rm May}\;{\rm and}\;{\rm October}} \hfill \cr {45\%} \hfill & {{\rm between}\;{\rm October}\;{\rm and}\;{\rm May}} \hfill \cr}} \right.$$
Table 4. Volume loss obtained for all cliffs over the study period
The bracketed May–October 2013 values for cliff 1 and cliff 2 are the estimates of Buri and others (Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016).
The uncertainty on the cliff area ((2) above) is negligible in comparison, so we represent the maximum uncertainty of volume calculations using Δbw.
5. RESULTS
We calculate a normalized volume loss, called hereafter average cliff backwasting rate, as the volume loss between t 1 and t 2 divided by the average of the cliff areas at t 1 and t 2, converted to cm w.e. assuming an ice density of 900 km m−3. Figure 7 and Table 4 show the backwasting rate calculated with the TIN method for all cliffs and seasons investigated. This normalization is useful for comparing cliff changes within different areas (Fig. 7a), but it may not be appropriate to normalize by the mean area when the area changes drastically between the two dates (for instance cliff 1 and cliff 4’ between May and October 2014; Fig. 7). Cliffs lose mass at much faster rates between May and October (3.8 ± 0.3 cm w.e. d−1 in 2013 and 3.1 ± 0.7 cm w.e. d−1 in 2014) than between October and May (1.0 ± 0.3 cm w.e. d−1; Table 4), in agreement with Steiner and others (Reference Steiner2015). We suspect that a large part of the measured melt between October and May happens just after observations in October or just before observations in April/May, corresponding to the end and start of the melt season in Langtang catchment (Steiner and others, Reference Steiner2015). The consistency in the backwasting rates between the two monsoon seasons suggests similar controls of backwasting. Cliffs with area spanning over two orders of magnitude lose mass at similar rates (Fig. 7). The reduced backwasting rate during winter is expected and can be explained by a reduction in the total energy reaching the cliff surface, especially the incoming shortwave and longwave radiation, which are the two main energy fluxes (Steiner and others, Reference Steiner2015). In winter, the solar elevation angle is lower, reducing incoming shortwave radiation over the entire glacier tongue, and especially affecting the cliff surfaces due to their northerly aspects and steep slopes. The incoming longwave radiation from the surrounding debris is also reduced due to the lower average temperature of the debris surfaces. In addition, refreezing may play a strong role in the reduction of melt in winter (Steiner and others, Reference Steiner2015).
Fig. 7. Area change (a) and backwasting rate (b) of all cliffs as a function of time. The backwasting rate is defined as the volume change between t 1 and t 2 divided by the number of days and the mean cliff area. The average melt rate (grey) is from Immerzeel and others (Reference Immerzeel2014). Note that the y-axis for the upper panel has a logarithmic scale. Shaded areas in both panels indicate the monsoon period (15 June–30 September).
The backwasting rates obtained using the TIN method are compared with the backwasting rates obtained using the dGPS method in Table 4. The data are sufficient to apply the latter only for cliff 1 and cliff 2 between May and October 2014. The cliff outlines were fully mapped for cliff 3, 4 and 4’ in May 2014 as well, but lack of time did not allow the survey to be repeated in October 2014. Cliff 5's bottom edge has never been accessible because of a large supraglacial lake. For cliff 1 the agreement is very good and the volume loss estimated with the dGPS method is <1.5% lower than the TIN method's result. On the other hand for cliff 2 the agreement is not as good and the dGPS method underestimates the volume loss by 9%. This is within the uncertainty of the TIN method and could be due either to the complex geometry of cliff 2 (which exhibits an overhanging face) or the incomplete dGPS outline of cliff 2 (Fig. 2), which leads to an underestimation of the volume loss. As they are, these results are insufficient to reach a conclusion about the dGPS method's accuracy.
6. DISCUSSION
6.1. Method's limitations and comparison with M3C2 method
The method developed in this paper has the potential to provide estimates of cliff backwasting that can be used to validate melt models and/or to provide a direct estimate of cliff contribution to total glacier mass balance. The method requires several steps of manual processing: delineating the cliff outlines, merging the different TINs and calculating the volume of the mesh. This was not a limiting factor in this study, but it will be useful to automate some of these steps to apply the method more systematically.
The TIN and dGPS methods are not suitable for estimating backwasting rates for cases, in which the cliff geometry changes significantly, as for cliff 1, which shrank between May and October 2014 (Fig. 7a). Despite the high confidence in the volume loss estimate, we can obtain a backwasting rate of 0.7 or 3.8 cm w.e. d−1 if we normalize by the cliff area of May 2014 or October 2014, respectively. A similar effect can be observed for cliff 4, which expanded between May and October 2014.
The method is also intrinsically limited by the accuracy of the topographic data and it can only be applied to calculate backwasting rates over periods long enough to measure changes that are significantly higher than the uncertainty of the method.
Another main limitation of the TIN and dGPS methods is the fact that they provide only estimates of total volume loss, and cannot be used to calculate the spatial distribution or patterns of cliff backwasting. A promising alternative method is the M3C2 algorithm, which was developed to measure distances between point clouds (Lague and others, Reference Lague, Brodu and Leroux2013). This method is based on calculation of local normals to the initial surface, which are then used to calculate distances from the other surface. This algorithm was tested without success on our data, which describe the geometry of the cliff surface alone. However, the test was unsuccessful for two reasons. First, as the backwasting is usually not perpendicular to the cliff face, mismatches occur between normals calculated at the cliff surface at t 1, which frequently miss the cliff surface at t 2. Based on our attempts with the algorithm, this results in a loss of up to 82% of points, with the highest efficacy when the backwasting is smallest (winter). Second, interpretation of spatial patterns of backwasting are subjected to high relative georeferencing errors – a slight lateral shift in the position of cliff geometry observed at t 2 may emphasize backwasting at one part of the cliff over another; the volume estimate will not be affected so heavily but the spatial patterns may be very sensitive to such a change. Finally, computation of volumetric change from M3C2 distances is not a trivial undertaking (Lague and others, Reference Lague, Brodu and Leroux2013). As a result, we believe that the method presented here has value as an alternate approach to assess volumetric change of ice cliffs.
However, algorithms for cloud comparison are fast developing and the M3C2 method is very promising for future analyses of the spatial patterns of backwasting, especially for short-return time intervals. Thus, we present here some recommendations for data collection targeting M3C2 use, and highlight the challenges that future applications will face. First, to maximize the likelihood of surface normals from t 1 intercepting the surface at t 2, we suggest a closed-volume approach to terrestrial SfM. This entails measuring the entire ice cliff depression, as well behind the cliffs (rather than the cliff geometry alone), such that the area of volumetric change is fully enclosed by the two point clouds. Second, the georeferencing needs to be more accurate than was performed for this study. Therefore, we recommend georeferencing of the TINs on artificial GCPs (i.e. painted rocks for instance) instead of natural features, for easier identification and stronger constraint in the SfM processing. Third, we recommend more frequent surveys of cliff geometry, studying changes over shorter time intervals than used for this study and maximizing the usable area of the cliff.
6.2. Comparison with a simplified approach
Both the TIN and dGPS methods require collection of substantial field data. It is beneficial to determine the minimum amount of input data required to obtain a satisfying estimate of cliff backwasting rate. To estimate ice cliff backwasting, Han and others (Reference Han, Wang, Wei and Liu2010) measured the distance between the cliff edge and a distinguishable fixed boulder at t 1, and measured it again at t 2 (1 month later in their case). We adapted their method and measured the horizontal distance between the cliff upper edge at t 1 and at t 2 on the ortho-photos. For comparison, we calculated the backwasting distance (i.e. volume loss divided by the cliff average area with no density correction (m)) taking into account the cliff geometry provided by our topographic dataset. Interpolating the two cliff edges to obtain the same number of points along each edge, we took the median distance between these pairs of points as the distance between the edges. After correcting for glacier flow and excluding cliff 4’ from the analysis because its shape changed too much, we compared the horizontal distance between the edges with our estimate of the perpendicular backwasting (called backwasting distance; Fig. 8 – left panel). The horizontal distance is higher than the backwasting distance (mean bias of 2.9 m). If we calculate the horizontal backwasting by projecting the backwasting distance (i.e. by dividing the backwasting distance by the sine of the slope), the mean bias is reduced (0.9 m) but the correlation is slightly worse (R 2 = 0.58 after projection versus 0.66 without correction; Fig. 8). Even if the agreement between the two methods is not fully satisfying, it has the potential to be used widely because cliff edges are visible on very high-resolution satellite imagery. However, a major drawback of this method is that it cannot provide estimates of the volume loss (because the cliff area and slope are a priori unknown) and is therefore of limited interest.
Fig. 8. Comparison of the cliff backwasting measured with our method (x-axis) and the horizontal displacement of the cliff edge (y-axis). The circles correspond to the period May–October 2013, the diamonds to the period October 2013–May 2014 and the squares to the period May–October 2014. The dashed line is the 1:1 line. The backwasting distance is calculated as the volume loss divided by the cliff area. Note that cliff 4’ is excluded from the analysis. The x error bars are calculated assuming 15% uncertainty in the volume loss results and the y error bars correspond to 1 m of error in the edge backwasting observation.
6.3. Comparison with model results
The present study derives mass losses from cliff backwasting that can be used to validate models. Buri and others (Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016) developed a distributed energy balance model that calculates melt for each grid cell of a high-resolution DEM of the cliff. They applied the model to cliff 1 and 2 of Lirung Glacier for 8 May–23 October 2013. The model was forced with data from an AWS located on-glacier. Buri and others (Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016) ran the model on a 20 cm grid extracted from an UAV DEM (Immerzeel and others, Reference Immerzeel2014; Buri and others, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016). Buri and others (Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016) calculated average melt rates of 3.6 and 3.1 cm w.e. d−1 for cliff 1 and 2 respectively, results consistent with ours: 3.8 ± 0.8 and 3.7 ± 0.7 cm w.e. d−1, respectively (Table 4). The volume losses calculated with the TIN method are slightly higher than the model estimates (even if the uncertainties on both are still large). This can be explained by the fact that Buri and others (Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016) calculated only the volume losses associated with melt due to energy exchange with the atmosphere, whereas other processes, such as the contribution to radiative fluxes of a lake or thermoerosion at the base of the cliff, may also contribute to cliff backwasting.
6.4. Insight into mechanisms of cliff area changes
The repeated high-resolution topographic surveys are very useful for documenting cliff evolution. The expansion of cliff 4, with the appearance of a new cliff on its western side (cliff 4’), or the shrinkage of cliff 1 are striking behaviours, in contrast to cliff 2 and 5, which maintain a constant shape (Figs 3, 9). The most interesting change observed in this study is the shrinkage of cliff 1, which nearly disappeared between May and October 2014 (Figs 3, 7, 9).
Fig. 9. Evolution of the cliffs between May 2013 and October 2014. The X, Y and Z coordinates are in meters with arbitrary origins. X and Y correspond to the easting and northing, respectively. All the TINs from which these figures are produced are shifted to the October 2014 base to correct for the ice flow.
A first possible mechanism responsible for the disappearance of cliff 1 may be the absence of a lake at its foot. The lake contributes to the ice cliff radiation budget and maintains the steep slope at the lower section of the cliff through thermoerosional undercutting (Miles and others, Reference Miles2016). Field evidence suggests that a pond was present during May–October 2013, which then disappeared, triggering the substantial changes in cliff 1 from May to October 2014.
A second possible mechanism may be the limiting volume of the topography behind the cliff. We analyzed the surrounding terrain within a 2 m buffer around the upper cliff edge. For cliff 1, we observed that a small ridge, parallel to the cliff edge and located ~1 m from the edge in the north-eastern side of the cliff disappeared between October 2013 and May 2014. It is likely that its disappearance limited the amount of ice available for melt.
These process-oriented hypotheses remain qualitative and speculative, but our method has produced one of the first quantitative datasets of cliff evolution over several seasons, showing interesting phenomena of growth and decay that should be incorporated into advanced numerical models of cliff evolution.
7. CONCLUSIONS
In this study we have developed a method based on high-resolution topographic data to assess volume loss due to ice cliff backwasting, and applied the method to selected ice cliffs of Lirung Glacier (Langtang Valley, Nepal). Our results confirm the significant contribution of ice cliffs to total mass loss, with backwasting rates of 3.5 ± 0.7 cm d−1 over the monsoon. This rate is 6 times higher than the average glacier ablation rate over the same period (Immerzeel and others, Reference Immerzeel2014), thus substantiating a potential explanation for the so-called debris-covered anomaly. This is, to our knowledge, the first accurate, high-resolution estimate of ice cliff backwasting from ground data, and is important as a confirmation of the results of the few existing modelling studies (Han and others, Reference Han, Wang, Wei and Liu2010; Reid and Brock, Reference Reid and Brock2014; Steiner and others, Reference Steiner2015; Buri and others, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016).
We have provided the first time series of cliff 3-D changes, encompassing a period of two monsoon seasons and one winter season. This time series shows geometric developments that are not uniform, with some cliffs exhibiting few changes and a striking persistence over the two monsoon seasons, while others grow or shrink substantially. These variable observations point to the heterogeneity of the processes controlling cliff evolution and dynamics, even on a relatively small glacier such as Lirung Glacier, and to the need for process-based modelling studies of cliff developments.
While there is some evidence that cliffs form preferentially on stagnant tongues, we still lack an understanding of their distribution, formation and future change. The datasets in this study can provide the necessary basis for numerical studies of cliff evolution.
In addition, the existing energy balance models of cliff backwasting are still affected by parametric uncertainties and errors due to the poorly understood variability of meteorological forcing at the cliff's local scale (Steiner and others, Reference Steiner2015; Buri and others, Reference Buri, Pellicciotti, Steiner, Miles and Immerzeel2016). Availability of accurate, high-resolution datasets of 3-D cliff backwasting such as presented here can help to significantly constrain these uncertainties and inform model development. While assessment of ice cliff volume changes with this method is not practical at the scale of large glaciers or basins, more datasets of the type presented in this paper can provide critical calibration and validation data for numerical model development, and be used to test simplified methods of cliff volume change estimates based on satellite images. These data are now easier to collect and process thanks to the recent advances in UAVs and photogrammetric software.
ACKNOWLEDGEMENTS
We thank Tek Rai, Ibai Rico, Anna Chesnokova and Bhanu Thakur for help with dGPS and photogrammetric data collection. All the participants in the four field campaigns made this work possible and are gratefully acknowledged, as all are the Nepali helpers who supported our work in a great way. We thank Joe Shea and the International Centre for Integrated Mountain Development (ICIMOD) for administrative and logistical support in Kathmandu. We thank Dimitri Lague and Matt Westoby for precious advice on M3C2 algorithm use. This study is funded by the Swiss National Science Foundation (SNF) project UNCOMUN (Understanding Contrasts in High Mountain Hydrology in Asia) and by the USAID High Mountain Glacier Watershed Programs Climber-Scientist (Grant No. CCRDCS0010), which are warmly acknowledged for their support. E. B. acknowledges support from the French Space Agency (CNES) through the TOSCA and ISIS programs. We thank Duncan Quincey and an anonymous reviewer who made very helpful comments on an earlier version of this manuscript, which greatly improved the final paper. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 676819).
Open access
