Using thermal UAV imagery to model distributed debris thicknesses and sub-debris melt rates on debris-covered glaciers

Abstract

Supraglacial debris cover regulates the melt rates of many glaciers in mountainous regions around the world, thereby modifying the availability and quality of downstream water resources. However, the influence of supraglacial debris is often poorly represented within glaciological models, due to the absence of a technique to provide high-precision, spatially continuous measurements of debris thickness. Here, we use high-resolution UAV-derived thermal imagery, in conjunction with local meteorological data, visible UAV imagery and vertically profiled debris temperature time series, to model the spatially distributed debris thickness across a portion of Llaca Glacier in the Cordillera Blanca of Peru. Based on our results, we simulate daily sub-debris melt rates over a 3-month period during 2019. We demonstrate that, by effectively calibrating the radiometric thermal imagery and accounting for temporal and spatial variations in meteorological variables during UAV surveys, thermal UAV data can be used to more precisely represent the highly heterogeneous patterns of debris thickness and sub-debris melt on debris-covered glaciers. Additionally, our results indicate a mean sub-debris melt rate nearly three times greater than the mean melt rate simulated from satellite-derived debris thicknesses, emphasising the importance of acquiring further high-precision debris thickness data for the purposes of investigating glacier-scale melt processes, calibrating regional melt models and improving the accuracy of runoff predictions.

1. Introduction

More than a quarter of the Earth's tropical glaciers are located in the Cordillera Blanca mountain range of Peru (RGI 6.0, 2017; Fig. 1). The meltwater from these Peruvian glaciers feeds into the Rio Santa river basin, providing water resources to ~250 000 people living in the Ancash region of Peru (Mark and others, 2010). The glaciers of the Cordillera Blanca have shown consistently negative mass balances and significant reductions in spatial coverage over the past decades. For example, Silverio and Jaquet (2017) reported that the total area covered by these glaciers shrunk by 46% between 1930 and 2017, while Veettil (2018) indicated an area loss of 33.5% between 1975 and 2016. Rabatel and others (2012) found an average mass balance of −0.76 m w.e. a⁻¹ between 1976 and 2010, while Seehaus and others (2019) reported a mass balance of −0.236 ± 0.042 m w.e. a⁻¹ between 2000 and 2016. Glacier retreat in the Cordillera Blanca is impacting agriculture and drinking water supplies in the region not only by modifying the quantity of water available, but also negatively impacting the water quality via acid rock drainage due to enhanced weathering of metal- and sulphide-rich bedrock (Fortner and others, 2011; Guittard and others, 2017; Mark and others, 2017). Many of the glaciers within this region are mantled with a layer of supraglacial debris (Seehaus and others, 2019), which is likely to be impacting both the retreat rate and the melt rate of these glaciers, thereby influencing downstream water toxification and long-term water resource depletion.

Fig. 1. Map of the study site location. (a) Location of the Ancash region within Peru (dark green shading). (b) Location of panel (c) within Ancash. The ice-covered areas within the Ancash region are shown by the white shaded areas in (b) and (c), while the rivers are shown by the dark blue lines. The coloured triangles in (c) show the locations of Llaca Glacier (red) and the Cuchillacocha weather station (installed by Bridgewater State University) where the meteorological data used within this study were collected (yellow).

A key factor controlling the melt rate of debris-mantled glaciers is the supraglacial debris thickness. Previous studies have shown that if the debris layer is thinner than a critical thickness, sub-debris melt rates are enhanced, while, if the debris layer is thicker than the critical thickness, sub-debris melt rates are reduced through insulation of the ice surface (e.g. Östrem, 1959; Nicholson and Benn, 2006; Vincent and others, 2016; Anderson and others, 2021). Furthermore, supraglacial ice cliffs and meltwater ponds, which are often abundant on debris-covered glaciers, create localised areas of enhanced ablation, further complicating the melt patterns on debris-covered glaciers (e.g. Sakai and others, 2000; Buri and others, 2016; Steiner and others, 2019). As the thickness of the debris layer and the presence of supraglacial features can change significantly over small spatial scales (e.g. Zhang and others, 2016; McCarthy and others, 2017; Nicholson and others, 2018), precisely mapping the debris thickness distribution is critical for effectively simulating their melt rates and meltwater contribution. Other debris characteristics, including moisture content, grain size and lithology can affect the relationship between debris thickness and sub-debris melt rate, for example, by altering the thermal conductivity of the supraglacial debris (e.g. Nicholson and Benn, 2012; Collier and others, 2014; Nakawo and Young, 1981). Therefore, it is also important to quantify the thermal conductivity of the debris layer in order to effectively simulate sub-debris glacial melt rates with a high degree of accuracy.

Debris thickness can be measured in situ by manual excavation through the debris layer to the debris-ice interface (e.g. Reid and others, 2012). However, such measurements are limited in scale, due to the challenges associated with accessing and navigating the surface of debris-covered glaciers. More recently, structure-from-motion (SfM) via terrestrial photogrammetry has been used to quantify the debris thickness exposed above ice cliffs (Nicholson and Mertes, 2017), while ground-penetrating radar has been used to quantify debris thicknesses over glacier surface transects (McCarthy and others, 2017). While these techniques have yielded greater spatial coverage compared to manual excavations, neither provide spatially complete, 3-D debris thickness observations. Additionally, since cliff-top debris thicknesses can differ considerably from surrounding debris thicknesses, the accuracy of debris thicknesses interpolated between observation sites is sometimes poor (McCarthy and others, 2017).

In order to provide spatially distributed estimates of debris thickness, several previous studies have used surface energy-balance modelling, combined with thermal satellite data, to derive the thermal resistance of the debris (e.g. Nakawo and Rana, 1999; Zhang and others, 2011). Since the thermal resistance is equal to the thermal conductivity divided by the debris thickness, in situ measurements of debris thermal conductivity can then been used, in conjunction with meteorological data, to model the thickness of the debris layer (e.g. Foster and others, 2012; Rounce and McKinney, 2014). Mihalcea and others (2008) used a different approach to model debris thickness from thermal satellite imagery, which involved finding the correlation between field-derived debris thickness and satellite-derived surface temperature and subsequently using this relationship to model glacier-wide debris thicknesses. Additionally, Herreid (2021) tested the use of an empirical approach using ground-based thermal imagery to model debris thickness, which further developed the approach of Mihalcea and others (2008) by using measured off-glacier surface temperatures as a reference in order to account for the effects of local radiative forcing on observed surface temperatures. While the methods of Mihalcea and others (2008) and Herreid (2021) offer reduced model complexity, neither account for the impact of spatial variations in meteorological conditions across the glacier surface on debris surface temperatures. For example, the quantity of incoming shortwave radiation that reaches the debris surface can vary considerably across relatively small spatial scales as a result of shading from both local glacier surface topography and surrounding mountain topography, leading to spatial inconsistencies in the relationship between surface temperature and debris thickness. Additionally, while the use of thermal satellite imagery enables debris thicknesses to be modelled across entire glaciers or regions, the resolution of satellite-derived thermal imagery is relatively coarse (>60 m). As a result, sub-pixel variations in debris thickness cannot be detected, while the presence of supraglacial ice ponds and ice cliffs can lead to underestimation of debris thickness values (Rounce and McKinney, 2014; Huang and others, 2017).

Imagery collected by uncrewed aerial vehicles (UAVs) offers significantly higher spatial resolution compared to satellite imagery. In recent years, a number of studies have used UAVs to collect high-resolution imagery in the visible spectrum (hereafter ‘visible imagery’) of debris-covered glaciers in the Himalaya, in order to investigate their surface characteristics (e.g. Immerzeel and others, 2014; Kraaijenbrink and others, 2016). In the Cordillera Blanca, visible UAV surveys of Llaca Glacier were conducted in 2014 and 2015 (Wigmore and Mark, 2017). Comparison of the data collected from these surveys showed spatially variable rates of ice loss, with the highest rates occurring where supraglacial ice cliffs and meltwater ponds were present (Wigmore and Mark, 2017). The use of UAVs to collect thermal imagery of a debris-covered glacier was demonstrated for the first time at Lirung Glacier in the Central Himalaya (Kraaijenbrink and others, 2018). The results showed high levels of spatial and temporal heterogeneity in the glacier's surface temperature, highlighting the potential drawbacks of using coarser-resolution thermal satellite data to model supraglacial debris thickness. The study also demonstrated that while UAVs facilitate the collection of high-resolution thermal imagery, it is important to account for factors such as surface emissivity variations and sensor bias to derive reliable absolute surface temperatures (Kraaijenbrink and others, 2018).

Here, we test the use of UAV-derived radiometrically calibrated thermal imagery, combined with local meteorological data, visible UAV imagery and thermal measurements taken within the debris layer, to produce centimetre-scale maps of distributed debris thickness for a portion of Llaca Glacier tongue (Fig. 1). Following calibration of the thermal imagery, we model the thermal conductivity of the debris layer using a time series of debris temperature measurements collected at varying depths within the debris layer. Spatially distributed debris thicknesses are estimated using a surface energy-balance modelling approach which accounts for the changes in meteorological conditions over the duration of the thermal UAV surveys. Using the high-resolution debris thickness maps produced in this study, daily spatially distributed melt rates are simulated over the duration of a 3-month period in 2019. The results are compared to melt rates simulated based on satellite-derived debris thicknesses in order to investigate the impact of debris thickness estimation accuracy on the simulated melt rates of debris-covered glaciers.

2. Methods

Figure 2 shows the workflow developed for this study, demonstrating schematically the links between the data acquisition, data processing and simulation steps of the methods, as described below.

Fig. 2. Workflow used in this study for simulating spatially distributed debris thickness and spatio-temporally distributed sub-debris melt rates.

2.1 Study site

Llaca Glacier is located in the central Cordillera Blanca, a 200 km-long mountain chain situated within the wider Peruvian Andes range (Fig. 1). Covering an area of ~5.1 km², the glacier extends from ~4460 to ~6090 m a.s.l. (RGI 6.0, 2017). The debris-covered tongue of the glacier has an area of ~0.22 km² and ranges in elevation from ~4460 to ~4620 m a.s.l. The supraglacial debris consists of a mixture of boulders, coarse gravel and fine dust. The meltwater from Llaca Glacier contributes to the supply of water for the Ancash region, which is inhabited by ~250 000 people. Llaca Glacier was selected as the site for this study due to its relative accessibility in comparison to other glaciers in the region, as well as due to the fact that Llaca Glacier has previously been surveyed in 2014 and 2015 (Wigmore and Mark, 2017).

2.2 UAV-based data collection

2.2.1 Thermal imagery acquisition

A standard DJI Phantom 4 UAV was fitted with a custom-built thermal camera system, comprising a FLIR Vue Pro R 640 (13 mm FOV) thermal camera and a U-BLOX GNSS GPS chip. This was used to collect ~15 000 radiometric thermal images across a total survey area of ~0.25 km² (Fig. 3a). The Phantom 4 was chosen because, unlike most commercially available drones, it is capable of flying at high altitudes of up to 6000 m a.s.l. Standard Phantom 4 propellers were used. The Vue Pro R camera was selected due to its ability to collect radiometrically calibrated thermal images at high thermal precision (30 mK/0.03°C). The built-in visible camera was removed from the Phantom 4 in order to reduce weight and allow greater flight times.

Fig. 3. UAV survey setup at Llaca Glacier. (a) Extents of the thermal and visible UAV surveys, and locations of the two UAV launch points and ground control points where GNSS data were collected and site where thermistors were installed within the debris layer. (b, c) Photographs of the materials used as ground control points for the thermal and visible UAV surveys respectively. (d) The custom-built UAV that was used to collect thermal imagery. (e) GNSS antenna setup for measuring the GPS position of each ground control point.

The UAV-based thermal imagery collection was conducted within two survey zones (Z _T1 and Z _T2) with differing launch point altitudes (LP₁: 4537 m a.s.l. and LP₂: 4576 m a.s.l.)(Fig. 3a), in order to ensure that the UAV maintained a safe altitude above the sloping glacier surface, since terrain correction was not used for the UAV flights. In total, four thermal UAV surveys (S _T1–S _T4) were conducted, each at different times of day on 18–19 August 2019 (Table 1). S _T1, S _T2 and S _T4 were launched from LP₁ and conducted within the bounds of Z _T1, while S _T3 was launched from LP₂ and conducted over the entire extent of Z _T2. S _T1, which covered an area of 94 000 m², was conducted between 16:25 and 17:20 h on 18 August 2019. S _T2 was conducted between 9:30 and 10:00 h on 19 August and covered an area of 87 000 m². The largest of the four surveys, S _T3, was conducted between 10:55 and 12:50 h on 19 August and covered an area of 137 000 m². The final survey, S _T4, was conducted between 14:25 and 15:45 h on 19 August and covered an area of 72 000 m².

Table 1. UAV survey information

The UAV was flown using an automated gridded flight plan, created using DroneDeploy flight planning software. As the option for terrain correction was not currently available with open-source flight planning software, the flight paths were along horizontal planes with a consistent altitude of 70 m relative to the launch point altitude. This flight altitude was chosen in order to provide a balance between obtaining high-resolution imagery (~5 cm) and providing coverage of a sufficiently large area (250 000 m² in total). The use of a consistent altitude relative to the launch altitude resulted in variations in the exact pixel spatial resolution and image overlap since the altitude above ground level (AGL) varied with surface topography. Since the surface elevation range of the complete survey area (~120 m) exceeded the average above-surface flight altitude (70 m), two separate launch points were used (Fig. 3a).

A flight speed of 7 m s⁻¹ and an image capture interval of 1 s were used, in order to provide a forward image overlap of 90%. The flight lines, which ran nearly perpendicular to the glacier flow direction, were spaced 7 m apart in order to provide an 80% lateral image overlap. During each of the four thermal surveys, the UAV was returned to its launch point multiple times for battery replacement. At an altitude close to sea level, the DJI Phantom 4 can fly for ~25 min between battery changes. However, due to the high altitude of Llaca Glacier (~4500 m a.s.l.), the air is considerably thinner and a significantly greater amount of power is required to create lift. Consequently, the average flight time between battery changes was roughly halved to ~12 min.

Many UAV-mountable thermal cameras, including the Vue Pro R 640 used in this study, use uncooled microbolometers, which are sensitive to changes in the temperature of the sensor, body and lens. While radiometric cameras apply corrections to account for these effects, Kelly and others (2019) highlight the need to allow time for the camera to stabilise after activation. For this reason, the camera was turned on ~20 min prior to the start of each survey, while a couple of extra flight lines were added to the start of each survey to allow the camera to adjust to meteorological conditions experienced during flight.

For calibration purposes, images of 40 × 40 cm anodised aluminium calibration targets (Fig. 3b) were collected with the Vue Pro R camera from an altitude of 10 m, at the beginning and end of every flight. Meanwhile, the temperatures of these panels were also recorded using an Apogee thermal infrared (TIR) radiometer for subsequent comparison against the UAV-acquired temperatures.

2.2.2 Visible imagery acquisition

Using a second DJI Phantom 4 UAV, with a built-in visible camera, 950 visible images of the glacier tongue were collected, covering a total survey area of ~325 000 m². Visible UAV data collection was also conducted using an automated gridded flight plan created using DroneDeploy flight planning software. An average flight altitude of 85 m was chosen, to allow the collection of high spatial resolution (2.5 cm) imagery, while providing coverage of a relatively large survey area. A flight speed of 5 m s⁻¹ and an image capture interval of 1 s were used in order to provide 90% forward overlap between images, while flight lines were spaced 45 m apart to allow 80% lateral image overlap.

Similar to the thermal UAV surveys, the UAV-based visible imagery collection was divided into two survey zones (Z _V1 and Z _V2), with corresponding launch points LP₁ and LP₂ respectively (Fig. 3a). In total, two visible UAV surveys (S _V1 and S _V2) were conducted (one for each of the two visible survey zones). S _V1 was conducted between 9:20 and 10:10 h on 21 August 2019 and S _V2 was conducted between10:50 and 12:20 h on the same day (Table 1). Due to the slightly lighter weight of the visible camera, in comparison to the thermal camera, a slightly longer flight time of ~15 min could be achieved between battery changes. At the beginning of S _V2, there was a technical camera error, which resulted in the camera changing from a nadir 0° angle to an oblique 90° angle, resulting in a small data gap within the visible imagery. This data gap did not impact the results as it was outside the area that debris thickness was modelled for.

2.3 Ground-based data collection

2.3.1 Ground control data acquisition for UAV surveys

In order to georeference the thermal and visible UAV imagery, two separate ground control surveys were conducted. The thermal ground control survey was carried out on 17 August 2019 (one day before the first thermal UAV survey) and the visible ground control survey was carried out on 20 August 2019 (one day before the visible UAV surveys). For each of the two ground control surveys, ground control point (GCP) targets were distributed across the UAV survey areas, with most of the points around the perimeters of the UAV survey areas, which were more accessible than the central survey areas (Fig. 3a). The GCP targets were fixed to flat surfaces using tape and rocks (Figs 3b, c).

For the thermal ground control survey, 20 GCP targets were assembled, each consisting of a 60 cm foam square with two triangles of insulated aluminium foil attached to the surface (Fig. 3b). Foam and aluminium were chosen due to their contrasting emissivity values of ~0.6 and ~0.1 respectively, making their central point clearly distinguishable from the UAV-mounted thermal camera. For the visible ground control survey, 22 GCP targets, each consisting of a 30 cm × 30 cm checkboard square (Fig. 3c), were set out across the glacier surface.

For each of the two ground control surveys, a Leica GNSS system was used to measure the position of each GCP target with high (sub-cm) spatial accuracy. A fixed-location GNSS reference station was set up in a flat area ~20 m in front of the glacier terminus, in order to collect continuous GPS measurements over the complete ~8 h duration of each ground control survey. Meanwhile, using a GNSS rover, the precise location of the centre of each GCP target was measured over a period of 5–10 min per target (Fig. 3e).

2.3.2 Vertical debris temperature profile measurements

In order to measure vertical changes in debris temperature within the debris layer, a vertical profile of five thermistors, each connected to a DataHog2 data logger, was installed within the debris layer near the western margin of the glacier (Fig. 3a). The thermistor probes were placed at depths of 5, 10, 20, 30 and 40 cm, with the 40 cm probe at the debris-ice interface. Once adjusted to local environmental conditions, the thermistors recorded temperatures at repeat intervals of 10 min between 17 August 2019 at 00:00 h and 19 August 2019 at 16:00 h. These measurements were used to model the thermal properties of the debris layer for integration within the surface energy-balance model (Fig. 2).

2.3.3 In situ surface temperature and debris thickness measurements

An Apogee TIR radiometer was used to collect a sequence of ground-based TIR measurements at 22 points across the glacier surface, with varying supraglacial debris thicknesses. At each measurement point, three emitted TIR measurements of the debris surface were taken. Subsequently, a pit was dug through the debris layer to the debris-ice interface and the depth of the debris layer was measured. These 22 measurements were taken in close succession over a total duration of 1 h 40 min (13:25–15:05 h) on 21 August, in order to minimise biases associated with temporal changes in meteorological conditions.

In order to validate the debris thickness model, an additional set of debris thickness measurements were taken, in conjunction with high-precision GPS positions measured with the Leica GNSS System. Unfortunately, due to a technical glitch with the pre-programmed UAV flights, several of these measurements were just outside the bounds of the thermal UAV survey. As a result, only three of the coupled GPS-debris thickness measurements could be used for validation of the debris thickness model.

2.4 UAV data processing

2.4.1 Producing surface temperature maps, DEMs and orthomosaics

To produce maps of surface temperature, the radiometric TIR images were processed using Pix4Dmapper software, which was selected due to its compatibility with the radiometric jpeg files collected by the Vue Pro R camera. To produce DEMs and orthomosaics, the visible images were processed using Agisoft Metashape Software. This software contains proprietary implementations of common SfM photogrammetric workflows, and includes feature recognition, image matching, bundle block adjustment, point cloud densification and ultimately the generation of high-resolution digital surface models and orthomosaics. To provide accurate georeferencing, the thermal and visible GCP targets were identified within the thermal and visible UAV images and linked to the known coordinates recorded during the thermal and visible ground control surveys. To account for the effects of emissivity on the amount of TIR energy emitted by the debris surface, an emissivity value of 0.94 was assumed when converting emitted TIR values measured by the Vue Pro R to surface temperatures (Salisbury and D'Aria, 1992).

2.4.2 Calibrating UAV-derived surface temperature maps

Images of a blackbody calibrator (a target object with an emissivity close to 1), captured in the lab using the same Vue Pro R camera that was used in the field, were used to calibrate the surface temperature maps to account for sensor bias (Fig. S4). These thermal images were captured for blackbody temperatures between 5 and 60°C, at 5°C intervals. The equation of the best-fit line between measured temperature and actual temperature was used to calibrate the surface temperature values collected by the thermal camera.

The surface temperatures derived from UAV-mounted TIR cameras can be influenced by atmospheric attenuation of thermal radiance (Maes and others, 2017). Since the UAV flights were conducted across horizontal planes with constant flight heights of 4607 m a.s.l. (S _T1, S _T2 and S _T4) and 4646 m a.s.l. (S _T3), the flight height AGL varied with surface topography. Consequently, the effect of atmospheric attenuation on measured surface temperatures is likely to have changed over the duration of the thermal UAV surveys. Since the flight height AGL is a function of elevation, a surface-altitude-dependent correction factor was applied to the thermal imagery in order to account for the effects of differential atmospheric attenuation on recorded surface temperatures. The surface-altitude-dependent correction factor was calculated based on differences between actual and recorded surface temperatures of exposed ice cliff surfaces, similar to the calibration approach used by Kraaijenbrink and others (2018). It was assumed that exposed areas of ice cliffs have a surface temperature of 0°C. This assumption was validated using spot measurements of surface temperature collected in the field with an Apogee TIR radiometer. Using a series of ice cliffs distributed from the lowermost to the uppermost part of each thermal UAV survey, the linear relation between the glacier surface elevation and the measured-actual ice cliff temperature difference was computed and subsequently used to correct the surface temperatures within the thermal orthomosaics (Fig. S4). Since the flight altitude AGL decreased continuously over the duration of each UAV survey (as the UAV gradually travelled up-glacier between sequential cross-sectional flight lines), it was assumed that this correction would also (at least partially) account for time-dependent sensor-related biases.

The accuracy of surface temperatures recorded by thermal cameras can be impacted by distortion caused by the lens optics, known as ‘vignetting’, where surface temperatures are slightly enhanced in the central region of each image and reduced in the outer portions of each image. It was assumed that, due to the continuously high overlap between subsequent images collected by the Vue Pro R camera, camera vignetting effects would be minimised by the averaging of temperature values during the image-stitching process and that any remaining vignetting effects, which may lower temperature values, were removed by the sensor bias correction.

2.5 Generating debris thickness maps

2.5.1 Estimating debris thermal properties

The effective conductivity of the debris was estimated at depths of 5, 10, 20, 30 and 40 cm within the debris layer, using the thermistor time series (Section 2.3.2). Following the methods of Conway and Rasmussen (2000), debris thermal diffusivity K (m² h⁻¹) was approximated as the gradient between the first derivative of debris temperature T (K) with respect to time t (hr) and the second derivative of debris temperature with respect to depth z (m):

(1)$$K = \displaystyle{{\dot{T}} \over {T^{\prime \prime}}} = \displaystyle{{\left({\displaystyle{{\partial T} \over {\partial t}}} \right)} \over {\left({\displaystyle{{\partial^2T} \over {\partial z^2}}} \right)}}.$$

This equation assumes that heat conduction is occurring solely vertically and that there are no significant heat sources or sinks within the debris layer. Using the approximated thermal diffusivity values, the effective thermal conductivity k _eff (W m⁻¹ K⁻¹) was computed at each depth within the debris layer, assuming a rock density ρ _rock of 2700 kg m⁻³, heat capacity c _rock of 750 J kg⁻¹ K⁻¹ (Clark, 1966) and porosity φ of 0.3 (Conway and Rasmussen, 2000):

(2)$$k_{{\rm eff}} = K\rho _{{\rm rock}}c_{{\rm rock}}( {1-\varphi} ).$$

The overall k _eff, which was used in the surface energy balance (SEB) model (discussed in Section 2.5.3), was calculated by treating the debris layer as a series of conductors corresponding to specific layers within the overall debris layer, each with different conductivities. These specific layers were: 0–5 cm depth (assigned the 5 cm modelled k _eff), 5–10 cm (assigned an average of the 5 and 10 cm k _eff values), 10–15 (assigned the 10 cm k _eff), 15–25 cm (assigned the 20 cm k _eff), 25–35 cm (assigned the 30 cm k _eff) and 35–40 cm (assigned the 40 cm k _eff). The overall k _eff was calculated as the arithmetic average of the k _eff values assigned to these layers, accounting for the relative depth of each layer. The arithmetic average was used in order to minimise skewing of the results due to a single layer with a very small or large k _eff.

In order to account for the non-linearity of the vertical temperature gradient, Rounce and McKinney (2014) introduced a non-linearity factor G _ratio, which is computed based on the difference in vertical temperature gradient between the top 10 cm of the debris layer and the total depth of the debris layer. However, as vertical temperature gradients are likely to vary between debris layers of differing depths, and since thermistor measurements were only available for one location on the glacier surface, we instead used an alternative approach. This approach involved introducing a tuning parameter X to account for the differences between the expected debris thickness (if the vertical temperature gradient were to be purely linear) d _e and the actual debris thickness measured in the field d _a:

(3)$$X = \displaystyle{{d_{\rm a}} \over {d_{\rm e}}} = \displaystyle{{d_{\rm a}} \over {\left({\displaystyle{{k_{{\rm eff}}( {T_{{\rm s}_{\rm g}}-T_{\rm i}} ) } \over {Q_{\rm c}}}} \right)}}, \;$$

where d _a values correspond to the 22 ground-based debris thickness measurements that were made in conjunction with ground-based radiometer-derived measurements of surface temperature $T_{{\rm s}_{\rm g}}$. d _e was modelled using surface energy-balance modelling (as in Section 2.5.3) at each of these 22 points based on $T_{{\rm s}_{\rm g}}$, alongside the meteorological data recorded at Cuchillacocha weather station at the time the measurements were made. Air temperature and incoming longwave radiation were corrected to account for altitude differences between the study site and the weather station, using a lapse rate of −6.61°C km⁻¹, derived from differences between values recorded at Cuchillacocha weather station (4630 m a.s.l.) and another weather station located further down the valley (3920 m a.s.l.). Since all of the in situ debris measurements of $T_{{\rm s}_{\rm g}}$ and d _a were made in non--shaded areas in the upper portion of the study area, with minimal variation in altitude (<20 m) between points, meteorological variables were assumed to be constant in space between the 22 measurement sites. X was calculated as the gradient of the linear relationship between d _e and d _a (Fig. S3). The y-intercept was not considered as it was calculated to be <0.01 (Fig. S3). The calculated X value of 2.21 was subsequently used within the surface energy-balance model to model debris thickness from UAV-derived surface temperatures (Section 2.5.3). This approach differs from the approach of Mihalcea and others (2008), which involved modelling debris thickness based on the empirical relationship between satellite-derived surface temperatures and manual in situ debris thickness measurements. Here, we have used a surface energy-balance model to estimate debris thickness and only used in situ surface temperature and debris thickness measurements to calibrate the model for the effects of non-steady temperature profiles through the debris layer.

2.5.2 Estimating the spatial and temporal distribution of meteorological variables

Since the thermal UAV surveys were conducted over periods of up to 2 h, it was necessary to account for the changing spatial distribution of incoming shortwave radiation (SW _in) over the duration of each survey. For this reason, the survey area was divided into cross-sectional segments corresponding to the areas surveyed during each 5 min period of the thermal UAV surveys and spatially distributed SW _in was modelled for each of these segments for the 5 min period within which each segment was surveyed. In order to do this, the 10 cm DEM of the glacier tongue (produced from the visible UAV imagery) was firstly joined with the ALOS 30 m DEM of the surrounding topography. The resulting joined DEM was used to model the spatial distribution of solar radiation (direct and diffuse) across the glacier surface, accounting for the effects of shading from both the glacier surface topography (e.g. supraglacial debris mounds and ice cliffs) and the surrounding mountain topography (similar to the approach used by Buri and others (2016) to model solar radiation distribution across supraglacial ice cliff surfaces). Solar radiation distribution was modelled for 5 min periods at 30 min intervals over the duration of each thermal survey. Through linear interpolation, a SW _in distribution map was produced for every 5 min period of each thermal UAV survey. In order to tune the modelled SW _in maps to the observations recorded at the weather station, each map was firstly divided by its maximum value to produce fractional SW _in maps for every 5 min period. The SW _in measurements recorded by Cuchillacocha weather station (Fig. 1) at 30 min intervals were subsequently linearly interpolated to each 5 min period within each survey and, based on the assumption that SW _in measured at the weather station was equal to the maximum radiation across Llaca Glacier tongue, each modelled fractional SW _in map for every 5 min period was multiplied by the corresponding weather-station-derived SW _in value for the same period. For each of the tuned SW _in maps produced, the cross-sectional area (corresponding to the 5 min flight time block for which SW _in was modelled) was extracted. Finally, all extracted cross-sectional segments were merged to produce a single map of SW _in, consisting of transverse swaths corresponding to each 5 min period of the thermal UAV survey.

Weather station observations of air temperature (T _air) and relative humidity (RH) were used to model the temporal variations in incoming longwave radiation, LW _in, over the duration of the thermal UAV surveys, using the Stephan–Boltzmann law:

(4)$$LW_{{\rm in}} = \varepsilon _{{\rm eff}}\sigma T_{{\rm eff}}^4 , \;$$

where $\varepsilon _{{\rm eff}}$ is the effective emissivity of the atmosphere, σ is the Stefan–Boltzmann constant (5.67 × 10⁻¹ W m⁻² K⁻⁴) and T _eff is the effective air temperature. T _eff is represented by T _air at screen level. As there were clear weather conditions with no clouds during the thermal UAV surveys, the clear sky emissivity ($\varepsilon _{{\rm clear}}$) was approximated using a parameterisation introduced by Dilley and O'Brien (1998) (Eqn 5). This approach has been found to provide the best parameterisation of LW _in over melting glaciers (Juszak and Pellicciotti, 2013).

(5)$$\varepsilon _{{\rm clear}} = \displaystyle{{a_{{\rm DO}} + b_{{\rm DO}}{\left({\displaystyle{{T_{{\rm air}}} \over {273.16}}} \right)}^ {\!\!6} + c_{{\rm DO}}\sqrt {\displaystyle{{4.65e_{\rm a}} \over {25T_{{\rm air}}}}} } \over {\sigma T_{{\rm air}}^4 }}, \;$$

where e _a is the atmospheric vapour pressure (Pa), which was approximated from the altitude-corrected T _air and relative humidity RH recorded at Cuchillacocha weather station (Fig. 1), using the Magnus formula (Bell, 1996). a _DO, b _DO and c _DO are parameters from Dilley and O'Brien (1998) which have fixed values of 59.38, 113.7 and 96.96, respectively. Using this method, modelled LW _in values were produced at 30 min intervals (corresponding to the frequency of weather station observations of T _air and RH), which were linearly interpolated to produce modelled LW _in values for every 5 min of each thermal UAV survey. Each of these LW _in values was assigned to each of the corresponding cross-sectional segments of the thermal survey areas associated with each 5 min flight time block of each survey. The cross-sectional segments produced were merged together to produce maps of LW _in which account for temporal variations in LW _in over the course of the thermal UAV surveys.

To account for altitudinal variations in T _air across the survey area, a spatially dependent altitudinal correction was applied, using the UAV-derived DEM of the glacier tongue and a mean lapse rate of −6.61°C km⁻¹, derived from differences between values recorded at Cuchillacocha weather station (4630 m a.s.l.) and another weather station located further down the valley (3920 m a.s.l.). To account for temporal variations in T _air and wind speed over the duration of the thermal UAV surveys, the weather station observations (30 min intervals) were again linearly interpolated to every 5 min flight time block of each thermal UAV survey. The resulting values were assigned to the cross-sectional segments corresponding to each 5 min flight time block of each thermal survey, to produce temporally corrected maps of T _air and wind speed for incorporation within the debris thickness model.

2.5.3 Producing modelled debris thickness maps

Using the calibrated UAV-derived surface temperatures, alongside the spatially and temporally distributed meteorological variables parameterised in Section 2.5.2, debris thickness maps were produced using surface energy-balance modelling, which has previously been used to model debris thickness from coarser-resolution thermal satellite imagery (e.g., Foster and others, 2012; Rounce and McKinney, 2014). As the surface temperature data collected between 11:55 and 12:50 h on 19 August (during survey S _T3) were most optimal for simulating debris thickness (as discussed further in Section 4.1), these data were used to model debris thickness across survey zone Z _T2 (Fig. 3a), which covers an area of ~137 000 m².

Firstly, the ground heat flux Q _c (W m⁻²) was calculated as:

(6)$$Q_{\rm c} = R_{\rm n} + H + LE, \;$$

where R _n is the net radiation flux (W m⁻²), H is the sensible heat flux (W m⁻²) and LE is the latent heat flux (W m⁻²). The net radiation flux was calculated as:

(7)$$R_{\rm n} = SW_{{\rm in}}( {1-\alpha } ) + \varepsilon ( {LW_{{\rm in}}-\sigma T_{\rm s}^4 } ) , \;$$

where SW _in is the incoming shortwave radiation (W m⁻²), α is the albedo (dimensionless), $\varepsilon$ is the emissivity (dimensionless), LW _in is the incoming longwave radiation (W m⁻²), σ is the Stephan–Boltzmann constant (5.67 × 10⁻⁸ W m⁻² K⁻⁴) and T _s is the surface temperature (K). Debris emissivity and albedo values of 0.94 and 0.3, respectively, were assumed (Salisbury and D'Aria, 1992; Nicholson and Benn, 2012).

The sensible heat flux was calculated as:

(8)$$H = \rho _{{\rm air}}\left({\displaystyle{P \over {P_0}}} \right)c_{{\rm air}}Au( {T_{{\rm air}}-T_{\rm s}} ) \;, \;$$

where ρ _air is the density of air (1.29 kg m⁻³), P is the atmospheric pressure (Pa), P ₀ is the atmospheric pressure at sea level (101 325 Pa), c _air is the specific heat capacity of air (1010 J kg⁻¹ K⁻¹), A is the transfer coefficient (dimensionless), u is the wind speed recorded at the weather station (m s⁻¹) and T _air is the air temperature. The atmospheric pressure was computed using the barometric pressure formula and the transfer coefficient was calculated as:

(9)$$A = \displaystyle{{k_{{\rm vk}}^2 } \over {{\rm ln}\left({\displaystyle{{z_{\rm h}} \over {z_0}}} \right){\rm ln}\left({\displaystyle{{z_{\rm h}} \over {z_0}}} \right)}}, \;$$

where k _vk is the von Kármán's constant (0.41), z _h is the height of meteorological measurements (2 m) and z ₀ is the surface roughness length, for which a value of 0.016 m was assumed (Rounce and McKinney, 2014).

The latent heat flux was assumed to be zero, based on the assumption that the debris was dry.

The debris thickness d was calculated as follows:

(10)$$d = X\displaystyle{{k_{{\rm eff}}( {T_{\rm s}-T_{\rm i}} ) } \over {Q_{\rm c}}}, \;$$

where T _i is the temperature at debris-ice interface (assumed to be 273.15 K based on the thermistor measurements).

The three in situ coupled GPS-debris-thickness measurements within the bounds of the thermal UAV survey were used as a guide to ensure that realistic debris thicknesses were being modelled. Negative modelled debris thickness values and values more than three median absolute deviations (MADs) from the mean were assigned as no data values.

To estimate the uncertainty associated with the mean debris thickness modelled across the study area, a sensitivity analysis was firstly conducted to determine the sensitivity of the model to input parameters (further details provided within the Supplementary Material). Based on these model sensitivities, in conjunction with estimated uncertainties associated with each input parameter, the overall error in mean modelled debris thickness was estimated through linear propagation:

(11)$$\sigma _{\rm d}^2 = \mathop \sum \limits_{i = 0}^n \displaystyle{{\delta d} \over {\delta y_i}}{\sigma_{y_{i}}} ^2, \;$$

where σ _d is the error in mean modelled debris thickness (m), δd/δy _i is the sensitivity of modelled mean debris thickness to changes in each input parameter y and σ _y is the estimated uncertainty associated with each input parameter, assumed independent from one another (further details provided within the Supplementary Material).

2.5.4 Omitting supraglacial ice cliffs and ponds from the model

Cliffs and ponds were not included in the model since the simulation approaches used to estimate debris thickness cannot be applied to areas beneath supraglacial ponds or on the surface of supraglacial ice cliffs. Cliffs were semi-automatically classified using the DEM and orthomosaic derived from the visible UAV survey. Firstly, areas with a surface slope of >40° were isolated. To eliminate false detection areas, which primarily occur along the edges of large boulders where the surface gradient is high, areas with a maximum inter-pixel difference (between central pixel and surrounding pixels) of more than 40, in the brightness of the greyscale orthomosaic, were removed. Interconnected areas of <1 m² were also removed, in order to eliminate any remaining small rocks from the areas classified as ice cliffs. Finally, any boulder edges which were not successfully eliminated during the previous steps were removed manually from the ice-cliff-classified areas. As supraglacial ponds were relatively rare in comparison to ice cliffs, these features were classified using manual digitisation. Once ice cliffs and ponds had been classified, the areas covered by these features were removed from the model.

2.6 Simulating daily sub-debris melt rates

Distributed daily sub-debris melt rates were simulated using the modelled debris thickness (derived from the surface energy-balance modelling approach), following the method of Nicholson and Benn (2006). This method uses mean daily meteorological data, assuming that the daily mean temperature profile through the debris layer is linear and that net changes in heat storage are negligible on diurnal timescales. This method was chosen as it significantly reduces the computational resources needed to run the model and has been demonstrated to yield reliable sub-debris ablation rate estimates despite model simplification (Nicholson and Benn, 2006). Firstly, the distributed average daily surface temperatures were solved for iteratively, based on Eqns 6–9.

Data from the Cuchillacocha weather station (Fig. 1) were used to calculate mean daily meteorological variables. Mean daily LW _in was calculated from T _air and RH (Eqns 4, 5), with a cloud correction (Unsworth and Monteith, 1975) being used in order to account for the effects of cloud cover on incoming longwave radiation:

(12)$$\varepsilon _{{\rm all}} = \varepsilon _{{\rm clear}}\cdot ( {1 + a\cdot n} ) + b \cdot n, \;$$

where $\varepsilon _{{\rm all}}$ is the all-sky emissivity, a and b are parameters from Unsworth and Monteith (1975) with fixed values of –0.84 and 0.84 respectively, and n is the cloud cover. Values of n recorded at the nearby Anta weather station (~15 km northwest of Llaca) since cloud cover values were not recorded at the Cuchillacocha weather station. The surface energy-balance model assumed that snow cover was absent and precipitation was negligible over the duration of the 90 d melt simulation period, since examination of all available visible satellite images collected during this period (43 images collected by Landsat 7–8 and Sentinel-2) showed that no snow cover was present and only 5.8 mm of precipitation was recorded at the nearest weather station over the duration of the melt simulation period.

The spatial variability of meteorological variables was accounted for using the following approaches. The spatial distribution of SW _in was modelled over a 24 h period every 10 d between the 5 July and 5 October. Through linear interpolation, maps of average SW _in distribution were produced for every day during this 3-month period. For each day, the corresponding map of SW _in distribution was divided by its mean value, before being multiplied by the mean SW _in recorded at the weather station on that day. To account for spatial variations in T _air, the UAV-derived DEM was used to produce a map of altitude-dependent T _air, using an environmental lapse rate of 6.5°C km⁻¹. u and LW _in were assumed to be spatially homogeneous.

Using the resulting surface temperature map, the ground heat flux was computed using Eqn 6 and the distributed daily melt rate M (m d⁻¹) was computed as:

(13)$$M = t\displaystyle{{Q_{\rm m}} \over {\rho _{{\rm ice}}L_{\rm f}}}, \;$$

Where t is time (seconds), Q _m is downward energy flux at the base of the debris layer (equal to Q _c), ρ _ice is the density of ice (assumed to be 900 kg m⁻³) and L _f is the latent heat of fusion (334 kJ kg⁻¹). Negative values of M were set to 0, while values more than three standard deviations from the median were set to no data values.

In order to investigate the impact of differing modelled debris thickness on simulated sub-debris melt rates, the model was run for a second time using debris thicknesses modelled from satellite data (Rounce and others, 2021) and the results were compared to the sub-debris melt rates modelled in this study.

In order to provide indirect validation for the sub-debris melt rates simulated in this study, the mean annual ablation rate was calculated independently, based on mass continuity, using UAV-derived glacier surface elevation changes and ice surface velocities. To calculate the mean surface elevation change across the study area, the UAV-derived DEM produced in this study (for the glacier surface in 2019) was differenced with another UAV-derived DEM of the glacier surface in 2014, produced by Wigmore and Mark (2017). Using the orthomosaics of the glacier surface produced both in this study (for 2019) and by Wigmore and Mark (2017) (for 2014), alongside modelled ice thicknesses produced by Farinotti and others (2019), the incoming and outgoing ice fluxes were computed for the survey area. Based on mass continuity, the mean ablation rate across the survey area between 2014 and 2019 was calculated.

To estimate the uncertainty associated with the mean sub-debris melt rate simulated across the study area for the 93 d model period, a second sensitivity analysis was conducted and overall uncertainty was estimated using the same approach used to estimate uncertainty in the mean modelled debris thickness (Eqn 11) (further details provided within the Supplementary Material).

3. Results

3.1 Vertical debris temperature profiles

The debris temperatures recorded by the thermistors at depths of 5, 10, 20, 30 and 40 cm (40 cm being the debris-ice interface) are shown in Figure 4. Temporal variations in debris temperature are greatest near the surface and reduce with increasing depth, with an average range in daily temperature of 17.9°C at 5 cm depth compared to 0.7°C at the debris-ice interface (40 cm depth). As depth in the debris layer increases, the times of the diurnal peaks and troughs in debris temperature are increasingly lagged. For example, there is a lag of 6.83 h between the average time of peak daily debris temperature at a depth of 5 cm (14:17 h) and at a depth of 40 cm (21:07 h) (Fig. 4a).

Fig. 4. Thermal properties of the supraglacial debris layer at Llaca Glacier, derived from thermistors located at depths of 5, 10, 20, 30 and 40 cm (with the debris-ice interface being at 40 cm depth) within the debris layer. Figure (a) shows the direct measurements that were recorded by each of the thermistors between 17 August 00:00 h and 19 August 16:00 h, with lines colour-coded according to the debris thicknesses shown in (b–f). Figures (b–f) show the relationship between the second derivative of debris temperature with respect to depth (d2T/dz2) and the first derivative of debris temperature with respect to time (dT/dt), derived from the time series recorded by each of the thermistors. The gradient of this relationship, which was used to approximate the thermal diffusivity, is shown for each of the thermistors in (b–f), along with the R ² value associated with each gradient.

3.2 Thermal conductivity

The effective thermal conductivity of the debris decreases by 0.963 W m⁻¹ K⁻¹ from near the surface of the debris layer (5 cm depth) to 30 cm depth (Table 2). However, the thermal conductivity appears to decrease below 30 cm, with a thermal conductivity 0.52 W m⁻¹ K⁻¹ lower at the ice-debris interface (40 cm depth) than at 30 cm depth. R ² associated with the gradients between ∂T/∂t and ∂²T/∂z ² are generally high (Figs 4c–f), except for at 5 cm depth, where the R ² value is considerably lower (Fig. 4b).

Table 2. Debris thermal properties on Llaca Glacier tongue

Thermal diffusivity and effective conductivity values, med from the thermistor measurements, are shown for each depth within the debris layer

3.3 Modelled debris thickness

Figure 5 shows the debris thickness modelled across Llaca Glacier using the workflow outlined in Figure 2. The results indicate a mean debris thickness of 0.18 ± 0.08 m, with a variance of 0.02 m, across the survey area. Debris thickness generally decreases up-glacier, with the lowest debris thicknesses being found in the uppermost portion of the survey area, where the debris layer was ~1–7 cm thick. Over >90% of the survey area (excluding supraglacial ice cliffs and ponds), the debris layer is thicker than 5 cm, indicating that a large proportion of the debris on Llaca tongue is likely to be thicker than the critical thickness required to reduce sub-debris melt rates through insulation of the ice surface (Östrem, 1959). Comparison against the three usable manually acquired debris thickness measurements, coupled with accurate GPS locations, within the survey area indicates good agreement between measured and modelled values, with differences of <5% between modelled and measured debris thickness values (Fig. 5b).

Fig. 5. Spatially distributed map of modelled debris thickness. Modelled debris thicknesses are shown in (a), with black triangles showing the locations of the in situ debris thickness measurements within the survey area. White areas show the presence of no data values, where modelled values were negative or more than three MADs outside the mean (discussed in 2.5.3). Grey areas show the presence of supraglacial ice cliffs and ponds, which were not included in the model (discussed in 2.5.4). RGB orthomosaic for the modelled area is shown in (b). A comparison between the modelled and measured debris thicknesses at these three sites is shown in (c), while (d) shows the spatial coverage of the debris thickness map and RGB orthomosaic shown in (a) and (b), respectively. This area corresponds to the thermal UAV survey area ZT2 (see Fig. 3).

3.4 Surface temperature

The results show that the spatial heterogeneity in surface temperatures is greatest during the middle of the day compared to early morning and late afternoon (Table S1). More specifically, the surface temperatures derived from the thermal imagery collected between 11:55 and 12:50 h (survey S _T3) show a variance of 35.8 K across the survey area, while the surface temperatures derived from the thermal surveys of 16:25–17:20 h (S _T1), 9:30–10:00 h (S _T2) and 14:25–15:45 h (S _T3) show lower variances of 14.2, 23.6 and 21.8°C, respectively, across their survey areas. Meanwhile, the surface temperatures from S _T3 have a std dev. of 6.0, while the surface temperatures from S _T1, S _T2 and S _T4 have lower std dev. of 3.8, 4.9 and 4.7°C, respectively. Note that some of the differences in variance may be partially attributable to different areas being covered by thermal surveys conducted at different times of day (as described in Section 2.2.1).

3.5 Simulated sub-debris melt rates

The results indicate a mean sub-debris melt rate of 0.70 ± 0.29 cm d⁻¹ across the total survey area, over the entire 3-month model period between 5 July and 5 October (Figs 6, 7), with high levels of spatial heterogeneity (mean melt rates range from 0.00 to 3.14 cm d⁻¹) across the study area (Figs 5d, 6). Sub-debris melt rates also generally increase up-glacier through the survey area (Fig. 6). Maximum sub-debris melt rates decrease as a function of debris thickness, with mean melt rates for the 93 d period of up to nearly 3 cm d⁻¹ for thin debris layers of a few mm, compared to mean melt rates of up to ~0.8 cm d⁻¹ where debris is 30 cm thick (Fig. 7c). Figure 7b shows that mean sub-daily melt rate for the period 5 July–October has a skewed distribution, with the most frequent mean melt rate being more than 0.5 cm d⁻¹.

Fig. 6. Spatially distributed simulated sub-debris melt rates on Llaca Glacier tongue. Maps of the mean simulated melt rates (across the area shown in Fig. 4c) are shown for three 31 d periods: 5 July–4 August (a), 5 August–4 September (b) and 5 September–5 October (c). These values were simulated backwards and forwards in time from the date of thermal UAV data collection, 19 July 2019. Black shaded areas show the presence of supraglacial ice cliffs and ponds (which were not included in the model) and white areas show the presence of no data values.

Fig. 7. Simulated sub-debris melt rates on Llaca Glacier tongue between 5 July and 5 October 2019. (a) Temporal variations (smoothed) in the average sub-debris daily melt rates (red line with std dev. shaded) and the 20 d moving average (yellow line) for the area shown in Figure 5c. (b) Probability density function (PDF) for the mean sub-debris melt rates of the period 5 July–5 October for the survey area shown in Figure 5c. (c) Mean simulated sub-debris melt rates for the same area and period shown as a function of modelled debris thickness. (d) Mean daily air temperature and incoming shortwave (SW) radiation recorded at Cuchillacocha weather station and the mean daily modelled incoming longwave (LW) radiation between 5 July and 5 October.

There are high levels of temporal variability in the sub-debris melt rates on Llaca Glacier tongue, with simulated mean daily melt rates for the whole survey area ranging from 0.00 to 1.78 cm d⁻¹ between 5 July and 5 October 2019 (Fig. 7a). The results indicate that sub-debris melt rates are generally increasing over the course of the 3-month simulating period (see moving average in Fig. 7a), with mean daily sub-debris melt rates of 0.53 cm between 5 July and 4 August, increasing to a mean daily melt rate of 0.82 cm between 5 September and 5 October (Fig. 6). Over the duration of the 3-month simulation period, the mean incoming SW and LW radiation fluxes are 195 and 207 W m⁻², respectively (Fig. 7d). The ranges in mean daily incoming SW and LW radiation are 194 (83–277) and 151 (170–321) W m⁻², respectively, over the 93 d period, while the range in mean daily air temperature is 3.8°C.

4. Discussion

4.1 Simulating debris thickness from thermal UAV imagery

The results of this study demonstrate that thermal UAV imagery can be used to effectively model spatially distributed supraglacial debris thicknesses. The results also demonstrated that, while high levels of precision can be gained from using such imagery, a number of calibrations and corrections are critical to ensure that (a) the thermal imagery is calibrated to account for biases associated with UAV-mounted thermal sensors, (b) the thermal imagery is corrected to account for the changing sensor-surface distance over the course of the thermal UAV flights (if terrain correction is not used), (c) the temporal changes in meteorological variables over the course of the thermal surveys are accounted for, and (d) the spatial variations in meteorological variables across the thermal survey area are accounted for. We recommend that future studies take the aforementioned steps in order to maximise the accuracy of debris thickness maps derived from thermal UAV imagery.

The results also show that thermal imagery acquired near the middle of the day is optimal for simulating debris thicknesses, due to (a) high spatial heterogeneity in surface temperatures, and (b) relatively low temporal variations in the modelled meteorological variables. As shown in Table S2, survey S _T3 yielded a variance in surface temperatures 50–150% greater, and a std dev. 20–60% greater, than the other three surveys. As a result, it is much easier to distinguish between different debris thicknesses using thermal imagery collected during S _T3, which was conducted during the middle of the day (10:55–12:50). In contrast, debris-thickness-driven differences in surface temperatures are less pronounced in the early morning, since debris has not yet heated up sufficiently, resulting in cooler and more homogeneous debris temperatures (further demonstrated in Fig. 4a). During the late afternoon, the debris has cooled significantly since the middle of the day (Fig. 4a), again partially obscuring some of the debris-thickness-driven differences in surface temperatures. Furthermore, during the early morning and late afternoon, temporal variations in meteorological variables are high (Figs S1a, b, d, e), making it difficult to account for variability in meteorological variables within the duration of the thermal UAV surveys, thereby impacting the accuracy of modelled debris thicknesses. We therefore recommend that future studies collect thermal imagery near the middle of the day (e.g. between ~11:00 and 13:00 h), when surface temperature variations are greatest and temporal meteorological variability is likely to be relatively low (Fig. S1).

The modelled debris thicknesses are likely to be the least reliable in areas where the debris layer is very thin, due to higher spatial variability in surface albedo (due to patchiness in the debris layer) and moisture within the debris (Nicholson and Benn, 2006; Fyffe and others, 2020). As a result, simulated melt rates are likely to decrease in reliability with decreasing debris thickness. The results of this study also suggest that the method of inverting debris thicknesses from surface temperature does not perform well where large boulders (>~1 m) are present. In these areas, debris thicknesses appear to be modelled as relatively thin debris areas, surrounded by borders of no data values, indicating that the surface heating of large boulders is likely to be independent from the thickness of the underlying debris.

The relationship between the second derivative of debris temperature with respect to depth (d²T/dz ²) and the first derivative of debris temperature with respect to time (dT/dt), which was used to estimate debris thermal conductivity, is generally strong (R ² > 0.75), with the exception of the uppermost portion of the debris layer (R ² value of 0.15 at 5 cm depth). A potential explanation for this lower R ² in the uppermost layer is that, since heat transfer is partly convective near the surface, convection can be considered as a heat source/sink outwith thermal conduction.

Based on the 22 coupled measurements of debris thickness and emitted TIR collected in the field, the tuning parameter × (which was used within the surface energy-balance model for estimating debris thickness) was calculated as 2.21. This value may vary considerably between glaciers and/or regions due to differences in factors such as lithology, which may result in differences in vertical temperature profiles within the debris layer. Therefore, we recommend that future studies simulating debris thickness from thermal UAV imagery collect site-specific in situ measurements of debris thickness and emitted TIR, to ensure that surface energy-balance models used are calibrated to account for these differences. Ideally, these measurements should cover a wide range of debris thicknesses, up to the thickness at which emitted TIR no longer changes with increasing debris thickness.

Figure 8 demonstrates the improved level of detail obtained using UAV-derived surface temperatures to model debris thickness, compared to the use of thermal imagery derived from satellites. Figures 8a and b show a comparison between the UAV-derived surface temperature data produced in this study and the highest-resolution thermal satellite imagery currently available (acquired by Landsat 7 at 60 m resolution and resampled to 30 m resolution). As demonstrated, there are large variations in debris thickness over relatively small spatial scales, which can only be distinguished using the thermal UAV imagery. This provides further evidence that glaciological models which use satellite-derived debris thicknesses as input data are likely to be affected by ‘inter-pixel mixing effects’, as described by Rounce and McKinney (2014). We compare our UAV-derived results (Fig. 8c) to a global debris thickness data modelled from satellite data (Fig. 8d), which was previously produced by Rounce and others (2021) using a surface temperature inversion method (based on Landsat 8 surface temperature data) in conjunction with a mass-continuity-based sub-debris melt inversion method. Our UAV results indicate a mean debris thickness of 0.18 ± 0.08 m across the survey area, which is 74% less than the mean debris thickness modelled by Rounce and others (2021) for the same area (0.70 m). This suggests that debris thicknesses on Llaca Glacier may be significantly overestimated by satellite-based models, with significant impacts on simulated sub-debris melt rates (as discussed further in Section 4.2). Further evidence is therefore required in order to confirm whether this could be indicative of a wider issue.

Fig. 8. Comparison between UAV and satellite-derived surface temperature information and modelled debris thicknesses. (a) Calibrated surface temperatures derived from thermal UAV imagery, acquired between 10:55 and 12:50 h on 19 August 2019. (b) Landsat Collection 2 Surface Temperature Product, generated from Landsat 7 imagery acquired at 15:03 h on 19 August 2019. (c) Debris thicknesses modelled in this study from thermal UAV imagery. Black shaded areas show the presence of ice cliffs and white areas show the presence of no data values. (d) Debris thicknesses modelled by Rounce and others (2021). The black line surrounding each of the four maps shows the same reference area on Llaca Glacier, as shown in Figure 5c.

4.2 Simulating sub-debris melt rates from UAV-derived debris thickness maps

This study demonstrates that by using high-resolution maps of debris thickness derived from thermal UAV data, greater-precision estimates of sub-debris melt rates can be produced. As shown in Figure 6, high levels of spatial heterogeneity in ice surface melt rates over sub-metre scales can be detected across the survey area on Llaca Glacier tongue. Inputting the derived debris thicknesses modelled by Rounce and others (2021) (which are nearly four times greater than our estimates) into our sub-debris melt model, we estimated a mean simulated melt rate of 0.25 cm d⁻¹, which is only about one-third of the melt rate simulated in our study (0.70 ± 0.29 cm d⁻¹). This indicates that since the debris thicknesses modelled from satellite data appear to be considerably overestimated for Llaca Glacier tongue, this would lead to significant underestimation of sub-debris melt rates as a consequence (by a factor of nearly 3 for this study). We therefore advocate for further high-precision studies of debris thickness, with better model validation, in order to better understand the uncertainties associated with satellite-derived debris thicknesses and more accurately calibrate regional-scale models of debris thickness and sub-debris melt rates.

The general increase in simulated sub-debris melt rate over the course of the 93 d simulation period (Fig. 7a) is likely to result from enhanced incoming LW radiation (Fig. 7d) resulting from an increase in cloud cover towards the end of the simulation period.

The mean total simulated sub-debris melt rate for the entire 93 d period between 5 July and 5 October 2019 was 0.65 m. For indirect validation of the simulated sub-debris melt rate, we calculated the mean annual ablation rate using UAV-derived surface elevation change and ice surface velocities, based on the mass continuity method, as described further in Section 2.6. Based on this approach, we calculated a mean ablation rate of 3.63 m a⁻¹ between 2014 and 2019. If minimal seasonal variations in ablation were to be assumed, this would equate to a mean total melt of 0.92 m across the survey area during the 93 d model period (5 July–5 October). A recent study found that at the neighbouring Cuchillacocha Glacier (which has similar debris cover characteristics to Llaca Glacier), although ablation occurs throughout the year, melt rates are generally greater during the wet season (October–March) than the dry season (April–September) (Fyffe and others, 2021), indicating that melt rates during the 93 d model period were likely to be lower than the annual average melt rate. Therefore, the mean total ablation is likely to have been <0.92 m during the model period, indicating that the melt rate simulated in this study (0.65 m over the 93 d model period) is a reasonable estimate.

A recent study indicated that seasonal melt patterns vary between two glaciers studied near Llaca Glacier (Fyffe and others, 2021), with one glacier (Shallap) having the greatest melt rates during the dry season (April–September) while another glacier (Cuchillacocha) has the greatest melt rates during the wet season (October–March). If wet season melting on Llaca Glacier tongue is similar to or greater than dry season melting, a lower melt rate than 0.65 m (simulated in this study) might be expected during the 93 d simulation period. However, given that ice flow velocity decreases down-glacier towards the terminus (Wigmore and Mark, 2017), it is likely that the consequent flux imbalance is partially offsetting melt-driven ice-surface lowering, indicating that 0.65 m of melt is a reasonable estimate for the simulation period.

To investigate further the effectiveness of the models used in this study, we compared the relationship between modelled debris thickness and modelled mean surface temperatures (produced within the sub-debris melt model) on 19 August (Fig. S2) to the relationship between in situ measurements of debris thickness and surface temperature collected on the same day (Fig. S5). In line with in situ measurements of surface temperature (Fig. S1), modelled surface temperatures show a logarithmic increase with increasing debris thickness, with surface temperatures beginning to stabilise when debris thickness exceeds ~30 cm (Figs S1, S2). As the modelled surface temperatures represent daily averages, while the in situ measurements were gathered at a specific time of day, the magnitude of modelled and observed surface temperatures cannot be compared. Therefore, in future, measurements of diurnal temperature variations in surface temperature at a range of different debris thicknesses would be valuable to better evaluate models similar to those used in this study.

4.3 Thermal properties of supraglacial debris

As shown by numerous previous studies (e.g. Conway and Rasmussen, 2000; Nicholson and Benn, 2012), our results further confirm that, during the daytime, the temperature of the supraglacial debris on Llaca Glacier is lowest at the debris ice interface and increases towards the surface of the debris layer (Fig. 4a). The increasingly lagged peaks in debris temperature from the surface to the base of the debris layer indicate that, with increasing depth within the debris layer, the delay time in the melt response to meteorological forcing also increases.

The depth-averaged thermal conductivity of the debris layer at Llaca Glacier (0.78 W m⁻² K⁻¹) appears to be slightly lower than the values generally reported by previous studies in the Himalaya. For example, at Ngozumpa Glacier in Nepal, Nicholson and Benn (2012) found thermal conductivities of 0.95 ± 0.10 and 1.29 ± 0.13 W m⁻² K⁻¹ for dry debris in winter and summer respectively. This lower thermal conductivity could potentially be attributed to greater debris porosity, which results in more air being trapped within the debris layer and less heat being transferred downwards through the debris layer (Juen and others, 2012).

As a result of its relatively low thermal conductivity, the debris layer on Llaca Glacier is likely to be providing a greater insulative effect on the ice below than on glaciers where the debris thermal conductivity is greater. This suggests that the debris is having a relatively high inhibiting effect on melt rates at Llaca Glacier, in comparison to some of the well-studied glaciers in the Himalaya. This could be indicative of a higher regional importance of debris cover in controlling glacier melt rates in the Cordillera Blanca. However, further evidence is required in order to better understand the role of debris cover in controlling glacial melt rates and downstream hydrology in the Ancash region of Peru.

4.4 Model sensitivity and limitations

The sensitivity analyses results (Tables S3, S4) indicate that the debris thickness model is most sensitive to incoming shortwave radiation and albedo, while the error analyses results (Tables S3, S4) show that uncertainties in albedo make the largest contribution (58%) towards overall uncertainty in the mean modelled debris thickness. In order to reduce the contribution of albedo uncertainties in the future, ground-based pyranometer measurements of albedo could be collected and surface classification of albedo could be divided into a greater number of categories (e.g. corresponding to different debris lithologies). Additionally, while incoming solar radiation recorded in the nearby valley is likely to be similar to Llaca Glacier (as there were cloud-free conditions during the thermal surveys), future studies could further minimise model uncertainties associated with incoming solar radiation and other meteorological variables by gathering meteorological data on site. The sensitivity analysis also shows that the debris thickness model is relatively sensitive to surface temperature (Table S3), with surface temperature uncertainties contributing towards 14% of the overall uncertainties in modelled mean debris thickness, emphasising the importance of accurately calibrating the surface temperatures derived from UAV-mounted thermal cameras. While various calibration procedures were performed to correct for the effects of sensor biases, sensor drift, atmospheric signal attenuation and surface emissivity variations, future efforts could be made to further improve these calibration procedures. For example, since the anodised aluminium panels used for quality control appeared to have relatively unstable surface temperatures, which varied considerably over short timescales and across their length, the use of calibration targets consisting of other materials and/or a portable calibration blackbody may help to improve calibration accuracy in the future.

Additionally, the sensitivity analyses indicate that the sub-debris melt model is most sensitive to albedo, incoming shortwave radiation and incoming longwave radiation (Table S4). Similar to the debris thickness model, uncertainties in the albedo make the largest contribution towards the overall uncertainty in mean simulated sub-debris melt rate, further emphasising the potential benefit of future studies collecting ground-based albedo measurements and increasing the number of surface classification categories used for assigning different surface albedo values. Additionally, the high sensitivity of the model to incoming shortwave radiation emphasises the importance of effectively accounting for the effects of shading from the surrounding topography when simulating the spatial distribution of solar radiation, highlighting the benefits of using a high-resolution UAV-derived DEM of the glacier surface to effectively account for the effects of local surface topography (e.g. debris mounds and ice cliffs) on shortwave radiation distribution.

The R ² value associated with the gradient between d²T/dz ² and dT/dt at the shallowest thermistor at 5 cm depth (0.15) is considerably lower than the R ² value associated with this gradient at the deeper thermistors (0.75–0.90). This is likely to result from the fact that temperatures recorded by shallow sensors within the debris layer are typically more noisy than temperatures recorded by deeper sensors (e.g. Collier and others, 2014). Consequently, thermal conductivity values associated with shallower depths within the debris layer are likely to be less reliable than those at greater depths. It is also possible that horizontal heat conduction could be contributing towards the temperatures observed by the thermistors within the debris layer, impacting the accuracy of the thermal conductivities calculated from the thermistor data. Furthermore, the thermal conductivity of the debris layer may vary spatially across the glacier surface (Laha and others, 2022). Nevertheless, the sensitivity analysis confirms that the debris thickness and sub-debris melt models are not highly sensitive to thermal conductivity (Tables S3, S4).

A key limitation associated with this study is the scarcity of debris thickness data suitable for validation (as discussed in Section 2.3.3), which results in some remaining uncertainty in the reliability of debris thicknesses modelled from thermal UAV data. In future, further thermal UAV surveys of debris-covered glaciers, coupled with a greater number of in situ debris thickness measurements within surveyed areas, would be highly beneficial for establishing this technique as a viable method for obtaining high-precision estimates of supraglacial debris thickness with improved levels of accuracy. While we do not have direct measurements of sub-debris melt rates, ice-surface-lowering rates were compared against simulated melt rates for indirect validation of the sub-debris melt model. In future, the collection of ablation-stake measurements at debris-thickness simulation sites would enable further validation of simulated sub-debris melt rates produced using UAV-derived modelled debris thicknesses.

5. Conclusions

This study has presented an approach for simulating high-resolution, spatially distributed supraglacial debris thicknesses and sub-debris melt rates from UAV-derived thermal imagery, in conjunction with local meteorological data, visible UAV imagery and vertically -profiled debris temperature measurements. We have demonstrated that by (a) effectively calibrating the radiometric thermal imagery, (b) accounting for the temporal variations in meteorological variables over the UAV survey duration, (c) estimating the spatial distribution of meteorological variables across the survey area, and (d) simulating the thermal conductivity of the debris layer, surface energy-balance modelling can be used to model effectively the debris thickness and sub-debris melt rates of debris-covered glaciers. We have also demonstrated that by obtaining high-resolution (10 cm) UAV imagery, as opposed to using coarser (>60 m) satellite imagery, the highly spatially heterogeneous debris thickness across Llaca Glacier tongue can be more precisely represented, facilitating higher-precision sub-debris melt simulation. Our findings have indicated that the mean debris thickness across the survey area on Llaca Glacier tongue is ~74% lower than the satellite-derived estimate, indicating that the accuracy of satellite-derived debris thicknesses are likely to be poor in places. Furthermore, sub-debris melt simulation has indicated that this overestimation of debris thicknesses would have resulted in a considerable underestimation of sub-debris melt rates across the survey area, with simulated melt rates being only about one-third of the melt rates simulated in this study from high-resolution UAV data. The reliability of debris thicknesses and sub-debris melt rates modelled in this study may decline with decreasing debris thickness, due to higher variability in surface albedo and greater moisture content. Our results also indicated that the debris layer on Llaca Glacier has an ~18–40% lower thermal conductivity compared to the debris on previously studied glaciers in the Himalaya, suggesting that the inhibiting effect of debris on melt rates may vary considerably between glaciers and/or regions. Overall, the results of this study emphasise the need for further high-precision UAV/ground-based studies of the thermal properties of supraglacial debris in the Cordillera Blanca, as well as in other mountain regions around the world, in order to better calibrate debris thicknesses within glaciological models and improve the accuracy of runoff predictions.