Thinning and surface mass balance patterns of two neighbouring debris-covered glaciers in the southeastern Tibetan Plateau

Debris-covered glaciers are a common feature of the mountain cryosphere in the southeastern Tibetan Plateau. A better understanding of these glaciers is necessary to reduce the uncertainties in regional water resource variability and to anticipate potential cryospheric risks. In this study, we quantified the seasonal thinning and surface mass balance patterns of two neighbouring debris-covered glaciers (23K Glacier and 24K Glacier) in the southeastern Tibetan Plateau with four unpiloted aerial vehicle surveys and in situ measurements. We observed that the thinning of 23K Glacier was ∼2–7 times greater than that of the 24K Glacier for annual and cold periods. The surface velocity of the 24K Glacier is higher than that of the 23K Glacier (∼5–6 times) for all periods. In contrast with the thinning patterns, the surface mass balance patterns of the two glaciers closely agree across the different periods. We found that the surface mass balance distribution strongly correlated with the spatial distribution of debris thickness for both glaciers. Ice cliffs and supraglacial ponds are prevalent on the surface of these glaciers (∼4.4–7.2±0.5 %) and enhance melt overall (enhancement factor: ∼2.5) but do not control the surface mass balance pattern of either glacier. This comparison study of two neighbouring glaciers confirms the significance of both glacier dynamics and debris thickness in controlling thinning and melt for the different debris-covered glaciers of the southeastern Tibetan Plateau in the context of climate change.

1 Introduction

Monsoon-influenced glaciers in the southeastern Tibetan Plateau have experienced more significant mass loss than those in most other regions of High Mountain Asia in the past 2 decades (Kääb et al., 2012; Yao et al., 2012; Brun et al., 2017; Shean et al., 2020; Hugonnet et al., 2021a). The percentage of debris-covered area as a proportion of the total glacierized area in the southeastern Tibetan Plateau was estimated at ∼17 %–19 %, which exceeds the percentage of debris-covered glacier area at the global scale (∼4.4 %–7.3 %; Scherler et al., 2018; Herreid and Pellicciotti, 2020). A better understanding of the evolution and mass balance patterns of debris-covered glaciers in the southeastern Tibetan Plateau is essential for constraining changes in regional water resources (Zhang et al., 2011; Neckel et al., 2017). Because proglacial lakes can develop on or in front of debris-covered glaciers (Wang et al., 2011; Allen et al., 2019; Racoviteanu et al., 2022), and because glacier thinning may affect slope stability (Kääb et al., 2021; An et al., 2022; Zhao et al., 2022), expanding this knowledge base may also inform the understanding of mountain geohazards.

The presence of debris can influence the response of glaciers to climate change. If the supraglacial debris is more than a few centimetres in thickness, it will provide a melt-buffering effect (Østrem, 1959; Nakawo et al., 1999; Nicholson and Benn, 2006; Reid and Brock, 2010; Anderson and Anderson, 2016; Yang et al., 2017). However, several satellite remote sensing studies have found similar thinning rates for debris-free and debris-covered glaciers (Kääb et al., 2012; Gardelle et al., 2013; Pellicciotti et al., 2015; Brun et al., 2019). Additionally, a higher thinning rate of debris-covered glaciers than debris-free glaciers was found in the Lahaul and Spiti region of the Indian Himalaya (Vincent et al., 2013), and several studies have also confirmed the strong surface thinning of debris-covered glaciers in the southeastern Tibetan Plateau (−0.52 to −0.83 m a⁻¹; Neckel et al., 2017; Ke et al., 2020). This phenomenon has been referred to as the “debris-cover anomaly” (Pellicciotti et al., 2015; Vincent et al., 2016). Ice cliffs and supraglacial ponds could partly explain this anomaly because they are directly exposed to incoming radiation and therefore act as melt “hotspots” (Sakai et al., 1998, 2002; Reid and Brock, 2014; Juen et al., 2014; Steiner et al., 2015; Buri et al., 2016; Miles et al., 2016, 2018; Buri et al., 2021). The area influenced by cliffs and ponds is characterized by high melt rates relative to the surrounding debris-covered area (not containing ice cliffs and ponds) based on the differencing of high-resolution digital elevation models (DEMs) and energy-balance modelling (Thompson et al., 2016; Buri et al., 2016, 2021; Brun et al., 2018; Miles et al., 2018, 2022; Mishra et al., 2022; Sato et al., 2021; Kneib et al., 2022). However, some studies have found that thicker debris cover has a larger effect on total thinning than enhanced ice ablation from ice cliffs and supraglacial ponds area (e.g. Hambrey et al., 2008; Vincent et al., 2016; Brun et al., 2018; Anderson et al., 2021a). Additionally, glacier dynamics play an essential role in debris-covered glacier elevation changes, and the rapid thinning of debris-covered glaciers is speculated to be partly caused by declining ice flow (Berthier and Vincent, 2012; Nuimura et al., 2017; Brun et al., 2018; Anderson et al., 2021a, b; Rounce et al., 2021). However, this has been evaluated with high accuracy at very few sites, none of which are in the southeast Tibetan Plateau.

These hypotheses (additional melt at areas influenced by cliffs or supraglacial ponds or extra thinning from reduced ice supply) therefore need to be supported with high-resolution data to account for the local effects of these melt hotspots. Since ice cliffs and supraglacial ponds tend to have a relatively small area and have a high interannual and interseasonal variability (Miles et al., 2017; Kneib et al., 2021), it is challenging to accurately track their evolution over time and quantify their contribution to melt (Mishra et al., 2022; Kneib et al., 2022). Because of these challenges, studies that have investigated their spatial distribution, dynamics, and melt contribution have typically utilized unpiloted aerial vehicles (UAVs) or very high-resolution optical satellite data (Immerzeel et al., 2014; Brun et al., 2016, 2018; Mölg et al., 2019; Anderson et al., 2021a; Kneib et al., 2021, 2022; Mishra et al., 2022; Sato et al., 2021), which can infer detailed observations of local processes (e.g. Westoby et al., 2020). Compared with the use of UAVs, in situ observations have the shortcoming of limited spatial representation, while remote sensing data have insufficient spatial and temporal resolution and are vulnerable to cloudy and rainy weather. Therefore, UAV technology is widely applied in glaciological studies, including in debris-covered glacier settings (Hugenholtz et al., 2013; Immerzeel et al., 2014; Kraaijenbrink et al., 2016; Wigmore and Mark, 2017; Fugazza et al., 2018; Rossini et al., 2018; Bash and Moorman, 2020; Westoby et al., 2020; Cao et al., 2021; Mishra et al., 2022; Xu et al., 2022). High-precision DEMs and orthophotos can be obtained relatively easily from UAV images processed using structure from motion (UAV-SfM) with multi-view stereo photogrammetry (Westoby et al., 2012; Benoit et al., 2019). However, few studies have estimated the seasonal or annual surface mass balance of debris-covered glaciers based on repeated UAV data.

Here, we systematically compare the glacier change patterns of two neighbouring debris-covered glaciers, 23K Glacier and 24K Glacier, in the southeastern Tibetan Plateau for the period of 13 August 2019 to 22 October 2020 using change detection applied to high-resolution repeated DEMs and orthoimages acquired via UAV-SfM surveys and in situ measurements. The glaciers are located in the same catchment and climatic setting, but the topography of the glaciers, as well as their dynamic behaviour and supraglacial debris thickness, differs considerably. The objective of this study is to explore the factors that control the interglacier variability in surface thinning and surface mass balance patterns of these two glaciers, with a view to advance the understanding of the key mechanisms that control debris-covered glacier changes in the southeastern Tibetan Plateau and to assess the different hypotheses at the glacier scale that contribute to the anomalous thinning of debris-covered glaciers.

2 Study area

The 23K (∼4 km²) and 24K (∼2 km²) glaciers are located in the southeastern Tibetan Plateau (∼29.77^∘ N, 95.70^∘ E; Fig. 1) and are mainly affected by two streams of humid air: the Bay of Bengal vortex (in spring) and the Indian summer monsoon system (in summer), respectively (Ye and Gao, 1979; Yang et al., 2013, 2016). Thus, the monthly precipitation distribution exhibits a double-peak pattern occurring in both spring and summer (Yang et al., 2013). This is significantly different from the glacier regime on the Tibetan Plateau, which is the so-called “summer accumulation” type (Fujita and Ageta, 2000; Maussion et al., 2014). Regional geodetic mass balance studies indicate that the magnitude of ice loss in the region has exceeded the average for High Mountain Asia (Kääb et al., 2012; Brun et al., 2017; Shean et al., 2020; Hugonnet et al., 2021a) in the past 20 years, also affecting debris-covered glaciers. The 23K and 24K glaciers are located on the northern slopes of the Gangrigabu mountains, ∼23 and 24 km, respectively, from Bomi City (Yang et al., 2017). The 23K Glacier spans an altitudinal range of 3760 to 5437 m a.s.l. and flows initially towards the southeast before turning to the northeast in the ablation zone (Fig. 1a). In contrast, the 24K Glacier flows to the northwest and spans an altitudinal range of 3900 to 5621 m a.s.l. Both glaciers are partly covered by a layer of rock debris (Fig. 1a and c), and the debris-covered area represents approximately 34 % and 41 % of the total area for the 23K Glacier and 24K Glacier, respectively, based on satellite data (Pléiades-1A false-colour image from 20 September 2021, 2 m resolution). The terminus retreat patterns of the two glaciers are also conspicuously different: the terminus of the 23K Glacier appears largely stagnant and is enclosed by a latero-terminal moraine complex, while the terminus of the 24K Glacier presents a large ice cliff, which is bounded by lateral moraines (Fig. 1c and d). Data from an automatic weather station (AWS, 3900 m a.s.l., running between June and September 2016) on the 24K Glacier indicate a warm and humid climate, with mean temperature and total precipitation reaching ∼9 ^∘C and ∼1700 mm, respectively (Yang et al., 2017; Fugger et al., 2022).

Figure 1 (a) Overview of the 23K Glacier and 24K Glacier basins including the UAV survey areas, accuracy test areas, in situ measurement locations (ground validation points, ground penetrating radar,debris thickness measurements), surface mass balance estimation zones (zones A–F), and glacier outlines (background image Pléiades-1A false-colour image from 20 September 2021). © CNES 2021, Distribution Airbus D&S. (b) Geographic location of the study area. (c–d) UAV-derived three-dimensional images of the two glaciers. (e–f) Spatial distribution of the surface mass balance estimation zones of the two glaciers (e: 23K Glacier; f: 24K Glacier), debris thickness measurement points (in June 2019), and distribution of ice cliffs and supraglacial ponds based on 2020 orthomosaics.

3 Data and methods

3.1 UAV flights and data processing

Optical imagery of the glacier surface and its immediate surroundings was acquired using UAVs during four field campaigns conducted between 13 August 2019 and 22 October 2020 (Table 1). We used an eBee Plus aircraft with built-in global navigation satellite system (GNSS) post-processed kinematic (PPK) functionality in October 2019, August 2020, and October 2020 surveys and a DJI Phantom 4 real-time kinematic (RTK) aircraft in an August 2019 survey to capture the high-resolution annual and seasonal patterns of glacier thinning and surface displacement.

The eBee Plus (Fig. S1a in the Supplement) is a fixed-wing UAV that has a 20-megapixel RGB digital compact camera (SenseFly S.O.D.A.). The flight management software eMotion3^® was used for flight planning. The DJI Phantom 4 RTK is a rotary-wing UAV (Fig. S1c) equipped with a 2.54 cm (1 in.), 20-megapixel CMOS camera. In this study, all longitudinal and lateral image overlaps were set to 65 % and 80 %, respectively. The flight lines for the eBee Plus and DJI Phantom 4 RTK aircraft maintained a relatively constant survey height above the glacier surface, which resulted in a constant ground resolution for each survey. A Huaxing A10 GNSS GPS was used as a static base station (Fig. S1b; fixed position for different surveys), and these data were attached to exchangeable image file format (EXIF) metadata of every geotagged image (Yang et al., 2020) and thereby integrated into a post-processed kinematic (PPK) correction workflow to improve the accuracy of the UAV-SfM reconstruction. The DJI Phantom 4 RTK aircraft was permanently connected to a GNSS receiver (D-RTK 2, Fig. S1d) so that each survey image already had its high-precision position information embedded (i.e. no post-processing required). The geotagged UAV images were used to create orthomosaics and DEMs using the SfM-based photogrammetric software Pix4Dmapper version 4.3.31.

To assess the final positional accuracy of the UAV-derived outputs, seven ground validation points were laid out in the vicinity of in situ measurements (Fig. 1a) on the surface of the 24K Glacier in October 2020 (synchronized with the October 2020 UAV survey). We used a Huaxing A10 differential Global Positioning System (dGPS) to measure the position (XYZ) of each ground validation point. The horizontal accuracy of the dGPS is ±8 mm, and the vertical accuracy is ±15 mm. By comparing the ground validation point measurements with the UAV-derived orthomosaics and DEMs, we obtain absolute XYZ accuracies for the October 2020 UAV-SfM survey product (Table 2). While ground validation points were not used for the other UAV-SfM surveys, the accuracies of all UAV-SfM products were indirectly assessed by comparing the horizontal (XY) and vertical errors (Z) between all orthomosaics and DEMs (Table 3). The horizontal (XY) error was estimated by measuring the displacements of 25 benchmark boulders located on stable off-glacier terrain. For the vertical (Z) error, we calculated the elevation difference over stable terrain, as outlined in Fig. 1a.

Table 1 UAV photogrammetric flights used for the acquisition of glacier images.

3.2 In situ measurements

To evaluate the role of debris on the surface mass balance patterns of both glaciers, debris thickness was manually measured by digging pits at selected sites in June 2019 (Fig. 1d and e). For the 23K Glacier, a total of 157 locations were measured within the elevation range of 3740–3880 m a.s.l. For the 24K Glacier, we measured a total of 349 points from the glacier terminus to 4250 m a.s.l. The uncertainty of the manual measurement of the debris thickness was assumed to be 2 cm.

In October 2019, three ice thickness cross-sections (two in the debris-covered area and one in the debris-free area, Figs. 1a and S2) were measured on the 24K Glacier using a Kentech ground penetrating radar (GPR) monopulse transmitter with 2.5 MHz antennas. The in situ ice thickness measurements were used to correct the 24K Glacier distributed ice thickness from Farinotti et al. (2019) using a linear regression (Fig. S2; Kneib et al., 2022).

3.3 Thinning patterns, glacier dynamics, and glacier driving stress

We compared the thinning and velocity patterns of the 23K and 24K glaciers on annual and seasonal timescales. The period between August 2019 and August 2020 (373 d, mean temperature at AWS: 1.66 ^∘C) was selected for the annual timescale analysis. The annual rates were adjusted according to the ratio of days (366 divided by 373) to make this result (13 August 2019–20 October 2020, 373 d) closer to that of the natural year (366 d). For the seasonal analysis, we refer to the period 12 October 2019–20 August 2020 (313 d, mean temperature: 0.42 ^∘C) as the “cold period” and the period 20 August 2020–22 October 2020 (63 d, mean temperature: 8.68 ^∘C) as the “warm period”. We applied the post-processed kinematic and real-time kinematic GPS correction technologies and ensured that the static base stations were fixed at the same location for different surveys, which led to DEMs with only very small minor offsets in XYZ (Yang et al., 2020). Therefore, we did not perform the co-registration of DEMs for the thinning calculation. For each period the thinning pattern was obtained by DEM differencing in ArcGIS 10.4.

A spatially distributed estimate of XY surface displacements was obtained for each period by applying a normalized cross-correlation algorithm to the multi-directional 0.15 m resolution DEM hillshades using ImGRAFT (Messerli and Grinsted, 2015). A search window of 10×10 pixels (1.5×1.5 m) was used to compute the magnitude and direction of the displacement vectors. Surface displacements greater than 30 m were considered as noise and were filtered out. Due to the high-resolution velocity data and the small number of gaps (<5 %), we interpolated the velocity values in the data gaps using nearest neighbour interpolation.

The glacier driving stress was calculated as (Cuffey and Paterson, 2010)

$$\begin{array}{}\text{(1)}& D_{\mathrm{stress}}={\mathit{\rho }}_{\mathrm{i}}\cdot g\cdot h\cdot \sin {\mathit{\alpha }}_{\mathrm{s}},\end{array}$$

where D_stress is the glacier driving stress (kPa), ρ_i is the density of ice (917 kg m⁻³), g is the gravitational acceleration (9.81 m s⁻²), h is the ice thickness (in m), and α_s is the glacier slope obtained from the AW3D30 (30 m resolution) DEM (Tadono et al., 2014) smoothed with a Gaussian filter (8 pixel window) to remove effects from local variations in surface topography (Brun et al., 2018).

3.4 Surface mass balance of the UAV survey area

To investigate the magnitude and distribution of the surface mass balance, each glacier was divided into six zones (A–F), and each zone was used as a separate section for the surface mass balance estimation. Each surface mass balance estimation zone was outlined manually, perpendicular to the main glacier flow line and with a similar area (Fig. 1a). To extract the melt contribution from the ice cliffs and ponds, we flow-corrected the DEMs following Brun et al. (2018). To perform the flow-correction, we used orthomosaics and DEMs from the August 2020 UAV-SfM surveys as a reference for each glacier, to which the August 2019, October 2019, and October 2020 UAV-derived data were flow-corrected. Corrections in XY were made using the ArcGIS georeferencing tool to manually track the surface flow for each given period by using surface tie points (mainly large boulders). For each correction period, the number of tie points was sufficiently large (>75), and the tie points were well distributed across the glacier surface to ensure a spatially representative correction, using a spline-based transformation. The resulting XY flow-corrected elevation change (dh_c, in m) is therefore equal to the sum of the surface lowering ($\dot{b}$, in m) from the surface mass balance (in m w.e.) by multiplying the conversion factor of 1.09 (the ratio of densities of water to ice) and the mean vertical displacement (ω, in m):

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{d}{h}_{\mathrm{c}}=\dot{b}+\mathit{\omega },\text{(3)}& {\displaystyle }\mathit{\omega }={\mathit{\omega }}_{\mathrm{s}}+{\mathit{\omega }}_{\mathrm{e}},\text{(4)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{s}}=u_{\mathrm{s}}\cdot \tan \left({\mathit{\alpha }}_{\mathrm{m}}\right),\text{(5)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{e}}={\displaystyle \frac{\mathrm{\Delta }q}{A}},\text{(6)}& {\displaystyle }\mathrm{\Delta }q=q_{n+\mathrm{1}}-q_n,\text{(7)}& {\displaystyle }q=\mathit{\mu }\cdot h_q\cdot v\cdot l,\end{array}$$

where ω_s (in m) corresponds to the elevation change resulting from the horizontal flow-correction of the DEMs, u_s (in m) is the mean horizontal surface displacement, and α_m (in ^∘) is the mean surface slope of a given zone. ω_e (in m) corresponds to the flux divergence. A (in m²) is the area of the zone, and Δq (in m³) is the ice flux difference at a given period. q_n+1 and q_n are the ice flux entering and leaving the estimation zone at a given period. q (in m³) is the ice flux through a glacier cross-section. μ is a coefficient to convert the surface velocity into a depth-averaged velocity which we assumed to be equal to 0.9 following Miles et al. (2018), h_q (in m) is the ice thickness for the corresponding cross-section, v (in m) is the surface displacement component normal to the cross-sections line, and l (in m) is the width of the cross-section. The following equation was used to evaluate the uncertainty of the surface mass balance (${\mathit{\sigma }}_{\dot{b}}$) for each zone:

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_{\dot{b}}=\sqrt{{\mathit{\sigma }}_{\mathrm{d}{h}_{\mathrm{c}}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathit{\omega }}^{\mathrm{2}}}.\end{array}$$

The uncertainties of dh_c and u_s (also for v) were obtained by calculating the mean surface elevation difference and displacement from the non-glacial test areas. They were determined to be 0.09 and 0.25 m, respectively, by averaging the values over all periods.

The equations below were applied in the uncertainty of the vertical component of velocity (σ_ω):

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathit{\sigma }}_{\mathit{\omega }}=\sqrt{{\mathit{\sigma }}_{u_{\mathrm{s}}}^{\mathrm{2}}+{\mathit{\sigma }}_{{\mathit{\alpha }}_{\mathrm{m}}}^{\mathrm{2}}+{\mathit{\sigma }}_{{\mathit{\omega }}_{\mathrm{e}}}^{\mathrm{2}}},\text{(10)}& {\displaystyle \frac{{\mathit{\sigma }}_{{\mathit{\omega }}_{\mathrm{e}}}}{{\mathit{\omega }}_{\mathrm{e}}}}=\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{\Delta }q}}{\mathrm{\Delta }q}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_A}{A}}\right)}^{\mathrm{2}}}.\end{array}$$

For the uncertainty associated with the slope correction, we assumed a 2^∘ uncertainty in the slope (α_m). We obtained the uncertainty of the zone area (A) using ±20 m buffers from the original outlines (Bolch et al., 2010; Miles et al., 2018). The uncertainty of the ice flux difference (σ_Δq) was calculated as

$$\begin{array}{}\text{(11)}& {\mathit{\sigma }}_{\mathrm{\Delta }q}=\sqrt{{\mathit{\sigma }}_{q_{n+\mathrm{1}}}^{\mathrm{2}}+{\mathit{\sigma }}_{q_n}^{\mathrm{2}}},\end{array}$$

where σ_q (m³) is the uncertainty of ice flux through a glacier cross-section, given by

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{{\mathit{\sigma }}_q}{q}}=\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_v}{v}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\mu }}}{\mathit{\mu }}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{h_q}}{h_q}}\right)}^{\mathrm{2}}},\end{array}$$

where the uncertainty of the ratio μ (column-averaged velocity divided by surface velocity) is assumed to be 0.1 (Cuffey and Paterson, 2010; Miles et al., 2018). The uncertainty in h_q is ∼10–35 m (26 %) for the 23K Glacier (Farinotti et al., 2019) and is assumed to be equal to 12 m for the corrected 24K Glacier (Fig. S2a) ice thickness.

3.5 Ice cliffs and supraglacial ponds outlines and melt contribution

In this study, the proportion of ice cliffs and supraglacial ponds area within each surface mass balance estimation zone was calculated to better understand the surface mass balance patterns on both glaciers. We manually extracted ice cliffs and supraglacial ponds outlines from the flow-corrected orthomosaics (August 2020 and October 2020, August 2019 and August 2020, October 2019 and August 2020). The outlines of the ice cliffs and supraglacial ponds area for each zone were combined into one polygon feature by taking the union of the areas of these features (Brun et al., 2018). In this way, we could effectively measure the surface elevation and mass balance changes at these locations, thereby providing us with an upper-bound estimate of the total contribution of ice cliffs and ponds to the surface mass balance (Kneib et al., 2022; Mishra et al., 2022). We conservatively assumed the uncertainty of the area density of ice cliffs and supraglacial ponds to be equal to 0.5 % (Brun et al., 2018; Steiner et al., 2019). We also defined the hotspots area enhancement factor as the ratio of melt rate of the ice cliffs and supraglacial ponds area to the sub-debris area melt rate (Miles et al., 2022).

4 Results

4.1 UAV product accuracy

The uncertainties in the UAV measurements (22 October 2020) reported as the mean absolute deviation in the X, Y, and Z directions are 0.14±0.11 m, 0.09±0.11 m, and 0.24±0.18 m, respectively (Table 2). The vertical uncertainty is twice as high as the horizontal uncertainty, which is in agreement with the findings of other studies (James et al., 2017; Li et al., 2019).

Table 2 XYZ geolocation accuracy of the UAV-SfM orthomosaic and DEM for the two glaciers based on the dGPS measurements of the seven ground validation points.

The average XY errors for the three survey pairs (August 2019 vs. October 2019, October 2019 vs. August 2020, and August 2020 vs. October 2020) were 0.07±0.06 and 0.08±0.06 m, respectively (Table 3). The relative vertical (Z) errors between the DEMs were ≤0.09 m for all pairs, which is appropriate for resolving fine-scale surface change for glaciological analysis.

Table 3 XYZ errors between different orthomosaics and DEMs for the two glaciers using the fixed benchmark boulders and validation areas.

4.2 Spatio-temporal variability of thinning patterns

The thinning during different periods (annual: August 2019–August 2020, cold: October 2019–August 2020, warm: August 2020–October 2020) for the two glaciers is shown in Fig. 2 and Table 4. The annual thinning rate of the 23K Glacier in the survey area was $-\mathrm{2.3}\pm \mathrm{0.1}$ m a⁻¹, whereas the thinning rate of the 24K Glacier was $-\mathrm{1.2}\pm \mathrm{0.1}$ m a⁻¹ (Fig. 2a and d). During the cold period (October 2019–August 2020), the average surface elevation changes over the 23K and 24K glacier survey areas were $-\mathrm{1.5}\pm \mathrm{0.1}$ and $-\mathrm{0.2}\pm \mathrm{0.1}$ m, with an average daily thinning of $-\mathrm{0.5}\pm \mathrm{0.03}$ and $-\mathrm{0.1}\pm \mathrm{0.03}$ cm d⁻¹, respectively (Fig. 2b and e). During the warm period (August 2020–October 2020), the magnitude of the thinning of both glaciers was very similar: $-\mathrm{0.7}\pm \mathrm{0.1}$ and $-\mathrm{1.0}\pm \mathrm{0.1}$ m, with an average daily thinning of $-\mathrm{1.2}\pm \mathrm{0.03}$ cm d⁻¹ for the 23K Glacier and $-\mathrm{1.6}\pm \mathrm{0.03}$ cm d⁻¹ for the 24K Glacier, respectively (Fig. 2c and f).

Table 4 Total thinning and its daily rate for 23K and 24K glaciers during the different periods.

Figure 2 Annual surface elevation changes between UAV-derived DEMs for August 2019–August 2020 (a, d) and thinning rates for the cold period (b, e) and the warm period (c, f). Upper panels refer to the 23K Glacier and bottom panels refer to the 24K Glacier.

In Fig. 3, we present the relationships between thinning rates and altitude for each glacier at 5 and 15 m intervals. The absolute thinning of the 23K Glacier increases with altitude based on the Mann–Kendall test (i.e. a negative gradient; Z value: −3.75 to −3.25) for every period of the analysis, while the thinning of the 24K Glacier decreases with altitude at the annual scale and in the cold period (positive gradient; Z value: +3.65 to +3.69). In particular, the thinning of the 24K Glacier in the warm period shows the opposite altitudinal trend to the annual and cold periods (i.e. consistent with the 23K Glacier; Z value: −0.77).

Figure 3 Annual average glacier surface elevation change rates within 5 m (23K Glacier) and 15 m (24K Glacier) elevation bands (dots) with the corresponding standard deviations (horizontal error bar) for August 2019–August 2020 (a, d), the cold period (b, e), and the warm period (c, f) across the monitoring area of the two glaciers. The shadowed red sections represent the terminal ice cliff on the 24K Glacier.

4.3 Glacier dynamics

Between August 2019 and August 2020, the mean surface velocity of the UAV survey domain was 1.7±0.2 m a⁻¹ for the 23K Glacier and 9.2±0.2 m a⁻¹ for the 24K Glacier (Fig. 4a and d). There is a zone of stagnation at the terminus of the 23K Glacier where the surface velocity is less than 2 m a⁻¹ (from terminus to 900 m up-glacier, Fig. 4a). During the cold period, the mean surface velocity (UAV survey domain) of the 23K Glacier was 0.4±0.03 cm d⁻¹, while it was 2.4±0.08 cm d⁻¹ for the 24K Glacier (Fig. 4b and e). During the warm period, the mean surface velocities of the UAV survey domains for the 23K Glacier and 24K Glacier were 0.6±0.14 and 3.0±0.40 cm d⁻¹, respectively (Fig. 4c and f). For the UAV survey domain, the mean driving stress of the 23K Glacier was 1.6×10⁵ Pa and that of the 24K Glacier was 2.1×10⁵ Pa, ∼30 % higher than that for the 23K Glacier (Fig. 4g). The estimated emergence velocities agreed with the patterns of surface velocities and driving stress, with the 24K Glacier having significantly higher emergence velocities than the 23K Glacier at any period (Fig. 4h). The annual emergence velocity of the 23K Glacier was 0.18±0.04 m a⁻¹, approximately 7.6 times lower than that of the 24K Glacier (1.36±0.14 m a⁻¹). The mean emergence velocity during the cold period for the 23K (resp. 24K) Glacier was 0.05±0.01 cm d⁻¹ (0.36±0.04 cm d⁻¹), and the emergence velocity during the warm period for the 23K (resp. 24K) Glacier was 0.11±0.05 cm d⁻¹ (0.29±0.06 cm d⁻¹). The 24K Glacier had consistently higher emergence velocity than the 23K Glacier during all periods.

Figure 4 Average surface velocity of 23K and 24K glaciers for August 2019–August 2020 (a, d), the cold period (b, e), and the warm period (e, f). Spatial distribution of driving stresses (August 2019–August 2020) of the two glaciers (g1) and central flowlines (a–a^′ and b–b^′) and the profile of driving stresses (g2). The background is Pléiades-1A false-colour image from 20 September 2021. © CNES 2021, Distribution Airbus D&S. Emergence velocity as a function of altitude for the two glaciers (h).

4.4 Surface mass balance patterns

At the annual scale, the mean surface mass balances for the 23K and 24K glacier survey areas were $-\mathrm{2.50}\pm \mathrm{0.11}$ m w.e. a⁻¹ ($-\mathrm{0.67}\pm \mathrm{0.03}$ cm w.e. d⁻¹) and $-\mathrm{2.76}\pm \mathrm{0.34}$ m w.e. a⁻¹ ($-\mathrm{0.74}\pm \mathrm{0.09}$ cm w.e. d⁻¹), respectively, and are therefore not significantly different (Table 5). During the cold period, the glacier mass balance was $-\mathrm{0.49}\pm \mathrm{0.03}$ cm w.e. d⁻¹ for the 23K Glacier and $-\mathrm{0.42}\pm \mathrm{0.09}$ cm w.e. d⁻¹ for the 24K Glacier (Fig. 5). In contrast, during the warm period, the mass balance of the 23K Glacier was $-\mathrm{1.25}\pm \mathrm{0.15}$ cm w.e. d⁻¹, while for the 24K Glacier it was $-\mathrm{1.89}\pm \mathrm{0.18}$ cm w.e. d⁻¹; the surface mass balance during the warm period for the 24K Glacier is thus ∼46 % larger than that for the 23K Glacier (Fig. 5). The surface mass balance values of each mass balance zone (A–F) for both glaciers exhibit a weak decreasing trend with altitude in all periods (Fig. S3).

To better evaluate the role of emergence velocity replenishment on thinning, we calculated the ratio of ice emergence to melt (Table 5). A higher value of this ratio indicates a greater influence of emergence velocity replenishment on thinning (a value of 1 indicates perfect compensation of ablation by ice flow). The ratios of annual ice emergence and annual melt for the 23K and 24K glaciers are 0.09 and 0.49, respectively. The ratio values for the 24K Glacier are always higher than those for the 23K Glacier, and this is especially evident in the non-ablation period.

Figure 5 Mean surface mass balance rates and their uncertainties for both UAV survey domains during the annual observation period (August 2019–August 2020), the cold period (October 2019–August 2020), and the warm period (August–October 2020).

Table 5 Ice emergence, melt, and their ratio in all periods for both glaciers.

The density of ice cliffs and supraglacial ponds for the 23K Glacier is 6.8 %–7.2 % and is 4.4 %–5.1 % for the 24K Glacier (Table 6). There are few ponds (∼5, area >100 m²) on the 23K Glacier, while there is no ponded area on the 24K Glacier. The average debris thicknesses in the UAV survey areas of the 23K Glacier and 24K Glacier are 47.1±2 cm and 24.2±2 cm, respectively. The debris thickness of both glaciers decreases with increasing altitude in the UAV survey area (23K Glacier: $\sim -\mathrm{57}$ cm 100 m⁻¹; 24K Glacier: $\sim -\mathrm{9}$ cm 100 m⁻¹). To disentangle the influence of the sub-debris area and of the ice cliffs and supraglacial ponds area on the surface mass balance pattern, respectively, we extracted the surface mass balance rates of the area of ice cliffs and supraglacial ponds and of the sub-debris area in both glaciers for all periods (Table 6). As shown in Table 6, the surface mass balance rates in the ice cliffs and supraglacial ponds area of both glaciers are higher than those in the sub-debris area during all periods (i.e. the ice cliffs and supraglacial ponds area have a higher melt efficiency). In this study, the enhancement factors for the two glaciers range from 1.6 to 4.4 during all periods. They are significantly higher in the warm period (23K: ∼4.4; 24K: ∼2.6) than at the annual scale (23K: ∼2.2; 24K: ∼1.6) or during the cold period (23K: ∼2.6; 24K: ∼1.7), indicating that the hotspots effect is more pronounced during the ablation period. By extracting the ablation contribution of the ice cliffs and supraglacial ponds (Fig. 6a and b), we found that the total melt from ice cliffs and supraglacial ponds area accounted for 31.5±2.2 % of the total ablation in the UAV survey area of the 23K Glacier (Fig. 6c) and 11.4±1.3 % for the 24K Glacier during the warm period (Fig. 6d).

Table 6 Area proportions of ice cliffs and supraglacial ponds, mean surface mass balance for ice cliffs and supraglacial ponds area, and sub-debris area, and their enhancement factors in all periods for both glaciers.

Figure 6 Spatial distribution of surface mass balance during the warm period for the UAV survey domains of the 23K Glacier and 24K Glacier (a, b) and the conceptual diagrams of vertical motion components for the 23K Glacier and 24K Glacier (c, d).

To better disentangle the effects of debris thickness, ice cliffs, and supraglacial ponds on glacier melt, we compared the mean debris thickness, the percentage of ice cliffs and supraglacial ponds, and the mean melt rate during the warm period in each of the six zones (A–F) of each glacier (Fig. 7). We assessed the correlations between the zonal surface mass balance in the warm period and the zonal debris thickness as well as the ice cliffs and supraglacial ponds density. For the 23K Glacier, the correlation coefficient r between debris thickness and surface mass balance during the warm period is 0.88 (p value is 0.02 at the 95 % confidence level), indicating that the debris thickness is highly correlated with melt. In contrast, the correlation coefficient between the percentage of ice cliff and supraglacial ponds area and the surface mass balance is −0.29 (p value is 0.58 at the 95 % confidence level), indicating the absence of correlation. During the warm period, the 24K Glacier also exhibits a strong relationship between debris thickness and melt, with an r of 0.82 (p value is 0.05 at the 95 % confidence level). The ice cliffs and supraglacial ponds area is weakly correlated with glacier melt ($r=-\mathrm{0.48}$, p value is 0.36 at the 95 % confidence level).

Figure 7 Daily surface mass balance during the warm period (solid black line), mean debris thickness (solid blue line), and mean percentage of ice cliffs and supraglacial ponds area (solid red line) for the individual zones of the two glaciers.

4.5 Glacier change in the early twenty-first century

The thinning of the 23K Glacier tongue has been greater than that of the 24K Glacier tongue since 2000 (Hugonnet et al., 2021a; Fig. 8a–d), which is generally in agreement with the thinning patterns we derived from the UAV surveys. Both glaciers (UAV survey domains) show an inverted relationship between thinning and altitude at all periods of the last 2 decades (Fig. 8e). Both glaciers experienced increases in thinning rates between the 2000–2004 and 2015–2019 periods (increase rate for 23K: +77 %; 24K: +100 %), although the increase in thinning rate was slightly lower for the 24K Glacier than for the 23K Glacier over the last decade (23K: +44 %; 24K: +20 %).

Figure 8 Annual surface elevation changes (Hugonnet et al., 2021a) for 2000–2004 (a), 2005–2009 (b), 2010–2014 (c), and 2015–2019 (d) for the 23K Glacier and 24K Glacier (dashed red rectangular box represents the UAV aerial survey area). Annual surface elevation changes and their uncertainties for both UAV survey domains at different periods (e).

Using the velocity data by Dehecq et al. (2015, 2019a, b), we explore changes in the glaciers' dynamic state by analysing the surface velocities over the last 2 decades (Fig. 9). The velocity magnitude and pattern obtained by Dehecq et al. (2015, 2019a, b) for the period 2013–2015 correspond well with our UAV results. The surface velocities of both glaciers (UAV survey domains) decreased significantly between the 1999–2003 and 2013–2015 periods. Between these periods, the surface velocities of the 23K and 24K glaciers decreased on average by 84 % and 54 %, respectively. However, the velocity magnitude and pattern for 2013–2015 correspond to our UAV results. The 24K Glacier tongue is still replenished by ice flux at present, which compensates for its higher melt due to thinner debris.

A prior study analysed the mean surface mass balance from 2000–2016 based on regionally available datasets (Miles et al., 2020, 2021), with values of $-\mathrm{1.5}\pm \mathrm{0.7}$ m w.e. a⁻¹ for the UAV survey domain of the 23K Glacier and $-\mathrm{1.6}\pm \mathrm{1.1}$ m w.e. a⁻¹ for the UAV survey domain of the 24K Glacier. These melt values for the period of 2000–2016 are a little less negative than our results for the period of 2019–2020, which might be attributed to warming or to changes in precipitation phase in the southeast Tibetan Plateau (e.g. Jouberton et al., 2022).

Figure 9 Annual surface velocities (Dehecq et al., 2015, 2019a, 2019b) for 1999–2003 (a) and 2013–2015 (b) for the 23K Glacier and 24K Glacier (dashed red rectangular boxes represent the UAV aerial survey areas). Annual surface velocities and their uncertainties for both UAV survey domains at different periods (c).

5 Discussion

5.1 Controls on the thinning patterns

The thinning patterns of the two glaciers are different at the annual scale and during the cold period, as indicated by their magnitude and spatial distribution (Fig. 2). Some studies have highlighted that large differences in thinning can be caused by different dynamic states (Fig. 4; Brun et al., 2018; Anderson et al., 2021a, b; Rounce et al., 2021). In this study, the emergence velocity replenishment is consistently higher on the 24K Glacier than on the 23K Glacier for all periods (Fig. 4h). The replenishment of ice into the ablation area by ice flow is crucial to the glaciers' long-term sustainability, and net annual ablation exceeds ice resupply for the study areas of both glaciers. We therefore assessed the ratio of ice emergence to melt as an indicator of the local balance between ablation and ice supply (a direct local metric of glacier health), and this study illustrates that glacier health can vary greatly even over small distances (Table 5). The 23K and 24K glaciers experience the same climatic forcing, but their distinct geometries lead to different expressions of ice dynamics. For instance, the mean slope of the 24K Glacier's UAV survey domain (∼10.2^∘) is consistently higher than that of the 23K Glacier (∼4.6^∘), which explains its higher driving stress (Fig. S4) and faster surface velocity. In addition, the different glacier geometries may lead to distinct dependence on avalanche and rockfall mass supply, which could also lead to marked differences in mass supply to the glacier terminus. Using the ratio of ice emergence to melt, we can identify, for each season, whether ice resupply or melt is the main factor leading to the thinning rates. In Table 5, we clearly see that at the 23K Glacier the ice supply is considerably smaller than the absolute value of melt for all periods, indicating that the melt is directly responsible for contemporary thinning patterns. The 24K Glacier also exhibits a strong imbalance between ice emergence and melt over the warm period (0.15), but ice emergence nearly compensates for melt during the cold period (0.87). Thus, the 24K Glacier exhibits a relatively healthier cold season and annual balance than the 23K Glacier. In contrast to the annual and cold periods, the thinning patterns of the two glaciers are similar during the warm period (i.e. the magnitude of thinning increases with altitude). In this period, the 24K Glacier ice emergence only represents a small fraction of the melt (Table 5), and the pattern of surface mass balance is clearly driven by the spatial distribution of debris thickness. Therefore, the thinning pattern of the 24K Glacier goes from being controlled primarily by ice dynamics (for the annual and cold period) to being controlled by debris thickness during the warm period.

Fu et al. (2022) also estimated thinning for a debris-covered glacier (Hailuogou) in the southeastern Tibetan Plateau. They revealed that the tongue area showed considerable thinning (−2.81 m) during the ablation period (June 2018–October 2018) and that its ice resupply was very weak, similar to what we observed for the 23K Glacier. He et al. (2023) also observed the debris-covered tongue of Zhuxi Glacier in the southeastern Tibetan Plateau by UAV for the period of 2020–2021. The lower tongue of this glacier is similarly characterized by high rates of thinning (>1 m a⁻¹) and slow movement (<7 m a⁻¹). Other studies concerning the mass balance of Himalayan debris-covered glaciers also report on glaciers that possess a high thinning rate ($\sim -\mathrm{0.9}$–1.8 m a⁻¹) and weaker dynamic state (Vincent et al., 2016; Nuimura et al., 2017; Brun et al., 2018; Rowan et al., 2021) and that possess thinning patterns that are consistent with those of the 23K Glacier. In summary, the thinning pattern for the 23K Glacier appears to conform to that of other Himalayan debris-covered glaciers, while that of the 24K Glacier is perhaps more anomalous due to its relatively healthy ice supply all the way to the terminus.

In situ observations show that the terminus type differs between each glacier (Figs. 1c, d and S5), which can have a strong influence on ice dynamics for debris-covered glaciers (Anderson and Anderson, 2016). The terminus of the 23K Glacier appears largely stagnant and is enclosed by a latero-terminal moraine complex, while the terminus of the 24K Glacier exhibits a large ice cliff bounded by lateral moraines. The mass of the 24K Glacier is lost due to progressive melt, and episodic “dry” calving promotes glacier retreat. Such terminal-cliff debris-covered glaciers are less well studied but have been noted in a variety of settings (e.g. Ferguson and Vieli, 2021). Terminus type can have an important influence on debris-covered glacier geometry (King et al., 2019) and thus may influence glacier geometric evolution in a warming climate. Our investigations suggest that the 24K Glacier exhibits a reduced climatic–geometric imbalance compared with the 23K Glacier, possibly because it has responded to climate warming with progressive retreat rather than thinning. In addition, the ablation area of the 23K Glacier is largely stagnant with a low thinning rate, lower mean surface velocities, and reduced driving stress, which favour the presence of supraglacial ponds that enable the persistence of relatively large pond-influenced cliffs (Quincey et al., 2007; Sakai and Fujita, 2010; Miles et al., 2017; Kneib et al., 2023). In contrast, due to a steeper longitudinal slope and stronger ice flux, the 24K Glacier is characterized by thinner debris, which, combined with the steeper gradient, allows for the development of supraglacial streams and corresponding stream-influenced cliffs (Mölg et al., 2020; Kneib et al., 2023).

5.2 Possible reasons for the surface mass balance patterns

The presence of a debris cover and its influence on ablation typically causes surface mass balance patterns for debris-covered glaciers to differ from those of largely debris-free glaciers due to the melt-buffering effect of supraglacial debris cover that exceeds a few centimetres in thickness (Østrem, 1959; Nakawo et al., 1999; Nicholson and Benn, 2006; Reid and Brock, 2010; Anderson and Anderson, 2016; Yang et al., 2017). Our results show inverted melt-season surface mass balance profiles for both the 23K and 24K glaciers. In addition, debris-covered glaciers tend to develop ice cliffs that enhance melt locally, even relative to clean ice (Sakai et al., 1998, 2002; Reid and Brock, 2014; Juen et al., 2014; Steiner et al., 2015; Buri et al., 2016, 2021; Miles et al., 2016, 2018, 2022; Kneib et al., 2022). Overall, the annual surface mass balances of the 23K Glacier and 24K Glacier are similar, but the rate of mass loss is higher for the 24K Glacier during the warm period due to its thinner debris cover. During all periods, the magnitude of ablation increases with elevation for both glaciers, with the same spatial distribution (Figs. 7, S3). We find that the correlation between the melt and the debris thickness distribution is strong for both glaciers (23K: r=0.88; 24K: r=0.82) during the warm period, while there is little correlation between the melt and the ice cliffs and supraglacial ponds area distribution (23K: $r=-\mathrm{0.29}$; 24K: $r=-\mathrm{0.48}$). Prior studies have also identified the effect of the debris thickness spatial distribution on the physical mechanisms of ablation (Mihalcea et al., 2008; Zhang et al., 2011; Reid et al., 2012; Juen et al., 2014; Zhang et al., 2016; Gibson et al., 2017; McCarthy et al., 2017).

The role of ice cliffs and supraglacial ponds area is not negligible, as ablation in the survey areas would be underestimated by 24.5±1.7 % (23K Glacier) and 7.0±0.7 % (24K Glacier), respectively, if the ice cliffs and supraglacial ponds area were not taken into account. Other studies have observed and simulated the contributions of ice cliffs to be 17 %–26 % of the total ablation for different glaciers, which is in agreement with our results (Brun et al., 2018; Anderson et al., 2021a; Buri et al., 2021). The combination of ice cliffs and supraglacial ponds densities (23K Glacier: 7.2±0.5 %, 24K Glacier: 4.4±0.5 %; Fig. 6a and b) and their enhancement factors (23K Glacier: 4.4, 24K Glacier: 2.6; Table 6) shapes the ablation contribution patterns of ice cliffs and supraglacial ponds together. Based on UAV and time-lapse camera observations, Kneib et al. (2022) carried out high-precision observations of a single ice cliff on the 24K Glacier and measured the daily cliff melt during the ablation period at 3.9–5.1 cm d⁻¹. In this study, the daily melt rate of the 24K Glacier ice cliffs and supraglacial ponds area during the warm period was estimated to be 4.4±0.2 cm w.e. d⁻¹ (∼4.9 cm d⁻¹), which is similar to the observed value for a single ice cliff. We note that in this study the outlines of ice cliffs and supraglacial ponds area as digitized for warm periods were based on the union of the cliffs and supraglacial ponds' outlines (Brun et al., 2018), which may lead to an overestimation of the extent of the ice cliffs and supraglacial ponds area (Kneib et al., 2022). While the ice cliffs and supraglacial ponds are local controls of melt patterns, debris thickness is the dominant control on the altitudinal surface mass balance pattern for both the 23K and 24K glaciers, similar to Kennicott Glacier in Alaska (Anderson et al., 2021a).

This study provides insight into the controlling role of debris thickness on the ablation patterns of debris-covered glaciers. In future research, it will be beneficial to improve our understanding of the responses of debris-covered glaciers to climate change by focusing on the debris supply and evacuation differences. Based on field photography (Fig. S6), we also found that paraglacial slope failure events have recently occurred in the catchment of the 23K Glacier and may be the result of complex interactions between geologic structure and stress-related slope response to glacier mass loss (Zhong et al., 2022). Such events may become an increasingly important component of the debris supply and transport cascade for these land systems, with implications for the future development of supraglacial debris cover, glacier mass balance, and flow dynamics. Large rock avalanches that emplace in supraglacial environments have led to melt suppression and glacier advance in other locations (e.g. McSaveney, 1975; Shugar and Clague, 2011), and if sufficiently thick, this event could thus temporarily rejuvenate the emergence velocity into the terminus area.

6 Conclusions

We used multi-temporal high-resolution UAV-SfM surveys combined with in situ observations to quantify the seasonal thinning and surface mass balance of two neighbouring, but contrasting, debris-covered glaciers in the southeastern Tibetan Plateau. The conclusions are summarized as follows.

1. The thinning patterns of the two glaciers display distinct characteristics at an annual scale. The annual thinning of the 23K Glacier UAV survey area is 1.9 times greater than that of the 24K Glacier. The magnitude of 23K Glacier thinning increases with altitude, similar to many other debris-covered glaciers, whereas the 24K Glacier shows the opposite pattern, except during the melt season. These contrasting patterns are mainly driven by the stronger dynamic state of the 24K Glacier, which has a much higher down valley emergence velocity replenishment. However, the thinning pattern of the 24K Glacier is similar to that of the 23K Glacier during the warm period (i.e. the magnitude of thinning increases with altitude). In this period, the thinning of the 24K Glacier is mainly controlled by melt.

2. The surface mass balance patterns of the two glaciers are generally in agreement. The magnitude of the surface mass balance on both glaciers increases with altitude in the ablation period, exhibiting a melt inversion, which is attributable to the debris thickness distribution (which decreases with altitude). Due to the low areal proportion of ice cliffs and supraglacial ponds (melt hotspots), sub-debris ablation accounts for the majority of the total ablation.

3. Both glaciers experience accelerated thinning and reduced flow in their ablation area. The 23K Glacier possesses a higher thinning rate and is in a weaker dynamic state than the 24K Glacier, a pattern that is confirmed in the longer term.

High-resolution and comparative observations provide a rare perspective into the controls of the thinning and melt of two debris-covered glaciers. We provide evidence that the rate of mass loss on such glaciers can be highly dependent on the dynamic state and that the relatively thin debris on these glaciers is the main control of the glacier surface mass balance patterns, with supraglacial cliffs and ponds locally enhancing mass loss.