1. Introduction
Rates of glacier mass loss have increased in the Himalaya over the past 50 years: from 1974 to 2000, regional mass balance was estimated at −0.25 ± 0.09 m w.e. a−1, this increased to −0.39 ± 0.12 m w.e. a−1 from 2000 to 2015 (King and others, Reference King, Bhattacharya, Bhambri and Bolch2019). These negative mass balance trends can be attributed to rising temperatures. The Himalaya have experienced faster temperature increases relative to surrounding low-lying regions, resulting in higher ablation rates that have not been offset by increased precipitation (Kraaijenbrink and others, Reference Kraaijenbrink, Bierkens, Lutz and Immerzeel2017). Furthermore, these rates of change are not uniform. In the Eastern Himalaya, mass loss is nearly twice that of the central Himalaya and slightly higher than in the Western Himalaya (Bolch and others, Reference Bolch2012; Kääb and others, Reference Kääb, Berthier, Nuth, Gardelle and Arnaud2012; Xiang and others, Reference Xiang2021; Zhang and others, Reference Zhang, Xu, Liu, Ding and Zhao2023). The rate of mass balance change is also highest in Bhutan; mass balance has become increasingly negative, from −0.20 ± 0.08 m w.e. a−1 (1974–2000), to −0.43 ± 0.12 m w.e. a−1 (2000–15), substantially lower than the Himalayan mean (King and others, Reference King, Bhattacharya, Bhambri and Bolch2019). Much of this variability can be explained by the increasing presence of proglacial lakes. 22% of glaciers in the Eastern Himalaya are lake-terminating, increasing from 14% in 2000, accounting for 29% of total mass loss (Komori, Reference Komori2008; King and others, Reference King, Dehecq, Quincey and Carrivick2018; Zhang and others, Reference Zhang, Xu, Liu, Ding and Zhao2023), with rates of lake growth projected to continue increasing in this region (Bajracharya and others, Reference Bajracharya, Maharjan and Shrestha2014; Maurer and others, Reference Maurer, Schaefer, Rupper and Corley2019; Zhang and others, Reference Zhang, Xu, Liu, Ding and Zhao2023).
Proglacial lake development results in accelerated mass loss through the initiation of calving-driven retreat and subaqueous melt (Carrivick and Tweed, Reference Carrivick and Tweed2013; Maurer and others, Reference Maurer, Rupper and Schaefer2016; Liu and others, Reference Liu2020; Watson and others, Reference Watson2020). As the lake level increases, ice at the terminus can become buoyant, permitting calving, which can result in much larger volumes of mass loss than could occur through surface ablation alone (Kirkbride and others, Reference Kirkbride1993; Diolaiuti and others, Reference Diolaiuti2006; Benn and others, Reference Benn, Warren and Mottram2007). Lake formation increases hydrostatic and subglacial water pressure, reducing bed friction and increasing basal sliding (Van der Veen, Reference Van der Veen2002; Carrivick and Tweed, Reference Carrivick and Tweed2013; Sutherland and others, Reference Sutherland2020), resulting in flotation and a substantial reduction in basal friction, which in turn can promote acceleration and thinning (Wagner and others, Reference Wagner, Payne, Cornford and Moon2016). This impact on the glacial flow regime has resulted in half of the lake-terminating glaciers in the region showing an accelerating trend towards their termini, increasing the rate at which ice is transferred from the accumulation zone to the ablation zone. This contrasts with the broader trend of Himalayan glacier deceleration driven by thinning at lower elevations (Dehecq and others, Reference Dehecq2019; Pronk and others, Reference Pronk, Bolch, King, Wouters and Benn2021).
However, the glacier response to lake formation and the subsequent rate of lake expansion in this region are highly variable, and the mechanisms driving this variability remain poorly understood. This is significant in the Eastern Himalaya, particularly in Bhutan, where glaciers not only provide water supply for hydropower and irrigation (Ministry of Economic Affairs, 2021), but also pose a risk to downstream communities from glacial hazards in the future, particularly from Glacial Lake Outburst Floods (GLOFs) (Komori, Reference Komori2008). Bhutan alone has recorded 18 GLOF events since the 1950s, and has some of the highest risk basins globally (Watanabe and Rothacher, Reference Watanabe and Rothacher1996; Iwata and others, Reference Iwata, Narama and Karma2002; Fujita and others, Reference Fujita2012; Sakai, Reference Sakai2012; Veettil and others, Reference Veettil2016; Rinzin and others, Reference Rinzin, Zhang and Wangchuk2021; Rinzin and others, Reference Rinzin2023; Taylor and others, Reference Taylor, Robinson, Dunning and Westoby2023).
Satellite-based remote sensing of ice velocities offers the most efficient path to the large-scale monitoring of surface velocities. However, the monsoonal climate of the Bhutanese Himalaya, with cloud cover in the summer months, makes the acquisition of satellite imagery challenging. Partial or full cloud coverage in the satellite images also limits the accuracy of feature tracking, and as a result, previous studies have omitted velocity measurements during the monsoon and, therefore, measurements have been seasonally restricted (Copland and others, Reference Copland2011; Sato and others, Reference Sato, Fujita, Inoue, Sakai and Karma2022; Thongley and others, Reference Thongley2023). Furthermore, existing velocity products such as ITS_LIVE, while effective for larger glaciers and ice caps, use feature tracking window sizes that are too big to accurately capture velocities for smaller mountain glaciers and lack the temporal resolution to measure short-term changes in surface velocity (Gardner and others, Reference Gardner2018; Thongley and others, Reference Thongley2023; Scoffield and others, Reference Scoffield2024).
The seasonality of ice dynamics in the region therefore remains poorly understood, particularly in the Lunana, Bhutan, where large-scale changes to glacier area have been observed to occur at a seasonal-subseasonal frequency. Thus, there is a pressing need to understand lake–glacier interactions to predict future glacier mass loss, lake formation and associated hazards. In this study we aim to: i) quantify large-scale recent changes in ice velocity in Lunana, Bhutan at monthly scale from 2017 to 2023 using high-resolution PlanetScope imagery, ii) examine the relationship between ice velocity, and the rapid calving retreat of Thorthormi that occurred over this study period, iii) investigate the relationship between changing ice dynamics and glacier thinning and the associated drivers behind the heterogeneous response of the four lake-terminating Lunana glaciers to lake development.
2. Study area
The study area comprises the four large south-facing lake-terminating glaciers in the Lunana region of Bhutan (Fig. 1). These glaciers originate from a high elevation icecap (6900 to 6100 m asl) forming part of Chomolhari Kang (28°9’59’’ N, 90°10’57’’ E). Bechung terminates at the lowest elevation (4406 m), followed by Raphstreng (4442 m), Thorthormi (4468 m) and Lugge (4518 m). Thorthormi was the largest by area in 2021, around 5.86 km2, followed by Bechung (4.71 km2), Lugge (3.48 km2) and Raphstreng (1.85 km2). Both Bechung and Raphstreng are detached from the icecap and receive much of their accumulation by avalanching. These two glaciers also have debris-covered terminus, with Thorthormi’s terminus being around 90% debris-covered (2021), the ablation area of Lugge was around 60% clean ice in 2021. Precipitation is highest in the summer monsoon, with around 70% of annual rainfall occurring between June and September, with around 900 mm annual total (Suzuki and others, Reference Suzuki2007). At Lugge, air temperatures have been shown to range from around −14 to 9°C (Suzuki and others, Reference Suzuki2007).
Lunana region in the context of Bhutan (a). Lunana overview, including the four study glaciers: Bechung, Raphstreng, Thorthormi and Lugge (b). Red and blue polygons indicate the lower and upper areas of interest for each glacier. (Basemap: PlanetScope. 2021).
Like elsewhere in the region, rates of glacier retreat and subsequent lake expansion for Bechung, Raphstreng, Thorthormi and Lugge have been heterogeneous. Lake development in the region began when Raphstreng Glacier transitioned from land to lake-terminating in 1967, retreating at 35 m a−1, accelerating to 60 m a−1 from 1988, until the glacier became decoupled from the lake in 1998 (Ageta and others, Reference Ageta2000). In contrast, Lugge underwent continuous retreat since the late 1970s following proglacial lake development, at an average rate of 50 m a−1, with similar retreat rates at Bechung Glacier (Tsutaki and others, Reference Tsutaki2019). This is contrasted with the recent rapid retreat of Thorthormi Glacier, and the resultant expansion of a new large, potentially dangerous lake that followed a period of apparent relative stability between 1990 and 2010. After the evolution of large lateral marginal lakes into a single proglacial lake in 2011, Thorthormi began an accelerated calving retreat, with an average annual retreat rate of 85 m a−1 until 2014, increasing to 250 m a−1 by 2020 (Sato and others, Reference Sato, Fujita, Inoue, Sakai and Karma2022), the next two years were followed by the break-up of the remaining terminus. Calving at this scale has not previously been observed in the Himalaya before (Sato and others, Reference Sato, Fujita, Inoue, Sakai and Karma2022), and may be indicative of dynamic thinning, a positive feedback processes that decreases a glacier’s effective pressure altering the longitudinal stress balance, resulting in extensional flow and increased rates of calving (Benn and others, Reference Benn, Warren and Mottram2007; Carrivick and Tweed, Reference Carrivick and Tweed2013; King and others, Reference King, Dehecq, Quincey and Carrivick2018). Previous work in Lunana has revealed substantially higher velocities than the regional mean of 27.7 m a−1 (Pronk and others, Reference Pronk, Bolch, King, Wouters and Benn2021), with velocities ranging from 160 to 280 m a−1 recorded at Thorthormi within 2000 m of the Terminus from 2016 to 2017, and around 50–75 m a−1 for the same region of Lugge (Sato and others, Reference Sato, Fujita, Inoue, Sakai and Karma2022).
3. Methodology
3.1. Imagery sources
We use the PlanetScope Dove constellation because of its daily/sub-daily return frequency. This enables us to maximise the number of clear summer days and hence the availability of imagery for feature tracking. This offers increased image availability during the monsoon relative to measurements made using platforms such as Sentinel-2 or Landsat 8, which have limited coverage between May and September (Fig. 2).
Total availability of cloud-free (>80% of pixels) PlanetScope Dove imagery between 2017 and 2023 for Lunana. Compared with Sentinel-2 and Landsat 8 for the same period.
Alongside its short return frequency, PlanetScope’s high spatial resolution (3 m) compared to other commonly used platforms (Table 1) enables higher resolution velocity mapping.
Summary of satellite capabilities used for observation in High Mountain Asia, including spatial resolution, sampling frequency, wavelength range and service durations.
Orthorectified PlanetScope imagery between October 2016 and March 2024 was downloaded from Planet Explorer (planet.com/explorer/). Images with >20% identified as cloud were disregarded, leaving 526 images for the study period (Figure S1). To avoid interference from shadow motion, images were acquired at similar times: 90.1% were imaged between 03.30 and 04.50 UTC (09.30 and 10.50 BTT). Initial testing by this study found that the default multiband PlanetScope images had the lowest stable ground uncertainty compared to velocimetry runs using individual Planet bands. Multiband image pairs also had consistently lower standard deviations and low percentage error compared to single-band image pairs (Figure S6). We therefore chose to use multiband composites comprised of either 4 band PS2 (PlanetScope Dove Classic) and PS2.SD (PlanetScope Dove-R) instruments) or 8 band PSB.SD (PlanetScope Super Dove) to calculate all velocimetry data. Multiband composites from PS2 and PS2.SD instruments return RGB and NIR bands, while PS2.SD instruments return RGB, Coastal Blue, Green i, Yellow, Red Edge and NIR (Planet, 2024a).
3.2. Processing
Imagery was clipped to include the full ice area for each glacier (Fig. 3). This also included off-ice areas of stable ground proximal to the glacier to which ice velocity was normalised in order to correct for any systematic georeferencing/orthorectification variability between image pairs (Van Wyk de Vries and Wickert, Reference Van Wyk de Vries and Wickert2021). Stable ground regions were delineated following basin-wide Sentinel-2 velocimetry runs over an extended time series from 2017 to 2023, and defined as regions with detected ground motion <5 m a−1. Prior to use, all images were checked for any shifts in areas of ground where no movement is expected in the short term (e.g. ridge crests and gullies), if any movement was found, these images were removed. Any images with shading of ice areas were removed following manual examination of each time-series, as the changing position of shadows between images can be misinterpreted as ice motion (Lemmens, Reference Lemmens1988). Only 37% of images had no detectable shadow, and the remainder had a mean shadow percentage of 1.53%.
The processing chain from PlanetScope Ortho Scenes to velocity outputs.
3.3. Feature tracking
Ice velocity was determined using a technique known as Digital Image Correlation (DIC). This method measures the displacement of ice across a glacier’s surface by comparing the position of identifiable features between two images taken at different times by tracking pixel brightness within a ‘cross-correlation window’ to give a value for the relative displacement (Kääb and Vollmer, Reference Kääb and Vollmer2000; Fahnestock and others, Reference Fahnestock2016; Davison and others, Reference Davison2020). Glacier-wide values for displacement are converted into surface velocity by dividing the displacement in meters by the separation period between image pairs. This value is converted into annual velocity in m a−1 (assuming 365 days in a year) (Kääb and Vollmer, Reference Kääb and Vollmer2000; Atkinson and Becker, Reference Atkinson and Becker2020). DIC requires trackable features to be present. Cloud, snow cover and shadowing can reduce the number of trackable features such as crevasses and surface debris. Likewise, snowline retreat and shadow position change can be detected as ‘ice motion’ complicating measurements (Rosenau and others, Reference Rosenau, Scheinert and Dietrich2015). To mitigate against these issues, this study employs the feature tracking toolbox ‘Glacial Image Velocimetry’ (GIV) (Van Wyk de Vries and Wickert, Reference Van Wyk de Vries and Wickert2021). The advantage of GIV, besides being accessible and simple to use, is the ‘temporal oversampling’ function, which determines monthly velocities by tracking features between both consecutive and non-consecutive image pairs in order to produce a more accurate average velocity map for each month. This is particularly useful for Planet imagery that is sporadically spaced, as it can account for the varying time difference between image pairs (Van Wyk de Vries and Wickert, Reference Van Wyk de Vries and Wickert2021). This function also helps to average out the effect of poor or variable orthorectification between Planet images, resulting in substantial noise reduction.
The maximum interval between image pairs was varied between glaciers in order to compensate for differing flow regimes to prevent trackable features from being lost on the fastest-flowing ice (Table 2). Over time, tracked features may move outside the cross-correlation window or simply change/break apart losing coherence and becoming untrackable. We therefore undertook a series of structured tests to establish the most appropriate max-min separation between image pairs for ice at different speeds. Maximum separations >100 days revealed increasing noise for ice velocities >280 m a−1. Maximum separations were therefore capped at 90 days to reduce feature loss and minimise the effect of changing seasonal illumination of features. This was reduced to 45 days for Thorthormi, due to the relatively higher ice velocity for this glacier, and rapid deformation and breakup of features at the steep icefall, resulting in increased loss of trackable features over time. The minimum period was maintained at 7 days for all glaciers in order to reduce the effect of inaccuracies in orthorectification relative to feature displacement.
Variable GIV input parameters for the four study glaciers. Maximum and minimum image separations and the defined oversampling range for each image pair.
The temporal oversampling value defines the maximum number of image pairs that each available image is compared against; i.e. for a value of 2, image 1 is compared against image 2 as well as image 3. This value was maintained as high as computationally possible for each glacier; this ranged from 10 to 15 image pairs (Table 2). A near-anisotropic orientation filter was applied to improve feature contrast, and images were further filtered through contrast-limited histogram equalisation, high-pass filtering to reduce the effect of variable haze and lighting between image pairs (Van Wyk de Vries and Wickert, Reference Van Wyk de Vries and Wickert2021; Agarwal and others, Reference Agarwal2023).
GIV uses multiple passes to capture displacement over three gradually reducing window sizes. Our target output resolution was set to 40 m. This resulted in an initial window size of 160 m, halving with each iteration, resulting in an intermediate window size of 80 m. While initial tests were carried out using minimum window sizes (12 m), this resulted in data loss for the fastest-flowing ice. We found that a larger window was more effective in capturing the faster ice flow associated with Thorthormi Glacier, and icefall regions on Lugge, with 40 m offering an optimal balance between capturing ice flow and maintaining a good image output resolution. This resolution is still significantly finer than other existing products, thus maximising the benefit of the high-resolution Planet imagery. Velocity outputs were then filtered via signal-to-noise ratio <3 and a peak ratio <1.2, following Agarwal and others (Reference Agarwal2023). GIV excludes values greater than 1.5 standard deviations away from the monthly mean, ensuring that only values below this threshold are maintained (Van Wyk de Vries and Wickert, Reference Van Wyk de Vries and Wickert2021). As GIV calculates monthly velocity averages, this is the frequency of velocity outputs used in this study.
3.4. Regions of interest
To compare ice velocity and surface elevation change over time, Bechung, Thorthormi and Lugge were divided into two regions of interest to separate the slower-moving low gradient terminus ice, found below the lowest icefall of each glacier, and the faster-flowing steeper ‘upper’ ice (Fig. 1) for analysis purposes. This was because ice velocity in these regions differs by an order of magnitude, and the ice velocity response of ice to seasonal drivers also varies. Due to its smaller area and consistent gradient, Raphstreng is defined in one region (Table 3). The lower region of Thorthormi was redefined for every month so that the lower boundary matched the location of the glacier terminus and excluded lake ice/icebergs to avoid spurious velocities. Ice velocities were extracted as the average value across each region of interest (ROI) for every month.
Surface area for the regions of interest across the four study glaciers, including the average slope for each region of interest. Thorthormi’s lower region of interest is described for 2021, as the region area is frequently varied throughout the study, given the rapid changes in the spatial extent of the terminus over the study period.
3.5. Estimating surface elevation change
We acquired two 2 m resolution digital elevation models (DEMs) from the Pléiades Glacier Observatory—for October 2017 and November 2022 (Berthier and others, Reference Berthier2024). These DEMs were generated using the Semi-Global Matching approach, following the DEM pipeline developed by Berthier and others (Reference Berthier2024), which employs the Ames Stereo Pipeline (Beyer and others, Reference Beyer, Alexandrov and McMichael2018). Both DEMs were coregistered to the Copernicus GLO-30 DEM (European Space Agency, 2024). Surface elevation change (dH) for the four Lunana glaciers was estimated using z-based pixel differencing. Berthier and others (Reference Berthier2024) measured an average uncertainty of 0.492 m for glacial area across the entire data set. We estimated an uncertainty of ±1.026 m by calculating the mean elevation difference across our stable ground area. The rates of thinning in the DEMs were compared with available directly measured data:
Measurements for surface elevation change between 2017 and 2022 for the terminus region of Lugge are comparable to DGPS measurements by Tsutaki and others (Reference Tsutaki2019), although higher than Sato and others (Reference Sato, Fujita, Inoue, Sakai and Karma2022) (Table 4).
Rate of surface elevation change compared with previous studies for Thorthormi and Lugge. Data up to the present study compiled by Sato and others (Reference Sato, Fujita, Inoue, Sakai and Karma2022).
3.6. Uncertainty in surface flow velocity
Following error quantification approaches for ice velocities in previous work (Millan and others, Reference Millan2019; Zhang and others, Reference Zhang, Xu, Liu, Ding and Zhao2023), we used the mean velocity across the stable-ground regions as an indirect reference for the relative uncertainty for each monthly velocity output. These regions are assumed to have close to zero displacement, relative to glacier ice, and provide a good reference for the relative uncertainty across the ROI being interrogated on each GIV run (Zhang and others, Reference Zhang, Xu, Liu, Ding and Zhao2023). PlanetScope monthly ice velocities were then clipped to these areas and averaged for each month to give a value for stable ground.
In the absence of ground-based velocity data, we compared the GIV velocity results for PlanetScope against the 2017–2018 global ice velocity dataset from Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) (this data is available only for Lugge Glacier). While the Millan dataset may not be as accurate as other products, such as ITS_LIVE (Gardner and others, Reference Gardner2018), the limited coverage of smaller glaciers in High Mountain Asia made ITS_LIVE unsuitable for comparison in Lunana. To facilitate direct comparison, the PlanetScope annual velocity means were resampled to 50 m using bicubic interpolation, aligning with the resolution of the Millan product (Figure S2).
In addition, a secondary Sentinel-2 velocity dataset was generated using GIV. Sentinel-2 imagery was acquired via Planet Explorer (planet.com/explorer/). Band 4 was used as this band was typically less saturated than other bands over snowy ground (Kääb and others, Reference Kääb2016). The values for oversampling and maximum/minimum image separations were maintained as the same for each glacier as for the Planet imagery (Table 2). The process of filtering images for quality was also the same as for Planet imagery. Data from Thorthormi’s upper region had a low signal-to-noise ratio, as snow and cloud cover hampered feature tracking, preventing comparison between Sentinel-2 and PlanetScope for this region. We compared velocities mean velocities across the upper and lower ROIs for Bechung and Lugge, the ROI for Raphstreng and the lower region of Thorthormi (a lack of cloud-free imagery and extensive snow cover resulted in the velocity data for the upper region of Thorthormi being unusable) (Figure S3).
Following a similar approach to Sato and others (Reference Sato, Fujita, Inoue, Sakai and Karma2022), we also conducted manual validation by measuring the horizontal displacement of 65 clearly identifiable crevasses and seracs between Sentinel-2 image pairs for 2019 (Figure S5). Measured displacement between image pairs was divided by the days of image pair separation time and multiplied by 365. Manual measurements were compared directly with PlanetScope GIV outputs for the same period for Bechung, Thorthormi and Lugge. Raphstreng was excluded due to the rapid loss of visible features between image pairs. Sentinel-2 was used, as its quality is well established and less prone to orthorectification issues, making pairwise comparison more reliable, while the resolution is high enough for large features to be visible.
4. Results
4.1. Spatial pattern of ice flow
There was a strong spatial dichotomy in flow regime between Thorthormi Glacier, Bechung, Raphstreng and Lugge (Fig. 4). At the latter three glaciers, velocities decelerated towards the terminus, with mean velocities in the terminus regions of Bechung and Lugge ROI (Fig. 1) ranging from 10 to 70 ± 8.5 to 21.9 m a−1. In contrast, the low gradient terminus at Thorthormi exhibited downstream ice acceleration, with ice velocities ranging from 200 to 400 ± 10.0 m a−1, an order of magnitude higher than its neighbouring glaciers. The steep upper regions of all glaciers were fast flowing, with Thorthormi and Lugge showing the fastest centreline velocities (200 to 420 ± 10.0 to 8.5 m a−1), whereas Bechung and Raphstreng were slightly slower (180 to 300 ± 21.9 to 10.4 m a−1).
Spatial distribution of mean ice velocities for the four study glaciers with average flow directions (2017 to 2023). Basemap; PlanetScope 2021.
4.2. Temporal evolution of ice velocity
The seasonal pattern of velocity variability was relatively consistent between the study glaciers, but the magnitude of this change varied substantially (Fig. 5). Ice velocities for the terminus regions of the four glaciers typically peaked around spring/early-summer, while upper regions showed acceleration into early autumn, with all regions showing a yearly minimum in the winter. Lugge had an average seasonal ice velocity range of 35.5 ± 8.5 m a−1 (63.2% of mean) for the lower region (Fig. 5), and a much lower 18.8 ± 8.5 m a−1 (8.3% of mean) for the upper region. In contrast, the lower region of Thorthormi had an average seasonal range of 144.6 ± 10.0 m a−1 (44.7% of mean). This contrasted with the upper region of Thorthormi, which again had a lower seasonal variability (55 ± 10.0 m a−1, 25.7% of mean), although substantially higher than Lugge’s. Raphstreng, despite its small size, showed a large mean seasonal range (65.4 ± 10.4 m a−1, 39.7% of mean), with ice velocities typically peaking around June–July, and reaching a seasonal low around March. The upper region of Bechung showed the highest seasonal variability, with a mean seasonal range of 66.5 ± 21.9 m a−1 (62% of mean). The velocity in the lower region of Bechung was also variable between months; however, the relative change in velocity was comparable with the range in stable ground motion each month, so no clear pattern can be determined.
Mean monthly ice velocities for the study glaciers across the regions of interest (red lower, blue upper): Bechung (a), Raphstreng (b), Thorthormi (c), Lugge (d). Stable ground uncertainty is indicated in grey along the bottom.
4.2.1. Temporal and spatial variability in seasonal ice velocity patterns from 2017 to 2023
Seasonal acceleration at Bechung was generally restricted to the upper region of the glacier, where ice acceleration occurred from spring (March–May) until autumn (September–November), followed by a rapid slowdown in winter (December–February) (Fig. 6a–c). The pattern of change for the lower region was not uniform and varied spatially across the terminus between all seasons. From 2017 to 2023, a diverging velocity trend emerged between the lower and upper regions. The lower region showed a mean increase in ice velocity of 7.5 ± 21.9 m a−1 (31.1% of mean), compared to a decrease of −9 ± 21.9 m a−1 (8.4% of mean) in mean ice velocity from 2017 to 2023 across the upper region.
Left: changes in average seasonal velocity from 2017 to 2023 for Bechung (a), Raphstreng (d), Thorthormi (f), Lugge (i). From the figure left to centre follows changes from winter to spring, spring to summer, summer to autumn, autumn to winter. Right: average monthly ice velocities 2017 to 2023 for the lower and upper regions of Bechung (b,c), Raphstreng (e), Thorthormi (g,h), Lugge (j,k).
Raphstreng showed a consistent annual pattern of summer acceleration and autumnal/winter deceleration (Fig. 6d and e). The lowest region of the glacier was the first to show a speed-up from winter to spring, by summer, this acceleration was glacier-wide, followed by a period of autumnal deceleration lasting into the winter. In the middle of the ROI, there was a small region where acceleration occurred from summer to autumn and autumn to winter, with deceleration in the following two seasons. This area corresponded with a point where the glacier became partially detached, and frequent serac avalanches occurred over a steep bedrock band. While there was a slight negative trend in ice velocity over the study period (−5.7 ± 10.4 m a−1, 3.5% of mean), there was little change in mean ice velocity from 2017 to 2023.
The lower region of Thorthormi exhibited continuous ice acceleration from winter into summer, with higher upstream ice velocities in early summer transferring ice downstream to the terminus (Fig. 6f–h). The seasonal acceleration in upstream ice velocity occurred later than at the terminus, with widespread ice acceleration occurring over the spring-summer period. This was followed by decelerating ice velocities in the highest regions (>6085 m a.s.l.) during autumn, although acceleration continued into the autumn 1.5 km downstream from the steep, heavily crevassed icefalls (5200 m a.s.l.) down to the lower headwall (4630 m a.s.l.).
From 2017 to 2023, there was a divergent trend in ice velocities between the upper and lower regions of Thorthormi (Figs. 5 and 6g,h). The upper region showed a decelerating trend, with a −53.7 ± 10.0 m a−1 (−23.7%) change in mean ice velocity from 2017 to 2023. In contrast, the terminus accelerated from 2017 to 2021, with maximum velocities increasing from 307 m a−1, peaking at 448.3 ± 10.0 m a−1 in June 2021. From 2017 to 2021, mean annual ice velocity in this region increased by 180.4 ± 10.0 m a−1 (74.8%); however, by 2023, mean ice velocities in this region decreased by 104.5 ± 10.0 m a−1 from the 2021 peak (−24.8%).
Lugge experienced a low magnitude seasonal acceleration relative to the other three glaciers. This peaked around late spring for all but the upper icefall (between 3.5 and 3.7 km from the terminus), with continued acceleration into the early autumn for the upper regions, but deceleration for the middle and lower parts of the glacier (Fig. 6j and k). Ice velocities increased slightly again in the autumn-winter transition, particularly in the middle region of the glacier. The inter-annual velocity trend for Lugge revealed little overall change from 2017 to 2023; however, as with Bechung, there was a slight decelerating trend for the upper region ± 8.5 m a−1, 5% of mean), in contrast to an accelerating trend for the lower terminus region (3.3 ± 8.5 m a−1, 5.9% of mean).
4.3. Surface elevation change from 2017 to 2022
Surface elevation change was greatest for the lower region of Lugge with a mean change of −24.1 ± 1.03 m between 2017 and 2022 (−4.82 ± 0.17 m a−1) (Fig. 7). The greatest change in surface elevation was seen at the transition between the lower region and the icefall of Thorthormi which showed a surface elevation change of more than −60 m (−12 ± 0.17 m a−1). This contrasted with a mean change of −11.6 ± 1.03 m (−2.32 ± 0.17 m a−1) in the lower region. The lower region of Bechung also showed a substantial reduction in surface elevation, with a mean change of −14.7 ± 1.03 m Raphstreng at a slightly higher elevation also showed a relatively consistent reduction in surface elevation with a mean difference of −9.4 ± 1.03 m. Surface elevation change was lower for the upper regions of all glaciers, with a −2.3 ± 1.03 m change for Lugge, and a −9.3 ± 1.03 m change for the upper region of Thorthormi. The upper region of Thorthormi, however, showed a sizeable variation in elevation in change (−68 to +30 ± 1.03 m). The upper region of Bechung presents a change in mean surface elevation of −6.6 ± 1.03 m.
Change in surface elevation (±1.03 m) (dH) between 2017 and 2022 for the four study glaciers.
5. Discussion
5.1. Heterogeneous glacier response to lake formation in Lunana
Our data showed that ice velocities for the four study glaciers are considerably higher than the Himalayan lake-terminating average of 18.8 m a−1, as well as the average for the Eastern Himalaya of 27.7 m a−1 (Pronk and others, Reference Pronk, Bolch, King, Wouters and Benn2021). Relative to previous work, we show that ice velocities for Thorthormi continued to increase, exceeding previous peak velocity measurements of 280 ± 10.0 m a−1, while Lugge showed little change (Sato and others, Reference Sato, Fujita, Inoue, Sakai and Karma2022). We also reveal a contrasting pattern of velocity change over the study period between Thorthormi and the other three Lunana glaciers, Bechung, Raphstreng and Lugge. Specifically, we identified increasing variability in the seasonal ice velocities of the lower region of Thorthormi, coupled with an overall acceleration trend at low elevations between 2017 and 2022. The other three glaciers did not show an accelerating trend during the study period, nor do we find a consistent pattern of seasonal variability from year to year. In the following, we examine the differing dynamic behaviour of Thorthormi Glacier, despite similar climate forcing, given the spatial proximity to the four Lunana glaciers.
The transition to much faster ice flow, exceeding 450 ± 10.0 m a−1 by 2021, followed a series of major break-up events of the terminus. Between 2017 and 2018, the ice in the lower region detached from the valley sides and began calving of large tabular icebergs (Fig. 8). This may have de-buttressed the terminus, leading to a reduction in lateral and frontal shear stresses due to a loss in lateral confinement contributing to acceleration (Fig. 8c and d, Figure S7) (Raymond, Reference Raymond1996; Schoof and others, Reference Schoof, Davis and Popa2017; Hill and others, Reference Hill, Gudmundsson, Carr, Stokes and King2020). In June 2019, a major breakup event coincided with a large landslide and GLOF event (National Centre for Hydrology and Meteorology and Namgay W, 2019; National Centre for Hydrology and Meteorology, 2020). This was followed by the re-advance of the lake calving terminus, with the winter lake freeze reducing ice motion in the lake, allowing ice to accumulate at the base of the icefall, allowing the lake calving terminus to extend to a size of 0.53 km2 (Fig. 8e). This feature was highly fractured (Fig. 8f), and high ice velocities suggest that it may have been lightly grounded or floating, and was likely to be comprised of floating sections of ice that reformed into a semi-cohesive feature at the base of the glacier/lakes headwall. We find no comparable examples of such a rapid lake calving terminus formation in the existing literature. The snout continued to extend down lake until the detachment of a large tabular iceberg in November 2021 (Fig. 8g). This was followed by little change over the following winter until 09/04/2022, when a large calving event led to the detachment of the remaining terminus ice (Fig. 8i); the white ice floating in the lake resulting from fracturing indicated the scale of this event. This detachment was followed by the continued breakup of the terminus ice, and evidence of calving of seracs from the glacier headwall was also observed (Fig. 8j).
Stages of the break-up of Thorthormi’s terminus from 2017 to 2023 (a–j), compared with ice mean velocity across the terminus for the same period (k). The width of each scene corresponds to a scale of 2 km.
The widespread flotation of a glacier terminus is rare in lacustrine environments because ice fronts typically become unstable upon reaching flotation (Warren and others, Reference Warren, Benn, Winchester and Harrison2001; Carrivick and others, Reference Carrivick, Tweed, Evans, Roberts and Russell2020). Consequently, observations of large, transitory floating termini in lake-terminating settings are uncommon. Where they have been documented, we find a number of cases with large tabular style calving, similar to that of Thorthormi. For instance, at the Mendenhall Glacier in Alaska, a floating terminus persisted for nearly two years before rapidly breaking into four large tabular icebergs during a period of accelerated ice flow (Boyce and others, Reference Boyce, Motyka and Truffer2007). Similarly, from 2012 to 2019, the San Quintín Glacier in Patagonia experienced a comparable retreat pattern. Bathymetry data at San Quintín Glacier revealed that this period of rapid calving coincided with the glacier retreating into an overdeepening, accompanied by a shift in surface velocities from deceleration to acceleration towards the terminus (Pętlicki and others, Reference Pętlicki, Rivera, Oberreuter, Uribe and Johannes2023). This retreat was also associated with the formation of longitudinal fractures and the calving of elongated tabular icebergs (Podgórski and Pętlicki, Reference Podgórski and Pętlicki2020). Rapid lacustrine calving has also been documented at Upsala Glacier, Patagonia, from a floating terminus, where high thinning rates drove periodic retreat (Naruse and Skvarca, Reference Naruse and Skvarca2000; Boyce and others, Reference Boyce, Motyka and Truffer2007; Moyer and others, Reference Moyer, Moore and Koppes2016, Minowa and others, Reference Minowa, Schaefer and Skvarca2023). Additionally, Holdsworth (Reference Holdsworth1973) observed a lake-terminating outlet of the Barnes Ice Cap calving buoyant ice blocks fracturing within 50 m of the terminus. The increasing frequency of large calving events at Thorthormi during the study period indicates that the terminus may have become progressively more buoyant, likely a result of progressive thinning of the terminus ice. This thinning would have also reduced the tensile and bending stresses of the terminus, resulting in higher strain rates, which would have led to higher ice velocities and acceleration (Sugiyama and others, Reference Sugiyama2011; Tsutaki and others, Reference Tsutaki2019).
The increasingly high magnitude, high frequency seasonal changes in ice velocity we see for the lower region of Thorthormi between 2020 and 2023 indicate that the ice became increasingly sensitive to changing driving stresses as the terminus approached buoyancy. Ice that is close to its critical flotation thickness is more sensitive to thinning than well-grounded ice. Thinning would drive increased buoyancy, lowering the overburden pressure and decreasing longitudinal stresses, resulting in acceleration (Warren and others, Reference Warren, Benn, Winchester and Harrison2001; Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017; Carrivick and others, Reference Carrivick, Tweed, Evans, Roberts and Russell2020; Kellerer-Pirklbauer and others, Reference Kellerer-Pirklbauer2021; Pronk and others, Reference Pronk, Bolch, King, Wouters and Benn2021). Importantly, we do not see an accelerating trend at the lower regions of Bechung and Lugge, or at Raphstreng in response to lake growth. This suggests that the termini of these glaciers are grounded and not close to buoyancy (Boyce and others, Reference Boyce, Motyka and Truffer2007; Wagner and others, Reference Wagner, Payne, Cornford and Moon2016). This suggests that the heterogeneous lake expansion may be primarily a response to substantially differing stress regimes at the terminus of Thorthormi and the other three glaciers.
5.2. Seasonal controls on ice dynamics in Lunana
While all four glaciers reveal a clear annual cycle in glacier speed, the timing of peak seasonal velocities for Bechung, Raphstreng and Lugge is also not coincident, and the seasonal velocity pattern for individual glaciers also varies between the years (Fig. 5). For example, ice velocities in the upper region of Bechung peak between September and November (Figs. 5a and 6a), where at Lugge this occurred between April and July (Figs. 5d and 6i). Meanwhile, Raphstreng’s ice velocities peaked in mid-summer (Figs. 5b and 6d) (June to August), and the upper region of Thorthormi undergoes two periods of annual acceleration, one in late spring/early summer (April to June), and another in autumn (September to November) (Figs. 5c and 6f). These variable trends in seasonal velocity change have been observed in the Himalaya generally (Zhang and others, Reference Zhang, Xu, Liu, Ding and Zhao2023), and in the Nyainqêntanglha mountains in Tibet, with some glaciers showing velocities peaking in the winter, and others in spring and summer (Zhang and others, Reference Zhang2023). We explore the relationship between glacier velocity in Lunana and its likely seasonal drivers, within the constraints of a lack of directly measured meteorological and hydrological data for Lunana.
5.2.1. Glacier velocity relationship with hydrological drivers
For all four glaciers, we observe a notable acceleration following the winter season (Figs. 5 and 6a,d,f,i). The initiation and magnitude of this acceleration are likely influenced by differences in each glacier’s response to changing basal water pressure, the structure of its basal till and the rate of subglacial drainage system closure due to creep (Hewitt, Reference Hewitt2013; Vincent and Moreau, Reference Vincent and Moreau2016; Armstrong and others, Reference Armstrong, Anderson and Fahnestock2017). At both Raphstreng Glacier and the upper and lower regions of Thorthormi, we observe a strong acceleration signal during the early melt season, between early spring and early summer (Figs. 5b,c and 6d,f). This is likely driven by increased surface melt overwhelming subglacial drainage systems that are adapted to winter conditions, such that as the input of meltwater increases, basal water pressures rise close to, or above the effective pressure, promoting enhanced basal sliding (Paterson, Reference Paterson1994; Sole and others, Reference Sole2011).
At Raphstreng and in the upper region of Thorthormi, the initial peak in velocity during the spring/early summer is followed by a gradual deceleration (Figs. 5b,c and 6d,f). This may be due to the evolution of a more efficient, channelised subglacial drainage system in response to the increased water flow. As this system becomes more developed, it reduces basal water pressure and thus slows down glacier movement. Additionally, Raphstreng, Thorthormi and Bechung all exhibit significant deceleration moving into the winter (Figs. 5a,b,c and 6a,d,f), likely driven by the reduction in surface meltwater input and the viscous creep closure of subglacial water channels (Anderson and others, Reference Anderson2004; Hewitt, Reference Hewitt2013).
5.2.2. Seasonal changes in driving stress
Seasonal acceleration may be driven in part by summer snowfall resulting in mass gain, and increased driving stresses (Jouberton and others, Reference Jouberton2022; Zhang and others, Reference Zhang, Xu, Liu, Ding and Zhao2023). The varying magnitude of the glacier response to a seasonal forcing is also a product of the relationship between the differing morphology of each glacier, such as varying overburden pressure (Zhang and others, Reference Zhang, Xu, Liu, Ding and Zhao2023). For instance, Lugge Glacier exhibits low seasonal velocity variability (Fig. 5), likely due to its location within a deep valley, resulting in higher resistive forces from the valley sides, with lateral resistive forces being shown to fluctuate less seasonally compared to basal resistive stress (Tripathi and others, Reference Tripathi, Singh, Rathore, Oza and Bahuguna2023).
In contrast, Bechung, Raphstreng and Thorthormi lack similar buttressing along the glacier margins, receiving less resistance from the valley sides (Fig. 5). As a result, the proportion of net drag from the margins is lower, making these glaciers more sensitive to changes in basal drag (Röthlisberger, Reference Röthlisberger1972; Tripathi and others, Reference Tripathi, Singh, Rathore, Oza and Bahuguna2023). This difference in resistive forces could explain the higher sensitivity of some of the Lunana glaciers to early-season warming, while others that show a pattern of seasonal acceleration later in the melt season may be more sensitive to peaks in monsoonal precipitation during July and September (Suzuki and others, Reference Suzuki2007; Naito and others, Reference Naito2012).
Monsoonal precipitation also increases accumulation from avalanche activity, which, in turn, raises driving stress and ice velocity (Bhushan and others, Reference Bhushan, Shean, Hu, Guillet and Rounce2024; Kneib and others, Reference Kneib2024). This could explain the late summer to early autumn peak in acceleration at the upper region of Bechung Glacier, which is primarily fed by avalanching (Figs. 5a and 6a), as mass increases following high rates of monsoonal avalanching. Similarly, mass input from avalanches off the overhanging ice cap at Raphstreng Glacier could account for years in which high ice velocities persist into late autumn. The small region of acceleration identified on Raphstreng that occurred from summer into early winter (Fig. 6d) is an area associated with frequent serac collapses and avalanching; the detected increase in ice velocity most likely corresponds to an increase in avalanching during late summer to early winter.
At Thorthormi, our data suggests that the winter freezing of surface lake ice could play an important role in controlling seasonal ice velocities across the terminus. Calving rates were much higher outside of the periods when the lake was frozen. In lake-terminating settings, the freezing of the lake and ice mélange can impart sufficient back stresses to have a stabilising effect on the calving front and lower region, contributing to a slowdown in ice loss as well as surface velocity (Geirsdóttir and others, Reference Geirsdóttir, Miller, Wattrus, Björnsson and Thors2008; Carrivick and others, Reference Carrivick, Tweed, Evans, Roberts and Russell2020). The widespread presence of mélange at Thorthormi may have also been partially responsible for its higher magnitude seasonal cycles in ice velocity relative to the other Lunana glaciers, as mélange has been shown to be an important factor in marine terminating settings in strengthening frozen sea ice, resulting in much stronger back stresses (Amundson and others, Reference Amundson, Fahnestock and Truffer2010; Greene and others, Reference Greene, Young, Gwyther, Galton-Fenzi and Blankenship2018), the same may be true in a lake environment. Measurements of surface lake ice throughout the winter, as well as modelling of lake ice buttressing potential, in a similar style to studies conducted at marine terminating glaciers (Greene and others, Reference Greene, Young, Gwyther, Galton-Fenzi and Blankenship2018; Parsons and others, Reference Parsons, Sun and Gudmundsson2024), would help to better constrain the importance of seasonal lake ice in controlling ice velocity in a high mountain setting.
We have compared several reanalysis datasets and found none to be representative of the few climate measurements available for Lunana. The nearest meteorological station with freely available data is in Gyatsa, around 79 km away from the Thorthormi glacier; this is therefore not utilised as it is also unlikely to be representative. Up-to-date meteorological data would be necessary to further constrain the environmental controls on ice dynamics in Lunana.
5.2.3. Velocity uncertainty
Ice velocities for Lugge, compared to the dataset produced by Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022), were on average 11% (32.11 ± 8.5 m a−1) faster for PlanetScope; however, in the slower-flowing terminus region, average ice velocities were 0.4% (−8.8 ± 8.5 m a−1) slower for PlanetScope (Figure S2). For our GIV-generated Sentinel-2 dataset, we found PlanetScope ice velocities to be 11% faster across all four glaciers. This ranged from 1.2% for the lower region of Bechung to 21.2% for Thorthormi (Figure S2). This compares to our manual feature tracking of Sentinel-2 image pairs, where between Bechung, Thorthormi and Lugge, we found PlanetScope velocities to be on average 2.89 m a−1 faster than measured velocities (r = 0.957) (Figure S5).
Mean detected stable ground motion from 2017 to 2023 was lowest for Lugge (8.5 m a−1), followed by Thorthormi (10.0 m a−1) and Raphstreng (10.4 m a−1). The stable ground velocities at Bechung Glacier were notably higher (21.9 m a−1). These ranges are comparable with other similar studies (Millan and others, Reference Millan2019; Tsutaki and others, Reference Tsutaki2019; Watson and others, Reference Watson2020; Tripathi and others, Reference Tripathi, Singh, Rathore, Oza and Bahuguna2023; Zhang and others, Reference Zhang, Xu, Liu, Ding and Zhao2023). Stable ground velocities for our Sentinel-2 dataset were higher than those for PlanetScope for Raphstreng and Lugge (12.3 m a−1–8.3 m a−1, respectively), and lower for Bechung and Thorthormi (8.2 m a−1–8.1 m a−1, respectively) (Figure S4). The difference in stable ground velocities for Bechung for our PlanetScope data (11.5 m a−1 higher) compared to Sentinel-2 is indicative of the variable orthorectification quality with the PlanetScope product, and while images were filtered prior to registration, variable levels of distortion between image scenes would have resulted in higher values for stable ground motion. It may be that a lack of ground control points proximal to Bechung contributed to the ground surrounding this glacier having higher levels of distortion (Planet Labs, 2019).
Ice velocities in this study consistently exceeded those derived from Sentinel-2; this is likely a result of PlanetScope’s shorter repeat imaging frequency, compared to Sentinel-2, allowing for improved tracking of the fastest-moving features. These features are typically located at rapidly deforming icefalls and therefore remain coherent and trackable for shorter periods. Higher frequency imaging increases the probability of tracking these features before they move beyond the search window between repeat intervals. Using higher frequency imaging, therefore, results in a higher value for mean velocity.
Environmental factors are also important. Changing snow cover between image pairs can produce false matches and therefore erroneous velocity values, while slight changes to illumination conditions can generate differently shaped shadows having a similar effect, partially over variably featured stable ground (Mueting and Bookhagen, Reference Mueting and Bookhagen2024). The rough, highly vegetated and boulder-strewn ground surrounding the study glaciers may have therefore contributed to the higher-than-average stable ground velocities for both the Sentinel-2 and PlanetScope datasets.
5.3. Surface elevation change
While the pattern of surface elevation change is relatively similar between all four glaciers in the basin, particularly in the upper regions, there are some notable disparities. The upper region of Thorthormi showed a considerable range in elevation change (−68 to + 30 ± 1.03 m) (Fig. 7). Much of this variability may be a result of changing positions of the large seracs present on the icefall over the sampling period. The large surface elevation change identified at the terminus headwall on Thorthormi (Fig. 7) may be indicative of a shift from compressive to extension flow in the terminus region, in part driving this thinning and movement of ice away from the headwall (Sato and others, Reference Sato, Fujita, Inoue, Sakai and Karma2022). We do not identify substantially higher thinning rates in the upper region of Thorthormi compared to Bechung, Raphstreng and Lugge, indicating the retreat of the terminus at Thorthormi has not resulted in notable extensional flow and accelerated thinning relative to neighbouring glaciers.
As observed in other regions of High Mountain Asia and Bhutan, slower ice velocities have been attributed to reduced driving stresses resulting from thinning (Dehecq and others, Reference Dehecq2019; Thongley and others, Reference Thongley2023). However, in Lunana, despite thinning across all four glaciers between 2017 and 2022, we did not find a net deceleration of glacier flow (Figs. 5–7). Although the lower regions of all glaciers show the most substantial thinning, we do not identify a large slowdown in ice velocity in these regions (Figs. 6c,h,k and 7). It is possible that ice velocities in the lower regions were sustained by increasing basal water flux and longer summer melt periods, driven by warmer mean annual air temperatures, offsetting the effect of reducing driving stresses (Vincent and Moreau, Reference Vincent and Moreau2016; Goosse and others, Reference Goosse2018; Dehecq and others, Reference Dehecq2019). If parts of Thorthormi’s terminus were close to a critical flotation point, seasonal increases in thinning through surface melt across the terminus region may have been sufficient to result in portions of the terminus becoming ungrounded, with the associated reduction in basal drag contributing to the observed seasonal acceleration. Such a relationship between seasonal thinning and ice acceleration has been observed at Kangerlussuaq and Helheim glaciers in Greenland (Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017).
Regardless, despite the high ice velocities, mean rates of thinning across Thorthormi’s terminus were not exceptionally high relative to the lower regions of its neighbouring glaciers (Fig. 7). The presence of very cold lake water, indirectly indicated by persistent icebergs and ice growlers, may have sustained the terminus by reducing rates of subaqueous melt and therefore thinning (Warren and others, Reference Warren, Benn, Winchester and Harrison2001). At the Yakutat Glacier in Alaska, flotation was sustained over longer periods (months to years) when the water temperature was close to freezing (Trüssel and others, Reference Trüssel, Motyka, Truffer and Larsen2013). The importance of a debris mantle on the terminus will have also lowered thinning rates by insulating the ice, lowering rates of surface melt (Rowan and others, Reference Rowan, Egholm, Quincey and Glasser2015), whereas Lugge had higher rates of thinning and at the terminus with a lower proportion of debris cover. Likewise, Bechung and Raphstreng both have significant debris coverage across their terminus and lower thinning rates, despite being at a lower elevation than Lugge. The debris mantle may have delayed buoyant calving at Thorthormi as sediment cover increases the total mass of the ice column, resulting in a lower freeboard height and grounding at a greater water depth (Watson and others, Reference Watson2020). Previous studies indicate that Thorthormi’s lake may be unusually deep for the region (National Centre for Hydrology and Meteorology and Namgay, 2019; Zhang and others, Reference Zhang, Xu, Liu, Ding and Zhao2023). It is possible that this depth accommodated the development of an exceptionally thick terminus. This thickness may have been sufficient to support the tensile strength required to resist buoyant forces and longitudinal strain associated with buoyancy (Boyce and others, Reference Boyce, Motyka and Truffer2007). Constraints on lake bathymetry would provide a better understanding of the changing ice thickness over the study period and, therefore, terminus conditions, with bed topography being necessary to quantify dynamic thinning and any potential retreat into an overdeepening.
The contrasting dynamic behaviour of Thorthormi Glacier compared to its neighbouring glaciers under similar environmental conditions highlights the importance of understanding the localised controls on lake-terminating glacier dynamics. Unlike Bechung, Raphstreng and Lugge glaciers, which show more stable or decelerating trends (Fig. 5), Thorthormi’s likely proximity to buoyancy made the terminus highly sensitive to short-term changes in driving stress, resulting in high magnitude changes in ice velocity. This highlights the role of subaqueous processes, such as thinning from melt and calving, which altered the glacier’s stress regime, making it more sensitive to dynamic changes. While lake freezing in winter may have provided temporary stabilisation, the progressive thinning and floating of Thorthormi’s terminus appear to have exacerbated its seasonal variability. In contrast, the grounded termini of the other Lunana glaciers exhibit less sensitivity to similar environmental changes. With lake-terminating glaciers becoming more numerous throughout the Himalaya (King and others, Reference King, Bhattacharya, Bhambri and Bolch2019), broader work should be undertaken to measure the changing pattern of ice velocity and surface elevation of different glaciers to lake formation in different regions of the Himalaya. This is crucial for predicting for assessing future risk from glacial hazards such as GLOFs, understanding water supply, as well as the wider scale implications of proglacial lake development on glacier dynamics throughout High Mountain Asia.
6. Conclusions
This study provides a comprehensive analysis of surface ice velocity and dynamic changes across four lake-terminating glaciers in the Lunana region of Bhutan between 2017 and 2023. Our findings highlight the contrasting behaviours of Thorthormi compared to Bechung, Raphstreng and Lugge, offering new insights into the role of proglacial lakes in glacier dynamics and retreat in the Himalaya.
We identified a substantial acceleration of Thorthormi Glacier, with annual peak velocities increasing from 307 ± 10.0 m a−1 in 2017 to 448 ± 10.0 m a−1 by 2021, driven by the breakup of its terminus and widespread flotation. We also observed an increasing magnitude in seasonal surface velocities variation for the lower region of Thorthormi over the study period, coinciding with thinning and calving of large tabular icebergs, suggesting widespread flotation of the terminus and retreat driven by dynamic thinning. In contrast, the other three glaciers did not exhibit the same accelerating trend, though Bechung and Raphstreng showed considerable seasonal variability in ice velocities, particularly in their upper regions, while Lugge exhibited much lower variability without such strong seasonality.
The study revealed that glacier thinning in the upper regions of all four glaciers corresponded with a slight deceleration, yet the lower regions of Bechung and Lugge experienced a minor increase in velocity despite high thinning rates. An increased meltwater flux may be compensating for the loss of overburden pressure caused by thinning. However, the relationship between meltwater and ice velocity is also complex and strongly influenced by the geometry and composition of the glacier bed and a glacier’s stress regime. Our findings underscore that while proglacial lake expansion has a critical impact on glacier retreat, glacier behaviour can vary substantially, and is highly dependent on specific terminus conditions.
This study is an important step in understanding how glacial outburst flood risk can increase over a short period of time, demonstrating the importance of routine surveying of lake-terminating glaciers that may be at risk of transitioning to a phase of rapid retreat. High-resolution mapping of ice velocity can provide a useful indicator of a changing flow regime that could be indicative of acceleration preceding rapid retreat. Given the increasing prevalence of the lake-terminating glaciers in the region, further work on the interactions between glacial dynamics, lake formation and GLOF susceptibility is crucial for predicting the broader impacts of climate change on Himalayan glaciers and understanding the potential changes to outburst flood risk over the coming century.
Supplementary material
The Supplementary Material for this paper can be found at https://doi.org/10.1017/jog.2025.10103.
Acknowledgements
This research has been supported by the National Environmental Research Council (NERC) funded ONE Planet Doctoral Training Partnership (NE/S007512/1), hosted jointly by Newcastle and Northumbria universities. The authors would like to thank Etienne Berthier for kindly producing the two DEMs from our Pleiades imagery, used to calculate surface elevation change. We would also like to thank the reviewers and the Scientific Editor for their helpful comments and contributions.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Open access
