Scholarly article on topic 'Seasonal surface velocities of a Himalayan glacier derived by automated correlation of unmanned aerial vehicle imagery'

Seasonal surface velocities of a Himalayan glacier derived by automated correlation of unmanned aerial vehicle imagery Academic research paper on "Earth and related environmental sciences"

Share paper
Academic journal
Annals of Glaciology

Academic research paper on topic "Seasonal surface velocities of a Himalayan glacier derived by automated correlation of unmanned aerial vehicle imagery"

Annals of Glaciology 57(71) 2016 doi: 10.3189/2016AoG71A072 103

© The Author(s) 2016. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons. org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.

Seasonal surface velocities of a Himalayan glacier derived by automated correlation of unmanned aerial vehicle imagery

Philip KRAAIJENBRINK,1 Sander W. MEIJER,1 Joseph M. SHEA,2 Francesca PELLICCIOTTI,3 Steven M. DE JONG,1 Walter W. IMMERZEEL1-2

1 Department of Physical Geography, Utrecht University, Utrecht, The Netherlands 2International Centre for Integrated Mountain Development, Kathmandu, Nepal 3Institute of Environmental Engineering, ETH Zürich, Zürich, Switzerland Correspondence: Philip Kraaijenbrink <>

ABSTRACT. Debris-covered glaciers play an important role in the high-altitude water cycle in the Himalaya, yet their dynamics are poorly understood, partly because of the difficult fieldwork conditions. In this study we therefore deploy an unmanned aerial vehicle (UAV) three times (May 2013, October 2013 and May 2014) over the debris-covered Lirung Glacier in Nepal. The acquired data are processed into orthomosaics and elevation models by a Structure from Motion workflow, and seasonal surface velocity is derived using frequency cross-correlation. In order to obtain optimal surface velocity products, the effects of different input data and correlator configurations are evaluated, which reveals that the orthomosaic as input paired with moderate correlator settings provides the best results. The glacier has considerable spatial and seasonal differences in surface velocity, with maximum summer and winter velocities 6 and 2.5 m a-1, respectively, in the upper part of the tongue, while the lower part is nearly stagnant. It is hypothesized that the higher velocities during summer are caused by basal sliding due to increased lubrication of the bed. We conclude that UAVs have great potential to quantify seasonal and annual variations in flow and can help to further our understanding of debris-covered glaciers.

KEYWORDS: debris-covered glaciers, glacier flow, glacier mapping, glaciological instruments and methods, remote sensing


Himalayan glaciers play a varying, but generally important, role in the water supply of many regions in Asia (Immerzeel and others, 2010; Kaser and others, 2010; Lutz and others, 2014). Most glaciers in High Mountain Asia are losing mass at rates similar to other regions in the world, except for the Karakoram mountain range, where there are indications of positive mass balances (Bolch and others, 2012; Gardelle and others, 2012). In the central Himalaya, for example, negative mass balances of — 0.26 ± 0.13 mw.e. a-1 for the Everest region and of —0.32 ± 0.13 mw.e. a-1 for west Nepal are reported for the period 1999-2011 (Gardelle and others, 2013), whereas for the Langtang catchment in central Nepal a mass balance of —0.33 ± 0.18 m w.e. a-1 is reported (Pellicciotti and others, 2015). These negative mass balances temporarily result in higher water availability, until the glaciers recede so far that absolute meltwater yield starts to decline (Immerzeel and others, 2013).

Around 10% of the Himalayan glacierized area is debris-covered (Bolch and others, 2012) and the debris-covered tongues are generally located at the lowest elevation. Most debris-covered tongues exhibit slower rates of retreat than debris-free glaciers, but they thin at substantial rates (Scherler and others, 2011). Theoretically, the debris, when thicker than a few centimetres, should insulate the ice from melt (Ostrem, 1959). However, recent work suggests that the debris-covered tongues lose mass at the same rates as debris-free glaciers (Kaab and others, 2012; Gardelle and others, 2013; Pellicciotti and others, 2015). The underlying reason may be the presence of supraglacial lakes and ice

cliffs that accelerate melt significantly (Sakai and others, 1998; Benn and others, 2012; Immerzeel and others, 2014). Little is known, however, about the behaviour and response of debris-covered glaciers, as they are generally inaccessible and the spatial and temporal resolution of satellite remote-sensing products limits our ability to understand the processes governing thinning.

The flow velocity and the associated mass turnover determine, to a large extent, the sensitivity of a glacier to climate change. Recent work has shown that many of the mountain glaciers are slowing considerably. Glaciers in the Pamir, for example, slowed by 43% between 2000 and 2010 (Heid and Kaab, 2012). A 70% reduction in flow velocity was reported for Yala Glacier in the Langtang catchment in Nepal between 1982 and 2009 (Sugiyama and others, 2013). In contrast, Karakoram glaciers again exhibit anomalous behaviour, as glaciers there generally have accelerated (Heid and Kaab, 2012). Recent work in the central Himalaya shows great variation in surface velocities. To the north side of the Himalayan arc, near Bhutan on the Tibetan Plateau, flow velocities of 100-200ma-1 are reported, whereas on the south side maximum flow velocities of a few tens of metres are reported (Kaab, 2005). This is confirmed for the south side of the Everest region, where flow velocities for debris-covered tongues vary from 0 to 37 ma-1 (Quincey and others, 2009a). Most of the glacier velocity studies in the Himalaya are based on optical, spaceborne satellite remote sensing and feature tracking. The first automated approach applied to glaciers was published >20 years ago (Scambos and others, 1992), and over the years the approach has

Fig. 1. Location of the study area (top left), view of Lirung Glacier from across valley (bottom left), and UAV-derived orthomosaic (middle) and DEM (right) for May 2014.

proved very powerful and reasonably accurate (Scherler and others, 2008; Copland and others, 2009; Quincey and others, 2009a,b).

Unmanned aerial vehicles (UAVs) have great potential in glaciology, in particular for debris-covered glaciers, as was shown in a recent study of Lirung Glacier in central Nepal (Immerzeel and others, 2014). The study revealed a highly heterogeneous pattern of mass loss on the debris-covered tongue over a single monsoon season, with a possibly important catalytic role for supraglacial lakes and ice cliffs. Additionally, the study showed that it is possible to determine the glacier's surface velocity and its general spatial pattern by manual digitization and interpolation of the displacements found between UAV image pairs. Such a digitization method, however, does not optimally use the full information content present in the UAV data, is subject to human error and is time-consuming. To overcome these issues, an automated feature-tracking approach could be applied to the high-resolution UAV imagery. This would result in surface velocity products with better accuracy and spatial resolution, that may achieve a level of detail that is currently unobtainable with spaceborne remote sensing.

In this study we derive surface velocities of Lirung Glacier in central Nepal for both the summer monsoon and the dry winter season by applying frequency cross-correlation algorithms (Leprince and others, 2007) on UAV-acquired imagery of May 2013, October 2013 and May 2014. Our objectives are twofold. First, we evaluate the effect of two high-resolution image products on the correlation output, i.e. digital elevation models (DEMs) and orthomosaics, and different settings for the cross-correlation algorithm. Based on the best-performing configuration we then assess differences between monsoon and winter velocities, and discuss implications for debris-covered glacier dynamics. The two 2013 datasets have already been presented (Immerzeel and others, 2014), but are reanalysed in this study, using the frequency cross-correlation technique to improve the detail of the surface velocity product and to increase its comparability with the 2014 data.


Lirung Glacier is located in Langtang National Park in the Nepalese Himalaya (Fig. 1). It is part of the Langtang catchment and it has a typical summer monsoon from June to September, in which most of the annual precipitation of ~800mm occurs. The glacier is characterized by a debris-covered tongue, 500 m wide and 3000 m long. The terminus of the glacier tongue lies at an elevation of ~4050 m above mean sea level (a.m.s.l.) and the remainder of the tongue slopes up to ~4350 m a.m.s.l.

The glacier tongue is detached from the steep accumulation slopes below Lantang Lirung Peak (7235 m a.m.s.l.) and it is currently fed only by avalanches and occasional snowfall on the tongue itself. The heterogeneous pattern of surface lowering found over the monsoon season was 1.09 m on average (Immerzeel and others, 2014), which is comparable to other debris-covered glaciers in these parts of the Himalaya (Bolch and others, 2011; Kaab and others, 2012). The ice cliffs present on the glacier appear more dynamic, with reported melt rates of up to ~8cmd-1 (Sakai and others, 1998; Buri and others, 2016; Steinerand others, in press). Surface velocities determined by manual feature tracking are ~2.5m over the monsoon season in the glacier's upstream area (Immerzeel and others, 2014), i.e. 5.8 ma-1, relatively low compared with other findings (Quincey and others, 2009b).


Lirung Glacier was surveyed by UAV three times, on 18 May 2013, 22 October 2013 and 1 May 2014. A Swinglet CAM UAV from the company senseFly (SenseFly, 2015) was used in 2013. In May 2014 an eBee from the same company was used. The months May and October were chosen as ideal survey and fieldwork conditions usually prevail, i.e. calm winds, moderate temperatures and little or no precipitation. Also, they are just before and after the monsoon, which is

Fig. 2. Overview of the 11 UAV flights over the three survey periods, their approximate ground coverage, positions of the gathered images selected for processing, locations of the GCP and locations of the tie points.

the most dynamic season when synchronous accumulation (at high altitudes) and ablation (on the tongue) processes prevail (Immerzeel and others, 2013). These months are therefore optimal for understanding the difference in behaviour of the glacier between the monsoon and the remainder of the year, which is generally colder and drier. The monsoon and dry seasons covered by the three datasets are hereafter referred to as, respectively, summer (May-October 2013; 157 days) and winter (October 2013-May 2014; 191 days).

To obtain the required imagery the UAV was deployed over the glacier in 11 separate flights (Fig. 2) over the course of the three survey dates. A total of 284, 307 and 314 usable JPEG images (May 2013, October 2013 and May 2014, respectively), with ground resolutions of 4-7cm and sufficient overlap of ~60% were acquired with a 16 megapixel consumer-grade digital camera, using a fixed focal length. Images that were either redundant, had too much motion blur or strong rolling shutter distortions were removed from the image set. Although the lossy compression associated with the JPEG format is not ideal for data analysis and consistent results, it is currently a limitation of the UAV system in use.

During the October 2013 field campaign a total of 19 ground-control points (GCPs) were collected on Lirung Glacier's lateral moraines, using differential GPS to geo-reference the imagery. During the other two campaigns no GCPs of sufficient quality were collected. It was therefore decided to tie the data together geodetically, by sampling 47 tie points from the October 2013 data, which were used as GCPs in the processing of the May 2013 and 2014 datasets (similar to the approach taken by Immerzeel and others, 2014).

UAV data processing

For each of the three dates, the UAV-acquired images were processed using a Structure from Motion (SfM) workflow (Westoby and others, 2012; Lucieer and others, 2013; Immerzeel and others, 2014). In the workflow, feature recognition and matching algorithms, together with multiview stereo techniques (Szeliski, 2010; Westoby and others, 2012), are applied to the overlapping input images to obtain per-image depth maps and camera orientations. This information is used to construct three-dimensional (3-D) point clouds that can be triangulated and interpolated into gridded DEMs and to stitch the input imagery into geometrically corrected image mosaics, called orthomosaics. By marking the measured GCPs and/or tie points on the input images during the SfM workflow, xyz-georeferencing of the output is obtained. In this study we use the SfM workflow as implemented in the software package Agisoft Photoscan Professional version 0.9.1 (Agisoft, 2013).

To obtain optimal results from the SfM workflow in Agisoft, each processing step was performed using high-quality settings. The 3-D point clouds were cleaned in a three-step iterative process, using the point reprojection error. High reprojection errors indicate poor localization accuracy of the corresponding point projections at the point-matching step and are also typical for false matches (Agisoft, 2013). Points with a reprojection error >1.5 pixels (i.e. ~10cm for most input images) were therefore removed at each iteration. After removal the point coordinates and camera calibrations were optimized each time by minimizing the sum of the reprojection error (Agisoft, 2013). If camera calibration estimates are inaccurate, the SfM matching algorithms can introduce a doming or bowl effect in the output 3-D model. This was counteracted using

spatially well-distributed GCPs and tie points during the optimizations. The output orthomosaics and DEMs, which have 0.1 m and 0.15-0.2 m resolutions, respectively, were all resampled to a 0.2 m resolution to reduce the effects of any remaining motion blur as well as JPEG artefacts for further processing.

During the May 2014 field campaign the UAV had difficulty acquiring GPS fixes, causing the UAV to skip a number of image captures. Consequently, a handful of tie points could not be placed during the SfM workflow and the quality of the output was reduced considerably. As a solution some off-glacier images of static areas from the October 2013 set (Fig. 2) were added to the May 2014 image set during processing, which allowed placement of all but two of the tie points.

An indication of the horizontal accuracy of the DEMs and the orthomosaics is obtained by measuring the difference between the GCP or tie-point coordinates and their positions on the output orthomosaics. The vertical accuracy is determined by calculating the differences between the GCP or tie-point elevations and the output DEMs, while correcting for the horizontal error. The number of GCPs collected in October 2013 was limited because of the inaccessible terrain and all points were required for the SfM processing. This led to the absence of redundant GCPs that could be used for independent accuracy checks.

Surface velocity determination

Cross-correlation feature tracking

COSI-Corr (Co-registration of Optically Sensed Images and Correlation) is a software tool developed to co-register pairs of satellite images, perform orthorectification and also subpixel automated image correlation (Leprince and others, 2007; Ayoub and others, 2009). Its correlation algorithms are used for the determination of surface velocity of glaciers, but until now they have only been applied to comparatively much coarser resolution satellite imagery (Leprince and others, 2008; Scherler and others, 2008; Herman and others, 2011). Here we apply COSI-Corr's correlation algorithms to the high-resolution UAV data.

The software provides two ways to correlate images, either statistical or frequential. Both act on a moving window level. It is advisable to use the frequential correlation method when performing feature tracking on optical images that are relatively noise-free and the statistical method when images have considerable amounts of noise or when image pairs have different content, such as when correlating an orthomosaic with an elevation model (Ayoub and others, 2009). As we have relatively noise-free data, frequency correlation is used. The correlator obtains image-to-image displacements by determining phase differences between Fourier transforms of the moving window of both images. It does this in a twofold process, first roughly at the pixel level and subsequently at a sub-pixel level (Leprince and others, 2007).

Multi-scale windows

Lirung Glacier, besides its general ice flow, has considerable temporal surface variations that are unrelated to the flow of the ice but are clearly noticeable in the high-resolution UAV data. Examples are the melt of ice cliffs, tumbling of boulders and collapse of debris slopes. Ideally these features are not detected by the correlation algorithm, as the aim is to extract surface velocities only. It was therefore decided to

use the COSI-Corr frequency correlator's multi-scale mode (Ayoub and others, 2009), as it has the potential to filter out these unwanted disturbances.

In the multi-scale mode, windows of decreasing sizes are correlated iteratively, using a preconfigured initial and final window size. The dominant displacement is first detected by correlation at the initial window scale. Increasingly smaller windows are then used while accounting for the dominant signals that were found. If a correlation at a current iteration deviates too much from the previous one, the iteration is stopped and the previous window's results are used.

The multi-scale mode decreases the amount of irregularly distributed noise in the output (Ayoub and others, 2009), i.e. small displacements present in the images that are unrelated to the dominant signal. The use of larger initial window sizes allows for the reduction of more noise. However, it is a trade-off, as too large initial window sizes may result in loss of detail that is relevant. Some of this detail can be retained by using a smaller final window size, but the use of smaller final windows introduces more uniformly distributed noise (Ayoub and others, 2009). A correct balance of the settings with respect to the input data is therefore of key importance to obtain the best results.

Input and setting assessment

Most studies use optical imagery as input data for automated feature tracking (Kaab, 2005; Scherler and others, 2008; Copland and others, 2009). However, a correlation algorithm has been applied successfully to a UAV-derived hillshade (Lucieer and others, 2013). In order to achieve optimal results we therefore first assess the use of three different input data types: the orthomosaic, a hillshade and the DEM processed by the Sobel edge-detection operator (Szeliski, 2010). COSI-Corr requires a single-band raster as input and it was decided to use the orthomosaic's red band, as its longer wavelength experiences less influence from atmospheric scattering (Lillesand and others, 2003). The hillshade was created using a solar azimuth of 120° and a zenith angle of 45°. The edge-detected DEM is added, as it accentuates the outlines of the boulders that are abundantly present on the glacier. It presents a strong contrast that may be picked up by the correlation algorithm. All assessments are performed using the summer dataset only, as the expected higher flow velocities in this season will allow for a better evaluation of possible differences in the correlation output. Initial and final window size settings are held equal for each input data type. Their optimal settings are determined by trial and error, keeping in mind the suggestion to work with window sizes that are at least five times the expected displacement (Leprince and others, 2007), which is ~3 m (15 pixels) in this case (Immerzeel and others, 2014).

To assess the effects of varying window sizes the frequency correlator is applied to the input dataset with the most satisfying correlation results by testing various combinations of initial and final window size. The start window sizes are varied so that they yield increasing levels of irregularly distributed noise reduction and detail retention with respect to the input data resolution. The final window sizes are always chosen to have a good balance, visually determined, between detail and noise. After the evaluation of the effects of different input data and correlation settings, the single optimal configuration found for the summer period is applied to the winter data as well. To be able to

Fig. 3. (a) Box plots of the horizontal errors between GCP (October 2013) and tie-point (May 2013 and 2014) coordinates and their positions on the orthomosaics. (b) Histogram of displacements at static off-glacier areas (0.18 km2), as calculated by frequency cross-correlation using W256 F64.

calculated by an ordinary kriging approach (Davis, 2002), applied to the values on the perimeter of the patches.

Correlation accuracy

A proper assessment of the accuracy of the frequency cross-correlation algorithms is difficult, as there are no quality reference data available. To estimate accuracy we assess whether the algorithm can reproduce surface velocities that are derived by a manual digitization. A comprehensive manual image matching for the summer period is performed by digitizing flow vectors on the image pair visually in a GIS. Only surface features are selected for matching that encounter no displacements that are unrelated to the flow of the ice, as determined by expert opinion. The differences between the two methods are assessed by plotting the digitized data against the correlation output value for each window setting, which is sampled from COSI-Corr's gridded velocity field output at the coordinates of the digitized vectors' origin. Linear regression models are fitted to the data to quantify the velocity differences.

To get another measure of the possible errors involved, horizontal displacements found by the correlator for a static off-glacier area of 0.18 km2 are evaluated and compared with the errors of the SfM output. Signal-to-noise ratios (SNRs) provided by COSI-Corr (Leprince and others, 2007) are evaluated for the three different window settings as another indicator of the quality of the frequency cross-correlation.

compare the surface velocities measured over summer and winter, the values are scaled to year-equivalent values (ma-1) throughout this paper.

COSI-Corr analyses the input data using a moving window that has the preconfigured initial window size. The moving window samples the input data at a configurable spatial interval, called the step size. This step size can be set to be both smaller and larger than the size of the moving window itself. Note that the chosen setting for the step size also predetermines the output pixel size of the velocity field, i.e. input resolution multiplied by the step size. For consistency, this parameter is held constant at a value that provides good output resolution while limiting noise in the output and the processing time required.

Post-correlation noise reduction

The use of optimal correlator settings will result in a good balance between noise and detail, but, to a certain extent, noise will persist after the cross-correlation procedure. This comprises both Gaussian noise and some remaining non-normally distributed noise. To further improve the final surface velocity products for summer and winter, two separate noise-reduction methods are applied to the velocity fields. First, Gaussian noise is targeted using non-local means filtering (Buades and Coll, 2005), as implemented in COSI-Corr (Ayoub and others, 2009). The algorithm is applied using moderate noise reduction settings that are able to reduce most noise without having too much of a smoothing effect, determined by visual inspection of the output. Patches of irregular noise are subsequently targeted, by removing velocity values above a threshold if the focal standard deviation is high. The threshold values for the velocity and the focal standard deviation are determined by trial and error. The noise is replaced by values that are


Figure 3a shows the horizontal errors of the UAV products obtained by SfM processing. It shows the differences found between the GCP or tie-point coordinates and their positions on the output orthomosaics for the three periods. Only the errors found for the May datasets are reflected in the accuracy of the derived surface velocities, as they are directly georeferenced to the October 2013 dataset using the tie points. The accuracy of the surface velocity products is not affected by the true geodetic accuracy of the data, which is indicated by the GCP errors for October 2013. The errors found for both May datasets have a similar distribution and range. About 75% of the tie points are located on the ortho-mosaic within 0.2 m of their original position, with only a few outliers that go up to 0.6 m. Errors found further away from the tie points on the off-glacier moraine area and on the glacier surface itself are assumed to be similar due to the high density of tie points used. The bulk of the vertical errors at the tie points are within 50 cm, and 75% are even within ^25 cm. The vertical errors, however, do not contribute much to the accuracy of the surface velocity product determined by feature tracking, as they have little to no influence on the orthomosaic, hillshade and edge-detected DEM.

Correlation assessment

Optical vs DEM derivatives

COSI-Corr's frequency correlator is applied to the UAV-derived orthomosaic, hillshade and edge-detected DEM using an initial window size of 128 pixels (px) and a final window size of 64 px (coded as W128 F64). It was found that a step size of 16 px provides a good balance between output detail, noise and processing times while working with the 0.2 m resampled input data. Note that the output

Fig. 4. Surface velocity results obtained by frequency cross-correlation of three different input data types for the summer period: the orthomosaic, hillshade and edge-detected DEM. Every input image product consisted of a 0.2 m resampled raster, and was processed using an initial window size of 128 px and a final window size of 64 px. Note that the vectors are not linearly scaled and that vectors with a magnitude of <1.4 ma-1 are not displayed.

resolution of COSI-Corr's velocity field is consequently 3.2 m. The resulting surface velocities are shown in Figure 4. The vectors plotted on the map denote the detected flow direction. Vectors that have a magnitude of less than the maximum horizontal error of ~0.6m (Fig. 3), i.e. 1.4 ma-1, are left out.

The general pattern of flow velocity and direction that is detected by the correlation algorithm is similar for each type of input data. A noticeable difference, however, is that irregularly distributed noise is abundant in the correlated hillshade and edge-detected DEM, while the orthomosaic reveals this type of noise almost exclusively at and around the ice cliffs. This higher noise abundance is possibly because the hillshade and edge-detected data both contain similar recurring patterns of crests and edges that result in mismatches of the correlator. The edge-detected DEM, while showing slightly less irregular noise than the hillshade, does experience more erratic variation in flow directions. This is probably because the edge-detection filter also enhances the pattern of tiny edges that results from triangulation of the point cloud, which spatially varies independently from image to image. Also, the very strong contrast and lack of clear gradients may play a role.

The noise near the ice cliffs in the correlated orthomosaic arises because the cliffs on Lirung Glacier are generally larger (Immerzeel and others, 2014) than the initial window size of 128 px, which equates to 25.6 m in case of the 0.2 m resampled input data. As a result, the melt of the ice cliff will be the local dominant signal picked up by the moving window. In terms of noise filtering, the opposite of the desired filter effect now occurs, as the cliff-unrelated displacement is filtered locally (e.g. the flow). Additionally, other mismatches might be introduced near the ice cliffs, as they are not merely displaced features, but represent an actual change in their appearance and shape due to melt (Immerzeel and others, 2014). Most other displacements that are unrelated to flow have been filtered out at W128

F64, except for some slope anomalies on the lateral moraine walls.

Window size assessment

The results show that use of the orthomosaic as input to the frequency correlator yields the best, noise-free output to be used for determination of the surface velocities. We therefore assess the effects of varying the window sizes by using this input data type.

As ice-cliff-related noise persists at a window size setting of W128 F64, it is chosen to assess changes in noise level and detail retention by using two larger initial windows, i.e. W256 F64 and W512 F128 (Fig. 5). To limit the amount of uniformly distributed noise in the output, a final window size of 128 px was chosen for the correlation with the initial window size of 512 px.

The pixel values for the surface velocity are very similar between the three settings. Excluding noise and outliers that are >6 m a-1, the averages found over the entire area for the different window settings (small to large) are 1.62, 1.59 and 1.57 ma-1. Furthermore, 90% of the pixel-by-pixel differences between W512 F128 and W128 F64 are within -0.56 and 0.29 ma-1 and 75% are even within —8.78 x 10~2 and 6.59 x 10_3ma-1. In terms of SNR, larger window sizes yield more pixels that are reported to have little-to-no correlation. The percentage of pixels reported to have a SNR of <0.75, i.e. little correlation quality, are 7.09, 11.15 and 13.80% (small to large windows).

As shown, larger initial window sizes are capable of reducing most cliff-related noise present in the output. However, they introduce sharper, unrealistic boundaries between areas with contrasting velocities. Additionally, much of the finer spatial variability that is present in the W128 F64 results is lost at W512 F128. To balance noise levels, artefact presence and measured correlation performance, the results from W256 F64 are chosen as the optimal configuration.

Fig. 5. Summer period frequency cross-correlation results for three different window size settings applied to the 0.2 m resampled orthomosaic. The initial and final window sizes used in each case are denoted by W and F respectively. Note that the vectors are not linearly scaled and that vectors with a magnitude of <1.4 ma-1 are not displayed.

To further improve the W256 F64 output, moderate nonlocal means filtering is effective regarding Gaussian noise. Almost all noise is removed, while the detail is largely retained. Replacement of correlated velocities with interpolated values proves to be effective in removing any remaining patches of irregular noise near the ice cliffs and terminus (Fig. 7). For the summer dataset, the use of a threshold of 7.5 ma-1 on the velocity values if the focal standard deviation of a 9 px by 9 px window is larger than the 95th percentile of all focal standard deviations shows good results. A threshold of 5 ma-1 while using the same settings for the focal standard deviation suffices for the winter velocity product.

Figure 3b shows a histogram of the displacements over a static off-glacier area of 0.18 km2, as calculated by COSI-Corr using W256 F64. The off-glacier displacements have a mean of 0.14 m, and 95% of the values are <0.34 m, which are acceptable errors. Note that the errors of the SfM output (Fig. 3a) are, besides the cross-correlation error, also reflected by these displacements and that the actual frequency cross-correlation output errors thus are probably smaller than the values shown in the histogram.

Digitized vs correlated flow

To evaluate and compare the results of the summer correlation output, visual interpretation and manual digitization were performed on the summer image pair by manually digitizing 459 vectors in a GIS. Note that this was performed on the original orthomosaic of 0.1 m resolution. A higher sampling density inareas with higher surface velocities was used to obtain more detail there, but generally the measurements are well distributed over the glacier's surface. Figure 6 shows a scatter plot of the correlation output for the three different window sizes against the manually digitized flow, as well as the spatial distribution of the point measurements. Extreme and unrealistic outliers that are due to noise in the correlation output are removed from the scatter plot, i.e. values >8 ma-1

(n = 6). The results of linear regression models that were fitted to the filtered results are shown in the inset table.

Reproduction of the manually digitized flow by the frequency correlator is very good and the overall flow pattern found is similar. Points scatter close to the 1:1 line with R2 values of 0.83 to 0.90 and with relatively small root-mean-square errors of ~0.6 m a-1 over the observed period. Mean velocity errors between the two methods are ~0.10-0.15 m. The slopes of the fitted model indicate a slight underestimation of surface velocity by the correlator, as compared with the manual digitization. The cause for this is unclear and it has not yet been possible to attribute this specifically to one of the methods.

The flow directions that follow from correlation and digitization have, similarly to the velocities, the same overall trend and only slight differences between the two methods are found locally. Compared with W256 F64, half of the digitized vectors have differences in bearing of <8.0° and 75% are within 19.6°.

Note that small differences in velocity and direction are expected here for two main reasons. Firstly, the digitized surface velocities at a point scale are compared with those that are measured for blocks of 16 px by a correlator that bases itself on windows of 128, 256 and 512 px. Secondly, manual digitization is not always completely accurate. Differences in lighting conditions can cause the small surface features on the glacier used for digitization to appear quite differently from image to image. It is estimated that the visual pixel-matching errors may be as large as 24 px on the 0.1 m resolution orthomosaic. This is equal to ~0.2-0.4m over the summer period, i.e. 0.46-0.93 m a-1.

Seasonal surface velocities

Our finding of different flow velocities in summer and winter is notable, as velocity measurements on debris-covered glaciers are rare (Quincey and others, 2009a). While it is possible that short-term and unmeasured variations in velocity may contribute to the overall

1-1-1-1-1-1-1-1- -1-1-1-1-1-

0 1 2 3 4 5 6 7012345 6<

Digitized surface velocity (m a 1)

Fig. 6. Manually digitized surface velocity measurements (n = 453) plotted against the three different frequency cross-correlation outputs for the summer period (left) and the locations of the measurements plotted over an interpolated surface obtained by ordinary kriging (right). The results of fitted linear models and mean velocity errors (MVE) are shown in the inset table.

differences, the seasonal patterns point towards differences in flow regimes. During the summer period, surface velocities of Lirung Glacier range from completely stagnant in the lower (southern) sections of the terminus to ~6ma-1 in the upper (northernmost) surveyed area (Fig. 7). Surface velocities decrease gradually down-glacier to ~2ma-1 at the junction between southeastward and southward flow vectors. Due to this velocity gradient, longitudinal ice

compression and a related emergence velocity are expected to occur here, which coincides with the reported elevation gain of +0.5 m over the summer (Immerzeel and others, 2014). We find summer velocities for Lirung Glacier that are considerably lower than those reported by Naito and others (1998). They state the glacier had moved between 2.8 and 7.5 m (~6.5 and 18.0 m a-1) in the middle part and between 1.9 and 2.5 m (~4.5 and 6.0 ma-1)

Summer Winter

1 I I I 1 I r

0 1 2 3 4 5 6 < Transverse distance (m)

Surface velocity (m a"1)

Fig. 7. Surface velocity and flow direction obtained by noise-filtered frequency cross-correlation (W256 F64) for the summer (left) and winter (middle) period. The plots on the right show transverse surface velocity profiles for both seasons taken at the three indicated locations.

in the lower area over the period June-October in the years 1994-96. In two decades the glacier thus seems to have reduced its flow velocity by about a factor of two, though short-term velocity variations could contribute to this difference.

The winter period shows a considerably different picture to that of the summer. Here the high velocities for the upper area are absent and have reduced to ~2-3 ma-1. Velocities on the lower (southward-flowing) portion of the terminus are similar to those found in the summer period. Combining and rescaling the summer and winter velocities, respectively valid for 157 and 191 days, to full-year values, the upper portion of the surveyed area has velocities of ~3.5 m a-1 and velocities at the transition zone between upper and lower regions are ~1.5 m a-1.

The distinct contrasts between (1) the surface velocities for the upper and lower portions of the terminus in the summer and (2) the summer and winter velocity fields, indicate the presence of two different dominant flow regimes (Copland and others, 2009). We hypothesize that the faster flow in summer is caused by basal-sliding-dominated processes, while the lower velocities found in winter and in the lower portions are mainly due to deformation. The large amounts of monsoon precipitation and the opening of sub- and englacial conduits due to rising temperatures (Benn and others, 2012) are likely to lubricate the base of the ice and introduce a basal-sliding component to the flow in the summer season. In areas where basal sliding dominates, the ice is expected to move in a block-like motion with relatively high and constant velocities in the centre and sharp lateral velocity gradients (Copland and others, 2009). Transverse velocity profiles over the glacier (Fig. 7) show that there indeed is a difference between the summer and winter flow in terms of lateral gradients. Especially near the northeastern ice boundary, the summer velocity profile is reduced by 3 ma-1 over a few tens of metres (profiles A and B). The lateral winter velocity variation generally shows a more parabolic pattern, as does the summer velocity at profile C. This is usually associated with more deformation-dominated flow (Copland and others, 2009). Why the basal sliding occurs only in the upstream area is not entirely clear, but it probably results from increased driving stresses caused by thicker ice, that is due to the regular avalanches and rockfall from the steep slopes of the Langtang Lirung Peak to the northwest. This difference in ice thickness may also play a role in the contrasting velocities found laterally, i.e. fast flow at the western ice boundary and limited flow on the other side of the tongue.

Glacier ice flow is a complex process and is governed by a wide range of processes and forces (Van der Veen, 2013). Nevertheless, an improved understanding of the ice thickness of Lirung Glacier and the local bedrock configuration underneath the ice will greatly contribute to a better understanding of the flow patterns found in this study. Furthermore, it would provide the ability to estimate mass turnover rates that are related to the flow velocities found. Ice thickness measurements of the glacier were performed ~15 years ago using radio-echo sounding techniques (Gades and others, 2000), but the quality, resolution, specific locations and time period of the measurements make them unsuitable in this case. A new survey of Lirung Glacier using ground-penetrating radar would help to fill some of the gaps raised by this study.

Value of UAV surveys

As the typical pixel size of spaceborne imagery that is suitable for glacier velocity monitoring is often considerably larger than the seasonal or even annual displacements of Himalayan debris-covered glaciers, data that span multiple years are generally required to extract meaningful velocity signals (Kaab, 2005; Scherler and others, 2008; Herman and others, 2011). This renders the quantification of seasonal variations in surface height change and velocities very difficult. Of course, this is even more of an issue when flow velocities are relatively low, such as for many debris-covered glaciers in the Himalaya (Quincey and others, 2009b). Although high temporal resolution can be achieved using in situ methods, they are unfeasible for high spatial resolution surveys of large glacier surfaces, as fieldwork on debris-covered glaciers is often difficult and time-consuming.

The use of UAVs allows high-resolution continuous values of the surface velocities of a single season to be obtained. The techniques used here would also allow for large-scale determination of interannual flow. This will improve our understanding of the relationship with local precipitation and temperature perturbations, which will eventually lead to the ability to provide better predictions of possible future changes in glacier volume under climate change scenarios. A deepened knowledge of the smaller-scale variations in flow, both spatially and temporally, also helps to unravel the bigger picture of heterogeneous mass wasting and distribution of surface features found on debris-covered glaciers (Immerzeel and others, 2014).


In this study, UAVs are used to acquire images of debris-covered Lirung Glacier for May and October 2013 and May 2014. The imagery is processed into orthomosaics and DEMs using a SfM workflow and georeferenced using GCPs and tie points. Displacements of the glacier surface are derived for both summer and winter using an automated frequency cross-correlation algorithm, which is tested for sensitivity to input datasets and parameters. From the study we draw the following conclusions.

Summer and winter surface velocities for Lirung Glacier are ~6 and 2.5 ma-1, respectively, in the upstream part of the tongue. In the bend and in the lower areas of the tongue both seasons show comparable slow flows of ~1.5-2 m a-1 and stagnancy. The differences in surface velocity and flow direction between the two seasons lead to the hypothesis that the fast flow in summer is caused by basal-sliding-dominated processes, while the lower velocities found in winter are mainly due to plastic deformation. Transverse velocity profiles over the glacier seem to confirm this hypothesis. For an improved understanding of the spatial surface velocity differences and flow patterns of Lirung Glacier found in this study it is important learn more about its ice thickness and bedrock configuration.

Frequency cross-correlation techniques applied to highresolution UAV imagery can determine surface velocities of a debris-covered glacier. Displacements unrelated to ice flow can largely be filtered out by the correlation method, and any remaining noise can be removed using post-correlation noise-reduction techniques. In comparison to a manual digitization technique, both methods have similar accuracies taking into account the associated errors. The continuous

surface output that the correlator yields provides more detail, however, and the method is less time-consuming.

It is found that using an orthomosaic as input for the correlation outperforms the use of a hillshade or an edge-detected DEM, in terms of irregularly distributed output noise. The use of different settings for the correlation algorithm does not alter surface velocities and flow directions significantly. There are, however, subtle differences present and small window sizes give better performance in terms of SNR, the retention of detail and the overall results in comparison with manual digitization. However, displacements such as the melt of ice cliffs are not filtered out; this requires the use of larger correlator windows which can result in loss of fine-scale detail. Optimal settings for the input resolution are found to be an initial window size of 256 px with a final window size of 64 px.

The use of UAV imagery and feature-tracking algorithms allows determination of high-resolution seasonal surface velocities, something not possible with most spaceborne remote-sensing techniques. Our approach yields insights into the smaller-scale temporal and spatial variations in glacier flow, and improves our understanding of heterogeneous mass wasting and surface features found on debris-covered glaciers.


This study was carried out with funding from the Climate Research and Information Services in South Asia project of the UK Department for International Development (DFID), the Innovation Fund of the International Centre for Integrated Mountain Development (ICIMOD), The Netherlands Organization for Scientific Research (NWO) and the European Institute of Innovation and Technology (EIT). We thank Rijan Kayastha and the Department of National Parks and Wildlife Conservation (Nepal) for facilitating our research permits. We also thank all the other individuals who helped with fieldwork or by providing equipment: Waqar Ali, Fionna Heuff, Martin Heynen, Sharad Joshi, Arthur Lutz, Pradeep Mool, Lene Petersen, Allen Pope, Arun Shrestha, Eduardo Soteras, Patrick Wagnon, Niko Wanders, Muhum-mad Atif Wazir and Simon Wicki. Finally, we thank Duncan Quincey and Adrian Luckman for their reviews of the manuscript and their constructive comments.


Agisoft (2013) PhotoScan Professional 0.9.1 user manual. Agisoft, St Petersburg

Ayoub F, Leprince S and Keene L (2009) User's guide to COSI-CORR co-registration of optically sensed images and correlation. California Institute of Technology, Pasadena, CA Benn D and 10 others (2012) Response of debris-covered glaciers in the Mount Everest region to recent warming, and implications for outburst flood hazards. Earth Sci. Rev., 114(1-2), 156-174 (doi: 10.1016/j.earscirev.2012.03.008) Bolch T, Pieczonka T and Benn D (2011) Multi-decadal mass loss of glaciers in the Everest area (Nepal Himalaya) derived from stereo imagery. Cryosphere, 5(2), 349-358 (doi: 10.5194/tc-5-349-2011)

Bolch T and 11 others (2012) The state and fate of Himalayan glaciers. Science, 336, 310-314 (doi: 10.1126/science.1215828) Buades A and Coll B (2005) A non-local algorithm for image denoising. Comput. Vision Patt., 2, 60-65 (doi: 10.1109/ CVPR.2005.38)

Buri P, Pellicciotti F, Steiner JF, Miles E and Immerzeel WW (2016) A grid-based model of backwasting of supraglacial ice cliffs over debris-covered glaciers. Ann. Glaciol., 71 Copland L and 8 others (2009) Glacier velocities across the central Karakoram. Ann. Glaciol., 50(52), 41 (doi: 10.3189/ 172756409789624229) Davis JC (2002) Statistics and data analysis in geology. John Wiley, New York

Gades A, Conway H, Nereson N, Naito N and Kadota T (2000) Radio echo-sounding through supraglacial debris on Lirung and Khumbu Glaciers, Nepal Himalayas. IAHS Publ. 264 (Workshop at Seattle 2000 - Debris-Covered Glaciers), 13-22 Gardelle J, Berthier E and Arnaud Y (2012) Slight mass gain of Karakoram glaciers in the early twenty-first century. Nature Geosci., 5(5), 322-325 (doi: 10.1038/ngeo1450) Gardelle J, Berthier E, Arnaud Y and Kââb A (2013) Region-wide glacier mass balances over the Pamir-Karakoram-Himalaya during 1999-2011. Cryosphere Discuss., 7(2), 975-1028 (doi: 10.5194/tc-7-1263-2013) Heid T and Kââb A (2012) Repeat optical satellite images reveal widespread and long term decrease in land-terminating glacier speeds. Cryosphere, 6(2), 467-478 (doi: 10.5194/tc-6-467-2012) Herman F, Anderson B and Leprince S (2011) Mountain glacier velocity variation during a retreat/advance cycle quantified using sub-pixel analysis of ASTER images. J. Glaciol., 57(202), 197-207 (doi: 10.3189/002214311796405942) Immerzeel W, Van Beek L and Bierkens M (2010) Climate change will affect the Asian water towers. Science, 328(5984), 1382-1385 (doi: 10.1126/science.1183188) Immerzeel W, Pellicciotti F and Bierkens M (2013) Rising river flows throughout the twenty-first century in two Himalayan glacierized watersheds. Nature Geosci., 6, 742-745 (doi: 10.1038/ngeo1896) Immerzeel W and 6 others (2014) High-resolution monitoring of Himalayan glacier dynamics using unmanned aerial vehicles. Remote Sens. Environ., 150, 93-103 (doi: 10.1016/j.rse. 2014.04.025)

Kââb A (2005) Combination of SRTM3 and repeat ASTER data for deriving alpine glacier flow velocities in the Bhutan Himalaya. Remote Sens. Environ., 94(4), 463-474 (doi: 10.1016/j. rse.2004.11.003) Kââb A, Berthier E, Nuth C, Gardelle J and Arnaud Y (2012) Contrasting patterns of early twenty-first-century glacier mass change in the Himalayas. Nature, 488(7412), 495-498 (doi: 10.1038/nature11324) Kaser G, Grosshauser M and Marzeion B (2010) Contribution potential of glaciers to water availability in different climate regimes. Proc. Natl Acad. Sci. USA, 107, 20223-20227 (doi: 10.1073/pnas.1008162107) Leprince S, Ayoub F and Avouac JP (2007) Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements. IEEE Trans. Geosci. Remote Sens., 45(6), 1529-1558 (doi: 10.1109/TGRS.2006.888937) Leprince S, Berthier E, Ayoub F, Delacourt C and Avouac JJP (2008) Monitoring earth surface dynamics with optical imagery. Eos, 89(1) (doi: 10.1029/2008E0010001) Lillesand TM, Kiefer RW and Chipman JW (2003) Remote sensing

and image interpretation. John Wiley & Sons, New York Lucieer A, Jong S and Turner D (2013) Mapping landslide displacements using Structure from Motion (SfM) and image correlation of multi-temporal UAV photography. Progr. Phys. Geogr., 38(1), 97-116 (doi: 10.1177/0309133313515293) Lutz AF, Immerzeel WW, Shrestha AB and Bierkens MFP (2014) Consistent increase in High Asia's runoff due to increasing glacier melt and precipitation. Nature Climate Change, 4(June), 1-6 (doi: 10.1038/nclimate2237) Naito N and 8 others (1998) Surface flow on the ablation area of the Lirung Glacier in Langtang Valley, Nepal Himalayas. Bull. Glacier Res., 16, 67-73

Ostrem G (1959) Ice melting under a thin layer of moraine, and the existence of ice cores in moraine ridges. Geogr. Ann., 41(4), 228-230

Pellicciotti F, Stephan C, Miles E, Herreid S, Immerzeel WW and Bolch T (2015) Mass-balance changes of the debris-covered glaciers in the Langtang Himal, Nepal, from 1974 to 1999. J. Glaciol., 61(226), 373-386 (doi: 10.3189/2015JoG13J237) Quincey DJ, Luckman A and Benn D (2009a) Quantification of Everest region glacier velocities between 1992 and 2002, using satellite radar interferometry and feature tracking. J. Glaciol., 55(192), 596-606 (doi: 10.3189/002214309789470987) Quincey DJ, Copland L, Mayer C, Bishop M, Luckman A and Belo M (2009b) Ice velocity and climate variations for Baltoro Glacier, Pakistan. J. Glaciol., 55(194), 1061-1071 (doi: 10.3189/002214309790794913) Sakai A, Nakawo M and Fujita K (1998) Melt rate of ice cliffs on the Lirung Glacier, Nepal Himalayas, 1996. Bull. Glacier Res., 16, 57-66

Scambos T, Dutkiewicz MJ, Wilson JC and Bindschadler RA (1992) Application of image cross-correlation to the measurement of glacier velocity using satellite image data. Remote Sens. Environ., 42(3), 177-186 (doi: 10.1016/0034-4257(92)90101-O)

Scherler D, Leprince S and Strecker MR (2008) Glacier-surface velocities in alpine terrain from optical satellite imagery:

accuracy improvement and quality assessment. Remote Sens. Environ., 112(10), 3806-3819 (doi: 10.1016/j.rse.2008.05.018) Scherler D, Bookhagen B and Strecker MR (2011) Spatially variable response of Himalayan glaciers to climate change affected by debris cover. Nature Geosci., 4(3), 156-159 (doi: 10.1038/ ngeo1068)

SenseFly (2015) Drones for professionals, mapping & photogram-metry, flight planning & control software. https://www.sensefly. com/home.html

Steiner JF, Pellicciotti F, Buri P, Miles ES, Immerzeel WW and Reid TD (in press) Modelling ice-cliff backwasting on a debris-covered glacier in the Nepalese Himalaya. J. Glaciol., 61(229) (doi: 10.3189/2015JoG14J194)) Sugiyama S, Fukui K, Fujita K, Tone K and Yamaguchi S (2013) Changes in ice thickness and flow velocity of Yala Glacier, Langtang Himal, Nepal, from 1982 to 2009. Ann. Glaciol., 54(64), 157-162 (doi: 10.3189/2013AoG64A111) Szeliski R (2010) Computer vision: algorithms and applications.

Springer, London Van der Veen C (2013) Fundamentals of glacier dynamics, 2nd edn.

CRC Press, Boca Raton, FL Westoby MJ, Brasington J, Glasser NF, Hambrey MJ and Reynolds JM (2012) 'Structure-from-Motion' photogrammetry: a low-cost, effective tool for geoscience applications. Geomorphology, 179, 300-314 (doi: 10.1016/j.geomorph.2012.08.021)