Articles | Volume 16, issue 6
The Cryosphere, 16, 2643–2654, 2022
The Cryosphere, 16, 2643–2654, 2022
Research article
01 Jul 2022
Research article | 01 Jul 2022

Land- to lake-terminating transition triggers dynamic thinning of a Bhutanese glacier

Land- to lake-terminating transition triggers dynamic thinning of a Bhutanese glacier
Yota Sato1, Koji Fujita1, Hiroshi Inoue2, Akiko Sakai1, and Karma3 Yota Sato et al.
  • 1Graduate School of Environmental Studies, Nagoya University, Nagoya, Japan
  • 2Multi-hazard Risk Assessment Research Division, National Research Institute for Earth Science and Disaster Resilience (NIED), Tsukuba, Japan
  • 3Cryosphere Service Division, National Centre for Hydrology and Meteorology (NCHM), Thimphu, Bhutan

Correspondence: Yota Sato (


There have been rapid increases in both the number and expansion of the proglacial lakes across High Mountain Asia. However, the relationship between proglacial lakes and glacier dynamics remains unclear in the Himalayan region. Here we present the surface elevation, flow-velocity changes, and proglacial lake expansion of Thorthormi and Lugge glaciers in the Lunana region, Bhutanese Himalaya, during the 2000–2018 period using photogrammetry and GPS survey data. The lake expansion and surface lowering rates and flow-velocity field of Lugge Glacier, a lake-terminating glacier, have remained approximately constant since 2000. Conversely, there have been accelerated proglacial lake expansion and a 2-fold increase in the thinning rate of Thorthormi Glacier since 2011, as well as a considerable speed-up in the flow-velocity field (>150 m a−1). We reveal that the lake formation and transition of Thorthormi Glacier from a land- to lake-terminating glacier have triggered glacier speed-up and rapid thinning via a positive (compressive) to negative (extensional) change in the emergence velocities. This study provides the first evidence of dynamic glacier changes that are associated with proglacial lake formation across the Himalayan region.

1 Introduction

A recent deglacial trend has been reported for numerous glaciers across High Mountain Asia (HMA; e.g. Brun et al., 2019; Maurer et al., 2019; Shean et al., 2020), with these glaciers exhibiting spatially heterogeneous thinning patterns (Bolch et al., 2012; Kääb et al., 2012; Brun et al., 2017). There has been a rapid increase in both the number and expansion of the proglacial lakes across HMA owing to this deglacial trend (Zhang et al., 2015; Nie et al., 2017; Shugar et al., 2020), which has been particularly pronounced across the eastern Himalayas (e.g. Gardelle et al., 2011; Chen et al., 2021). Proglacial lakes form via the coalescence of supraglacial lakes near the glacier terminus (Quincey et al., 2007); their formation suggests the final phase of retreat for these contracting glaciers (Sakai and Fujita, 2010; Benn et al., 2012). The increasing number and evolution of proglacial lakes have led to a rise in the hazardous potential of glacial lake outburst floods (GLOFs). GLOF hazards can be triggered by either unstable terminal moraines or snow/rock avalanches (e.g. Fujita et al., 2008; Westoby et al., 2014), and they can cause significant damage to hydropower stations, bridges, and buildings that exist downstream of proglacial lakes (Richardson and Reynolds, 2000).

Proglacial lake formation accelerates glacier mass loss via thermal undercutting and calving at the glacier terminus (e.g. Benn et al., 2007; Sakai et al., 2009). Previous studies have analysed the interaction between proglacial lakes and glacier dynamics using in situ measurements and remote-sensing methods across HMA (e.g. King et al., 2018; Haritashya et al., 2018; Wei et al., 2021). Recent high-resolution satellite and aerial photogrammetry techniques have led to improved glacier and proglacial lake studies. For example, Watson et al. (2020) acquired in situ measurements and unmanned aerial vehicle (UAV) photogrammetry across Thulagi Glacier in the Nepal Himalaya, and they estimated the calving volume at the terminus based on iceberg size. Furthermore, previous studies have also reported the retreat and thinning of lake-terminating glaciers in their catchments to a broad regional scale (e.g. Song et al., 2017; Zhang et al., 2019; Maurer et al., 2019). King et al. (2019) reported that the mass loss of lake-terminating glaciers was greater than that of land-terminating glaciers across broad Himalayan regions, with an observed increase in mass loss after 2000. Pronk et al. (2021) analysed the surface flow velocities of more than 300 glaciers in the Himalayan region and determined that the flow velocities of lake-terminating glaciers were twice as high as those of land-terminating glaciers on average. The existence of a proglacial lake might be a factor enhancing the glacier flow velocity, retreat, and thinning of HMA glaciers. The response of lake- and land-terminating glaciers can fluctuate with different patterns even if they are located near each other and/or exist under similar climatic conditions (Liu et al., 2020). Therefore, advancing our understanding of lake-terminating glacier fluctuations is essential for making robust future predictions of the HMA glacier response.

Numerous proglacial lakes have an exceptionally high-risk potential for GLOFs throughout the Bhutan Himalaya (e.g. Fujita et al., 2013; Zheng et al., 2021), and lake expansion appears to continue unabated (Ageta et al., 2000; Komori, 2008). Previous lake-terminating glacier studies have been conducted across the Bhutan Himalaya using either in situ measurements or satellite remote-sensing methods to assess their dynamics and evolutions (e.g. Suzuki et al., 2007; Fujita et al., 2008). Tsutaki et al. (2019) revealed contrasting fluctuations between two neighbouring glaciers in the Lunana region using in situ GPS measurements, satellite remote-sensing data, and numerical modelling. They reported a greater thinning rate along lake-terminating Lugge Glacier than along land-terminating Thorthormi Glacier during the 2004–2011 period, which was attributed to their contrasting terminus conditions. They also projected that the thinning rate and flow speed of Thorthormi Glacier could be accelerated if the current land terminus changed to a lake terminus. The terminus of Thorthormi Glacier is now detached from the terminal moraine and has evolved into a lake-terminating glacier. The associated changes in glacier dynamics due to proglacial lake formation have been studied worldwide (e.g. Boyce et al., 2007; Tsutaki et al., 2013); however, no such study has been undertaken in the Himalayan region to date. Therefore, this study aims to (1) update the fluctuations of two glaciers that have been affected by proglacial lakes in the Bhutan Himalaya and (2) evaluate the changes in glacier dynamics associated with the transition from land- to lake-terminating conditions. We analysed past in situ measurements and satellite and airborne remote-sensing datasets to achieve this goal.

2 Study site

Our target glaciers, Thorthormi and Lugge, are located in the Lunana region of northern Bhutan (Fig. 1; 28.06 N, 90.18 E). Thorthormi Glacier covers 11.6 km2 based on the GAMDAM (Glacier Area Mapping for Discharge from the Asian Mountains) glacier inventory (Nuimura et al., 2015; Sakai, 2019) and the 2017 terminus position. Its elevation range spans 4400–6900 m above sea level (a.s.l.). Thorthormi Glacier had been in contact with the terminal moraine before 2011 and then detached from the terminal moraine and transitioned into a lake-terminating glacier (Tsutaki et al., 2019). Lugge Glacier is located to the east of Thorthormi Glacier and covers an area of 10.0 km2 based on its 2017 terminus position, with its elevation range spanning 4500–6900 m a.s.l. Lugge Glacier Lake has expanded since the 1960s (Komori, 2008), with a maximum lake depth of 126 m reported in 2002 (Yamada et al., 2004). This lake caused an outburst flood in October 1994 and damaged the downstream areas (Fujita et al., 2008; Maurer et al., 2020). Both glaciers are debris-covered and have been reported to be experiencing long-term mass-loss and thinning trends (Bajracharya et al., 2014; Maurer et al., 2016).

Figure 1Details of the study site. (a) Location of the Lunana region (inset) and helicopter photogrammetry (HP) orthoimage of Thorthormi and Lugge glaciers (acquired on 24 March 2018). (b) Surface elevation map generated from the HP-DEM using ground control points (GCPs) for terrain data processing (open circles) and 2011 GPS tracks (dots). (c, d) Aerial photographs of Thorthormi and Lugge glaciers. Red arrows in panel (a) indicate the directions from which the aerial photographs were taken. The dashed box in (b) shows the domain of (a). The green GPS track in (b) was not used for the DEM accuracy check or elevation change analysis.


3 Observations and analysis methods

3.1 DGPS and aerial photogrammetry survey

We used the global positioning system (GPS) dataset of Tsutaki et al. (2019), who conducted a kinematic survey with a differential GPS (DGPS; GEM-1, GNSS Technologies) across the on- and off-glacier terrains during the 19–22 September 2011 field campaign. The base station for this survey was installed to the west of Thorthormi Glacial Lake (Fig. 1a). These GPS data points were used to validate the satellite and photogrammetry digital elevation models (DEMs) and compute the surface elevation changes (Sect. 3.2 and 3.3).

We conducted a helicopter photogrammetry survey on 24 March 2018. Four action cameras (GoPro HERO5 Black) were attached to the skids of a helicopter and acquired 4000×3000 pixel images in 1 s shooting mode. The shutter speed, focal length, and ISO were fixed at 1/1250 s, 28.3 mm, and 100, respectively. We obtained ∼7500 photos in total and cropped each photograph by preserving the central 2500×2500 pixel area to eliminate the “fisheye effect” of the GoPro camera lens (Girod et al., 2017). We finally employed 3560 images based on the image quality and glacier coverage. These images were processed in Agisoft Metashape Professional Edition 1.7.1 (Agisoft LLC), and the sky view was masked for the terrain data processing.

3.2 Terrain data processing

We extracted ground control points (GCPs) from a Pléiades panchromatic orthoimage (0.5 m resolution), which was acquired on 7 November 2017 (Berthier et al., 2014), and its DEM (2.0 m resolution) for the photogrammetry terrain data processing. We first generated a GPS-derived DEM (GPS-DEM) to assess the vertical accuracy of the Pléiades-derived DEM (PL-DEM). The 2011 GPS data points (UTM Zone 46N, WGS84) were interpolated using the inverse distance weighted method and then exported to the same grid size as the PL-DEM in ArcGIS (Tshering and Fujita, 2016; Sato et al., 2021). We employed the standard deviation (SD; σ) of the elevation difference between the PL- and GPS-DEMs on the off-glacier stable terrain as the vertical accuracy of the PL-DEM. We did not use the grid cells with steep slopes (>30; Fujita et al., 2008; Nuimura et al., 2012). We then eliminated the validation points that were greater than ±3σ from the mean elevation difference as extreme outliers (Mertes et al., 2017). Berthier et al. (2014) reported that the vertical accuracy of the PL-DEM was improved by shifting the DEM horizontally. We therefore shifted the PL-DEM by ±2 pixels (±4 m) in the northing and easting directions, computed the elevation difference against the GPS-DEM, and confirmed that there was no improvement in the vertical accuracy. Finally, the PL-DEM vertical bias (mean elevation difference: MED) was assessed for 12 009 grid cells, yielding a mean value of 0.26±3.86 m (MED ± SD). We extracted the GCP coordinates from the orthoimage and bias-corrected PL-DEM. Specific topographic features (e.g. boulders, river bending points, and dense vegetation spots) on the stable ground were used as GCPs for the photogrammetry terrain data processing.

We used the structure from motion (SfM) software in Agisoft Metashape to generate orthoimages and a DEM from the helicopter photogrammetry data (hereafter HP-ortho and HP-DEM, respectively). We overlaid the 77 GCPs that were extracted from the Pléiades products onto the helicopter photogrammetry images (Fig. 1b) and generated both the HP-ortho and HP-DEM at a 0.5 m resolution (Fig. 1a and b). We employed the same approach used in the PL-DEM evaluation to evaluate the vertical bias and accuracy of the HP-DEM by re-generating a 0.5 m resolution GPS-DEM. The vertical accuracy of the HP-DEM was 0.25±3.70 m (N=25 474 GPS-DEM grid cells; Fig. S1a in the Supplement); we also applied an elevation change correction (Sect. 3.3) to correct for the vertical bias of the HP-DEM.

3.3 Changes in the glaciers and glacial lakes

We calculated the surface elevation change rates (dh/dt) by comparing the 2011 GPS-DEM and 2018 HP-DEM (both at 0.5 m resolution). We used 9491 and 15 604 grid cells to calculate dh/dt for Thorthormi and Lugge glaciers, respectively (red tracks in Fig. 1b). We then compared our results with previous elevation change studies. Tsutaki et al. (2019) calculated dh/dt from the overlapping 2004 and 2011 DGPS data; they also computed the spatial distribution of dh/dt for the same period using Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER)-derived DEMs. We also employed the long-term dh/dt data from Brun et al. (2017) and Maurer et al. (2019) to assess the thinning trends of these two glaciers. Brun et al. (2017) computed dh/dt for the 2000–2016 period over the Hindu Kush Himalaya region using ASTER-based DEMs, and Maurer et al. (2019) calculated dh/dt for the 1975–2000 and 2000–2015 periods using satellite-based DEMs. These datasets (hereafter RS-based dh/dt) are provided as 30 m resolution raster data. We extracted the dh/dt data from our DEMs at the same positions for a comparative analysis.

We computed the surface velocity field using the ImGRAFT (Image GeoRectification and Feature Tracking) open-source feature tracking toolbox in MATLAB (Messerli and Grinsted, 2015). The normalized cross-correlation algorithm (NCC; Heid and Kääb, 2012) in the feature tracking toolbox (Templatematch) identifies the displacement patterns of the glacier surface features and computes their magnitude from a pair of images. We selected a Sentinel-2 image pair that was acquired on 16 November 2016 and 11 November 2017 (post-monsoon seasons). After visual trial and error, we chose a 20×20 template size (200×200 m) and 75×75 search window size (750×750 m) to compute the surface feature displacements and calculate the annual flow velocity. We set the correlation value for image matching (r>0.6) and signal to noise ratio (SNR <0.7) and eliminated the low-quality pixels, all of which served as a confidence level threshold for successful image matching. We then estimated the uncertainty of the glacier surface velocity and corrected systematic error by checking the stable-ground (off-glacier) displacement (e.g. Liu et al., 2020). We calculated the stable-ground surface displacement (slopes <20; Quincey et al., 2009) and set the corrected median values of Vx (east–west component) and Vy (north–south component) to zero. The flow speed V (m a−1) is calculated as follows:

(1) V = V x 2 + V y 2 .

The mean and median stable-ground V values were 2.2 and 1.6 m a−1, respectively, after the displacement correction. The velocity profiles were extracted along the glacier central flowlines every 10 m, and the pixel values where the flow directions differed by >90 from the flowlines were eliminated. We also eliminated the velocity data along the upper section of Lugge Glacier (>5100 m elevation) because of its heavy snow cover, which can cause incorrect image matching via feature tracking (e.g. Nuimura et al., 2017). We extracted the surface velocity from the regional velocity product derived from the ITS_LIVE (Inter-Mission Time Series of Land Ice Velocity and Elevation) project (Gardner et al., 2019), which covered the entire HMA region. The ITS_LIVE velocity product is generated from the Landsat series with the auto-RIFT feature tracking processing chain yielding a 240 m spatial resolution (Gardner et al., 2018). We extracted the annual velocity data along the central glacier flowlines from the ITS_LIVE product (2010–2018) to trace the temporal changes in flow velocity; the ITS_LIVE velocities possess mean uncertainties of 1.0 m a−1 (maximum of 3.1 m a−1) and 0.6 m a−1 (maximum of 6.0 m a−1) along Thorthormi and Lugge glaciers, respectively. We also employed the velocity data produced by Tsutaki et al. (2019), which were calculated annually from the ASTER-derived optical satellite images at 15 m resolution during the 2002–2011 period.

We delineated the glacial lake area from Landsat 7 and 8 (ETM+: Enhanced Thematic Mapper Plus; OLI: Operational Land Imager) images with false-colour image composites that were acquired between November 2012 and November 2018 (30 m resolution). We then combined the proglacial lake polygons before 2012 (Tsutaki et al., 2019) and traced the annual lake area changes for the entire 18-year study period. The total lake area uncertainties were estimated to be ±0.14 and ±0.08 km2 for Thorthormi and Lugge glaciers, respectively, depending on the user-induced error and satellite image resolution (Paul et al., 2013). We removed the DEM and velocity data where the glacier surface turned into the lake surface in successive images. We compared the recent lake expansion rates with a glacial lake inventory (High Mountain Asia Glacier-lake inventory: Hi-MAG; Chen et al., 2021), which was generated for the entire Himalayan region using data from the 2008–2017 period. We finally chose the 2011 and 2017 proglacial-type lakes (n=832) and calculated the expansion rates between 2011 and 2017 in the eastern Himalaya region (including the Lunana region).

3.4 Emergence velocity and ice flotation of Thorthormi Glacier

We calculated the emergence velocities along Thorthormi Glacier to evaluate the change in glacier dynamics since its detachment from the terminal moraine. We estimated the emergence velocity of a given section (Ve, m a−1) from the ice fluxes along the upper and lower boundaries of the section as follows (e.g. Nuimura et al., 2011; Vincent et al., 2016; Brun et al., 2018):

(2) V e = q in - q out W d x ,

where qin and qout are the ice fluxes (m3 a−1) along the upper and lower boundaries, respectively, and W and dx are an averaged glacier width (m) and length (200 m in this study) for the analysed Ve section, respectively. The ice flux q (qin or qout, m3 a−1) is calculated as follows:

(3) q = W h V c ,

where W, h, and Vc are the glacier width (m), ice thickness (m), and depth-averaged ice velocity (m a−1), respectively. We then applied a simplified assumption that the glacier width is constant such that Eq. (2) can be rewritten as follows:

(4) V e = h up V c , up - h low V c , low d x ,

where hup/low are the ice thicknesses (m) and Vc,up/low the depth-averaged ice velocities (m a−1) along the upper/lower boundaries. We assumed that both the glacier thickness and width were constant in the transverse and longitudinal directions, respectively, to calculate the emergence velocities along the central flowline. The ice thickness h along the central flowline was calculated from the HP-DEM-derived glacier surface elevation and estimated bedrock elevation in Tsutaki et al. (2019). Tsutaki et al. (2019) estimated the glacier-bed topography using the Farinotti et al. (2009) ice thickness model and tuning a model parameter based on the observed lake depth (ice thickness). Tsutaki et al. (2019) simulated that the basal velocity reaches 90 % of the surface velocity of Thorthormi Glacier. We therefore calculated the emergence velocity using two assumptions regarding the surface velocity: the depth-averaged velocity is 90 % of the surface velocity based on Tsutaki et al. (2019), and 100 % is assumed for a floating condition after terminus detachment. We then used the surface velocity component of the same vector in the central flowline direction. The sections without flow velocities (2520–3020 m from the 2002 terminus) were linearly interpolated using the surface velocities of the surrounding upglacier and downglacier sections. We calculated the emergence velocity for a 200 m section by shifting the section in 50 m increments and obtained a mean emergence velocity around the current terminus (2400–3500 m from the 2002 terminus). We also calculated Ve in 2011, when Thorthormi Glacier was still a land-terminating glacier, to compare the land- and lake-terminating conditions. We estimated the ice thickness from the glacier surface elevation of the ASTER-derived DEM acquired on 6 April 2011 and the glacier-bed elevation along the central flowline. The depth-averaged ice velocities were calculated from the surface velocities in Tsutaki et al. (2019) (Sect. 3.3).

We evaluated the ice flotation potential of the Thorthormi Glacier terminus based on the ice flotation thickness (hf, m), which was calculated as follows (Boyce et al., 2007; Watson et al., 2020):

(5) h f = ρ w ρ i d ,

where ρw is the density of water (1000 kg m3), ρi is the density of ice (917 kg m3) (e.g. Boyce et al., 2007; Carrivick et al., 2017), and d is lake depth (m). We then defined an index of potential ice flotation as follows:

(6) P f = h f h × 100 ,

where Pf is the potential ice flotation (%). The glacier can attain flotation when the glacier ice thickness reaches hf such that Pf is ≥100 %. We extracted the lake surface elevation (4415 m a.s.l.) from the HP-DEM-derived lake perimeter and estimated the lake depth from the glacier-bed elevation (Tsutaki et al., 2019). We then calculated Pf in 100 m intervals in 2011 and 2018 along the glacier central flowline in the terminus region (up to 3500 m from the 2002 terminus).

Figure 2Temporal variations in the spatial extents of (a) Thorthormi and (b) Lugge proglacial lakes. (c) Cumulative changes in lake area for Thorthormi and Lugge glaciers relative to 2000. The background images in (a) and (b) are Sentinel-2 satellite images that were acquired on 11 November 2017. The 2000–2011 lake outlines are from Tsutaki et al. (2019). The dA/dt values in (c) are the 2000–2011 (upper left) and 2011–2018 (lower right) lake expansion rates.

4 Results

4.1 Lake expansion

We traced the lake expansion for the 2000–2018 period (Fig. 2a and b). The proglacial lake areas at the termini of Thorthormi and Lugge glaciers were 3.05 and 1.58 km2 in 2018, an increase of 2.01 (193 %) and 0.48 km2 (44 %) from the 2000 lake areas, respectively (Fig. 2c). Both lakes have expanded throughout the study period, and the lake expansion rates (dA/dt) were calculated via a linear regression of the cumulative areas during the 2000–2011 and 2011–2018 periods (Fig. 2c). Lugge Glacial Lake steadily expanded during the 2000–2018 period, with 0.03 and 0.02 km2 a−1 observed before and after 2011, respectively. However, there has been accelerated expansion of Thorthormi Glacial Lake since 2011, with 0.07 km2 a−1 observed before 2011 and 0.13 km2 a−1 observed after 2011. A comparison of these observations with the Hi-MAG data (Chen et al., 2021) indicates that the expansion rates are in the upper 2.5 % (Thorthormi) and 10 % (Lugge) of the observed proglacial lakes across the eastern Himalayas.

Table 1Comparison of the measured surface elevation change rates of Thorthormi and Lugge glaciers from various studies.

Download Print Version | Download XLSX

4.2 Thinning rates

The dh/dt values of both glaciers were calculated from the 2011 GPS-DEM and 2018 HP-DEM, with the 2002 terminus position used as the base position for the comparison. We also extracted the calculated dh/dt values from previous studies that had focused on different time periods (Fig. 3 and Table 1). The thinning rate of Lugge Glacier was more than 3 times greater than that of Thorthormi Glacier for the 2004–2011 period, when Thorthormi Glacier was a land-terminating glacier, and was then comparable (−2.78 m a−1) to that of Lugge Glacier (−2.87 m a−1) for the recent 2011–2018 period, when Thorthormi Glacier had evolved into a lake-terminating glacier (Fig. S1). There was a 2-fold increase and 0.61-fold decrease in the thinning rates of Thorthormi and Lugge glaciers between the 2004–2011 analysis by Tsutaki et al. (2019) and our presented 2011–2018 analysis. The RS-based dh/dt values for the 2000–2016 period (−3.81 to −3.50 m a−1; Brun et al., 2017; Maurer et al., 2019) are similar to the Tsutaki et al. (2019) values (Table 1). The RS-based dh/dt values of Thorthormi and Lugge glaciers are −0.16 and −1.20 m a−1 for the 1975–2000 period, respectively, which suggests that the lower section of Thorthormi Glacier experiences minimal thinning before 2000 (Maurer et al., 2019). The spatial distribution of dh/dt along Thorthormi Glacier exhibited a decreasing trend in the upglacier direction during the 2011–2018 period, whereas the dh/dt values during the 2004–2011 period were almost constant across the same region. The thinning rate was >4 m a−1 (dh/dt of less than −4 m a−1) in the upglacier area (>2500 m from the 2002 terminus) during this later period. The dh/dt profiles obtained in previous studies do not reveal such a remarkable trend; however, similar dh/dt plots are independent of the distance from the terminus (Fig. 3a). The results of this study reveal a large spatial variability compared with the RS-based dh/dt distributions of previous studies (Brun et al., 2017; Maurer et al., 2019), which is likely due to the differences in the spatial resolution of the data. Tsutaki et al. (2019) reported that Lugge Glacier has a heavily crevassed, bumpy surface. We consider that the high-resolution photogrammetry data (0.5 m) identified the displacements due to these steep surface features.

Figure 3Surface elevation change rates (dh/dt) along (a) Thorthormi and (b) Lugge glaciers (based on the distance from the 2002 glacier terminus). Each panel shows the elevation change rates for 1975–2000 and 2000–2016 (Maurer et al., 2019; JM19), 2000–2016 (Brun et al., 2017; FB17), 2004–2011 (Tsutaki et al., 2019; ST19), and 2011–2018 (this study). HP-DEM and GPS-DEM are resampled to a 30 m resolution for comparison with the dh/dt datasets from previous studies.


4.3 Flow velocity

We calculated the surface velocity field between November 2016 and November 2017 and extracted the velocities along the central flowline (Figs. 4 and S2). We also plotted the ITS_LIVE product (2010–2018) and Tsutaki et al. (2019) velocity profiles. We found fast flow velocities (>200 m a−1) from the 2017 terminus to the middle part (∼4200 m from 2002 terminus) of Thorthormi Glacier. The ITS_LIVE velocity also exhibited a similar flow-velocity magnitude near the 2017 terminus; however, the ITS_LIVE velocities decreased more rapidly in the upglacier direction than our calculated results. We are able to confirm the large displacement of a surface feature (∼200 m a−1) from the Sentinel-2 satellite images at ∼3500–4000 m from the 2002 terminus (Fig. S3), which suggests that the flow-velocity profile for Thorthormi Glacier that is calculated in this study should be more reliable than the automatically derived ITS_LIVE flow-velocity profile. A comparison of our velocity profile (2016–2017) with the 2002–2010 velocity profile (Tsutaki et al., 2019) reveals a substantial 2–4-fold increase at ∼2400–4000 m from the 2002 terminus. The ITS_LIVE flow-velocity profiles indicate that until 2016 the Thorthormi Glacier shows similar flow velocity profiles to the 2002–2010 average, and remarkable increases are observed in 2017 and 2018 (Fig. 4a). Tsutaki et al. (2019) projected the flow velocities of Thorthormi Glacier under the assumption of a “lake-terminating condition” (dashed line in Fig. 4a). Although the magnitude is ∼100 m a−1 less than that in this study, the increasing flow velocities toward the calving front are similar to the trend in the recent velocity profile.

Figure 4Central flowline velocities of (a) Thorthormi and (b) Lugge glaciers. Dashed vertical lines indicate the glacier terminus positions in 2002, 2011, and 2017. The red lines represent the flow velocity that was calculated in this study (2016–2017). Interannual velocities (2010–2018) are extracted from the ITS_LIVE velocity product. The black lines and grey shaded regions represent the mean and standard deviation of the flow velocities for the 2002–2010 period, respectively (Tsutaki et al., 2019). The thick dashed line in (a) denotes the simulated surface velocities with a lake-terminating assumption (Tsutaki et al., 2019).


Conversely, the calculated surface velocities of Lugge Glacier are ∼50 m a−1 up to ∼2000 m from the 2002 terminus in this study (2016–2017). There also appears to be a gradual decrease up to ∼2700 m from the 2002 terminus, and it possesses a similar velocity magnitude/trend to the 2002–2010 mean velocity calculated by Tsutaki et al. (2019; Fig. 4b). The ITS_LIVE velocity profiles (2010–2018) show flow velocities of <5 m a−1 for the entire glacier, which is probably due to the coarser resolution (240 m) of the velocity field compared with that in this study (10 m) and Tsutaki et al. (2019; 15 m). Although the terminus position of Lugge Glacier has retreated almost 1 km since 2002, the mean velocity profile appears to have remained persistent between the 2000–2010 (Tsutaki et al., 2019) and 2016–2017 observation periods.

Table 2Comparison of the emergence velocity of Thorthormi Glacier in 2011 and 2017. The mean values are calculated for the 2400–3500 m section from the 2002 terminus (Fig. S4). Two basal sliding conditions are assumed, whereby depth-averaged velocity equals either 90 % or 100 % of the surface velocity.

Download Print Version | Download XLSX

4.4 Ice emergence velocity and flotation potential of Thorthormi Glacier

We calculated the emergence velocity of the Thorthormi Glacier terminus under the assumption that the depth-averaged velocity equals either 90 % or 100 % of the surface velocity. The resultant Ve values are -0.69±11.65 and -0.77±12.94 m a−1, respectively (2400–3500 m from the 2002 terminus; Fig. S4 and Table 2), although there are large variations depending on the computational area (Fig. S4). The land-terminating condition yielded Ve values of 5.20±3.78 and 5.78±4.20 m a−1, respectively, for the above-mentioned depth-averaged velocity assumption (Table 2). These results suggest that the mean Ve has decreased and become negative after transitioning to a lake-terminating glacier. Tsutaki et al. (2019) also estimated Ve via numerical modelling of the lake- and land-terminating conditions, yielding 2.2±1.9 and 6.3±2.2 m a−1, respectively.

We also estimated the potential ice flotation index in the terminus area of Thorthormi Glacier (up to 3500 m from the 2002 terminus). The mean Pf values for 2011 (land-terminating) and 2018 (lake-terminating) are 86 % and 97 %, respectively, with this increase attributed to the surface lowering of the terminus area during the 2011–2018 period. As a result of surface lowering, some parts of the ice in the terminus area reached ice flotation thickness (Pf>100 %) by 2018.

5 Discussion

5.1 Contrasting temporal changes in the glacier regimes

Thorthormi and Lugge glaciers possess contrasting dh/dt trends and flow velocities even though they are adjacent to each other. The dh/dt trend and the flow-velocity magnitude and spatial distribution of Lugge Glacier are approximately constant between the 2004–2011 and 2011–2018 periods and 2002–2010 and 2016–2017 periods, respectively (Figs. 3b and 4b). Conversely, remarkable increases in the thinning rate and flow velocity are observed across Thorthormi Glacier over the same study periods (Figs. 3a and 4a). Such a drastic velocity increase within a decade has not been reported in the Himalayas, although the multi-decadal acceleration of glacier thinning and deceleration of glacier flow have been reported (Dehecq et al., 2019; Maurer et al., 2019).

Tsutaki et al. (2019) performed finite-element simulations of present (land-terminating) and future (lake-terminating) Thorthormi Glacier dynamics. Their simulations reproduced the flow velocities for the land-terminating condition with a small root-mean-square error (<10 m a−1) using satellite-based flow velocities. Their future prediction for a lake-terminating condition, which suggested an increase in flow velocity, is inconsistent with our 2017 velocity analysis (Fig. 4a). However, Tsutaki et al. (2019) highlighted that changes to the sliding coefficient and ice thickness parameters could alter the flow velocity significantly as their sensitivity tests demonstrated that the simulated flow velocity increased (decreased) by 33 % (51 %) if the sliding coefficient and ice thickness changed by +30 % (−30 %) for the land-terminating condition of Thorthormi Glacier. Therefore, the difference between the observed and simulated velocities is likely due to the uncertainties in the sliding coefficient, ice thickness, and state of the terminus position. Despite this underestimation, Tsutaki et al. (2019) have reasonably demonstrated the change in sliding conditions associated with the transition from land- to lake-terminating conditions.

5.2 Dynamic thinning triggered by terminus detachment

Lateral lakes had formed on both sides of the Thorthormi Glacier terminus several years before 2011 (Fig. 2a and c). The ice thickness was near flotation (Pf>85 %) such that the flow velocities could accelerate near the terminus. However, the flow velocities decreased toward the terminus, producing a longitudinal compression field and subsequent surface lowering that might have been less than that of lake-terminating Lugge Glacier (Tsutaki et al., 2019). The longitudinal stress field regime changed from compressional to extensional after the terminus detached from the terminal moraine and transitioned to a lake-terminating condition, and its flow increased owing to efficient basal sliding. Although the satellite imagery shows that the glacier terminus detached from the terminal moraine in 2011 (Fig. 2a), the rapid increase in flow velocity since 2017 suggests that the glacier terminus was in contact with the terminal moraine underwater until 2016 and detached between 2016 and 2017 (Fig. 4a). Furthermore, the lakes that formed on both sides of the glacier terminus may also have reduced the lateral resistive stresses that prevented glacier flow (e.g. Adhikari and Marshall, 2012). These factors might have led to the observed dramatic increase in flow velocities (Fig. 4a). Such an increase in flow velocities due to proglacial lake formation has been observed in other regions (e.g. Boyce et al., 2007; Tsutaki et al., 2011; Sakakibara and Sugiyama, 2014); however, this is the first observation of such a phenomenon in the Himalayan region.

The rapid increase in flow velocities may have enhanced the ice flux towards the glacier terminus due to the longitudinal strain. Positive emergence velocities are distributed up to 2400–3500 m from the 2002 terminus for the 2011 land-terminating condition (Sect. 4.4). However, Ve decreased in 2017 and became negative due to the increase in flow velocities toward the terminus. The unchanged thinning rate and velocity regime of Lugge Glacier (Figs. 2b and 3b and Table 1) suggest that any recent climatic changes in the Lunana region could not have yielded a significant increase in surface ablation. Therefore, the increased thinning rate of Thorthormi Glacier can be largely attributed to the decrease in Ve. However, the decrease in Ve (approximately −6 m a−1; Table 2) seems to be too large to account for the increased thinning rate from 2004–2011 to 2011–2018 (−1.38 m a−1; Table 1). We quantified the change in emergence velocity by hypothesizing that this change occurred in the last 2 years (2017 and 2018), during which time the surface velocity accelerated (Fig. 4a). The weighted average of the emergence velocity for the 2011–2018 period (Ve,2011–2018, m a−1) is described as follows:

(7) V e , 2011 2018 = Δ t land V e , land + Δ t lake V e , lake Δ t land + Δ t lake ,

where Ve and Δt are the emergence velocity and duration of the land- or lake-terminating conditions, respectively. We obtained a time-weighted mean emergence velocity of 3.52 m a−1 for the 2011–2018 period based on emergence velocities of Ve,land=5.20 and Ve,lake=-0.69 m a−1 (assuming a depth-averaged velocity that is 90 % of the surface velocity; Table 2) and periods of Δtland=5 and Δtlake=2. This means that Ve decreased by −1.68 m a−1 around 2011, which is consistent with the 2004–2011 to 2011–2018 change in dh/dt of −1.38 m a−1. The highly variable Ve profiles suggest that there are large uncertainties in the estimates (Fig. S4 and Table 2); however, our first-order evaluation can explain the cause of the drastic change in the thinning rate of Thorthormi Glacier.

We simply assumed a constant glacier width to calculate Ve along the central flowline. However, this glacier terrain tends to widen in the downglacier direction, yielding an extensional velocity regime. The lateral proglacial lakes on both sides of the terminus before it transitioned to a lake-terminating condition may have further contributed to a more negative Ve than that estimated along the central flowline. Despite these favourable conditions to enhance dynamic thinning, the surface lowering of Thorthormi Glacier has likely been suppressed by the compressive flow regime of the land-terminating condition. The transition to a lake-terminating condition should have caused a 2-fold increase in the thinning rate during such a short period (Fig. 2a and Table 1).

These above-mentioned mechanisms might cause a positive feedback between glacier thinning and the increase in flow velocity by enhancing each other. Therefore, increased glacier thinning and surface velocity speed-up will continue along Thorthormi Glacier in the future. The dynamic thinning of lake-terminating glaciers has been discussed in other HMA regions (e.g. Nuimura et al., 2012; King et al., 2018; Liu et al., 2020). However, our study is the first reported observation of the dynamic changes during the transition from land- to lake-terminating conditions, which have led to the enhanced thinning of a Himalayan glacier.

This study employed a modelled ice thickness (lake depth) that was tuned using point measurement data (Tsutaki et al., 2019) to estimate the dynamics of Thorthormi Glacier. Previous studies have suggested that the surface flow velocity of lake-terminating glaciers is sensitive to the terminus ice thickness and lake water depth (Benn et al., 2007; Pronk et al., 2021). Therefore, constraints on the lake bathymetry may allow us to better understand past and current terminus conditions and quantify the dynamic thinning process.

6 Conclusions

We presented the surface elevation and velocity changes and proglacial lake expansion of lake-terminating Thorthormi and Lugge glaciers in the Lunana region, Bhutanese Himalaya. We analysed satellite and photogrammetry data and compared our results with those in previous studies to reveal the recent glacier and proglacial lake changes of Thorthormi Glacier, which are associated with the transition from land- to lake-terminating conditions. Whilst the lake expansion and surface lowering rates of Luge Glacier have been approximately constant since 2000, those of Thorthormi Glacier have exhibited a continued increase after the terminus reached flotation and detached from the terminal moraine. There has been a 2-fold increase in the thinning rate of Thorthormi Glacier since this transition to lake-terminating conditions. The flow-velocity field of Thorthormi Glacier has also sped up considerably (>150 m a−1), whereas that of Lugge Glacier has remained unchanged. We estimate that the rapid thinning and increased flow-velocity field of Thorthormi Glacier were due to this transition to lake-terminating conditions. This study provides the first evidence of the dynamic glacier changes associated with proglacial lake formation in the Himalayan region and will contribute to advancing our understanding of the dynamics of lake-terminating glaciers, as well as their potential evolution in the future.

Data availability

The Landsat 7 ETM+, Landsat 8 OLI, and Sentinel-2 satellite data are distributed by the United States Geological Survey (, USGS, 2021). The ASTER-DEM data are distributed by the National Institute of Advanced Industrial Science and Technology (, AIST, 2021).


The supplement related to this article is available online at:

Author contributions

KF designed the study. HI and K conducted the photogrammetry survey. YS processed the photogrammetry data and analysed the satellite data. YS, KF, and AS wrote the manuscript. All of the authors contributed to the discussion.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank Jamyan Chopel and Shiro Ohmi for supporting the aerial photogrammetry survey. We are indebted to Shun Tsutaki and Takayuki Nuimura for providing their data and supporting our data analysis. We thank Etienne Berthier for providing the Pléiades satellite data. We appreciate Homa Kheyrollah Pour and two anonymous referees for their insightful and constructive comments.

Review statement

This paper was edited by Homa Kheyrollah Pour and reviewed by two anonymous referees.


Adhikari, S. and Marshall, S. J.: Parameterization of lateral drag in flowline models of glacier dynamics, J. Glaciol., 58, 1119–1132,, 2012. 

Ageta, Y., Iwata, S., Yabuki, H., Naito, N., Sakai, A., Narama, C., and Karma: Expansion of glacier lakes in recent decades in the Bhutan Himalayas, IAHS-AISH Publ., 264, 165–175, 2000. 

Bajracharya, S. R., Maharjan, S. B., and Shrestha, F.: The status and decadal change of glaciers in Bhutan from the 1980s to 2010 based on satellite data, Ann. Glaciol., 55, 159–166,, 2014. 

Benn, D. I., Warren, C. R., and Mottram, R. H.: Calving processes and the dynamics of calving glaciers, Earth-Sci. Rev., 82, 143–179,, 2007. 

Benn, D. I., Bolch, T., Hands, K., Gulley, J., Luckman, A., Nicholson, L. I., Quincey, D., Thompson, S., Toumi, R., and Wiseman, S.: Response of debris-covered glaciers in the Mount Everest region to recent warming, and implications for outburst flood hazards, Earth-Sci. Rev., 114, 156–174,, 2012. 

Berthier, E., Vincent, C., Magnússon, E., Gunnlaugsson, Á. Þ., Pitte, P., Le Meur, E., Masiokas, M., Ruiz, L., Pálsson, F., Belart, J. M. C., and Wagnon, P.: Glacier topography and elevation changes derived from Pléiades sub-meter stereo images, The Cryosphere, 8, 2275–2291,, 2014. 

Bolch, T., Kulkarni, A., Kääb, A., Huggel, C., Paul, F., Cogley, J. G., Frey, H., Kargel, J. S., Fujita, K., Scheel, M., Bajracharya, S., and Stoffel, M.: The state and fate of Himalayan glaciers, Science, 336, 310–314,, 2012. 

Boyce, E. S., Motyka, R. J., and Truffer, M.: Flotation and retreat of a lake-calving terminus, Mendenhall Glacier, southeast Alaska, USA, J. Glaciol., 53, 211–224,, 2007. 

Brun, F., Berthier, E., Wagnon, P., Kääb, A., and Treichler, D.: A spatially resolved estimate of High Mountain Asia glacier mass balances, 2000–2016, Nat. Geosci., 10, 668–673,, 2017. 

Brun, F., Wagnon, P., Berthier, E., Shea, J. M., Immerzeel, W. W., Kraaijenbrink, P. D. A., Vincent, C., Reverchon, C., Shrestha, D., and Arnaud, Y.: Ice cliff contribution to the tongue-wide ablation of Changri Nup Glacier, Nepal, central Himalaya, The Cryosphere, 12, 3439–3457,, 2018. 

Brun, F., Wagnon, P., Berthier, E., Jomelli, V., Maharjan, S. B., Shrestha, F., and Kraaijenbrink, P. D. A.: Heterogeneous Influence of Glacier Morphology on the Mass Balance Variability in High Mountain Asia, J. Geophys. Res.-Earth, 124, 1331–1345,, 2019. 

Carrivick, J. L., Tweed, F. S., Ng, F., Quincey, D. J., Mallalieu, J., Ingeman-Nielsen, T., Mikkelsen, A. B., Palmer, S. J., Yde, J. C., Homer, R., Russell, A. J., and Hubbard, A.: Ice-Dammed Lake Drainage Evolution at Russell Glacier, West Greenland, Front. Earth Sci., 5, 100,, 2017. 

Chen, F., Zhang, M., Guo, H., Allen, S., Kargel, J. S., Haritashya, U. K., and Watson, C. S.: Annual 30 m dataset for glacial lakes in High Mountain Asia from 2008 to 2017, Earth Syst. Sci. Data, 13, 741–766,, 2021. 

Dehecq, A., Gourmelen, N., Gardner, A. S., Brun, F., Goldberg, D., Nienow, P. W., Berthier, E., Vincent, C., Wagnon, P., and Trouvé, E.: Twenty-first century glacier slowdown driven by mass loss in High Mountain Asia, Nat. Geosci., 12, 22–27,, 2019. 

Farinotti, D., Huss, M., Bauder, A., Funk, M., and Truffer, M.: A method to estimate the ice volume and ice-thickness distribution of alpine glaciers, J. Glaciol., 55, 422–430,, 2009. 

Fujita, K., Suzuki, R., Nuimura, T., and Sakai, A.: Performance of ASTER and SRTM DEMs, and their potential for assessing glacial lakes in the Lunana region, Bhutan Himalaya, J. Glaciol., 54, 220–228,, 2008. 

Fujita, K., Sakai, A., Takenaka, S., Nuimura, T., Surazakov, A. B., Sawagaki, T., and Yamanokuchi, T.: Potential flood volume of Himalayan glacial lakes, Nat. Hazards Earth Syst. Sci., 13, 1827–1839,, 2013. 

Gardelle, J., Arnaud, Y., and Berthier, E.: Contrasted evolution of glacial lakes along the Hindu Kush Himalaya mountain range between 1990 and 2009, Global Planet. Change, 75, 47–55,, 2011. 

Gardner, A. S., Moholdt, G., Scambos, T., Fahnstock, M., Ligtenberg, S., van den Broeke, M., and Nilsson, J.: Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years, The Cryosphere, 12, 521–547,, 2018. 

Gardner, A. S., Fahnestock, M. A., and Scambos, T. A.: ITS_LIVE Regional Glacier and Ice Sheet Surface Velocities, National Snow and Ice Data Center,, 2019. 

Girod, L., Nuth, C., Kääb, A., Etzelmüller, B., and Kohler, J.: Terrain changes from images acquired on opportunistic flights by SfM photogrammetry, The Cryosphere, 11, 827–840,, 2017. 

Haritashya, U. K., Kargel, J. S., Shugar, D. H., Leonard, G. J., Strattman, K., Watson, C. S., Shean, D., Harrison, S., Mandli, K. T., and Regmi, D.: Evolution and Controls of Large Glacial Lakes in the Nepal Himalaya, Remote Sensing, 10, 798,, 2018. 

Heid, T. and Kääb, A.: Evaluation of existing image matching methods for deriving glacier surface displacements globally from optical satellite imagery, Remote Sens. Environ., 118, 339–355,, 2012. 

Kääb, A., Berthier, E., Nuth, C., Gardelle, J., and Arnaud, Y.: Contrasting patterns of early twenty-first-century glacier mass change in the Himalayas, Nature, 488, 495–498,, 2012. 

King, O., Dehecq, A., Quincey, D., and Carrivick, J.: Contrasting geometric and dynamic evolution of lake and land-terminating glaciers in the central Himalaya, Global Planet. Change, 167, 46–60,, 2018. 

King, O., Bhattacharya, A., Bhambri, R., and Bolch, T.: Glacial lakes exacerbate Himalayan glacier mass loss, Sci. Rep.-UK, 9, 18145,, 2019. 

Komori, J.: Recent expansions of glacial lakes in the Bhutan Himalayas, Quatern. Int., 184, 177–186,, 2008. 

Liu, Q., Mayer, C., Wang, X., Nie, Y., Wu, K., Wei, J., and Liu, S.: Interannual flow dynamics driven by frontal retreat of a lake-terminating glacier in the Chinese Central Himalaya, Earth Planet. Sc. Lett., 546, 116450,, 2020. 

Maurer, J. M., Rupper, S. B., and Schaefer, J. M.: Quantifying ice loss in the eastern Himalayas since 1974 using declassified spy satellite imagery, The Cryosphere, 10, 2203–2215,, 2016. 

Maurer, J. M., Schaefer, J. M., Rupper, S., and Corley, A.: Acceleration of ice loss across the Himalayas over the past 40 years, Sci. Adv., 5, eaav7266,, 2019. 

Maurer, J. M., Schaefer, J. M., Russell, J. B., Rupper, S., Wangdi, N., Putnam, A. E., and Young, N.: Seismic observations, numerical modeling, and geomorphic analysis of a glacier lake outburst flood in the Himalayas, Sci. Adv., 6, eaav3645,, 2020. 

Messerli, A. and Grinsted, A.: Image georectification and feature tracking toolbox: ImGRAFT, Geosci. Instrum. Method. Data Syst., 4, 23–34,, 2015. 

Mertes, J. R., Gulley, J. D., Benn, D. I., Thompson, S. S., and Nicholson, L. I.: Using structure-from-motion to create glacier DEMs and orthoimagery from historical terrestrial and oblique aerial imagery: SfM on Differing Historical Glacier Imagery Sets, Earth Surf. Proc. Land., 42, 2350–2364,, 2017. 

National Institute of Advanced Industrial Science and Technology (AIST): MADAS (METI AIST Data Archive System),, last access: 19 October 2021. 

Nie, Y., Sheng, Y., Liu, Q., Liu, L., Liu, S., Zhang, Y., and Song, C.: A regional-scale assessment of Himalayan glacial lake changes using satellite observations from 1990 to 2015, Remote Sens. Environ., 189, 1–13,, 2017. 

Nuimura, T., Fujita, K., Fukui, K., Asahi, K., Aryal, R., and Ageta, Y.: Temporal Changes in Elevation of the Debris-Covered Ablation Area of Khumbu Glacier in the Nepal Himalaya since 1978, Arct. Antarct. Alp. Res., 43, 246–255,, 2011. 

Nuimura, T., Fujita, K., Yamaguchi, S., and Sharma, R. R.: Elevation changes of glaciers revealed by multitemporal digital elevation models calibrated by GPS survey in the Khumbu region, Nepal Himalaya, 1992–2008, J. Glaciol., 58, 648–656,, 2012. 

Nuimura, T., Sakai, A., Taniguchi, K., Nagai, H., Lamsal, D., Tsutaki, S., Kozawa, A., Hoshina, Y., Takenaka, S., Omiya, S., Tsunematsu, K., Tshering, P., and Fujita, K.: The GAMDAM glacier inventory: a quality-controlled inventory of Asian glaciers, The Cryosphere, 9, 849–864,, 2015. 

Nuimura, T., Fujita, K., and Sakai, A.: Downwasting of the debris-covered area of Lirung Glacier in Langtang Valley, Nepal Himalaya, from 1974 to 2010, Quatern. Int., 455, 93–101,, 2017. 

Paul, F., Barrand, N. E., Baumann, S., Berthier, E., Bolch, T., Casey, K., Frey, H., Joshi, S. P., Konovalov, V., Le Bris, R., Mölg, N., Nosenko, G., Nuth, C., Pope, A., Racoviteanu, A., Rastner, P., Raup, B., Scharrer, K., Steffen, S., and Winsvold, S.: On the accuracy of glacier outlines derived from remote-sensing data, Ann. Glaciol., 54, 171–182,, 2013. 

Pronk, J. B., Bolch, T., King, O., Wouters, B., and Benn, D. I.: Contrasting surface velocities between lake- and land-terminating glaciers in the Himalayan region, The Cryosphere, 15, 5577–5599,, 2021. 

Quincey, D. J., Richardson, S. D., Luckman, A., Lucas, R. M., Reynolds, J. M., Hambrey, M. J., and Glasser, N. F.: Early recognition of glacial lake hazards in the Himalaya using remote sensing datasets, Global Planet. Change, 56, 137–152,, 2007. 

Quincey, D. J., Luckman, A., and Benn, D.: Quantification of Everest region glacier velocities between 1992 and 2002, using satellite radar interferometry and feature tracking, J. Glaciol., 55, 596–606,, 2009. 

Richardson, S. D. and Reynolds, J. M.: An overview of glacial hazards in the Himalayas, Quatern. Int., 65–66, 31–47,, 2000. 

Sakai, A.: Brief communication: Updated GAMDAM glacier inventory over high-mountain Asia , The Cryosphere, 13, 2043–2049,, 2019. 

Sakai, A. and Fujita, K.: Formation conditions of supraglacial lakes on debris-covered glaciers in the Himalaya, J. Glaciol., 56, 177–181,, 2010. 

Sakai, A., Nishimura, K., Kadota, T., and Takeuchi, N.: Onset of calving at supraglacial lakes on debris-covered glaciers of the Nepal Himalaya, J. Glaciol., 55, 909–917,, 2009. 

Sakakibara, D. and Sugiyama, S.: Ice-front variations and speed changes of calving glaciers in the Southern Patagonia Icefield from 1984 to 2011: calving glaciers in southern Patagonia, J. Geophys. Res.-Earth, 119, 2541–2554,, 2014. 

Sato, Y., Fujita, K., Inoue, H., Sunako, S., Sakai, A., Tsushima, A., Podolskiy, E., Kayastha, R., and Kayastha, R.: Ice cliff dynamics of debris-covered Trakarding Glacier in the Rolwaling region, Nepal Himalaya, Front. Earth Sci., 9, 623623,, 2021. 

Shean, D. E., Bhushan, S., Montesano, P., Rounce, D. R., Arendt, A., and Osmanoglu, B.: A Systematic, Regional Assessment of High Mountain Asia Glacier Mass Balance, Front. Earth Sci., 7, 363,, 2020. 

Shugar, D. H., Burr, A., Haritashya, U. K., Kargel, J. S., Watson, C. S., Kennedy, M. C., Bevington, A. R., Betts, R. A., Harrison, S., and Strattman, K.: Rapid worldwide growth of glacial lakes since 1990, Nat. Clim. Change, 10, 939–945,, 2020. 

Song, C., Sheng, Y., Wang, J., Ke, L., Madson, A., and Nie, Y.: Heterogeneous glacial lake changes and links of lake expansions to the rapid thinning of adjacent glacier termini in the Himalayas, Geomorphology, 280, 30–38,, 2017. 

Suzuki, R., Fujita, K., and Ageta, Y.: Spatial distribution of thermal properties on debris-covered glaciers in the Himalayas derived from ASTER data, Bull. Glaciol. Res., 24, 13–22, 2007. 

Tshering, P. and Fujita, K.: First in situ record of decadal glacier mass balance (2003–2014) from the Bhutan Himalaya, Ann. Glaciol., 57, 289–294,, 2016. 

Tsutaki, S., Nishimura, D., Yoshizawa, T., and Sugiyama, S.: Changes in glacier dynamics under the influence of proglacial lake formation in Rhonegletscher, Switzerland, Ann. Glaciol., 52, 31–36,, 2011. 

Tsutaki, S., Sugiyama, S., Nishimura, D., and Funk, M.: Acceleration and flotation of a glacier terminus during formation of a proglacial lake in Rhonegletscher, Switzerland, J. Glaciol., 59, 559–570,, 2013. 

Tsutaki, S., Fujita, K., Nuimura, T., Sakai, A., Sugiyama, S., Komori, J., and Tshering, P.: Contrasting thinning patterns between lake- and land-terminating glaciers in the Bhutanese Himalaya, The Cryosphere, 13, 2733–2750,, 2019. 

United States Geological Survey (USGS): Earthexplorer,, last access: 19 October 2021.  

Vincent, C., Wagnon, P., Shea, J. M., Immerzeel, W. W., Kraaijenbrink, P., Shrestha, D., Soruco, A., Arnaud, Y., Brun, F., Berthier, E., and Sherpa, S. F.: Reduced melt on debris-covered glaciers: investigations from Changri Nup Glacier, Nepal, The Cryosphere, 10, 1845–1858,, 2016. 

Watson, C. S., Kargel, J. S., Shugar, D. H., Haritashya, U. K., Schiassi, E., and Furfaro, R.: Mass Loss From Calving in Himalayan Proglacial Lakes, Front. Earth Sci., 7, 342,, 2020. 

Wei, J., Liu, S., Wang, X., Zhang, Y., Jiang, Z., Wu, K., Zhang, Z., and Zhang, T.: Longbasaba Glacier recession and contribution to its proglacial lake volume between 1988 and 2018, J. Glaciol., 67, 473–484,, 2021. 

Westoby, M. J., Glasser, N. F., Brasington, J., Hambrey, M. J., Quincey, D. J., and Reynolds, J. M.: Modelling outburst floods from moraine-dammed glacial lakes, Earth-Sci. Rev., 134, 137–159,, 2014. 

Yamada, T., Naito, N., Kohshima, S., Fushimi, H., Nakazawa, F., Segawa, T., Uetake, J., Suzuki, R., Sato, N., Karma, Chhetri, I. K., Gyenden, L., Yabuki, H., and Chikita, K.: Outline of 2002 – research activities on glaciers and glacier lakes in Lunana region, Bhutan Himalaya, Bull. Glaciol. Res., 21, 79–90, 2004. 

Zhang, G., Yao, T., Xie, H., Wang, W., and Yang, W.: An inventory of glacial lakes in the Third Pole region and their changes in response to global warming, Global Planet. Change, 131, 148–157,, 2015. 

Zhang, G., Bolch, T., Allen, S., Linsbauer, A., Chen, W., and Wang, W.: Glacial lake evolution and glacier–lake interactions in the Poiqu River basin, central Himalaya, 1964–2017, J. Glaciol., 65, 347–365,, 2019. 

Zheng, G., Allen, S. K., Bao, A., Ballesteros-Cánovas, J. A., Huss, M., Zhang, G., Li, J., Yuan, Y., Jiang, L., Yu, T., Chen, W., and Stoffel, M.: Increasing risk of glacial lake outburst floods from future Third Pole deglaciation, Nat. Clim. Change, 11, 411–417,, 2021. 

Short summary
We investigate fluctuations in Bhutanese lake-terminating glaciers focusing on the dynamics change before and after proglacial lake formation at Thorthormi Glacier (TG) based on photogrammetry, satellite, and GPS surveys. The thinning rate of TG became double compared to before proglacial lake formation, and the flow velocity has also sped up considerably. Those changes would be due to the reduction in longitudinal ice compression by the detachment of the glacier terminus from the end moraine.