Journal topic
The Cryosphere, 14, 549–563, 2020
https://doi.org/10.5194/tc-14-549-2020
The Cryosphere, 14, 549–563, 2020
https://doi.org/10.5194/tc-14-549-2020

Research article 12 Feb 2020

Research article | 12 Feb 2020

# Surface melt and the importance of water flow – an analysis based on high-resolution unmanned aerial vehicle (UAV) data for an Arctic glacier

Surface melt and the importance of water flow – an analysis based on high-resolution unmanned aerial vehicle (UAV) data for an Arctic glacier
Eleanor A. Bash and Brian J. Moorman Eleanor A. Bash and Brian J. Moorman
• Department of Geography, University of Calgary, Calgary, Alberta, Canada

Correspondence: Eleanor A. Bash (eleanor.bash@gmail.com)

Abstract

Models of glacier surface melt are commonly used in studies of glacier mass balance and runoff; however, with limited data available, most models are validated based on ablation stakes and data from automatic weather stations (AWSs). The technological advances of unmanned aerial vehicles (UAVs) and structure from motion (SfM) have made it possible to measure glacier surface melt in detail over larger portions of a glacier. In this study, we use melt measured using SfM processing of UAV imagery to assess the performance of an energy balance (EB) and enhanced temperature index (ETI) melt model in two dimensions. Imagery collected over a portion of the ablation zone of Fountain Glacier, Nunavut, on 21, 23, and 24 July 2016 was previously used to determine distributed surface melt. An AWS on the glacier provides some measured inputs for both models as well as an additional check on model performance. Modelled incoming solar radiation and albedo derived from UAV imagery are also used as inputs for both models, which were used to estimate melt from 21 to 24 July 2016. Both models estimate total melt at the AWS within 16 % of observations (4 % for ETI). Across the study area the median model error, calculated as the difference between modelled and measured melt (EB =−0.064 m, ETI =−0.050 m), is within the uncertainty of the measurements. The errors in both models were strongly correlated to the density of water flow features on the glacier surface. The relation between water flow and model error suggests that energy from surface water flow contributes significantly to surface melt on Fountain Glacier. Deep surface streams with highly asymmetrical banks are observed on Fountain Glacier, but the processes leading to their formation are missing in the model assessed here. The failure of the model to capture flow-induced melt would lead to significant underestimation of surface melt should the model be used to project future change.

1 Introduction

The Canadian Arctic Archipelago holds approximately 14 % of the world's area of glacier ice outside the major ice sheets, and rates of glacier melt in the region have increased since the late 1990s . show that recent melt rates on Canadian Arctic ice caps are the highest in 4000 years, while found that glaciers of the southern Canadian Arctic Archipelago contributed over 16 % of 2003–2010 sea-level rise. Future projections indicate that the glaciers from across Arctic Canada will continue to be the largest mountain glacier contributors to sea-level rise through the end of this century . Given the importance of Canadian Arctic glaciers to global sea-level rise and the rapid change in melt rates observed in recent years, it is critical to better understand the processes contributing to melt rates on these glaciers.

Direct measurements of melt rates on Arctic glaciers are scarce; only five glaciers from the Canadian Arctic have current data available online (WGMS2018). Where in situ data are collected, they are often based on ablation stake networks and data from automatic weather stations (AWSs; e.g. Bash et al.2018; WGMS2018). These in situ measurements must be extrapolated to provide estimates of melt at other locations on a glacier or at other glaciers where no measurements exist. Several kinds of models are commonly used to extrapolate glacier surface melt measurements, including the temperature index (TI; e.g. Hock1999, 2005), enhanced temperature index (ETI; e.g. Irvine-Fynn et al.2014; Bash and Marshall2014), and energy balance (EB; e.g. MacDougall and Flowers2011; Shaw et al.2016). The TI or ETI models are often preferred for regional-scale models because of their lower computational needs and the ease with which input variables can be estimated. Energy balance models require more inputs and as a result are preferred for in-depth studies of individual glaciers.

A lack of data for Arctic glaciers also makes it difficult to validate model results. Datasets are often split into training and validation periods to assess model performance; however, these measurements represent only a few locations on the glacier surface . Where outlet streams are measured, studies use total stream discharge and aggregated surface melt to validate results . Validation with total melt provides insight into the average performance of these models over an entire glacier, but neither this method nor the validation based on point data provides insight into model performance across a glacier surface.

Technological advances have allowed for increasingly detailed change detection from imagery obtained by satellites, unmanned aerial vehicles (UAVs), and terrestrial photography. Structure from motion (SfM) is widely used across a range of disciplines for reconstructing topography using imagery from the aforementioned sources . Employing SfM, use multiple reconstructions of a debris-covered glacier to measure seasonal ice cliff retreat from terrestrial photographs, including spatial variation in rates across individual cliff faces. Using similar methods and UAV imagery, measured spatially variable melt rates in the ablation zone of an Arctic glacier over 4 d. This new capacity for measuring change in high temporal and spatial detail provides opportunities for examining spatial patterns in ways that were previously not possible.

The methods of surveying that are needed for detailed surface change measurements using SfM are time-consuming . These surveys involve multiple site visits, measurement of control points using differential GPS or total stations, and imagery collection . Given the effort involved in these data collection campaigns, it is not feasible at this time to extend these methods to regional-scale glacier change detection. Modelling efforts will continue to play an important role in understanding glacier melt in the future, and thus it is critical to improve the models employed.

The aim of this study is to use melt measurements available in high spatial and temporal resolution to assess the performance of an energy balance model and an enhanced temperature index model across the glacier surface. We do this using previously published surface melt data which were measured using UAV surveys in the ablation zone of Fountain Glacier, Nunavut .

2 Methods

## 2.1 Study area

Fountain Glacier, located on Bylot Island, Nunavut (Fig. 1), has been studied in detail since 1991 . The glacier stretches across 16 km, from higher elevations in the Byam Martin Mountain Range to its terminus roughly 10 km inland from the coast. The focus area of this study, in the lower ablation zone, spans a narrow elevation range (250–400 m) and terminates in a cliff face with two calving fronts. The lower portion of the glacier faces east, with several supraglacial streams of varying sizes . The largest of these streams form deeply incised asymmetrical canyons, with a steeper north-facing valley wall and a more gently sloping south-facing valley wall (Fig. 1).

Figure 1Location of Fountain Glacier on Bylot Island, Nunavut, in the Canadian Arctic. (a) The automatic weather station (AWS) installed on a tripod, measuring incoming and outgoing shortwave radiation and temperature. An SR50 measuring surface position was installed on a separate pole drilled into the glacier (foreground; photo credit: E. A. Bash). (b) A 2011 orthomosaic of the lower ablation zone of Fountain Glacier . The boundary of the study area is shown, which extends to the 2016 glacier boundary on the north side. The AWS location is also indicated (green triangle). Two canyons formed by stream incision are indicated by red arrows. (c) An aerial view of the largest incised canyon on Fountain Glacier. The photo is taken facing east, and the steep north-facing canyon walls can be seen in contrast to the more gently sloping south-facing walls (photo credit: B. J. Moorman).

During July 2016 aerial surveys were conducted with a UAV to reconstruct 0.185 km2 of the glacier surface multiple times . measured surface lowering between 21 and 24 July using point cloud differencing across the study area indicated in Fig. 1. Concurrently, an AWS (Fig. 1) recorded surface melt with a Campbell Scientific SR50 sonic ranger as well as temperature, incoming and outgoing shortwave radiation, relative humidity, wind speed, and direction. assessed the surface lowering measured through point cloud differencing against 17 ablation stakes, as well as the AWS data, and found that the measured surface lowering agreed with other melt measurements throughout the study area. In this study we assume that surface lowering derived from point clouds is equivalent to surface melt. The uncertainty of this surface lowering was 0.048 m.

measured an average daily melt rate of 0.060 m d−1 across the study area between 21 and 24 July 2016. Between 2010 and 2011, measured melt rates of 0.030–0.055 m d−1 across the ablation zone. also showed, however, that melt was highly variable across the study area during July 2016, ranging from 0.010 to 0.110 m d−1. The authors found that these differences in melt were not tied to the elevation or aspect of the glacier surface.

## 2.2 Description of data

All AWS measurements were taken at 2 min intervals and recorded as hourly averages from 1 to 31 July 2016 (Fig. 2). The SR50 was installed on a pole drilled into the glacier surface immediately next to the AWS on 13 July and measured hourly average surface height. The pole holding the SR50 began to tilt due to melt out of the pole during the afternoon of 27 July and was reinstalled in the morning on 28 July. The data from that time period were removed from the time series for model validation.

Figure 2Hourly averaged melt model input variables, measured at the AWS for the period 1–31 July 2016. Variables include near-surface air temperature (a), incoming (SWin) and outgoing (SWout) shortwave radiation (b), albedo (c), wind speed (d), and relative humidity (e). Inputs which are estimated for distributed models are shown in red for incoming shortwave radiation (a) and albedo (b). The albedo derived from UAV imagery is scaled using the average value at the AWS for 21–24 July so that it remains constant throughout the model. The dates of the UAV study are shown in grey.

Models run at high temporal resolution can have significant errors which stem from SR50 readings . These errors stem from the uncertainty in SR50 readings due to the instrument accuracy (0.01 m) as well as uneven topography underneath the instrument. In addition to hourly surface height, the SR50 records a standard deviation of all measurements in the hourly average (typically 30 measurements). SR50 values where the standard deviation was greater than 0.01 m were removed, expected hourly melt rates are lower than the instrument accuracy, and standard deviations greater than that are assumed to be related to noise. After removing those points, a 5 h rolling average was calculated to smooth the data and fill in removed values. The 5 h window was used to assure that all windows had multiple measurements to average after some measurements had been removed. SR50 values which differed significantly from others in the 5 h window were also removed (defined as more than the standard deviation of hourly change in the 5 h window) under an assumption that it was physically unlikely for melt rates to change dramatically over short periods. Melt was then calculated in 12 h periods (00:00–12:00 LT) using the remaining values.

Processing of UAV imagery was done in Agisoft PhotoScan and is described in detail by . High-accuracy settings were used to retain detail in the final point cloud. During processing, point clouds of the glacier surface at midday were produced for 21 and 24 July 2016. The point clouds were used to create orthomosaics and digital surface models (DSMs) for each day with horizontal resolutions of 0.10 m.

The multiscale model to model cloud comparison (M3C2) was used to calculate melt across the study area at 0.02 m horizontal resolution , with vertical uncertainty of 0.048 m based on the root-mean-square error (RMSE) of surface lowering compared to ablation stakes. The M3C2 algorithm calculates change normal to the local surface, which was desirable for measuring detailed change at such fine resolution. After differencing the point clouds were filtered to remove large outliers and differences calculated from fewer than five points, as suggested by . This filtering resulted in variable density of points in the final point cloud. To align with DSMs and orthomosaics of the study area, for this study melt data were gridded to 0.10 m resolution. Melt measurements were averaged within each grid cell, and empty cells were filled using bilinear interpolation.

## 2.3 Model formulation

### 2.3.1 Energy balance model

An EB model accounts for all energy inputs and outputs at the glacier surface to estimate energy available for melt or freezing. As a representation of surface energy exchange, these models are often assumed to be the best estimators of surface melt (Hock2005). We employ an energy balance model described by as a check on both SfM-derived change and the ETI model described below:

$\begin{array}{}\text{(1)}& {Q}_{\mathrm{N}}={Q}_{\mathrm{SW}}^{↓}-{Q}_{\mathrm{SW}}^{↑}+{Q}_{\mathrm{L}}^{↓}-{Q}_{\mathrm{L}}^{↑}+{Q}_{\mathrm{H}}+{Q}_{\mathrm{E}},\end{array}$

where QN is the net energy flux at the surface and ${Q}_{\mathrm{SW}}^{↓}$, ${Q}_{\mathrm{SW}}^{↑}$, ${Q}_{\mathrm{L}}^{↓}$, ${Q}_{\mathrm{L}}^{↑}$, QH, and QE represent incoming and outgoing shortwave radiation, incoming and outgoing longwave radiation, and sensible and latent heat flux, respectively. The energy fluxes are all calculated in megajoules per square meter per hour (MJ m−2 h−1). For consistency with AWS measurements, and energy fluxes are positive when they contribute energy to the glacier surface. When QN>0 melt ($\stackrel{\mathrm{^}}{M}$; m ice) is calculated by

$\begin{array}{}\text{(2)}& \stackrel{\mathrm{^}}{M}=\frac{{Q}_{\mathrm{N}}}{{\mathit{\rho }}_{\mathrm{ice}}{L}_{\mathrm{f}}},\end{array}$

where ρice is the density of glacier ice, assumed to be 900 kg m3 , and Lf is the latent heat of fusion. Longwave radiation (QL) is parameterised based on the absolute temperature (T) and emissivity (ε):

$\begin{array}{}\text{(3)}& {Q}_{\mathrm{L}}=\mathit{\epsilon }\mathit{\delta }{T}^{\mathrm{4}},\end{array}$

where δ is the Stefan–Boltzmann constant. For ${Q}_{\mathrm{L}}^{↑}$, the calculation is straightforward if the surface is assumed to be melting (273.15 K); the surface emissivity is then 1.0. The atmospheric emissivity (εa) must be estimated in order to obtain ${Q}_{\mathrm{L}}^{↓}$. Here we use a parameterisation from based on vapour pressure (ev) and relative humidity (H), which was found to be robust when transferred directly without local training:

$\begin{array}{}\text{(4)}& {\mathit{\epsilon }}_{\mathrm{a}}=\mathrm{0.445}+\mathrm{0.0055}H+\mathrm{0.0052}{e}_{\mathrm{v}}.\end{array}$

The sensible and latent heat fluxes are calculated by Eqs. (5) and (7), respectively:

$\begin{array}{}\text{(5)}& {Q}_{\mathrm{H}}={\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{p}}{k}^{\mathrm{2}}v\left[\frac{{T}_{\mathrm{pa}}-{T}_{\mathrm{ps}}}{\mathrm{ln}\left(z/{z}_{\mathrm{0}}\right)\mathrm{ln}\left(z/{z}_{\mathrm{0}\mathrm{H}}\right)}\right],\end{array}$

where ρa is the air density, calculated from the near-surface air temperature and average air pressure at the AWS (Paws=969), cp is the heat capacity of air, v is the measured wind speed, and k=0.4 is von Karman's constant. The potential temperature is calculated for the near-surface air (Tpa) and ice surface by (Tps)

$\begin{array}{}\text{(6)}& {T}_{\mathrm{p}}=T{\frac{{P}_{\mathrm{ref}}}{{P}_{\mathrm{aws}}}}^{R/{c}_{\mathrm{p}}},\end{array}$

where Pref is a reference pressure (1000 mbar), R=288.5, and the surface is assumed to be at the melting point.

$\begin{array}{}\text{(7)}& {Q}_{\mathrm{E}}={\mathit{\rho }}_{\mathrm{a}}{L}_{\mathrm{v}}{k}^{\mathrm{2}}v\left[\frac{{q}_{\mathrm{a}}-{q}_{\mathrm{s}}}{\mathrm{ln}\left(z/{z}_{\mathrm{0}}\right)\mathrm{ln}\left(z/{z}_{\mathrm{0}\mathrm{E}}\right)}\right],\end{array}$

where Lv is the latent heat of vaporisation, qa is the specific humidity, based on T and H, and qs is the specific humidity calculated from the saturation vapour pressure at the melting surface. z0, z0H, and z0E are the roughness length scales for turbulent exchange of momentum, heat, and moisture. The height of all AWS measurements is z=1.5 m, and z0=0.005 m was used as a tuning parameter at the AWS. Both Eqs. (5) and (7) neglect a correction for atmospheric stability, which has been shown in previous studies to provide better results for these parameters .

A 1-D formulation of this model (EBaws) runs using measured inputs from the AWS. A distributed version was also developed (EBdist) using modelled incoming radiation at the AWS location and albedo derived from the UAV imagery to determine SWnet (both described below; Figs. 2b and c and 3). Near-surface air temperature across the study area was modelled using the uncorrelated AWS temperature and modelled radiation (Sect. 2.3.2). Measured values from the AWS are used for all other variables in Eqs. (3)–(7).

Distributed surface albedo

Figure 3Albedo was derived from UAV imagery using average values from RGB bands and scaling the lower and upper bounds to match the albedo of bare rocks (based on ) and the average albedo measured on Fountain Glacier during the study period. (a) 21 July albedo values. (b) 23 July albedo values. (c) A hillshade model of the study area based on the DSM from 21 July.

Distributed incoming radiation was modelled by modifying the hourly measured incoming radiation at the AWS by the slope (S) and aspect (A) of each grid cell in the 21 July DSM (Fig. 3c). The terrain correction was based on :

$\begin{array}{}\text{(8)}& {\mathrm{SW}}_{\mathrm{in}}={Q}_{\mathrm{SW}}^{↓}\cdot \mathrm{cos}\mathit{\psi }\cdot \mathrm{cos}\left(\mathit{\omega }-A\right)\cdot \mathrm{sin}S\cdot \mathrm{cos}S,\end{array}$

where ψ is the solar altitude and ω the solar azimuth. This correction alone causes radiation to drop to zero overnight, which is not realistic during Arctic summers. Although direct radiation may drop to zero at some cells when sun angles are low overnight (solar elevation dips to 3.6), diffuse radiation is at measurable levels (Fig. 2b). A diffuse correction (SWdiff; based on ) was added to incoming radiation after correction for the incidence angle to estimate the total incoming radiation (Fig. 2b):

$\begin{array}{}\text{(9)}& {\mathrm{SW}}_{\mathrm{diff}}=\mathrm{16}{\mathit{\psi }}^{\mathrm{0.5}}-\mathrm{0.4}\mathit{\psi }.\end{array}$

Lower peak radiation in the model (Fig. 2b) is due to a slight difference in slope and aspect between the radiometer (which was level) and the grid cell of the DSM containing the AWS (slope = 5; aspect = 177). Estimated radiation at each cell was then multiplied by (1−α) to estimate SWnet at each cell.

Albedo was estimated across the study area using the orthomosaics from 21 and 23 July. Each orthomosaic contains the digital number recorded by the camera in the red, green, and blue (RGB) channels. The digital number cannot be directly used to calculate surface reflectance without conversion equations proprietary to the camera manufacturer. However, several studies have used other approaches to derive surface reflectance from digital numbers . use the sum of RGB digital number values as a proxy for albedo, similar to the approach we employ here.

Values for RGB channels were averaged and the total range scaled to match the value at the AWS location (Fig. 3). Scaling in this way allows the values to be directly input into models, as opposed to providing a proxy for albedo, as was done in . The average measured albedo value for the study period (0.31; 21–24 July 2016) was used to fix the rescaled value for both 21 and 23 July at the AWS location. The lower end of the range was based on albedo values of cryoconite holes reported by Ryan et al. (0.1; 2017); an assumption was made that cryoconite holes would have similar albedo to debris on the glacier surface. The scaling of RGB values was used to create two gridded albedo products that could be directly input into the melt model (Fig. 3); on 21 July albedo was used as a model input up until 12:00 LT on 23 July, while the 23 July albedo was used from 12:00 LT on 23 July onward.

Mean albedo over the grid was 0.35 for 21 July and 0.33 for 23 July, and the median difference between the two was 0.020 (Fig. 3). Albedo measured at the AWS on 21 and 23 July at 12:00 LT (the time of image acquisition) was 0.315 and 0.309. The more pronounced difference in the gridded albedo likely reflects darker imagery in the 23 July mosaic in addition to a true lowering of the surface albedo over the time period. To test the importance of this lower albedo value, total absorbed radiation at the AWS for July 2016 was adjusted to reflect a lower albedo, beyond what was measured at the AWS. This lowering of 0.014 could produce approximately 0.028 m of additional melt over the course of the month.

### 2.3.2 Enhanced temperature index model

This study assesses the performance of an enhanced temperature index model because these models are commonly used for regional-scale estimations of glacier melt. The model we use is formulated after , which uses absorbed radiation (SWnet; MJ h−1) and temperature (C) in a linear regression model to estimate melt ($\stackrel{\mathrm{^}}{M}$; m). The formulation of the model is unique in that it controls for correlation between independent variables which can make model results unstable when not addressed. The uncorrelated temperature (Tresidual) is calculated by

$\begin{array}{}\text{(10)}& & {T}_{\mathrm{SW}}=\mathit{\eta }+\mathit{\beta }{\mathrm{SW}}_{\mathrm{net}},\text{(11)}& & {T}_{\mathrm{residual}}=T-{T}_{\mathrm{SW}},\end{array}$

where η and β are coefficients fit using a linear regression. Tresidual is then used in the following equation to determine melt in metres of surface lowering consistent with SR50 and UAV measurements:

$\begin{array}{}\text{(12)}& \stackrel{\mathrm{^}}{M}=\left\{\begin{array}{ll}\mathrm{TF}\cdot {T}_{\mathrm{residual}}+\mathrm{RF}\cdot {\mathrm{SW}}_{\mathrm{net}}& :T>{T}_{\mathrm{T}},\\ \mathrm{0}& :T\le {T}_{\mathrm{T}},\end{array}\right\\end{array}$

where SWnet can be calculated by the difference between incoming and outgoing radiation or by multiplying incoming radiation by (1−α), where α is the surface albedo. The coefficients TF and RF were fit using a linear regression. Threshold temperatures (TT) ranging from 1 to 2 C have been used in previous studies . Here we use TT=0C. Although we have not performed any testing of this threshold, it is unlikely that it would affect model outcomes in this study, as the air temperature rarely fell below 2 C during the month of July 2016 and never during the period where distributed melt information is available.

Cumulative positive Tresidual and SWnet were calculated for 12 h periods beginning on 14 July at 12:00 LT. Data from 14 to 20 July were used as a training period to determine the model coefficients (η, β, TF, and RF) using multiple linear regression. The fitted model coefficients are shown in Table 1. Model coefficients compare well to those reported by , a study which also modelled an Arctic glacier, although they report a wide range of values for RF depending on the training period.

Table 1Model coefficients obtained using linear regression modelling for 14–20 July 2016. Units for TF are meters per hour per degree Celsius (m h−1C−1), and units for RF are meters per hour per megajoule (m h−1 MJ−1).

A 1-D formulation of the model was run at the AWS site (ETIaws) using measured SWnet and near-surface air temperature. A distributed model (ETIdist) was developed by estimating model variables across the study area at each grid cell. Given the relatively small size of the study area, Tresidual was assumed to remain constant across the area, calculated by Eqs. (10) and (11) with modelled incoming radiation and measured temperature at the AWS site. Equations (10) and (11) can be used to back-calculate temperature when Tresidual and SWnet are known; this allowed modelled temperature to vary across the grid. SWnet was calculated from modelled distributed radiation and albedo, described above.

## 2.4 Model performance

The four models (EBaws, ETIaws, EBdist, and ETIdist) were run for the period 00:00 LT on 1 July to 24:00 LT on 21 July 2016, allowing for comparison between models over a longer time period than when melt measurements were available. Outputs from 00:00 LT on 14 July to 24:0 LT on 31 July were compared directly to melt measured with the SR50. To compare EBdist and ETIdist estimates, melt values were extracted from the grid cell containing the AWS. All models were assessed using the mean and standard deviation (σ) of the residuals ($\stackrel{\mathrm{^}}{M}-M$).

Total melt estimates for 21–24 July from both distributed models were compared to the gridded melt measurements across the study area for the same period. Model error was calculated at each grid cell in the study area; negative values indicate underestimation and positive values indicate overestimation of the model as compared to the UAV-derived data ($\stackrel{\mathrm{^}}{M}-M$). The median model error (ME) and normalised median absolute deviation (NMAD) were used to describe the model error, following . suggest using these metrics to describe errors when the error distributions are non-normal, which is typical for the large datasets produced with SfM . In the case of non-normal distributions, standard descriptive statistics (based on data means and standard deviations) are disproportionately influenced by outliers, which the ME and NMAD are not . The median and range of measured and modelled melt were also calculated for comparison.

In addition to model performance, the model error across the study area was used to investigate factors influencing model performance. We examined model errors in relation to surface characteristics which are known to influence the energy available for melt (Hock2005). The Pearson correlation coefficient is often used to assess variable relationships:

$\begin{array}{}\text{(13)}& R=\frac{\sum _{i=\mathrm{1}}^{n}\left({x}_{i}-\stackrel{\mathrm{‾}}{x}\right)\left({y}_{i}-\stackrel{\mathrm{‾}}{y}\right)}{{\mathit{\sigma }}_{x}{\mathit{\sigma }}_{y}},\end{array}$

where $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{y}$ are the variable means, and σx and σy are the variable standard deviations. As noted above, however, in the case that the mean and standard deviation are poor descriptors of a population, an alternative measure of relationship is necessary, using robust estimators for the location and variance . Based on , we use the median correlation coefficient to assess the relationship between model error and potential explanatory variables:

$\begin{array}{}\text{(14)}& {R}_{\mathrm{r}}=\frac{\sum _{i=\mathrm{1}}^{n}\frac{\mathrm{1}}{n}\left({x}_{i}-\mathrm{med}\left(x\right)\right)\left({y}_{i}-\mathrm{med}\left(y\right)\right)}{{\mathrm{NMAD}}_{x}{\mathrm{NMAD}}_{y}}.\end{array}$

In the case of large samples such as those we have here ($n>\mathrm{1.8}×{\mathrm{10}}^{\mathrm{7}}$), p values are not a suitable measure for correlation significance ; instead we evaluate the strength of the correlation based on its value. Robust correlations were computed between model error and terrain variables (slope and aspect), albedo, and an estimate of surface water flow.

Figure 4(a) Modelled 12 h melt values at the AWS site on 1–31 July 2016. (b) Cumulative measured and modelled melt at the AWS on 14–31 July 2016. SR50 measurements were only available for the second half of July, and a gap exists in the record due to melt out of the pole holding the SR50. Melt was extrapolated during this time for better visual comparison, and extrapolation was based on average melt rates for the day prior to and following the gap. The dates of the UAV study are shown in grey.

Water flow direction and potential upstream catchment were quantified using the Hydrology toolset in ArcGIS 10. Assuming that water is produced at every grid cell, flow paths were calculated from the 21 July DSM based on potential upstream accumulation. Cells with more than 1500 upstream cells were converted into a set of linear features, and the number of these features within 20 m of each grid cell was calculated (W20 m). The 1500-cell threshold represents 15 m2 of upstream drainage area and was chosen based on the assumption that 15 m2 is likely to provide enough meltwater to start flow. The resulting layer provided a representation of areas with significant surface water flow, including streams and thin sheet flows.

3 Results

## 3.1 Inter-model comparison at AWS

The four models were compared at the AWS site for 1–31 July 2016. All four models capture daily patterns of melt, with lower melt totals on average in the second half of the month, which is consistent with lower solar radiation recorded at the AWS (Figs. 4a and 2b). The distributed models both consistently estimate lower melt totals than their point model counterparts, with differences more pronounced when melt values are high. This is most likely related to the lower shortwave radiation in the distributed model, which is in total 17 % less than measured for the study period (Fig. 2b). The EB models produce lower melt estimates than the ETI models when solar radiation is low. During periods of high incoming solar radiation, EBaws and ETIaws produce similar melt.

The total melt estimated during the month of July was lower in the two distributed models, EBdist=1.223 m and ETIdist=1.394 m, compared to the point models, EBaws=1.436 m and ETIaws=1.508 m. The EB models estimate less melt for the month than the ETI counterparts.

## 3.2 Model performance at AWS site

Model estimates from 14 to 31 July 2016 were compared to melt measured at the AWS site with the SR50. Total melt at the SR50 during the study period was 0.186±0.03 m, while the imagery-derived measurement for the same period was 0.274±0.048 m. After tuning z0, the mean residual of 12 h model estimates at the AWS was 0.000 m for EBaws and −0.003 m for EBdist. The fitted ETIaws also had a mean residual of 0.000 m, while that of ETIdist was −0.001 m. Three of the four models underestimate melt at the site; ETIaws overestimates melt by 0.016 m (Table 2; Fig. 4b). Both distributed models produce lower melt estimates than their point model counterparts. The largest underestimate was from EBdist. The variability in modelled melt is notably lower than measured melt in all cases.

Table 2Descriptive statistics for measured melt, modelled melt, and model residuals at the AWS site from 14 to 31 July 2016. All values were calculated based on 12 h melt estimates and are reported in metres. We denote the standard deviation of all 12 h melt values with σ and the standard deviation of all residuals in 12 h model estimates with σRes.

## 3.3 Model performance over the study area

Across the study area, median measured melt was 0.184±0.048 m, and the NMAD was 0.026 m (Fig. 5b). Median melt across the area was lower in both EBdist and ETIdist, as was the variability (Fig. 5a and c). The total range of measured melt was 0.305 m, while the range of EBdist is 0.123 and ETIdist is 0.083.

Figure 5Distribution of measured and modelled melt across the study area between 21 July (12:00 LT) and 24 July (12:00 LT). (a) Melt estimated using EBdist applied across the study area. (b) Melt measured through differencing surfaces reconstructed from UAV imagery. (c) Melt estimated using ETIdist applied across the study area. (d, e) Model error, calculated by the difference between modelled and measured melt. Cells with an absolute value lower than the measurement uncertainty (0.048 m) are transparent to emphasise areas where the model residuals are high.

The ME and NMAD value of EBdist is $-\mathrm{0.064}±\mathrm{0.022}$, while that of ETIdist is $-\mathrm{0.050}±\mathrm{0.023}$. The ME of ETIdist is of a similar magnitude to the measurement uncertainty (0.048 m), indicating general agreement with measurements. When areas of error above the measurement uncertainty are highlighted (Fig. 5d and e), ETIdist shows clustered patterns representing 53 % of the total study area. Errors greater than twice the uncertainty (0.096 m) account for 5 % of the study area. EBdist errors greater than 0.048 m cover 79 % of the study area, while those greater than 0.096 m cover 11 %.

Figure 6Model error with magnitude greater than 0.048 m underlain by flow path features with potential upstream accumulation greater than 15 m2 for (a) EBdist and (b) ETIdist. Density of flow accumulation features within a 20 m radius of each cell (c).

Model uncertainty can be derived from a number of metrics, including variability in the model residuals, represented by the NMAD or standard deviation. Based on assessment at the AWS site, uncertainty in 12 h estimates from EBdist and ETIdist are 2σRes=0.018 and 0.022 m, respectively. use validation at an AWS and incorporate uncertainties related to modelling input variables to estimate an uncertainty of 15 %. In this study, we also have the benefit of examining distributed results across the study area through comparison with UAV-derived melt. The variability in the distributed models is characterised by the NMAD, taking $\mathrm{2}\cdot \mathrm{NMAD}=\mathrm{0.028}$ m for EBdist and $\mathrm{2}\cdot \mathrm{NMAD}=\mathrm{0.044}$ m for ETIdist; these uncertainties are 4 %–6 % of total estimated melt. The total uncertainty in the model errors can be derived using these uncertainties in combination with the uncertainty in measured surface change (0.048 m): $\sqrt{\left(\mathrm{2}\cdot \mathrm{NMAD}{\right)}^{\mathrm{2}}+{\mathrm{0.048}}^{\mathrm{2}}}$. This gives total uncertainties of EBdist=0.056 m and ETIdist=0.065. In either case, the errors which are strongly correlated with W20 m are higher than the combined uncertainty of the measurement and model, particularly for ETIdist.

## 3.4 Relationship to surface characteristics

Moderate correlation was found between model error and albedo, aspect, and slope for both EBdist and ETIdist (Table 3). Through field observation, we were able to note a visual relationship between surface water features and high error during analysis (Fig. 6a and c). This association was confirmed through correlation with W20 m, Rr=0.470 for EBdist and Rr=0.462 for ETIdist.

Table 3Robust correlation statistics (Rr) for surface characteristics and distributed models. Strong correlations are highlighted in bold.

Although most of the errors are negative, positive errors (i.e. overestimation) were also found. These are primarily associated with boulders and other large stationary objects on the glacier surface.

4 Discussion

Model results at the AWS indicate that both tuned point models perform very well. Performance of the distributed models is poorer, with EBdist performing worse than ETIdist, despite the physical basis of the energy balance model. This is likely due to a combination of factors, including limited data inputs available, and need for further refinement of energy exchange parameterisations. Lower modelled melt across the study area can be partly explained by lower modelled radiation, which was noted at the AWS site. However, the discrepancy between measured and modelled melt is greater across the grid than at the AWS site, indicating that the difference is not due to solar radiation alone. Previous work has shown that energy balance model performance is sensitive to variations in roughness length, wind speed, and albedo . Of these components, roughness length is the only EB model input used here that is not measured. An analysis of the model sensitivity to z0 revealed a 10 % decrease in total estimated melt (1–31 July) for z0=0.001 m, compared to the reference run with the tuned z0. Similarly, z0=0.01 m resulted in a 7 % increase in estimated melt. The deviation from the reference model run was caused by a 33 % decrease and 22 % increase in both QH and QE. Despite this sensitivity, z0=0.005 provided the best fit for the training data.

Variability in both EB and ETI models is linked to solar radiation cycles. In addition to diurnal cycles in SWin due to zenith angle variations, topography blocks direct solar radiation for a time during the night and causes stronger diurnal variation in SWin. The strength of the signal from the diurnal cycle is much greater in the melt produced by EB models than the ETI models (Fig. 4a). As noted above, on clear-sky days with relatively low temperatures, the EB and ETI models produce similar results (1–3 July). Under the opposite conditions, low SWin and high temperatures (e.g. 14–15 and 28–29 July), the EB and ETI models produce the most dissimilar results. This is due to the relationship between temperature and radiation in the ETI model, where at low temperatures and high SWin, Tresidual is very small and melt is driven primarily by radiation. At high temperatures and low SWin, Tresidual is relatively larger and drives melt. In contrast, melt in the EB model is always primarily driven by solar radiation. Examining 12 h averages of SWnet, LWnet, and Qh+QE in the EB model outputs, there is a strong correlation between SWnet and available melt energy (R=0.90). The correlation between melt energy and LWnet, as well as Qh+QE, is much lower (R=0.25 and $R=-\mathrm{0.19}$, respectively).

This contrast between the EB and ETI models and the local tuning of the ETI model allows ETIdist to better estimate melt recorded at the AWS. For comparison to other studies, σRes from ETIdist was converted into an hourly rate (ETIdist=0.0005 m h−1), which is similar to other ETI models in the literature . report a standard error of 0.00018–0.00053 m h−1, and the range in error stems from model runs with different training datasets. Hourly ETI model error reported by is as high as 0.0055 m w.e. h−1 compared to an energy balance model as reference data. report a standard error of 0.00031 m h−1, when the ETI model is transferred in time. Overall these results at the AWS indicate that, in the absence of additional validation data, it would be reasonable to employ this model to estimate melt across the glacier surface.

Hock (1999) found that when TF is calculated based on melt and temperature records, it is highly variable. The inclusion of shortwave radiation in ETI models is an attempt to account for this variability. The low variability seen in the results of both ETIaws and ETIdist when compared to AWS measurements is an indication of potential problems in the ETI models. Similar variability was noted in other studies . However, Hock (1999), , and all cite R2 values of 0.60, 0.86, and 0.80 when comparing total modelled melt to measured discharge. These previous results show that although the ETI models capture total melt, melt rates at individual points are captured less effectively. Similarly, we found across the study area that median modelled and measured melt agreed, but model results showed variability an order of magnitude lower than measurements (Fig. 5b and c). The availability of melt measurements at 10 cm horizontal resolution shed light on the areas where modelled and measured melt were (or were not) in agreement, which has not been shown before. These measurements also allowed us to examine possible explanations for the model errors, including the aspect and surface water flow.

It is worth noting that surface lowering measured through point cloud differencing is not always equivalent to surface melt. Vertical ice motion in the ablation zone offsets melt in measurements of surface change. However, in the case of Fountain Glacier, estimate surface uplift to be 0.2 m over the course of the ablation season, or 0.002 m d−1. Given the magnitude of measured surface elevation change over the 3 d period and the correspondence of those measurements to melt measured at ablation stakes, we assume that vertical uplift is a negligible contribution and that the point cloud differences represent surface melt. Additionally, each point cloud was independently georeferenced each day, removing complications from horizontal flow.

measured higher melt rates in active supraglacial streams than on surrounding ice. This higher rate of melt leads to streams downcutting into the glacier surface over time, as long as meltwater routing paths remain constant. also noted the difference in melt rate between an active stream and the nearby channel which the stream formerly occupied. During preliminary analysis of distributed model results, we observed the occurrence of high model errors in areas with visual evidence of surface water flow, including thin sheet flow and streams. Further investigation revealed a correlation between model errors and the density of linear flow features (designated through DSM analysis and representing streams as well as thin sheet flow). The northern portion of the study area has a predominantly northern aspect, which exhibits high density in W20 m but low model error. This pattern may result from lower melt on the northern aspect, leading to less water production and flow than what is suggested by inferred flow paths in W20 m. The relation between model underestimation and surface water features suggests that the water flowing on the glacier surface contributes a measurable amount to surface melt. This is not only consistent with observed patterns of downcutting streams and higher melt in active stream channels but also shows the importance of water in areas with thin sheet flow over the glacier surface.

note the role of kinetic energy in development of surface weathering crusts and the complexity of near-surface glacier hydrology, which has only recently come to light. Temperatures measured in a supraglacial stream on Fountain Glacier in 2014 ranged from 0 to 0.1 C . These temperatures alone are not sufficient to explain the magnitude of model errors in the vicinity of water features, suggesting that other forms of energy transfer play an important role. Idealised crevasse bottom incision equations presented by are driven by discharge rates, suggesting a role in flow-related incision for both heat advection as well as energy transfer from kinetic and potential energy. note that the relative importance of frictional heat contributing to incision rates in non-temperate glaciers remains unresolved. The results of this study suggest that further investigation into energy transfer (heat, kinetic, potential, and frictional) and albedo effects from flowing water on non-temperate glaciers is necessary.

The coincidence of high model error with water feature density was noted in areas with thin sheet flow in addition to the locations of supraglacial streams. The contribution of melt energy to the glacier surface from less concentrated water flow found here has not been examined in previous studies, and these results add new insight into the emerging picture of complex surface hydrology on non-temperate glaciers, as suggested by . The results of this study suggest that on Fountain Glacier the relative importance of energy flux from water flow is high enough to cause significant errors in the distributed model.

We believe that the relative importance of factors other than solar radiation which are driving melt on Fountain Glacier will have cumulative effects if ETIdist or EBdist were employed for long-term modelling. Although this study looks at a short time frame of only 3 d, weather conditions during the study are similar to those found throughout the ablation season on Fountain Glacier and are likely to be representative of average conditions throughout the season. The high melt rates measured in areas of surface water flow, which are absent in both distributed models, are enough to change the terrain characteristics of the glacier over time (as seen in the formation of deeply incised canyons on Fountain Glacier). The changes in terrain feed back into melt dynamics related to slope and aspect; the inability of both models to capture the highest melt rates will slow these feedback mechanisms compared to observed evolution of the glacier surface . show that migration of small surface streams on Fountain Glacier on annual timescales slows or prevents development of deeply incised streams. Similarly, the migration of thin sheet flow over the course of a season is likely to prevent deep heterogeneity across the glacier surface due to this type of water flow. However, the additional melt caused by this type of flow remains unaccounted for in modelling results.

The availability of high-resolution melt measurements now allows for analysis of model performance that is not possible with other methods of measuring melt (i.e. ablation stakes or satellite measurements). Building on previous studies of distributed energy balance models, high-resolution distributed melt data may be used to develop simplifications for glaciers similar to Fountain Glacier, which are more appropriate than primarily radiation driven models. Deeply incised canyons such as those found on Fountain Glacier have been noted on 83 other glaciers across the Arctic, including 10 on Bylot Island (Sarah St. Germain, personal communication, 2019). Other studies of glacier terrain characteristics have found that gradual shifts in glacier aspect due to differential melt have important impacts on total glacier mass balance . This reinforces the importance of developing a model which captures melt driven by both radiation and water flow.

5 Conclusions

We present the first study using high-resolution measurements derived from SfM and UAV imagery to assess performance of distributed energy balance and enhanced temperature index models; additionally we present a new method of deriving albedo from the UAV imagery that can be used directly in the grid-based model. The model developed in this study was compared directly to AWS measurements and to distributed measurements on a cell-by-cell basis. The availability of high-resolution data revealed patterns in model performance that could not be found with other traditional methods of measuring glacier melt (i.e. ablation stakes or satellite measurements).

The ETIdist exhibits similar performance at the AWS to other studies implementing ETI models. In the absence of high-resolution distributed data, this model might be applied with confidence across the glacier surface. The bulk performance of the model across the study area is also similar to distributed melt measurements, with a ME similar to the measurement uncertainty. However, the range in melt from both models is an order of magnitude lower than that seen in distributed melt measurements. The lower range in modelled melt, combined with spatial patterns found in model error, allowed for investigation of sources of model error.

Significant correlation between surface water flow and model error highlighted the important role of water flow in melt dynamics on Fountain Glacier. Energy transfer from water flow in supraglacial streams is known to cause stream incision and was previously noted by on Fountain Glacier. The role of thin sheet flow on the glacier surface has not been studied previously, although investigate near-surface hydrology and the role of kinetic energy in weather crust development. Our results corroborate the important role of kinetic energy in surface hydrology and suggest that on Fountain Glacier energy transfer from water flow is a significant driver of melt. Further investigation is necessary to determine the most important contributors of additional energy to the glacier surface from water flow, including the potential role of the lower water albedo when compared to ice.

Given the importance of water flow in the development of deep canyons on Fountain Glacier, a model which accounts for water-related energy inputs is necessary to effectively capture long-term surface evolution on the glacier. A model which better reflects the drivers of melt on Fountain Glacier, and potentially similar glaciers across the Arctic, will provide insight into changing melt dynamics in these environments. In future work we hope to extend the time series of melt measurements from UAV imagery and use the distributed dataset to derive a new melt model which captures the melt drivers highlighted in this study.

Data availability
Data availability.

Model outputs analyzed in this study are available at: https://doi.org/10.5683/SP2/S3V7G2 (Bash2019). Additional inquiries can be directed to the corresponding author.

Author contributions
Author contributions.

The study was conceptualised by EAB, with supervision from BJM. Field data collection was designed and performed by EAB, with supervision from BJM. EAB developed and implemented models, analysed all outputs, and prepared the paper. BJM edited the paper and provided feedback on analysis.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

Field support for this work was provided by the Polar Continental Shelf Project and Parks Canada. Eleanor A. Bash was supported by the University of Calgary's Eyes High Recruitment Scholarship. We would like to thank Allison Gunther for fieldwork and technical assistance and Mustafizur Rahman and Shawn Marshall for coding assistance. In addition, we would like to thank to thank the editor, Daniel Farinotti; three anonymous reviewers; and Shawn Marshall for valuable feedback on this paper.

Financial support
Financial support.

This research has been supported by the Natural Sciences and Engineering Research Council of Canada, Northern Scientific Training Program, and the Arctic Institute of North America.

Review statement
Review statement.

This paper was edited by Daniel Farinotti and reviewed by three anonymous referees.

References

Arnold, N. S., Willis, I. C., Sharp, M. J., Richard, K. S., and Lawson, W. J.: A distributed surface energy-balance model for a small valley glacier. I. Development and testing for Haut Glacier d'Arolla, Valais, Switzerland, J. Glaciol., 42, 77–89, 1996. a

Arnold, N. S., Rees, W. G., Hodson, A. J., and Kohler, J.: Topographic controls on the surface energy balance of a high Arctic valley glacier, J. Geophys. Res.-Earth, 111, F02011, https://doi.org/10.1029/2005JF000426, 2006. a, b

Avanzi, F., Bianchi, A., Cina, A., De Michele, C., Maschio, P., Pagliari, D., Passoni, D., Pinto, L., Piras, M., and Rossi, L.: Centimetric accuracy in snow depth using unmanned aerial system photogrammetry and a multistation, Remote Sensing, 10, 765, https://doi.org/10.3390/rs10050765, 2018. a

Bash, E., Moorman, B., and Gunther, A.: Detecting Short-Term Surface Melt on an Arctic Glacier Using UAV Surveys, Remote Sensing, 10, 1547, https://doi.org/10.3390/rs10101547, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q

Bash, E. A.: Supplementary data for “Surface melt and the importance of water flow – an analysis based on high-resolution unmanned aerial vehicle (UAV) data for an Arctic glacier”, Scholars Portal Dataverse, V1, https://doi.org/10.5683/SP2/S3V7G2, 2019. a

Bash, E. A. and Marshall, S. J.: Estimation of glacial melt contribution to the Bow River, Alberta, Canada, using a radiation-temperature melt model, Ann. Glaciol., 55, 138–152, 2014. a, b, c, d, e, f, g

Bugler, J.: The determination of hourly insolation on an inclined plane using a diffuse irradiance model based on hourly measured global horizontal insolation, Solar Energy, 19, 477–491, 1977. a

Cook, K. L.: An evaluation of the effectiveness of low-cost UAVs and structure from motion for geomorphic change detection, Geomorphology, 278, 195–208, 2017. a

Corripio, J. G.: Snow surface albedo estimation using terrestrial photography, Int. J. Remote Sens., 25, 5705–5729, 2004. a

Cuffey, K. M. and Paterson, W. S. B.: Physics of Glaciers, Academic Press, Burlington, MA, 2010. a

Ebrahimi, S. and Marshall, S. J.: Parameterization of incoming longwave radiation at glacier sites in the Canadian Rocky Mountains, J. Geophys. Res.-Atmos., 120, 12536–12556, 2015. a

Ebrahimi, S. and Marshall, S. J.: Surface energy balance sensitivity to meteorological variability on Haig Glacier, Canadian Rocky Mountains, The Cryosphere, 10, 2799–2819, https://doi.org/10.5194/tc-10-2799-2016, 2016. a, b

Fisher, D., Zheng, J., Burgess, D., Zdanowicz, C., Kinnard, C., Sharp, M. J., and Bourgeois, J.: Recent melt rates of Canadian arctic ice caps are the highest in four millennia, Global Planet. Change, 85, 3–7, 2012. a

Fitzpatrick, N., Radić, V., and Menounos, B.: Surface Energy Balance Closure and Turbulent Flux Parameterization on a Mid-Latitude Mountain Glacier, Purcell Mountains, Canada, Front. Earth Sci., 5, 67 https://doi.org/10.3389/feart.2017.00067, 2017. a

Fountain, A. G. and Walder, J. S.: Water flow through temperate glaciers, Rev. Geophys., 36, 299–328, 1998. a

Gabbi, J., Carenzo, M., Pellicciotti, F., Bauder, A., and Funk, M.: A comparison of empirical and physically based glacier surface melt models for long-term simulations of glacier response, J. Glaciol., 60, 1140–1154, 2014. a

Gardner, A., Moholdt, G., Arendt, A., and Wouters, B.: Accelerated contributions of Canada's Baffin and Bylot Island glaciers to sea level rise over the past half century, The Cryosphere, 6, 1103–1125, https://doi.org/10.5194/tc-6-1103-2012, 2012. a, b

Gardner, A. S., Moholdt, G., Wouters, B., Wolken, G. J., Burgess, D. O., Sharp, M. J., Cogley, J. G., Braun, C., and Labine, C.: Sharply increased mass loss from glaciers and ice caps in the Canadian Arctic Archipelago, Nature, 473, 357–360, 2011. a

Goswami, D. Y., Kreith, F., and Kreider, J. F.: Principles of solar engineering, CRC Press, Boca Raton, FL, 2000. a

Hock, R.: A distributed temperature-index ice- and snowmelt model including potential direct solar radiation, J. Glaciol., 45, 101–111, 1999. a, b, c, d, e

Hock, R.: Glacier melt: a review of processes and their modelling, Prog. in Phys. Geogr., 29, 362–391, 2005. a, b, c

Hock, R. and Holmgren, B.: A distributed surface energy-balance model for complex topography and its application to Storglaciären, Sweden, J. Glaciol., 51, 25–36, 2005. a, b

Höhle, J. and Höhle, M.: Accuracy assessment of digital elevation models by means of robust statistical methods, ISPRS J. Photogram. Remote Sens., 64, 398–406, 2009. a, b

Irvine-Fynn, T., Barrand, N., Porter, P., Hodson, A., and Murray, T.: Recent High-Arctic glacial sediment redistribution: A process perspective using airborne lidar, Geomorphology, 125, 27–39, 2011. a

Irvine-Fynn, T. D., Hanna, E., Barrand, N., Porter, P., Kohler, J., and Hodson, A.: Examination of a physically based, high-resolution, distributed Arctic temperature-index melt model, on Midtre Lovénbreen, Svalbard, Hydrol. Process. 28, 134–149, 2014. a, b, c, d, e, f, g, h

James, M. R. and Robson, S.: Sequential digital elevation models of active lava flows from ground-based stereo time-lapse imagery, ISPRS J. Photogram. Remote Sens., 97, 160–170, 2014. a

Lague, D., Brodu, N., and Leroux, J.: Accurate 3D comparison of complex topography with terrestrial laser scanner: application to the Rangitikei canyon (N–Z), ISPRS J. Photogram. Remote Sens., 82, 10–26, 2013. a, b

Lenaerts, J. T., van Angelen, J. H., van den Broeke, M. R., Gardner, A. S., Wouters, B., and van Meijgaard, E.: Irreversible mass loss of Canadian Arctic Archipelago glaciers, Geophys. Res. Lett., 40, 870–874, 2013. a

Lovitt, J., Rahman, M. M., Saraswati, S., McDermid, G. J., Strack, M., and Xu, B.: UAV Remote Sensing Can Reveal the Effects of Low-Impact Seismic Lines on Surface Morphology, Hydrology, and Methane (CH4) Release in a Boreal Treed Bog, J. Geophys. Res.-Biogeo., 123, 1117–1129, 2018. a

MacDougall, A. H. and Flowers, G. E.: Spatial and temporal transferability of a distributed energy-balance glacier melt model, J. Climate, 24, 1480–1498, 2011. a, b

Mair, D., Willis, I., Hubbard, B., Fischer, U., Nienow, P., and Hubbard, A.: Hydrological controls on patterns of surface, internal and basal velocities during three “spring events”: Haut Glacier d'Arolla, Switzerland, J. Glaciol., 49, 555–5671, 2003. a

Maronna, R. A., Martin, R. D., Yohai, V. J., and Salibián-Barrera, M.: Robust statistics: theory and methods, John Wiley & Sons, West Sussex, England, 2006. a, b

Matthews, T., Hodgkins, R., Wilby, R. L., Guðmundsson, S., Pálsson, F., Björnsson, H., and Carr, S.: Conditioning temperature-index model parameters on synoptic weather types for glacier melt simulations, Hydrol. Process., 29, 1027–1045, 2015. a, b, c

Montgomery, D. C. and Runger, G. C.: Applied statistics and probability for engineers, John Wiley & Sons, Hoboken, NJ, 2007. a

Moorman, B. J.: Glacier-permafrost hydrology interactions, Bylot Island, Canada, in: Proceedings of the 8th International Conference on Permafrost, 20–25 July 2003, Zurich, Switzerland, 783–788, 2003. a

Noël, B., van de Berg, W. J., Lhermitte, S., Wouters, B., Schaffer, N., and van den Broeke, M. R.: Six decades of glacial mass loss in the Canadian Arctic Archipelago, J. Geophys. Res.-Earth, 123, 1430–1449, 2018. a

Pellicciotti, F., Brock, B., Strasser, U., Burlando, P., Funk, M., and Corripio, J.: An enhanced temperature-index glacier melt model including the shortwave radiation balance: development and testing for Haut Glacier d'Arolla, Switzerland, J. Glaciol., 51, 573–587, 2005. a, b, c, d, e, f

Radić, V. and Hock, R.: Regionally differentiated contribution of mountain glaciers and ice caps to future sea-level rise, Nat. Geosci., 4, 91–94, 2011. a

Rippin, D. M., Pomfret, A., and King, N.: High resolution mapping of supra-glacial drainage pathways reveals link between micro-channel drainage density, surface roughness and surface reflectance, Earth Surf. Proc. Land., 40, 1279–1290, 2015. a, b, c

Ryan, J. C., Hubbard, A., Box, J. ., Brough, S., Cameron, K., Cook, J. M., Cooper, M., Doyle, S. H., Edwards, A., Holt, T., Irvine-Fynn, T., Jones, C., Pitcher, L. H., Rennermalm, A. K., Smith, L. C., Stibal, M., and Snooke, N.: Derivation of High Spatial Resolution Albedo from UAV Digital Imagery: Application over the Greenland Ice Sheet, Front. Earth Sci., 5, 40, https://doi.org/10.3389/feart.2017.00040, 2017. a, b, c

Shaw, T. E., Brock, B. W., Fyffe, C. L., Pellicciotti, F., Rutter, N., and Diotri, F.: Air temperature distribution and energy-balance modelling of a debris-covered glacier, J. Glaciol., 62, 185–198, 2016. a

Shevlyakov, G. and Smirnov, P.: Robust estimation of the correlation coefficient: An attempt of survey, Aust. J. Stat., 40, 147–156, 2011. a, b

Stevens, I. T., Irvine-Fynn, T. D., Porter, P. R., Cook, J. M., Edwards, A., Smart, M., Moorman, B. J., Hodson, A. J., and Mitchell, A. C.: Near-surface hydraulic conductivity of northern hemisphere glaciers, Hydrol. Process., 32, 850–865, 2018. a, b, c

St. Germain, S. and Moorman, B.: The development of a pulsating supraglacial stream, Ann. Glaciol., 57, 31–38, 2016. a, b

St. Germain, S. L. and Moorman, B. J.: Long-term observations of supraglacial streams on an Arctic glacier, J. Glaciol., 65, 900–911, 2019. a, b

Wake, L. and Marshall, S.: Assessment of current methods of positive degree-day calculation using in situ observations from glaciated regions, J. Glaciol., 61, 329–344, 2015. a

Watson, C. S., Quincey, D. J., Smith, M. W., Carrivick, J. L., Rowan, A. V., and James, M. R.: Quantifying ice cliff evolution with multi-temporal point clouds on the debris-covered Khumbu Glacier, Nepal, J. Glaciol., 63, 823–837, 2017. a, b, c

WGMS: Fluctuations of Glaciers Database, https://doi.org/10.5904/wgms-fog-2018-06, 2018.  a, b

Whitehead, K., Moorman, B. J., and Hugenholtz, C. H.: Brief Communication: Low-cost, on-demand aerial photogrammetry for glaciological measurement, The Cryosphere, 7, 1879–1884, https://doi.org/10.5194/tc-7-1879-2013, 2013. a

Whitehead, K., Moorman, B., and Wainstein, P.: Measuring daily surface elevation and velocity variations across a polythermal arctic glacier using ground-based photogrammetry, J. Glaciol., 60, 1208–1220, 2014. a, b, c