Monitoring Tropical Debris Covered Glacier Dynamics from High Resolution Unmanned Aerial Vehicle Photogrammetry , Cordillera Blanca , Peru

The glaciers of the Cordillera Blanca Peru are rapidly retreating and thinning as a result of climate change, altering the timing, quantity and quality of water available to downstream users. Furthermore, increases in the number and size of proglacial lakes associated with these melting glaciers is increasing potential exposure to glacier lake outburst floods (GLOFs). Understanding how these glaciers are changing and their connection to proglacial lake systems is thus of critical importance. Most satellite data 30 are too coarse for studying small mountain glaciers and are often affected by cloud cover, while traditional airborne photogrammetry and LiDAR are costly. Recent developments have made Unmanned Aerial Vehicles (UAVs) a viable and potentially transformative method for studying glacier change at high spatial resolution, on demand and at relatively low cost. Using a custom designed hexacopter built for high altitude (4000-6000masl) operation we completed repeat aerial surveys (2014 35 and 2015) of the debris covered Llaca glacier tongue and proglacial lake system. High resolution orthomosaics (5cm) and digital elevation models (DEMs) (10cm) were produced and their accuracy assessed. Analysis of these datasets reveals highly heterogeneous patterns of glacier change. The most rapid areas of ice loss were associated with exposed ice cliffs and melt water ponds on the glacier surface. Considerable subsidence and low surface velocities were also measured on the sediments within the pro-glacial lake, indicating the presence of extensive regions of buried ice and continued connection to the glacier tongue. Only 40 limited horizontal retreat of the glacier tongue was observed, indicating that simple measurements of changes in aerial extent are inadequate for understanding actual changes in glacier ice quantity.


Introduction
The Peruvian Andes hold around 70 % of the world's tropical ice.Within Peru the Cordillera Blanca is the most glacierised mountain range with around 40 % of Peru's glacier area (ANA, 2014).These glaciers are rapidly retreating primarily as a result of warming temperatures (Vuille et al., 2008(Vuille et al., , 2003)).Most recently Burns and Nolin (2014) calculated a 25 % reduction in glacier area from 1987 to 2010.Similarly, Racoviteanu et al. (2008) calculated a 22.4 % glacier area reduction from 1970 to 2003, during which time the mean glacier terminus elevation increased by 113 m.Glacier changes are altering downstream hydrology, changing water supply in terms of quantity, quality and timing (Baraer et al., 2009(Baraer et al., , 2012;;Chevallier et al., 2011;Condom et al., 2011;Juen et al., 2007;Kaser et al., 2003;Rabatel et al., 2013).These changes have profound impacts on downstream communities, livelihoods and ecology (Bury et al., 2013(Bury et al., , 2011;;Carey et al., 2014;Mark et al., 2010;Postigo et al., 2008;Young, 2009Young, , 2014)).Additionally, rapidly melting glaciers are increasing the number and volume of unstable proglacial lakes that present a serious natural hazard through potential glacier lake outburst floods (GLOFs) (Carey, 2008;Carey et al., 2012;Huggel et al., 2002;Lliboutry et al., 1977; O. Wigmore and B. Mark: Monitoring tropical debris-covered glacier dynamics Portocarrero, 2014).Understanding how these glaciers are changing, their connection to proglacial lake systems and the role of glacier surface heterogeneity in controlling melt processes is critical to inform future management and adaptation strategies.Unfortunately, current understanding has been limited due to the spatial and temporal resolution of available datasets and difficulties of accessibility to these remote mountain locations.
Processes of glacier melt and the formation of glacier lakes is controlled by a multitude of factors including local energy balance, topography, geology, aspect, valley dimensions, moraine stability, debris thickness and the presence of surface features such as melt ponds and ice cliffs (Benn et al., 2012;Buri et al., 2016;Harrison et al., 2006;Huggel et al., 2002;Immerzeel et al., 2014;Lejeune et al., 2013;Portocarrero, 2014;Reid et al., 2012;Reid and Brock, 2010;Thompson et al., 2012).The impact of these variables has been studied extensively in other areas, particularly the Himalaya, where large debris-covered glacier tongues are common (Benn et al., 2012;Buri et al., 2016;Immerzeel et al., 2014;Scherler et al., 2011;Shroder et al., 2000).However, few detailed studies of these processes have focussed on the Cordillera Blanca.Remote sensing (satellite and airborne) has been widely used to quantify aerial changes in glacier extent for the Cordillera Blanca and is routinely applied in national-level inventories (ANA, 2014;Burns and Nolin, 2014).However, detailed studies of glacier melt processes and volumetric changes have been limited in the Cordillera Blanca due in part to the extreme topographic variation, generally small glacier size (typically less than 0.5 km 2 ) and paucity of suitable datasets.Mark and Seltzer (2005) integrated digital elevation model (DEM) differencing and climate variables to investigate the role of aspect and elevation on glacier ablation.Huh et al. (2017) combined airborne lidar and photogrammetry to quantify changes in glacier volume for a handful of valleys within the Cordillera Blanca, but large errors exist for the earlier date photogrammetric DEMs.Most recently Emmer et al. (2015) combined a single DEM with high-resolution satellite imagery and field investigations to understand the geomorphic evolution and measure surface velocities of the debris-covered Jatunraju Glacier.Most research in the Cordillera Blanca has addressed the entire range with a focus on glacier response to larger-scale climatic perturbations (Rabatel et al., 2013;Vuille et al., 2008) and their role in the hydrologic budget of the Rio Santa basin (Baraer et al., 2012;Juen et al., 2007;Kaser et al., 2010).A number of studies have investigated mass balance response of individual glaciers, though these have generally aggregated mass balance at the glacier scale (Hastenrath and Ames, 1995a, b;Kaser et al., 1990;Winkler et al., 2009).Few studies have documented glacier surface and volume changes at fine resolution in the Cordillera Blanca.
There are three primary methods for quantifying glacier change: glaciologic, hydrologic and geodetic (Cuffey and Paterson, 2010;Ostrem and Brugman, 1993).Only geodetic methods can provide detailed spatially continuous measurements of glacier surface changes.Geodetic methods rely on repeat remote sensing observations (satellite, airborne and terrestrial) to map surface characteristics and measure displacement velocities and vertical changes through differencing of DEMs.Traditional methods of aerial remote sensing (satellite, airplane/helicopter) are constrained by three competing primary factors: spatial resolution, temporal resolution and cost.To make a gain in one of these factors there is a simultaneous loss in one or both of the other factors.Additionally, in mountain regions cloud cover is a notable issue as is safety for airborne platforms.However, recent developments in military, remote control, hobby and opensource hardware have made it viable and potentially transformative to consider new small (1-5 kg), relatively inexpensive unmanned aerial vehicle (UAV) platforms for mapping glacier dynamics.UAV platforms include, helicopters, multirotors, fixed wings and paragliders; they are typically small (1-5 kg) and relatively cheap, ranging in cost from USD 400 to 20 000 (Eisenbeiß, 2009;Eisenbeiss and Sauerbier, 2011;Friedli, 2013).These platforms may be operated by remote control, first person or GPS-enabled autopilot systems.They can be fitted with a diversity of sensors for collecting observations of the earth surface and atmosphere, including multispectral cameras, meteorology sensors and lidar units (Colomina and Molina, 2014;Hardin and Jensen, 2011;Nex and Remondino, 2014).Due to the low flight elevations (50-1000 m a.g.l.), very high-resolution imagery (< 5 cm pixels) can be collected from below the cloud level with limited risk to ground personnel or property.In addition, these ondemand platforms can collect data at whatever temporal resolution is desired with no major additional acquisition costs.
In addition to these developments in UAV hardware, software developments have also occurred.Easy-to-use, off-theshelf software is now capable of processing large UAV image collections to derive orthomosaics and dense photogrammetric point clouds for the generation of DEMs.By integrating these data with differential GPS surveyed ground control points (GCPs), it is possible to generate DEMs with a ground resolution of under 5 cm and orthomosaics of around 2 cm resolution (Harwin and Lucieer, 2012;Hugenholtz et al., 2013).Accuracy of these datasets has been demonstrated to be very high in comparison to survey data and often equal or better than lidar (Harwin and Lucieer, 2012;Hugenholtz et al., 2013).
The use of UAVs in earth science research is undergoing rapid expansion with numerous researchers exploiting the financial and spatiotemporal resolution benefits of UAV technology (Ahmad, 2011;Colomina and Molina, 2014;D'Oleire-Oltmanns et al., 2012;Hardin and Jensen, 2011;Koh and Wich, 2012;Vivoni et al., 2014;Watts et al., 2012).A few projects have utilised UAV platforms in mountainous glaciated regions (Bhardwaj et al., 2016).Friedli (2013) has investigated their potential for glacier monitoring on the Rhône Glacier, while Whitehead et al. (2013) conducted a detailed glacier study in the Canadian Arctic, but this was at low elevation.At higher elevations Immerzeel et al. (2014) and Kraaijenbrink et al. (2016) successfully used a small fixed-wing platform to map glacier changes and surface velocities at 5000 m in the Himalaya.Recently researchers have begun to explore the use of UAVs within the Cordillera Blanca, however, to our knowledge no published studies to date have deployed UAVs to study glacier dynamics within the Cordillera Blanca or the Tropical Andes.

Objectives
There are three principle objectives of this study: 1. Successfully deploy a multirotor platform suitable for mapping glacier changes and capable of operating at ∼5000 m a.s.l. in mountainous terrain.
2. Measure glacier volume changes, surface velocities and proglacial lake changes at high resolution using decimetre DEMs and centimetre orthomosaics and provide a quantitative assessment of their accuracy.
3. Use these datasets to investigate processes and patterns of glacier change over a debris-covered glacier tongue and proglacial lake system in the Cordillera Blanca.

Study area
The Cordillera Blanca is located in the central Andes of Peru (Fig. 1).The range extends 180 km (north-south), encompassing numerous peaks over 6000 m, including Peru's highest, Huascaran, at 6768 m.The region contains more than 700 individual glaciers, which terminate at over 4500 m a.s.l.The climate is dominated by a strong wet-dry seasonal signal.The wet season corresponds to the austral summer (October-May), when roughly 80 % of the 800-1200 mm annual precipitation falls (Bury et al., 2011;Mark et al., 2010).
As with other tropical mountain regions, diurnal temperature variation greatly exceeds seasonal temperature variability (Kaser, 2001).Maximum glacier ablation occurs at the same time as maximum accumulation during the wet season months (Kaser, 2001).The focus of this study is Llaca Glacier, which is located in the central Cordillera Blanca above the city of Huaraz (Fig. 1).The Llaca Valley is relatively narrow, constricting to around 500 m at its narrowest.The head of the valley opens into a large cirque over 3 km wide with headwalls over 5500 m that include the summits of Ranrapalca (6162 m) and Vallunaraju (5686 m).Due to this large accumulation zone and narrow valley, Llaca Glacier extends to a considerably lower elevation (4500 m a.s.l.) than many other glaciers in the region that typically terminate between 4700 and 5100 m a.s.l.; for example, Vallunaraju shares a divide with Llaca and terminates at ∼ 4800 m (300 m higher).The glacier snout is covered in a thick (> 1 m) layer of debris that derives from the collapse of steep lateral moraines onto the glacier surface.This debris insulates the ice and enables the glacier to extend to this lower terminal elevation.The lower debris-covered tongue remains connected to the upper glacier and terminates in a moraine dammed and partially ice-cored proglacial lake that is approximately 1 km long, 0.2 km wide, up to 17 m deep (in 2004) and stores an estimated 274 300 m 3 of liquid water (Portocarrero, 2014).The moraine dam underwent engineered stabilisation and lowering (completed in 1977) as a component of the broader campaign carried out by the Huaraz-based glaciology office throughout the Cordillera Blanca to mitigate the risk of GLOFs that have been responsible for a number of major natural disasters and large loss of life in the region (Carey, 2008;Carey et al., 2012;Lliboutry et al., 1977;Portocarrero, 2014).Lake stabilisation efforts included constructing a reinforced concrete and steel dam within the terminal moraine and artificially lowering the lake level.Despite these efforts, a significant GLOF risk is still present due to avalanche risk from hanging glaciers above, tectonic activity and the location of the stream channel draining the lake that passes through the northern edges of the city of Huaraz to join the Rio Santa (Emmer and Vilímek, 2014;Portocarrero, 2014).

UAV platform design and specifications
A hexa-multirotor UAV was custom designed and built for this project (Fig. 2).The UAV is capable of operating at elevations between 3500 and 6000 m a.s.l.The 1 m (diameter) frame utilises a lightweight carbon fibre construction and is fitted with high-speed motors and 33 to 38 cm (13 × 5.5 to 15 × 5.5 inch) carbon fibre propellers (depending on flight altitude).Total weight excluding cameras is ∼ 2 kg s.Low weight is essential for operation at this altitude and eases transport when hiking/climbing to the glaciers.The maximum flight time for the platform at 5000 m a.s.l. is around 12-15 min on a single 10 000 MAH 4S lithium polymer battery, giving a range of 1.5-2.5 km (depending on above ground flight altitude) from launch at flight speeds acceptable for aerial mapping (∼ 8 m s −1 ).Autonomous navigation is managed by a 3DR Pixhawk flight controller.Compared to fixed-wing platforms, a multirotor platform sacrifices flight time in exchange for greater navigational precision, vertical takeoff and landing (essential over debris-covered glacier surfaces and around water bodies) and typically higherquality images due to greater platform stability.Communication with the UAV is maintained with an extended-range (2.5 km) 2.4 GHz radio controller (FrSky Taranis LR) and a long-range (> 40 km) ∼ 915 MHz (RFD900+) telemetry link between the autopilot and ground station.Ground station control was managed by a field tablet running APM Mission  Planner for Android (2014) and APM Mission Planner for Windows (2015).
The UAV was fitted with a single Canon S110 Powershot camera.The camera was installed in nadir view with no gimbal stabilisation.This camera allows full manual control of exposure settings including ISO, F-stop and shutter speed.The best settings are an ISO ∼ 100-200 (to reduce image graininess), F-stop > 4 (to improve depth of field) and shutter speed < 1/1000 s (to reduce motion blur).Camera control was managed by the Canon Hack Development Kit (CHDK).CHDK reboots the camera using firmware stored on the SD card and allows the use of custom scripts.In this case the kap_uav.luascript (CHDK, 2016) was used.This script allows the specification of the camera parameters above and automatic image capture at set intervals.An interval of 3 s was selected to provide > 90 % image overlap, with additional redundancy for blurred image removal.

Ground control
Before each survey a comprehensive Global Navigation Satellite System (GNSS) survey was conducted.This included installation and survey of highly visible targets for use as GCPs in georectification of the photogrammetric point cloud and collection of "check points" for accuracy assessment of derived DEMs (Figs. 2 and 3, Table 1).Target installation points were distributed fairly evenly across the survey area making sure to include high and low points within the survey region.Access to the upper portion of the study area is very difficult and targets could not be installed here.Targets were anchored to the ground using metal stakes and rocks as appropriate to ensure they would not move before flights.Selection of check points was semi-random, making sure to include high and low points within the land surface.In 2014, points were taken atop and between individual rocks as well as over more uniform terrain.In 2015, check points were placed primarily in areas that exhibit relative terrain homogeneity within a roughly 50 cm radius.
The 2014 base station data (for rover position correction) was collected with a Trimble 5700 receiver with a Zephyr geodetic antenna provided by UNAVCO that was installed on an apartment roof in Huaraz (baseline ∼ 14 km) for almost 7 days, collecting at 1 Hz (Table 2).In 2015, base station data were collected with a Topcon Hiper-SR integrated receiver/antenna installed at a site on the southern end of the moraine in Llaca Valley (baseline ∼ 1 km) that was occupied for 6.5 h collecting at 1 Hz.Base station positions were resolved using the Natural Resources Canada (NRCAN) Precise Point Positioning online service.Rover positions (targets and check points) were occupied for 5 min at 1 Hz frequency following a stop-go post-processing workflow: observations were collected using a Topcon GRS1 L1 receiver with a PG-A1 external antenna (Table 2).Rover observations were post-processed using Topcon Magnet TM Office Tools (version 2.7.1) software with base station coordinates adjusted as per NRCAN positioning results.

Aerial survey flight paths
Flight planning was completed in the office using mission planner software and Google Earth.Because the UAV does not have collision avoidance, it was critical to carefully plan a flight path that avoids topographic obstacles such as moraine walls and cliffs, especially considering that for much of the survey the UAV would be flying beyond line of sight.Much of Google Earth's topographic data relies on SRTM90; this dataset has too coarse a resolution and is often out of date (collected in 2000) for near surface navigation in the dynamic and heterogeneous terrain of the Cordillera Blanca.Additionally, limited sky view and possible dropped satellites present real navigational threats for the UAV itself.Therefore, an above ground level (a.g.l.) of ∼ 100 m was selected, as this would provide sufficient ground resolution (∼ 3 cm s) while also keeping the UAV well above the moraine walls (Fig. 3, Table 3).To maintain this a.g.l.altitude, the flight plan required the UAV to first ascend to 100 m, then gain an additional ∼ 200 m as it moved up the valley over the glacier.This also maintained a uniform ground resolution.Five flight lines were constructed from the launch point below the calving face to provide ∼ 60 % sidelap, which is necessary for robust DEM recovery (Fig. 3).The maximum length of an individual flight line was 1.2 km from the take off point, resulting in a total (return) flight length of roughly 2.5 km.Total flight time was around 12 min per flight travelling at ∼ 8 m s −1 .Rapid battery drain was experienced in the take-off portion and when gaining altitude up-valley, which limited the range and thus the total area surveyed.

Structure from motion (SfM) workflow -DEM and orthomosaic generation
UAV images were processed using the SfM workflow as implemented in Agisoft PhotoScan Professional Edition (Version 1.1) (Agisoft, 2016).This is proprietary software; thus the specific algorithms used are not available.However, the sequential steps that were followed in this study are outlined below and are similar across all commercial and open-source software packages.Further details on SfM can be found in Fonstad et al. (2013), Lowe (2004), Snavely et al. (2008), Szeliski (2010) and Verhoeven (2011).
1.A feature recognition algorithm is used to identify unique tie points across the image dataset, e.g.scale invariant feature transform (SIFT) (Lowe, 2004).
2. A bundle block adjustment is performed to reconstruct camera positions and construct a sparse point cloud.
3. Ground control targets are identified within the scene and positioned on each individual photo in which they are visible.Placement accuracy of individual markers is within ∼ 3 cm (width of the X on targets).For this study, DEMs were generated at 10 cm pixel resolution and RGB orthomosaics were created at 5 cm pixel resolution.

Accuracy assessment
Accuracy of the DEMs is estimated in two ways.Firstly, the photogrammetric software provides a number of different error statistics including tie point matching error and GCP placement error (difference between actual coordinates and estimated position from SfM processing) (Table 4).These parameters describe how well the point cloud and camera positions fit the in-scene ground targets.Error is provided as RMSE, x − y − z differences and pixel matching errors.These errors are reduced within the SfM workflow by optimising the point cloud while minimising potential position variance from control points.Outliers are selected and removed through a few iterations to reduce error estimates.Secondly, the surveyed GPS elevation at the check point locations is compared against the elevation extracted from the DEM surface to provide a more precise estimate of absolute elevation error for each surface.A third method for estimation of error is to compare difference in the DEM elevations over regions that have experienced no change.Unfortunately, there were no regions within the study area that could be identified unequivocally as not having undergone elevation change over the study period.This method was therefore not used.

Elevation change, surface velocity and orthomosaic comparison
Analysis of glacier changes was conducted in ESRI Ar-cMap 10.2.Volumetric change was calculated for the entire survey area and the debris-covered glacier section separately.This was done by clipping the DEMs to the glacier boundary and subtracting the 2014 DEM from the 2015 DEM (i.e.negative values equal elevation lowering and/or ice loss).Surface velocity was determined by manual feature tracking of 72 points on the glacier tongue to produce velocity vectors.Following vector creation, a velocity surface was interpolated for the glacier tongue using the spline interpolation tool in ArcMap.To investigate melt patterns across the glacier surface, horizontal movement must be removed so that the same area of the glacier surface is compared.To do this, the velocity vectors above were used to orthorectify the 2015 orthomosaic to the same location as the 2014 orthomosaic.Vectors of ice cliff movement were then calculated between the two dates.Within the lake section of the survey area, elevation changes of water bodies were measured, ice-cored regions were investigated and limited feature tracking was used to estimate horizontal movement.

GNSS and DEM accuracies
Base station positions are located with an estimated combined (horizontal and vertical) error of 0.5 cm (2014) and 7 cm (2015) (Table 2).Positional errors of the GCPs and check points are all estimated to be under 2.5 cm (2014) and 0.6 cm (2015) (Table 2).Providing an estimated combined positional error of 3 cm (2014) and 8 cm (2015).Thus maximum expected errors in positions compared between the two dates are ∼ 11 cm.The 2014 rover positions are less accurate due to the long (15 km) baseline, but the longer base station occupation interval means they are overall better constrained.Rover positions in 2015 are more accurate due to a shorter (1 km) baseline but the base station position is less accurate due to a much shorter occupation interval (6.5 h vs. almost 7 days) (Table 2).This could not be avoided as the base station could not be left in place overnight due to security concerns.Estimated tie point (image) matching errors are summarised in Table 4, with a mean value of 0.63 pixels (at 3.5 cm pixel size; Table 3) for 2014 and 0.43 pixels (at 3.6 cm pixel size; Table 3) for 2015.Errors in both cases are sub pixel at the minimum pixel size; these results are improved when the DEM and orthomosaics are aggregated for analysis to 10 and 5 cm respectively.Additional error estimates are provided as the measured difference (in metres and pixels) between GCP coordinates and estimated GCP position within the SfM model (Table 4).These errors are sub-centimetre for both years: 0.002 m in 2014 and 0.003 m in 2015.The error statistics provided above however tend to overestimate DEM accuracy.DEM accuracy should therefore be assessed by either comparing differences over no change areas or comparing to additional survey points not used in generating the DEM.
To gain a better understanding of the absolute accuracy of the two DEM surfaces, elevations from the DEM were compared against those from check points surveyed in their respective years (Figs. 3 and 4).The results of this analysis are summarised in Figs. 5 and 6.The mean difference between DEM and GNSS check points in 2014 is −0.037 m with a standard deviation (SD) of 0.191 m; in 2015 the mean difference is 0.001 m with a SD of 0.044 m.Negative values indicate underestimation of DEM surface with respect to surveyed positions and thus imply flattening of the DEM surface.Most errors are less than ±5 cm and are the combined result of errors in image matching/SfM processing, GNSS positional errors and DEM smoothing of small surface features.The 2014 data have a much larger error range than in 2015, but this is deemed to be the result of experimental error as opposed to true inaccuracies in the DEM surface.In 2014, check points were collected in areas of non-uniform topography, i.e. atop rocks, between cracks; these minor topographic features are smoothed out of the DEM resulting in a larger difference between GNSS and DEM elevations.The 2015 check points were located in areas with only minor topographical variation within a roughly 0.5 m radius, thereby minimising the impacts described above and resulting in higher estimated accuracy.

Measured changes over the survey area
Changes in surface elevation were highly variable across the survey area, including areas of both loss and gain (Fig. 7).Over the glacier tongue (blue boundary Fig. 7), the www.the-cryosphere.net/11/2463/2017/The Cryosphere, 11, 2463-2480, 2017   There is a roughly linear decrease in velocity down-glacier due to the relatively uniform slope.However, minor variations are evident (Fig. 7c), corresponding to variations in gradient, with steeper slopes resulting in increased flow velocity.Higher ice velocities also coincide with regions experiencing reduced rates of ice loss and/or some limited gains in elevation.Similarly, gentler slopes coincide with higher rates of ice loss.
The largest elevation change of 18 m of ice loss was measured at the calving face (Figs. 9 and 10).The greatest measured increase in surface elevation is the result of a large boulder moving down-glacier between the survey dates (Fig. 9), other areas of elevation increase are also primarily associated with the downward movement of large surface features.Apart from the calving face, the most dramatic areas of surface ablation are associated with exposed ice cliffs and melt ponds (Figs. 9 and 11).Figure 11 shows the horizontal movement and expansion of these ice cliffs and the retreat of the calving face (note: the 2015 image has been georectified to remove any horizontal glacier movement based on calculated velocity vectors; see Fig. 7).Rates of horizontal retreat measured for these features range from 2 to 25 m yr −1 with larger cliffs retreating more rapidly.Ap-pearance, expansion and disappearance of melt ponds is also visible in the imagery.
Figure 12 shows the mean and SD of elevation changes measured within 10 m altitude bands over the glacier tongue.Mean elevation change for the lowest-altitude band (4480-4490 m) is the highest at −10.4 m due to collapse of the calving face.Variability (SD) in mean ice loss is also highest for this section.Mean ice loss is greatest for elevations below 4530 m, or roughly the lowest half of the glacier tongue, and generally decreases moving up the glacier tongue.Figure 12 also shows there is some connection between highest mean mass loss and the slope of the glacier (contour spacing).Ice emergence as a result of higher ice velocities over slightly steeper sections may be responsible for the observed slight increases in mean elevation (ice gain) for these steeper slopes (Figs. 7 and 12).
Considerable changes were also measured within the proglacial lake system below the calving face.The lake system can be divided into two parts.The lower section comprises the main (larger) lake and is the section for which the glaciology office has completed bathymetric surveys and volume assessments (Portocarrero, 2014).This section was not surveyed and likely experiences limited change as its water level is controlled by the dam culvert.The upper section between the lake and the calving face includes numerous smaller ponds interspersed with debris piles.Within this upper lake section, numerous changes were observed over the survey period.Noticeable changes in the arrangement and size of the surface ponds occurred, with new ponds forming and others disappearing between the survey dates (Fig. 13).The elevation of these ponds varies greatly, with some sitting up to 4 m above the main lake elevation of 4480 m a.s.l.The largest pond is 15 000 m 2 and lies ∼ 2.5 m above the main lake.Additionally, two locations experienced large elevation changes with surface lowering of 5-10 m recorded over the study period (Fig. 13), suggesting the melting of buried blocks of remnant glacial ice within the sediments.Low velocities of approximately 2.5 m yr −1 were also mea-sured close to the calving face decreasing to ∼ 0.5 m yr −1 at the southern edge of the surveyed area.

Glacier and proglacial lake dynamics
Surface melt patterns over the glacier tongue were highly varied across small spatial scales, though generally ice loss was higher at lower elevations, presumably as a result of slightly higher temperatures (Fig. 12).The largest amount of ice loss occurred at the calving face where a major collapse took place between the two survey dates (Figs. 8  and 9).From speaking with local climbing guides much of this loss occurred over a roughly 2-week period in March-April 2015.Considerable downwasting (roughly 1 to 8 m) occurred in the region behind the calving face, likely re-sulting from the collapse of the calving face that had previously provided a buttressing effect for the ice behind.No notable changes in glacier area or terminus position were measured between the two survey dates despite the large reduction in ice volume (156 000 m 3 of ice loss).This finding illustrates the importance of measuring volume changes and not just areal changes when quantifying glacier melt.Elsewhere, the highest melt rates were associated with ice cliffs and surface ponds (Fig. 11).Exposed ice cliffs lacking the insulation of thick debris cover are exposed directly to intense incoming solar radiation at this high-altitude, lowlatitude setting and thus experience rapid melting and horizontal recession (Fig. 11) (Buri et al., 2016;Sakai et al., 1998Sakai et al., , 2000)).At the Lirung glaciers debris-covered tongue in Nepal, Sakai et al. (1998) observed an ice cliff horizontal retreat rate of 15-44 m yr −1 (4-12 cm day −1 ) and more recently Buri et al. (2016) measured a rate of 11-18 m yr −1 (3-5 cm day −1 ) for the same glacier.This is comparable to the 2-25 m yr −1 observed at Llaca Glacier.Surface melt ponds have a lower albedo than exposed ice, which enhances melt processes and increases melt pond size (Sakai et al., 2000).It is probable that the two features are closely related: a surface melt pond may develop at a low point on the glacier surface where surface meltwater collects; this raises surface albedo and enhances melt, further expanding the pond.As the pond expands and deepens, ice cliffs develop around it; then as they steepen, they lose their debris cover and begin to recede rapidly across the exposed glacier surface.Also, a number of surface melt ponds disappeared between the two survey dates, indicating potential subglacial drainage pathways through englacial conduits (Benn et al., 2012).Immerzeel et al. (2014) identified similar features through analysis of UAV datasets collected over a debris-covered glacier in the Himalaya and associated their presence with ice cliffs.The fine scale of these surface features could not be resolved using medium-and coarse-resolution satellite imagery and would be extremely labour intensive to observe with in situ measurements.
Apart from areas of ice cliffs and melt ponds, the glacier is covered in a thick debris layer.Going back to Ostrem (1959), numerous studies have investigated the role of debris cover in either enhancing ablation through reduced albedo when debris cover is thin or reducing ablation when debris cover is sufficiently thick to have an insulative effect (Juen et al., 2014;Pratap et al., 2015;Shahi et al., 2015).No systematic measurements of debris thickness and distribution were collected for this study.Nevertheless, field observations suggest it was generally in excess of 1 m thick and composed of mixed sediment sizes, ranging from silts and fine sands through boulders over 1 m in diameter.Thus, Llaca's debris cover is probably sufficiently thick that it primarily provides an insulative effect for the ice beneath.The insulative effect of the debris cover is at least in part responsible for its low terminus elevation of 4500 m, which is 200-600 m lower than many of the other glaciers in the region.Addi-  tional factors probably include the large accumulation zone within a high-elevation cirque over 6000 m high; a narrow valley width that restricts the glacier tongue width, increasing length and thickness; and a greater degree of surface shading than in some of the wider valleys.Debris thickness, or the lack there of, thus has a primary control on observed patterns of ablation rates over the Llaca Glacier tongue.Where debris cover is absent (ice cliffs) or melt ponds are present ablation rates are dramatically increased, while in areas of thick debris cover insulation of the ice surface dominates, reducing melt rates.The interplay of these processes across the glacier tongue is responsible for the high degree of spatial heterogeneity observed in glacier melt patterns and the complex surface topography.
Many sections of the glacier tongue show a gain in surface elevation between 2014 and 2015 (Fig. 7).These are not zones of net accumulation but are better explained by the highly irregular nature of the glacier surface.The glacier surface has numerous surface irregularities including topographic undulations, cliffs and large boulders.As the glacier moves downslope, these surface features move through the scene; thus elevation comparisons at an individual point are not actually comparing the same point on the glacier surface (which has moved) but instead compare the same position in real world coordinates.Therefore, if a large boulder or surface feature is displaced down-slope between surveys, it will be measured in the resulting DEM as an increase in surface elevation at its new position.The inverse effect (elevation loss) can occur as large depressions move through the scene.
Examples of this can be seen in Fig. 9f, where a large boulder roughly 10 m high moves down slope, producing a measured elevation decrease at is location in 2014 and a measured elevation increase at its location in 2015.Care must therefore be taken when comparing measurements of vertical change at a specific location.This is not an issue faced by non-geodetic methods such as ablation stakes, which move with the glacier and thus provide a measure of actual ablation at that position on the glacier surface.Vertical emergence may also contribute to the observed increases in ice elevation.There appear to be some minor increases in measured elevation and/or lower ablation in both areas of steeper slope (Fig. 12) and faster velocities (Fig. 7), suggesting some limited emergence could be occurring in these regions, presumably driven by bed topography.However, any emergence is likely to be limited as the glacier bed in the study section is likely smooth and of relatively uniform gradient.In addition, the glacier in this section undergoes no notable direction changes that could produce large emergence velocities.
Very few measurements of glacier surface velocities are reported for the Cordillera Blanca.Emmer et al. (2015) reported maximum values of 11.2 m yr −1 for a similar debriscovered section at the much steeper Jatunraju Glacier, which is less than half what is presented here.However, the relative size of the accumulation area at Llaca is much larger than for Jatunraju.It is also possible that the glacier-bedrock interface is somewhat lubricated through contact with the proglacial lake, potentially increasing horizontal velocity.The observed decrease in velocity down-glacier is essentially linear.Given the fairly constant surface gradient over the glacier tongue (Fig. 7), this velocity decrease is likely due to a reduction in ice mass closer to the glacier terminus resulting from the higher melt rates occurring over the lower sections of the glacier tongue.
The upper section of the proglacial lake was surprisingly dynamic, with a number of major changes occurring between the two survey periods, including the expansion and disappearance of ponds and large areas of subsidence (5-10 m yr −1 ).These changes suggest that sediments within the lake system are interspersed with buried chunks of remnant glacier ice that is slowly melting.As these ice remnants melt, surface and subsurface hydrologic connections can open and close, impacting the size and distribution of surface ponds and changing the surface topography.The variation in pond elevation, offset up to 4 m above the main lake, indicates that they are distinct hydrologic features, separate from one another and from the main lake, as opposed to comprising a single connected lake system.Drainage of these ponds is likely rapid, caused by surface and subsurface channels opening as buried ice bodies melt or through the impacts of waves triggered by talus falls from the lateral moraine and calving of the glacier tongue.Low surface velocities were also measured within this section ranging from 2.5 m yr −1 below the calving face to 0.5 m yr −1 at the terminal end of the survey area where the main lake starts.These low velocities suggest that the sediments below the calving face are still being actively pushed down valley by the glacier.It is also possible that some direct ice connection to the main glacier remains, buried below the surface sediments.
Measured changes within the proglacial lake sediments illustrate that this system is actively evolving as a result of continued ice melt.While no change in the glacier terminus was measured in this study, 18 m of downwasting occurred at the calving face; as this downwasting continues, the proglacial lake volume will increase and the surface area will eventually expand north-eastward.The Huaraz glaciology office estimated the main lake to hold ∼ 274 300 m 3 of water in 2004 (Portocarrero, 2014).From an inspection of the surrounding topography and comparison with other proglacial lakes in the Cordillera Blanca, it appears the lake area could expand some 600 to 1000 m further to the north-east with progressive melt, potentially doubling in area.This larger lake will be in closer proximity to the steep headwalls of the valley, which are probable locations for avalanche and rockfall debris that could potentially trigger hazardous GLOFs.The assessment of Emmer and Vilímek (2014) of the GLOF risk at Llaca defines this as the most likely risk scenario for this lake.Presently, the risk of a large GLOF is limited because any major rockfall and/or avalanche debris would land atop the glacier tongue.However, rapid drainage of perched meltwater ponds (the largest of which is 15 000 m 2 and sits ∼ 2.5 m above main lake) within this upper lake section could produce rapid fluctuations in water inputs to the main lake.Large rapid inflows to the main lake likewise have the potential to impact downstream infrastructure such as the moraine dam, power generator and water diversions in the valley below.

Benefits and challenges of UAVs for glacier mapping
UAVs are capable of providing centimetre-resolution imagery that is able to provide unique insights into spatiotemporal dynamics of the cryosphere.The ease and relatively low cost of collecting highly accurate repeat DEMs makes UAVs a particularly suitable system for the geodetic assessment of glacier mass balance changes.UAVs are especially suited to the study of mountain glaciers where steep topography, frequent cloud cover, surface heterogeneity and generally smaller glacier size (as compared to ice sheets) can limit the use of satellite imagery.Furthermore, UAVs are able to access and accurately survey potentially dangerous areas of the glacier surface (e.g.avalanche paths).In this study we exploited the high spatial resolution of UAV image data to analyse annual surface changes, although UAVs are also suited for conducting rapid return surveys over shorter temporal intervals.It is entirely feasible to complete weekly, daily or even diurnal surveys of the glacier surface, which could provide unique insights into glacier processes that are nearly impossible to obtain with other technologies.Similarly, multispectral, hyperspectral and thermal imaging sensors could be deployed to provide further insights into debris thickness distribution (Schauwecker et al., 2015), surface and melt pond albedo, and surface energy processes (Aubry-Wake et al., 2015).The application of UAVs to glaciology is just beginning, and the potential gains are great.Still, there are a few significant limitations to the methods described in this paper that are addressed below.
To generate high accuracy DEMs, in-scene ground control targets must be installed and surveyed.To minimise distortions in the SfM processing workflow these GCPs should be well distributed across the glacier surface, particularly at the survey margins where "fishbowling" of the DEM can occur due to the propagation of SfM modelling errors.However, installing GCPs across a glacier is extremely labour intensive and in many cases highly dangerous or impossible without considerable mountaineering experience.A possible solution to this is through the integration of real-time kinematic (RTK) GNSS positioning on the UAV platform.If the camera positions are accurately known in real world coordinates (including latitude, longitude elevation, pitch, roll and yaw), then GCPs are unnecessary for derivation of the DEM through SfM processing.These parameters can be added to image metadata by integrating the sensor package and autopilot/GNSS equipment.However, since the UAV is moving quickly (∼ 8 m s −1 ), even a slight delay between camera trigger time and the RTK position can produce significant errors.These may not be an issue for single surveys of the surface but could critically limit accuracy when repeat DEMs are to be used for assessment of surface change.Furthermore, the addition of RTK GNSS to the UAV increases cost and more importantly weight, which reduces flight time.Despite these potential issues, we believe the integration of RTK position-ing on board the UAV platform is a necessary addition for the operational deployment of UAVs for glacier monitoring.Thankfully, RTK systems are rapidly reducing in price, and the entry of numerous low-cost and UAV-specific options to market in recent years will facilitate greater availability of these systems.Another way to reduce DEM error is through survey design.Firstly, overlap and sidelap can be increased, though the latter increases flight line density and thus total survey time.Secondly, a crosshatch survey grid pattern can be employed, which may improve SfM processing accuracy and reduce distortions at the survey periphery (James and Robson, 2014).Nevertheless, deploying a crosshatch grid doubles survey time, which may be unfeasible.
An additional constraint is the maximum possible survey area that is fundamentally limited by the total flight time permitted given power drain and UAV design.Lower air density at higher elevations significantly reduces propeller generated lift and negatively impacts maximum flight time.For example, the multirotor platform used in this study is capable of approximately 25 min of flight time at sea level but only 12 min of flight time at 5000 m.With only 12 min of flight time, the maximum range for this platform is around 1.5 to 2.5 km from the launch point while flying at a suitable survey speed (around 8 m s −1 ) that further depends on local wind conditions and flight altitude above ground.Fixed-wing platforms are capable of considerably longer flight times (> 1 h), but multirotors are easier to fly, can land in a smaller location, can better navigate variable terrain features and are generally more stable in variable winds.Due to the complex terrain of mountain environments featuring often rough landing conditions and frequently strong, gusty winds, multirotors may be the preferred and safer option for many studies.However, the trade-offs between maximum survey area and local conditions should be carefully considered before deciding on a suitable platform.

Conclusions
This study deploys a custom designed hexacopter UAV to map dynamics of a debris-covered glacier tongue in the Cordillera Blanca, Peru.To our knowledge, it is the first time a UAV has been deployed for this purpose in the tropical Andes and is the highest altitude deployment of a multirotor UAV for mapping purposes in the current literature.We completed repeat aerial surveys (23 July 2014 and 28 July 2015) of the glacier tongue and proglacial lake system.Using a SfM workflow, we generated highly accurate 10 cm DEMs and 5 cm orthomosaics.Analysis of these data reveals heterogeneous changes in the glacier surface with some areas of the glacier losing as much as 18 m of vertical elevation due to collapse of the calving face.Average glacier downwasting was −0.75 m over the study period, equating to roughly 156 000 m 3 of total glacier ice loss.The most rapid areas of ice loss were associated with exposed ice cliffs and melt-water ponds on the glacier surface.Surface velocities were recovered from manual feature tracking and ranged from 27 m yr −1 in the upper surveyed section to around 4 m yr −1 at the calving face.Surface lowering and low surface velocities were also measured over the sediments within the proglacial lake, indicating that extensive regions of glacier ice remain buried within the lake sediments that are still connected horizontally to the glacier tongue.Only limited horizontal retreat of the glacier tongue was recorded, affirming that monitoring efforts that deploy only measurements of area changes are unable to document actual changes in glacier ice mass and melt volume.Continued downwasting of the glacier surface and melting of buried ice within the proglacial lake sediments will increase the size of the proglacial lake over time.As the lake boundary migrates towards the valley headwall, GLOF risks associated with avalanches from the cirque above are likely to increase.
This study achieves two primary goals.Firstly, it clearly illustrates the viability and suitability of UAV technology for studying high-resolution changes in mountain glaciers and the feasibility of operating multirotors at high elevation (> 5000 m).Secondly, it provides insights in to the processes that control glacier melt over a debris-covered tropical glacier.Our findings regarding the importance of ice cliffs and surface melt ponds in debris-covered glacier melt processes corroborate findings for debris-covered glaciers in the Cordillera Blanca (Emmer et al., 2015) and the Himalaya (Buri et al., 2016;Immerzeel et al., 2014).The high spatial resolution of UAV-derived datasets facilitates a much better understanding of surface heterogeneity and spatial variability in debris-covered glacier melt processes and changes in proglacial lake systems.UAVs offer a low-cost, on-demand alternative to existing imaging technologies; they can provide unique insights into dynamic earth system processes like glaciers and are ideal for small-scale studies in complex mountain regions.
Data availability.Due to the large size of these imagery, DEM and associated datasets, they are available via FTP upon request from the author at oliver.wigmore@colorado.edu or at www. oliverwigmore.com.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Llaca Glacier location map within the Cordillera Blanca.Showing UAV launch point and survey area boundaries.UAV survey boundary (yellow line) is maximum extent of study area covered by both 2014 and 2015 surveys.The glacier survey boundary (blue line) is the maximum extent of the glacier tongue covered by both 2014 and 2015 surveys.

Figure 3 .
Figure 3. Flight coverage areas, GCPs and check point locations for 2014 and 2015 UAV survey campaigns.Background imagery from Worldview 13 July 2011.

Figure 5 .
Figure 5. Histogram plot of elevation differences (errors) between check points and GCP surveyed elevations and DEM elevations for 2014 and 2015 respectively.Note: negative values indicate underestimation of DEM surface with respect to surveyed positions and thus imply flattening of the DEM surface.Most errors are less than ±5 cm and all within ±10 cm except for 2014 check points.Large 2014 check point errors are assumed to be the result of poor location selection.

Figure 6 .
Figure 6.Box and whisker plot of elevation differences (errors) between check points and GCP surveyed elevations and DEM elevations for 2014 and 2015 respectively.Thick lines represent median value, boxes show upper and lower quartiles and whiskers indicate minimum and maximum values, and circles are outliers.Note: negative values indicate underestimation of DEM surface with respect to surveyed positions and thus imply flattening of the DEM surface.Most errors are less than ±5 cm and all within ±10 cm except for 2014 check points.Large 2014 check point errors are assumed to be the result of poor location selection.

Figure 8 .
Figure 8. Histogram of measured elevation changes over glacier tongue (blue boundary in Fig. 7) from 23 July 2014 to 28 July 2015.Calculated from subtraction of 2014 DEM from 2015 DEM (i.e.negative values are ice loss) at 10 cm DEM pixel size.

Figure 9 .
Figure 9. Locations of maximum and minimum change in ice elevation on Llaca Glacier tongue.Panels (a) and (b) show 2014 orthomosaic; panels (c) and (d) show 2015 orthomosaic; and panels (e) and (f) show elevation change where 2014 elevation is subtracted from 2015 elevation (i.e.negative value is ice loss).Maximum ice loss of 18 m occurred at the calving face (pink circle lower left panel).Maximum gain of 11 m is recorded in upper right section where a large boulder (roughly 10 m tall) has moved through the scene (visible by comparison of 2014 and 2015 orthomosaics).Blue line is glacier tongue boundary; red arrows are velocity measurements.Note: elevation change values are truncated to ±10 m for display purposes.Location of panels is depicted in Fig. 7 inset boxes: 3 -left panels; 1 -right panels.

Figure 10 .
Figure 10.Collapse of the calving face.

Figure 11 .
Figure 11.Rapid ice loss around exposed ice cliffs and surface melt ponds.(a) Images are from upper right of surveyed area; (b) images are behind calving face.Yellow arrows indicate horizontal movement of feature from 2014 to 2015 of ice cliff position: 2-25 m yr −1 .More rapid retreat of ice cliffs is observed in the lower glacier, where cliffs are larger and temperatures are warmer.Note: 2015 images have been shifted (georectified) to 2014 positions based on velocity vectors (Fig.7); i.e. horizontal glacier movement is removed so that the same location relative to glacier surface is compared.Location of panels is depicted in Fig.7inset boxes: 2 -upper panels; 3 -lower panels.

Figure 12 .
Figure 12.Changes in glacier tongue elevation by 10 m elevation band.(a) Mean elevation change within a 10 m elevation band.(b) SD of elevation change within a 10 m elevation band.

Figure 13 .
Figure 13.Locations of major elevation change (ice loss) within lake sediments (pink circles) indicate that buried glacier ice is still present.Low surface velocities were also measured here (∼2.5 m yr −1 below calving face in left panels and ∼0.5 m yr −1 at end of surveyed region in right panels), indicating these sediments are still connected to the main glacier tongue.Panels (a and b) show 2014 orthomosaic; panels (c and d) show 2015 orthomosaic; panels (e and f) show elevation change where 2014 elevation is subtracted from 2015 elevation (i.e.negative value is ice loss).Blue line is glacier tongue boundary.Note: elevation change values are truncated to ±10 m for display purposes.Location of panels is depicted in Fig. 7 inset boxes: 4 -left panels; 5 -right panels.

Table 1 .
Llaca Glacier GNSS ground control survey and hardware specifications.

Table 2 .
Estimated GNSS survey positional errors for use in structure from motion ground control and absolute DEM error analysis.Note: total error estimates include base station and rover error estimates.
5.A dense point cloud is generated from the sparse point cloud model.6.A 3-D model, raster DEM and orthomosaics are output for analysis.

Table 4 .
Estimated pixel matching and DEM reconstruction errors from SfM processing workflow.