Detecting Short-Term Surface Melt on an Arctic Glacier Using UAV Surveys

Abstract

Current understanding of glacier mass balance changes under changing climate is limited by scarcity of in situ measurements in both time and space, as well as resolution of remote sensing products. Recent innovations in unmanned aerial vehicles (UAVs), as well as structure-from-motion photogrammetry (SfM), have led to increased use of digital imagery to derive topographic data in great detail in many fields, including glaciology. This study tested the capability of UAV surveys to detect surface changes over glacier ice during a three-day period in July 2016. Three UAV imaging missions were conducted during this time over 0.185 km${}^2$ of the ablation area of Fountain Glacier, NU. These were processed with the SfM algorithms in Agisoft Photoscan Professional and overall accuracies of the resulting point clouds ranged from 0.030 to 0.043 m. The high accuracy of point clouds achieved here is primarily a result of a small ground sampling distance (0.018 m), and is also influenced by GPS precision. Glacier surface change was measured through differencing of point clouds and change was compared to ablation stake measurements. Surface change measured with the UAV-SfM method agreed with the coincident ablation stake measurements in most instances, with RMSE values of 0.033, 0.028, and 0.042 m for one-, two-, and three-day periods, respectively. Total specific melt over the study area measured with the UAV was 0.170 m water equivalent (w.e.), while interpolation of ablation measurements resulted in 0.144 m w.e. Using UAVs to measure small changes in glacier surfaces will allow for new investigations of distribution of mass balance measurements.

1. Introduction

The negative trend in glacier mass balance is particularly pronounced in the Canadian Arctic, where rapid glacier recession has been documented in the early 21st century [1,2,3,4,5]. Understanding future glacier mass changes in this difficult to access region depends on accurate monitoring and modelling of glacier mass balance. Until recently, research in this field has relied on sparse data points and significant generalizations [6], but new technologies, such as unmanned aerial vehicles (UAVs) and structure-from-motion photogrammetry (SfM), are becoming widely used tools for improving understanding of glacier change [7,8].

Advances in UAV technology has enabled their widespread use in scientific research [9]. They are a particularly attractive tool for glaciological research due to their ability to collect data over wider areas and inaccessible regions of glaciers, with relatively low costs and efforts. As a result, there have been several recent glaciological studies relying on UAV data [7,10,11,12,13,14,15,16,17,18]. Most commonly, these studies are using imagery in the visual spectrum from UAV-mounted digital cameras and SfM to recreate glacier surfaces and track motion. Glaciological studies using UAVs have tracked glacier motion (e.g., [7,10,14]); measured seasonal or multi-year melt (e.g., [7,10,13]); analyzed calving and crevasse dynamics (e.g., [12,15,16]); and mapped surface characteristics such as drainage networks [11], albedo [17], debris cover [18], and cryoconite holes [19].

With quickly advancing technology, the capacity of UAVs as a tool for investigating scientific questions is always growing. Thus far, glaciological studies using UAVs have focused on phenomena that happen on relatively large scales. To the best of our knowledge, Cook et al. [19] is the only example in glaciology which tests the current capacity to detect small (<0.10 m) glacier surface features using UAVs, but no examples of detecting change on the same spatial scale are currently published. A better understanding of small scale processes will improve models and predictions of future glacier change. UAVs are a tool that can provide new insights in these, in ways previously available methods may not be able to.

This study aimed to measure daily glacier surface ablation using a UAV. We achieved this through a case study over a small area in the ablation zone of Fountain Glacier in the Canadian Arctic. There were three main objectives to this study:

• Reconstruct glacier surfaces using UAV imagery and SfM, and compare those surfaces to differential GPS surveys. While other studies have done similar analyses, our methods show significant improvement in measured accuracies over previously published results in glaciological studies.
• Measure glacier surface change over three days through differencing the reconstructed glacier surfaces and compare those measurements to in situ measurements, which has not been previously done in glaciology to our knowledge.
• Calculate total melt measured over study area during the three-day period using both the UAV measurement and ablation stake measurements for comparison.

Conducting the study within the ablation zone avoids the potential influence of snow density changes and firn compaction on surface elevation changes, although we discuss the implications of other influences on glacier surface elevation in Section 5. While we are not measuring mass balance directly in this case study, the methods we describe below will help to improve distributed mass balance measurements through the use of UAVs. Measuring surface melt with UAVs could also provide significant insight into the heterogeneous nature of ice melt on shorter time-scales. This study provides similar results to those found by Nield et al. [20] with terrestrial laser scanning and adds a temporal dimension to detailed glacier surface mapping studies in the existing literature (e.g., [11]).

2. Field Data Collection

2.1. Study Region

Fountain Glacier is located on Bylot Island in Sirmilik National Park, Nunavut, Canada (Figure 1A). It is approximately 1 km wide and 16 km long, stretching from an icefield in the centre of the island, to a dry calving front 9 km from the coast. Fountain Glacier has been part of a larger study program since 1991, which has included investigations in dynamics, hydrology and weather patterns [21,22,23,24]. From 2008 to 2011, terrestrial and aerial imagery were used to look at seasonal and annual movement, as well as melt patterns in the ablation area. Using time lapse photography and traditional photogrammetry, surface lowering in the ablation area was estimated at 2.5 m and 2.9 m in the 2009 and 2010 ablation seasons, respectively [24]. The same study determined horizontal flow rates in fall 2008 to be 6–12 mm day${}^{-1}$. Whitehead et al. [7] used aerial surveys and SfM in 2010 and 2011 to estimate the distribution of annual surface lowering and flow in the ablation area. They found the glacier surface lowered 1.5–2.5 m on average, but with significant heterogeneity. Horizontal annual flow reached up to 8 ma${}^{-1}$ 1300 m from the glacier terminus [7]. In June 2015 an Automatic Weather Station (AWS), equipped with an ultra sonic ranger (SR50) and temperature sensor, was set up in the ablation area of the glacier and an array of 17 ablation stakes were installed in a 0.185 km${}^2$ area to monitor surface melt (Figure 1C). Metal ablation stakes were used in order to last through multiple seasons in the ice. Between June 2015 and June 2016 the weather station recorded an average temperature of $-11.6$ ${}^{\circ }\mathrm{C}$, and 0.30 m of accumulated snow at the end of winter (31 May 2016, after which air temperatures were consistently above freezing).

Fountain Glacier is a polythermal glacier [23]. The bare ice surface in July and August is characterized by several small supraglacial streams, two deeply incised canyons, and a series of thrust planes that are exposed in the lower ablation zone. The ice surface has relatively little debris cover, but scattered rock fall is present in greater concentration on the northern side of the ablation area. Finer debris accumulates along thrust planes, and when pressurized water flow is present in the englacial hydrological system can accumulate more thickly where water flows out from between thrust sheets.

2.2. UAV Description

The UAV used in this study was a MikroKopter Inc. electric rotary hexacopter equipped with a Panasonic Lumix DMC-GF1 camera with a 12 Megapixel Live MOS sensor and Lumix G 20 mm lens (Figure 1B). The onboard inertial measurement unit, GPS and autopilot enables the unit to travel autonomously along a preprogrammed route using Mikrokopter software. The four-cell, 8000 mAh battery provided a flight time of ~17 min.

2.3. Study Area and Ablation Measurements

A study area was established in the lower ablation area of the glacier in July 2016, which encompassed the ablation stake locations and the AWS at the western edge. The area ran roughly 600 m from the northern edge of the glacier south to the edge of a deeply incised supraglacial stream, and 400 m west beginning near the steep frontal cliff (Figure 1C). This area included a steep ramp on the northern edge of the glacier, as well as a gradual slope near the western edge of the area, with elevation ranging from 309 to 360 m. In addition, the area contained numerous supraglacial streams.

The ablation stake network established in 2015 was maintained during 2016, following standard practice, with seventeen 50 mm diameter aluminum poles set to depths of >1.5 m in the ice. Prior to each UAV survey, ablation stake heights were measured to the nearest 0.005 m and the surface position at the base of the stake was measured with the differential GPS. Ablation stakes generally influence ice melt in their immediate vicinity (0.05–0.10 m). To avoid this area of influence, measurements were taken relative to the ice surface outside this zone by levelling a ruler at the top of the stake and measuring 0.10 m away from the stake.

It is difficult to assess the uncertainty in a single ablation stake measurement, because that uncertainty is a combination of human influence and surface variability in the immediate vicinity of the stake [6]. The true variability in melt, however, is relatively small when influencing factors (such as slope, aspect, albedo, and elevation) are constant. Braithwaite et al. [25] reported daily average differences of 10 stakes within 100 m${}^2$ as 5–10 mm water equivalent (w.e.) d${}^{-1}$. We adopted the upper end of this range for our analysis, given that measurements were taken with 0.005 m precision and surface change is calculated from the difference in two measurements.

The tripod mounted AWS was situated near the western edge of the study area. A Campbell Scientific SR50 ultrasonic ranger was also set up on a separate pole drilled into the ice to record average hourly surface position based on 2 min readings. Due to surface roughness at the meter-scale, it is difficult to situate the SR50 over a uniform surface. The ice surface in the field of view of the instrument sloped south east at 12°. Given the uneven surface and the reported accuracy of the instrument (0.01 m), we estimated an uncertainty associated with SR50 measurements of ±0.03 m. A 3-h moving average was used to smooth the SR50 data for comparison to UAV-SfM results in the discussion below.

2.4. UAV Survey Design

Three UAV surveys were conducted, on 21, 23 and 24 July. Each survey used four flights to cover the study area, and was conducted between 13:00 and 16:00 local time. The flights were planned using the MikroKopter flight planning software to maintain at least 60% overlap between images and a flying altitude of approximately 90 m above the ice surface. This resulted in a ground sampling distance (GSD) of 0.018 m. Waypoints were adjusted to follow glacier topography based on a 2011 DEM with 10 m resolution. The UAV flew autonomously between waypoints generated in the flight planning software and uploaded to the onboard flight controller before each flight. During the flights, the camera was triggered every 3 s as the UAV travelled at 5 m s${}^{-1}$ and all images were stored in RAW format to aid later processing. The camera was set with an F-stop of f/4, ISO-400, and an exposure time of 1/1600 s. Each survey produced approximately 1000 images in total between take-off and landing; all these images were used in SfM processing.

A grid of 26 targets was laid out in the study area on 20 July. Each one was a 0.25 m square sheet of corrugated plastic and was fixed to the surface with a 0.20 m metal stake drilled into the ice (Figure 1D). The targets were painted with alternating black and coloured quadrants to aid in locating the target centre during processing. Prior to each UAV survey, the 26 targets were repositioned to be level with the ice surface and surveyed with a Trimble R8 differential GPS. The R8 base station was mounted on the AWS.

3. Data Processing and Calculations

The imagery and GPS data collected during each survey were processed to produce a final orthophoto and 3D ice surface points of the study area for each date. These steps included post-processing the location of the GPS base station using Natural Resource Canada’s Precise Point Positioning (PPP) software; normalizing image exposure in Adobe Lightroom; calculating locations, producing point clouds and orthophotos using Agisoft Photoscan Professional (referred to hereafter as Photoscan); and calculating errors in the final 3D ice surface points using the LAScontrol tool in LAStools and ArcGIS 10.4. These steps are described in the following section.

3.1. GPS

The PPP application provides a precise correction to the base station location using the positions of the receiver logged throughout each survey and GPS satellite positions during the same period [26]. After processing base station data from July 2016 using the PPP tool, the base location was accurate to within 0.001 m and 0.004 m in the horizontal and vertical, respectively. The average precision of the target measurements was 0.016 m (horizontal) and 0.036 m (vertical).

3.2. Imagery

Some studies edit imagery to overcome variations in natural lighting due to clouds, which can complicate image alignment and increase uncertainty in surface reconstruction [27,28]. Capturing the original images in RAW format allowed for corrections to be made to exposure levels without losing detail necessary to optimize image alignment quality. Adobe Lightroom Software was used to correct dark images resulting from clouds and to export all RAW images into TIFF format, which was used in the structure from motion process. Through visual comparison between non-cloudy and adjusted images, the “auto tone” adjustment was found to produce the best match to unaltered images. The “auto-tone” function stretches the histogram of each color channel in the RAW image to match a black-white continuum, alleviating contrast problems in dark images. This adjustment was performed on darker RAW images from all UAV surveys in this study before they were exported to TIFF format.

3.3. Structure-from-Motion

Photoscan was used to align images into a mosaic, reconstruct UAV locations for each image (no GPS data from the UAV was used) and generate surface reconstructions. Photoscan uses SfM techniques described by Snavely et al. [29]. The steps and settings used to reconstruct the glacier surface for this study are outlined below.

The TIFF images from each UAV survey were loaded into a single Photoscan project as separate sub-sections to aid standardizing processing. The photos of each sub-section were aligned using the high accuracy setting, with a 500,000 key point limit and no limit to tie points. In particular, not limiting tie points aids in maximizing pixel matches for low-contrast surfaces, such as snow or ice. This alignment generates a sparse point cloud. The sparse point cloud is a reconstruction of tie point locations in a relative coordinate system.

There was some noise in the sparse point cloud due to reconstruction errors. Using the Gradual Selection tool, all points with a reprojection error greater than 0.5 pixels were removed. After using the Gradual Selection tool, obvious outliers were removed manually; these are primarily isolated points that have an elevation different than the surrounding points greater than 5 m.

After filtering the sparse point cloud, 13 targets were located in the imagery to serve as Ground Control Points (GCPs). The remaining 13 target positions were withheld to test the accuracy of the final reconstructions as ground validation points (GVPs). Because of the overlap of images, GCPs appear in anywhere from 10 to 30 images. Each GCP was manually marked in at least six images, selecting those where the centre of the target was easiest to identify. The maximum reconstruction error at all GCPs was under 0.26 pixels for each survey. The camera positions were optimized after GCP coordinates were loaded.

After optimization, the 3D location of each image pixel was calculated, resulting in a dense point cloud. The dense point cloud was constructed using the high quality and aggressive depth filtering settings. The high quality reconstruction balances reconstruction accuracy with computation time, while aggressive depth filtering removes noise on a relatively smooth surface. Photoscan reports accuracy statistics for all reconstructions, described below in Section 4.1. Finally, a single mosaiced orthophoto was constructed with a pixel size of 0.05 m.

3.4. Point Cloud Accuracy

The accuracy of each dense point cloud was assessed using the GVPs reserved from Photoscan processing. Vertical accuracy was assessed using the LAScontrol tool from LAStools. This tool calculates differences in z-coordinates between a set of control points and a dense point cloud at the x-y location of a control point [30]. The horizontal accuracy was assessed in ArcMap using the orthophoto. The orthophoto generated for each day was loaded into ArcMap along with the control point locations from the same day. Horizontal offsets were measured manually using the measuring tool, at a scale of 1:5. This scale allows for the best identification of the target centre. The RMSE of vertical and horizontal offsets were calculated separately, as well as the total RMSE of both (combined horizontal and vertical). Additionally, Photoscan reports reprojection errors, as well as errors based on control points located in the imagery, which are discussed below.

3.5. Change Calculation

3.5.1. Multiscale Model to Model Cloud Comparison Algorithm

Change between UAV surveys was calculated using the Multiscale Model to Model Cloud Comparison (M3C2) tool in the software CloudCompare [31]. The M3C2 algorithm computes the distance, positive or negative, from a reference point cloud (C${}_1$) to a target point cloud (C${}_2$) at points within the reference point cloud. In addition to the distance, the algorithm calculates a 95% confidence interval, which is used to determine whether the calculated change is significant. Performing the analysis on the point cloud directly removes the influencing factors of transforming the data into a digital surface model, which could otherwise have a large impact on the outcome of the study given the magnitude of changes we measured.

Lague et al. [31] described the M3C2 algorithm in detail; here, we provide only a brief explanation. The calculation described below can be made at every point in C${}_1$ or at a subset of points in C${}_1$. We chose to calculate differences at every point in C${}_1$, due to the the small changes we were attempting to detect. The algorithm first computes a normal direction to C${}_1$ based on points within a radius (D/2) of the calculation point (P; Figure 2A). A cylinder of diameter d is projected onto both clouds along the normal direction and the average position of the points within the cylinder are calculated for both C${}_1$ and C${}_2$ (Figure 2B). The distance between the two average positions is the distance between the clouds at P (M3C2${}_{\mathrm{dist}}$). A positive distance indicates C${}_2$ is above C${}_1$, while a negative distance indicates C${}_2$ is below.

The selection of D and d is dependent on the properties of the clouds themselves. The normal direction will vary with D due to the local roughness of the point cloud, which is a result of the nature of the physical surface as well as noise due to SfM calculations (Figure 2C). If D is too small, M3C2${}_{\mathrm{dist}}$ can be overestimated. It is recommended that D should be >20–25 times the local roughness (R) in order to ensure the M3C2${}_{\mathrm{dist}}$ is not influenced by local surface roughness. Lague et al. [31] suggested that d should be chosen such that C${}_1$ and C${}_2$ have more than four points within the cylinder.

The confidence interval (C${}_{95\%}$) is calculated based on the local roughness of each point cloud within the cylinder ($\sigma$d), the number of points within the cylinder in each point cloud (n), and a registration error (discussed below). The selection of d can influence local roughness and thus uncertainty, particularly where surfaces are locally steep (Figure 2D).

$$C_{95\%}\left(P\right)=\pm 1.96\left(\sqrt{\frac{{\sigma }_1{\left(d\right)}^2}{n_1}+\frac{{\sigma }_2{\left(d\right)}^2}{n_2}}+reg\right)$$

(1)

Lague et al. [31] found that Equation (1) provides a good estimate of confidence as long as n${}_1$ and n${}_2$ > 4. The registration error is set by the user and is related to instrument noise or errors in calculations. It can be thought of as the uncertainty of the point cloud itself.

3.5.2. Parameter Choices

The parameters D and d were chosen based on roughness calculations of the 21 July point cloud at a variety of scales (

Table 1. Roughness calculations of 21 July dense point cloud.

Kernel (m)	$\overline{\mathit{R}}$ (m)	$\mathit{\sigma }\mathit{R}$ (m)	$\frac{\mathit{D}}{\overline{\mathit{R}}}$	$\frac{\mathit{D}}{\overline{\mathit{R}}+\mathit{\sigma }\mathit{R}}$
0.10	0.003	0.003	29.16	15.380
0.15	0.004	0.004	34.10	17.005
0.20	0.006	0.006	35.30	17.018
0.25	0.007	0.008	35.32	17.016
0.30	0.009	0.009	35.08	16.828
0.40	0.017	0.013	34.54	16.508
0.50	0.015	0.016	34.23	16.339

). The roughness tool in CloudCompare was used with the kernel size set to various possible values of D. In each case, D was greater than 25 times the mean roughness ($\overline{R}$). The final value of D was chosen to maximize the ratio of D to the mean roughness plus one standard deviation ($\frac{D}{\overline{R}+\sigma R}$). The mean surface density of points was then calculated for the same kernel size. A value of d was chosen such that at least five points would lie within the cylinder based on the density less one standard deviation. The final values for D and d were 0.2 m and 0.1 m, respectively. With a value of D = 0.2 m, the condition suggested by Lague et al. [31] was satisfied for 99.2% of core points.

The registration error is related to errors in SfM reconstruction and GPS measurement uncertainty. This stems from manual identification of ground control points, as well as small stochastic differences in camera alignment in Photoscan. The registration error was assigned based on the mean of the total (vertical and horizontal) RMSE (described in Section 3.4) from both C${}_1$ and C${}_2$.

These parameters were used to calculate differences between point clouds for 21–23 July, 23–24 July, and 21–24 July. In each case, the reference point cloud was the earlier time period.

The M3C2 tool output includes M3C2${}_{\mathrm{dist}}$, as well as n${}_1$, n${}_2$, C${}_{95\%}$, and a significance value. The significance value is 1 where $|{\mathrm{M}3\mathrm{C}2}_{\mathrm{dist}}|>{\mathrm{C}}_{95\%}$ or otherwise 0. These outputs were used to filter outliers and noise out of the difference clouds. First, all points where n${}_1$ or n${}_2$ < 5.0 were removed. Then, a mean and standard deviation were calculated for the remaining M3C2${}_{\mathrm{dist}}$ values. Points where M3C2${}_{\mathrm{dist}}$ were more than two standard deviations from the mean were also removed. The remaining points were included in the analysis below.

The filtering process of 21–23 July resulted in removal of 53% of points (Figure 3A); 117.5 million points were retained (Figure 3B). Of the remaining points, 61% had a significance of 1 (M3C2${}_{\mathrm{dist}}>$ C${}_{95\%}$). The points that were filtered out for being more than two standard deviations from the mean were clustered on the edges. Similar results were seen after filtering the other differenced point clouds, with the exception that 97% of points had a significance of 1 for the 21–24 July point cloud. The average C${}_{95\%}$ for all three time windows was also similar (0.082, 0.075, 0.070).

3.5.3. Accuracy of Calculated Change

To assess the calculated change as a measure of surface ablation, M3C2${}_{\mathrm{dist}}$ values were compared to measured ablation with the SR50 and each ablation stake. Values of M3C2${}_{\mathrm{dist}}$ were averaged within a 0.25 m radius around each ablation stake, as well as the SR50, for each time period. The uncertainty of average M3C2${}_{\mathrm{dist}}$ values is reported as the RMSE of all ablation stake locations for a given time window.

4. Results

4.1. Point Cloud Accuracy

The 21 July point cloud with RGB colouring can be seen in Figure 4A. A profile view in Figure 4B shows the three surface reconstructions across a supraglacial stream; the three discrete surfaces indicate that the method employed here retrieved locally robust surface reconstructions across UAV surveys. Accuracy assessments of the point clouds from each time window agree with qualitative assessments of the profiles (

Table 2. Reported and calculated errors of point clouds and combined error of point cloud pairs. Reprojection error and GCP RMSE are calculated within Photoscan, calculation of GVP RMSE is described in Section 3.4. The registration error is the mean of the total RMSE from two point clouds.

	21 July	23 July	24 July
Reprojection error (pixels)	0.666	0.556	0.639
Reported GCP RMSE (m)
XY	0.037	0.029	0.032
Z	0.015	0.031	0.016
Total	0.040	0.042	0.036
GVP RMSE (m)
XY	0.038	0.038	0.033
Z	0.036	0.047	0.027
Total	0.037	0.043	0.030
Registration Error (m)	21–23 July	23–24 July	21–24 July
	0.040	0.036	0.033

). The mean x-y offsets between surveyed GVPs and GVP locations in orthomosaics was 0.033 m for 21–24 July, and 0.038 m for both 21–23 July and 23–24 July. These values are similar to reported values in Photoscan. Offsets in the z-direction are higher than reported by Photoscan, but similar to x-y offsets. The vertical offset for 21–23 July is 0.036 m, for 23–24 July is 0.043 m, and for 21–24 July is 0.030 m.

4.2. Accuracy of Calculated Change

Average changes calculated at ablation stake locations are shown in Figure 5 and

Table 3. Summary of surface change (m) from point cloud differences and in situ measurements at ablation stakes. Residuals are the differences between M3C2 and stake measurements at each stake location.

	M3C2 Results		Ablation Stakes		Residuals
	Mean M3C2 ${}_{\mathrm{dist}}$	Std Dev	Mean	Std Dev	Mean	RMSE
21–23 July	$-0.108$	$0.026$	$-0.105$	0.017	$-0.004$	0.034
23–24 July	$-0.082$	$0.032$	$-0.056$	0.014	0.026	0.044
21–24 July	$-0.180$	$0.028$	$-0.160$	0.046	0.030	0.048

. In most cases, change measured with ablation stakes and M3C2${}_{\mathrm{dist}}$ agree within the uncertainty ranges for the all three time periods. For 21–23 July, the mean residual between ablation stakes and M3C2${}_{\mathrm{dist}}$ was $-0.004$ m, while the RMSE was 0.034 m. For 23–24 July, the mean residual between ablation stakes and M3C2${}_{\mathrm{dist}}$ was 0.026 m, while the RMSE was 0.044 m. For 21–24 July, the mean residual between ablation stakes and M3C2${}_{\mathrm{dist}}$ was 0.030 m, while the RMSE was 0.048 m.

Calculations of M3C2${}_{\mathrm{dist}}$ within 0.25 m of the SR50 location were averaged for comparison with measured surface change during each time period. The average M3C2${}_{\mathrm{dist}}$ was $-0.051$± 0.034 m for 21–23 July, $-0.178$± 0.044 m for 23–24 July, and $-0.237$± 0.048 m for 21–24 July. Surface ablation at the SR50 (calculated using a 3 h average around noon for each day—the approximate time of GPS surveys) was $-0.088\pm 0.03$ m between 21 and 23 July, $-0.059\pm 0.03$ m between 21 and 24 July, and $-0.147\pm 0.03$ m between 21 and 24 July.

4.3. Total Melt

Total specific melt over the study area was calculated through an inverse distance weighted (IDW) interpolation of ablation stake data, as well as using the M3C2${}_{\mathrm{dist}}$ values (

Table 4. Total specific melt in m w.e. calculated for the study area using the UAV-SfM measurements and ablation stake measurements. Ice density was assumed to be 917 kg m${}^{-3}$.

	UAV-SfM	Ablation Stake
21–23 July	0.092	0.094
23–24 July	0.079	0.050
21–24 July	0.170	0.144

). It is common practice to assign melt measured at ablation stakes based on elevation bands of 50–100 m when interpolating ablation stake data over an entire glacier. Given the small study area and lack of elevation range between stakes, a simple distance interpolation method was preferred here. IDW is a common and straightforward local spatial interpolation method based on the principle of spatial autocorrelation of sample data points by distance. IDW assumes a smooth underlying surface and is considered suitable for quick interpolation from sparse, non-clustered data points [32]. The primary limitation of IDW is that there is no assessment of error in the resulting surface. As with any interpolation, IDW can potentially over- or underestimate values that are influenced by factors other than distance. However, exploratory analysis of potential contributors to melt (i.e., elevation, aspect, etc.) revealed no underlying trends, suggesting that more complex geostatistical interpolation techniques are not appropriate here.

The total specific melt calculated with each method are internally consistent (i.e., melt for 21–23 July and 23–24 July sum to the calculated melt for 21–24 July). The totals are also very similar between measurement methods for all time periods. The greatest difference is seen between measurements for 23–24 July, where the ablation stake and UAV-SfM methods differ by 0.029 m w.e. over the study area.

5. Discussion

Previous studies using UAVs and SfM in glaciology report vertical accuracies of 0.1 m to 0.3 m [7,10,33]. The accuracies of dense point clouds constructed in this study are on the order of a few centimetres. Using the M3C2 algorithm, we were able to measure glacier surface change over one-, two-, and three-day time windows. RMSE values were smaller than detection thresholds reported in previous work (e.g., [13,34,35,36]). It should be noted that the methods presented here differ from previous studies, particularly in the GSD achieved with the camera set-up and flying altitude. The goal in this study was to achieve high vertical accuracy, where other studies may have prioritized horizontal accuracy or areal coverage of imagery.

The change measured using the UAV-SfM method agrees with change measured at ablation stakes in the majority of instances. The locations where ablation stake and M3C2 measurements differ was investigated using a correlation coefficient, including analysis of slope, aspect, spatial distribution across the area, but no clear pattern emerged. In addition, albedo was assessed qualitatively at each ablation stake using the orthomosaic, but again no clear relationship was present. This leads us to believe that differences are primarily a result of the SfM reconstruction process and uncertainties in GCP positioning. The magnitude of difference between SR50 measurements and M3C2 measurements is likely higher (compared to ablation stake locations) due to poor ground control placement, which leaves the surface in the vicinity of the AWS poorly constrained [37].

Ablation is not the only process that can affect surface positions at the temporal and spatial scales examined in this study. Glacier dynamics and hydrologic changes can also raise and lower glacier surfaces on short time scales [24,38]. Whitehead et al. [24] measured flow and surface elevation changes on Fountain Glacier using ground-based photogrammetry. They estimated that in the summer of 2009 and 2010 ablation outpaced dynamic uplift approximately 5:1, with ice flow from up glacier raising the ice surface level in the terminus region by 0.2 m over the ablation season. If the rate of dynamic uplift was similar in 2016, this would amount to approximately 2 mm of surface uplift daily and potentially causing a small underestimation of ablation using the methods described here.

In 2009, uplift of the surface of Fountain Glacier associated with the presence of pressurized subglacial water was measured in both February and July [38]. Uplift was tied to a large lake drainage up glacier and measured at 0.05 to 0.30 m over a period of six days in July, with a steady lowering after the event over 2–3 days. In 2016, we noted the presence of pressurized water spilling onto the glacier surface between 13 and 20 July in two places indicating the presence of pressurized subglacial water prior to the study period discussed in this study. While uplift could not be measured for this period, there is no evidence for a drainage event of the same magnitude as that observed in 2009. Assuming the pressurized water was from a smaller magnitude event, and less uplift occurred prior to 20 July, settling of the surface may have contributed 1–2 cm of surface lowering daily to M3C2 measurements.

Total melt calculated from ablation stake measurements and M3C2 measurements agree remarkably well. The internal consistency between time periods (mentioned above) is expected with ablation stake measurements, given the total melt each day is calculated as a difference from the day before. However, the three surfaces measured with UAV-SfM are independent of each other and the internal consistency is interpreted as an indication the measurements are accurate. If the 23–24 July measurement was high (which could be suspected given the uncertainties in surfaces and the higher value compared to ablation stakes), the sum of one-day and two-day totals should be greater than the measured melt for the three-day period. The differences between methods (<0.03 m) are less than the RMSE values of the UAV-SfM method.

Given the accuracy of measurements from UAV imagery, it is possible to examine spatial patterns of ablation over the study area. Over a one-day period, ablation was relatively uniform throughout the centre of the study area, with higher ablation to the west (which was an east facing slope; Figure 6A). For example, patterns begin to emerge over the two-day time window, where differences between north- and south-facing stream banks can be seen (with high ablation on the south facing banks; Figure 6B). In addition, in the two-day period, enhanced ablation was observed along thrust planes (which run north-south in Figure 6) where debris was concentrated on the surface, decreasing albedo and enhancing melt. Over the three-day period, higher ablation can be seen overall (Figure 6C). The differences between north- and south-facing banks was even more apparent, as was higher ablation along thrust planes. Similarly, the two- and three-day periods show enhanced melt in two additional regions with higher amounts of supraglacial debris. During the full three-day period, the difference between an active (higher ablation) and abandoned (lower ablation) channel can also be seen clearly in the southwest of the study area. The lowering of the active stream was primarily due to thermal erosion of running water and the presence of water in the stream channels adds uncertainty to ablation measurements with the UAV-SfM methods. However, streams account for 0.35% of the study area and as such contribute less uncertainty to the melt totals than other factors discussed.

In general, the strongest influence on point cloud accuracies come from the GSD and GPS unit accuracy [39]. Lower flying heights and larger camera sensors both lead to smaller GSD, which can improve automatic tie-point location in the SfM process and decreases errors in manually locating GCP target centres [40]. In this study, we have prioritized a small GSD in our methodology to achieve the high accuracy reported. To carry a camera with a larger sensor, however, UAVs use more power. Ground coverage of imagery from UAV flights is dependent on both flying height and battery life, where higher flights and longer battery life can lead to greater spatial coverage of a single flight. In the case where small GSD is important, multiple flights are needed to cover the same area as a single flight at higher altitude.

The similarity of point cloud accuracies and RMSE values of M3C2 measurements at ablation stakes suggests that point cloud accuracy is a controlling factor for the accuracy of change measurements, and as a result controls the magnitude of measurable change with this method. In short, the point cloud accuracy is effectively a detection threshold which stems from GPS accuracy and GSD, which both influence GCP placement in the SfM process. Image texture plays an important role in SfM reconstruction accuracy [41] and the bare ice surface in this study provides more texture than snow covered glacier surfaces, thus improving the reconstruction process. The detection threshold will determine at what time scales ablation can be confidently detected using UAV data. During this study period, average daily ablation measured at stakes was 0.048 m, which allowed for change measurements over a one-day period. When applying these methods to other locations and time periods, the scale of ablation being measured must be balanced with survey set-up to achieve the proper accuracy in resulting point clouds.

6. Conclusions

Unmanned aerial vehicle surveys were conducted to reconstruct a glacier surface, with accuracies higher than previously reported in the literature allowing for small changes to be measured. Surface changes were calculated by differencing point clouds from repeat surveys. Average M3C2${}_{\mathrm{dist}}$ at 17 ablation stake locations was $-0.108$ m, $-0.082$ m, and $-0.180$ m for the 21–23 July, 23–24 July, and 21–24 July time windows. UAV-SfM derived surface lowering and ablation stake measurements agree in most instances, with RMSE values of 0.034 m, 0.044 m, and 0.048 m for 21–23 July, 23–24 July, and 21–24 July. Total melt for the study area was also calculated using each method and total values are both internally consistent and in agreement with each other.

The RMSE at ablation stake locations was on the same order as accuracies of the point clouds for all three periods, indicating photogrammetric reconstruction is the main control on the accuracy of the method. The achievable accuracy of point cloud reconstruction must therefore be balanced with desired outcomes in future studies.

High point cloud accuracies were achieved through a combination of high-precision GPS surveys of ground control points and small ground sampling distances (0.018 m). The acquisition of the data with such high accuracy is time consuming, however, due to the ground coverage of a single flight at low flying altitude and the necessity of GPS surveys. Thus, it is important to balance the detail obtained in UAV surveys with the time involved in the data collection.