The purpose of this study was to determine the glacier displacement, velocity, and thickness of seven major glaciers of Bhutan and predicted potential glacier lake formation site with its depth. We named the glacier identification (ID) number 1–7 for seven glaciers. From the study, the glacier velocity between the central trunk and snout saw rapid fluctuations in 1976–1978 with an average uncertainty velocity of ± 27.10 m/year and a decreasing velocity trend. The year 2013–2014 has the lowest uncertainty in glacier velocity, with a value of ±1.24 m/year. The glacier velocity progressively increases from the snout to the main central trunk for all the glaciers with a value of 0 to 98.63 m/year. The glacier with the highest average velocity is glacier ID 5, which has a velocity of 25.58 m/year. From 2000 and 2022, all of the glaciers’ thicknesses significantly decreased from 0 to 468.2 m. The thickness of glacier ID 6 was lowered by −192.3 ± 1.99 m, making it the highest among the seven glaciers. In the future, a glacier lake is predicted to form at the base of each glacier. Glacier ID 6 is predicted to form the largest lake with a surface area of 2.572 km2 and a depth of 208.5 m.

  • The paper is related to glacier dynamics in Bhutan.

  • Research related to glacier dynamics was never published until now.

  • This paper presents glacier displacement, velocity, and thickness change over the year.

  • This paper also presents predicted potential glacier lake formation sites with its possible depth.

  • Seven largest glaciers of Bhutan were considered as a study area for this study.

Glaciers are a piece of huge ice moving over land mass due to their weight and steep underlying terrain (Shreve 1985). Glaciers are normally found in the high mountains and polar regions of the world. The glaciers are retreating or advancing owing to ablations and accumulation (Hooker & Fitzharris 1999). A gradual retreat may result in the disappearance of the glacier in the future (Paul et al. 2007). According to numerous studies, the majority of the world's glaciers are significantly shrinking as a result of the increase in average global temperature (Dyurgerov & Meier 2000). The melting of the glaciers might lead to environmental and economic problems at local, regional, continental, and even global scales (Haeberli & Weingartner 2020).

Glaciers hold 70% of the freshwater worldwide (Hasnain 1999). The Himalayan Mountain is also known as the third pole and is the home of many glaciers. Many rivers in South Asian countries are fed by snow and glaciers in the Himalayan region (Khajuria et al. 2022). The rivers that originated from the glaciers are crucial for economic growth (Lal 2000). Nonetheless, the glaciers are undergoing several changes, including recession, the formation of glacier lakes, a decrease in thickness, and changes to the glacial mass balance (Wang et al. 2013). Similarly, global warming is causing the Himalayan glaciers to retreat at an alarming rate. Global warming-induced glacier retreat may result in dwindling water supplies and an increase in the danger of glacier lake outburst flood (Hasnain 1999). Kumar et al. (2022) stated a 2 °C rise in temperature led to a 53% increase in catchment discharge in the Himalayan region. The knowledge of the dimensions of Himalayan glaciers and their behavior in response to climate change is still limited, due to their remoteness, inaccessibility, and harsh topography (Momblanch et al. 2019). The advancement of remote sensing, digital elevation models, and geographic information systems (GIS) provides an effective approach to assessing and quantifying glacier changes. Remote sensing also plays an instrumental role in tracking glacier geometry, study of glacial lake formation, determination of equilibrium line altitude (ELA) (Garg et al. 2022), monitoring annual mass balance changes, and estimating various movement rates (Racoviteanu et al. 2008).

Glaciers in Bhutan constitute about 1.64% (630 sq. km) of the total area (NCHM 2019). As per the Bhutan Glacier Inventory 2018, there are 700 glaciers in Bhutan (NCHM 2019). The Bhutanese national newspaper Kuensel reported that 1,146.16 million tons of glacier ice are lost every year (Dema 2022). The south-facing glaciers encounter increasing retreats every year compared to north-facing glaciers (Iwata 2010). Bhutanese glaciers are very important for economic development such as hydropower generation, irrigation, and agriculture. Concerns regarding the melting of glaciers in Bhutan have been expressed by climate experts and decision-makers. However, despite these concerns, there is a lack of concrete evidence regarding glacial recession. This knowledge gap underscores the need for comprehensive research on glacier dynamics to address the lack of empirical data. While the significance of glaciers in Bhutan has been acknowledged through verbal discussions, there is a notable scarcity of publications focusing on the dynamics of glaciers, such as glacier velocity, thickness, and potential glacier lakes. Only a few studies exist that cover glacier velocity in Bhutan, including Luggye glaciers (16 m/year) in 1988–1993 (Ageta & Iwata 1998), Tarina Glacier (30–35 m/year) in 1976–1998 (Ageta & Iwata 1998), and Gangju La glacier (17 m/year) in 2013–2014 (Tshering & Fujita 2016). However, no prior research has specifically investigated the dynamics of glaciers in Bhutan. Furthermore, the challenging terrain characterized by ruggedness, high elevation, large area, and severe weather conditions makes traditional field-based techniques for studying glacier dynamics impractical. In light of these challenges, the use of satellite data and GIS is presented as the most feasible and practical method for studying such complex terrains. To address these gaps and challenges, this study aims to calculate glacier displacement and velocity over the period from 1976 to 2022, analyze variations in glacier thickness between 2000 and 2022, and identify potential sites for glacier lake formation using remote sensing data combined with GIS techniques. The research aims to provide valuable insights for policymakers, enabling them to make informed decisions regarding glacier-related concerns. Additionally, the study result is expected to serve as a warning system for downstream communities who may be affected by glacier-related hazards.

As shown in Figure 1, the study involves extracting the study area, calculating displacement, velocity, thickness, and the topography of the glacier bed, as well as identifying prospective lake formation areas.
Figure 1

Overview of workflow.

Figure 1

Overview of workflow.

Close modal

Study area

The study area (Figure 2) was extracted using the 2010 glacier inventory of Bhutan prepared by ICIMOD. The study area consists of seven large glaciers of Bhutan. Since the glacier name is not available in any database or research papers, we have assigned glacier ID 1 to glacier ID 7 for seven major glaciers. Detailed information on individual glaciers used for the study is mentioned in Table 1.
Table 1

Detail of the study area

Glacier IDLocation
GLIMS_IDBasinElevation (m)
Slope (degree)
Area (km2)
LatitudeLongitudeMinimumMaximumMinimumMaximum
28.133 89.996 G089996E28133N Phochu 4,386 7,084 69 28.07 
28.125 90.161 G090161E28125N Phochu 4,085 6,663 74 36.55 
28.116 90.081 G090081E28116N Phochu 4,155 6,485 61 10.56 
28.005 90.476 G090476E28005N Mangdechu 4,679 6,706 54 14.71 
28.022 90.42 G090420E28022N Mangdechu 4,943 7,231 69 20.252 
28.031 90.521 G090521E28031N Chamkhar 4,700 6,171 63 15.58 
28.026 90.786 G090786E28026N Chamkhar 4,864 5,886 56 5.84 
Glacier IDLocation
GLIMS_IDBasinElevation (m)
Slope (degree)
Area (km2)
LatitudeLongitudeMinimumMaximumMinimumMaximum
28.133 89.996 G089996E28133N Phochu 4,386 7,084 69 28.07 
28.125 90.161 G090161E28125N Phochu 4,085 6,663 74 36.55 
28.116 90.081 G090081E28116N Phochu 4,155 6,485 61 10.56 
28.005 90.476 G090476E28005N Mangdechu 4,679 6,706 54 14.71 
28.022 90.42 G090420E28022N Mangdechu 4,943 7,231 69 20.252 
28.031 90.521 G090521E28031N Chamkhar 4,700 6,171 63 15.58 
28.026 90.786 G090786E28026N Chamkhar 4,864 5,886 56 5.84 
Figure 2

Study area overlapped with Google Earth image.

Figure 2

Study area overlapped with Google Earth image.

Close modal

Datasets used

The Landsat datasets were downloaded from USGS earth explorer for Landsat 2 MSS, Landsat 5 TM, Landsat 7 ETM, and Landsat 8 OLI from 1976 to 2022. September to November are the optimum months to use Landsat images for glacier studies since there is little cloud cover, less seasonal snow cover, and the glaciers are fully exposed to their actual positions. The best images that are available for this study span from September 30 to December 17 as listed in Table 2. The Landsat images were used to derive glacier displacement and velocity while SRTM DEM and Copernicus DEM were used to derive glacier thickness in conjunction with the derived velocity. ICIMOD Glacier Inventory 2010 shapefile was downloaded from ICIMOD regional database (http://rds.icimod.org/Home/). The glacier inventory was used to extract the study area.

Table 2

Dataset list of datasets used for estimation of the glacier velocity

SensorDateRowPathBand usedResolution (m)
Landsat 8 OLI 19-11-2021 41 138 Pan 15 
Landsat 8 OLI 21-10-2022 Pan 15 
Landsat 8 OLI 12-10-2013 Pan 15 
Landsat 8 OLI 21-10-2014 Pan 15 
Landsat 7 ETM + 30-09-2000 Pan 15 
Landsat 7 ETM + 07-11-2002 Pan 15 
Landsat 5 TM 09-11-1994 Red 30 
Landsat 5 TM 30-11-1996 Red 30 
Landsat 5 TM 11-11-1989 Red 30 
Landsat 5 TM 14-11-1990 Red 30 
Landsat 2 MSS 17-12-1976 Red 60 
Landsat 2 MSS 07-12-1978 Red 60 
SensorDateRowPathBand usedResolution (m)
Landsat 8 OLI 19-11-2021 41 138 Pan 15 
Landsat 8 OLI 21-10-2022 Pan 15 
Landsat 8 OLI 12-10-2013 Pan 15 
Landsat 8 OLI 21-10-2014 Pan 15 
Landsat 7 ETM + 30-09-2000 Pan 15 
Landsat 7 ETM + 07-11-2002 Pan 15 
Landsat 5 TM 09-11-1994 Red 30 
Landsat 5 TM 30-11-1996 Red 30 
Landsat 5 TM 11-11-1989 Red 30 
Landsat 5 TM 14-11-1990 Red 30 
Landsat 2 MSS 17-12-1976 Red 60 
Landsat 2 MSS 07-12-1978 Red 60 

Estimation of glacier displacement and velocity using an optical feature tracking method

Glacier surface displacement was calculated using the sub-pixel correlation of a pair of Landsat images on a sliding window of a phase plane. The pair of Landsat images consist of pre- and post-event images with low snow and cloud coverage. The only limitation of the feature tracking method is that it used optical satellite images and the cloud coverage may result in misinterpretation of the images. Therefore, this study used Landsat images with cloud coverage less than 10%. The correlation of Landsat images is estimated using the Co-registration of Optically Sensed Images and Correlation (COSI-Corr). The COSI-Corr is freely available from http://www.tectonics.caltech.edu/ and it is patched in ENVI software. In the COSI-CORR, the two images are iteratively cross-correlated with an initial window size of 64 and a final window size of 32 pixels. A total of four iterations were conducted. Three outputs were obtained, namely: east-west displacement (E-W), north-south displacement (N-S), and signal-to-noise ratio (SNR). As per Gopika et al. (2021), SNR is used to estimate the correlation accuracy, and all the pixels with SNR values less than 0.9 are considered erroneous. Therefore, SNR values less than 0.9 are discarded in this study. The remaining pixels can be used for the calculation of horizontal and vertical displacements.

The resultant displacement is calculated using Equation (1):
(1)
The resultant displacement is calculated using displacement produced from north-south and east-west. The velocity of the glacier per year is calculated using Equation (2):
(2)
where velocity is in m/year, displacement is the average shift of the glacier in meters, and time is the number of days between two images acquired.

Estimation of glacier thickness and its bed topography

To estimate the glacier thickness, SRTM DEM of 2000 and Copernicus DEM of 2019 were used. The DEMs were reclassified into 100 m contour intervals and converted to contour lines with a vertical contour interval of 100 m. Glacier boundary and moraine are also required to find the glacier thickness. The flowlines were digitized along the central trunks of the glaciers. All the datasets were projected to UTM45N. The projected datasets were inserted into HIGHTHIM python package to estimate the glacier thickness.

The glacier thickness was estimated based on the user manual of HIGHTHIM as published by Kulkarni et al. (2018), and it is calculated using Equation (3):
(3)
where H is the thickness of the glacier, Us is the surface velocity, A is the creep parameter with its value 3.24 × 10−24 Pa−3 S−1, f is the shape factor and its value is 0.8. The shape factor is the ratio between the driving stress and basal stress along a glacier, ρ is the ice density with its constant value of 900 kg m−3, and g is the acceleration due to gravity which is 9.8 m s−2. Slope (α) is calculated from 100 m contour derived from SRTM DEM for the year 2000 and Copernicus DEM for the year 2019. Ice thickness is calculated along the central flow line for each sectional area between successive 100 m contours derived from SRTM DEM and Copernicus DEM.
As per Kulkarni et al. (2018), the glacier thickness was derived from the equation of laminar flow of glacier ice and basal shear stress which is expressed as:
(4)
where Us and Ub are surface and basal velocities, respectively. H is the ice thickness (m), A is the creep parameter, Glen's flow law exponent (n), is assumed to be 3. The basal stress is modeled as:
(5)
where ρ is the ice density, g is the acceleration due to gravity, and f is a shape factor.

Uncertainties of glacier velocity and thickness

There will be an error in the velocity and thickness of the glaciers due to low-contrast images, snow, and cloud coverage using optical satellite images. To calculate the velocity error, the stable area of the glaciers was used to evaluate the error with a minimum terrain slope with a zero or known velocity and stable ground (Bolch et al. 2011). Similarly, to calculate glacier thickness error, 150 stable points on SRTM DEM and Corpinicus DEM were used. The SRTM DEM which was launched in 2000 and the Copernicus DEM of 2019 was used to find an error in glacier thickness. As per Bolch et al. (2011), the law of propagation (e) is used to find the uncertainty of velocity and thickness as follows:
(6)
(7)
where STDDEV is the standard deviation of the points over the glacier-free region, MED is the mean velocity or thickness difference of stable point, SE is the standard error, and n is the number of points used in the stable areas.

Glacier displacement and velocity

The glacier displacement and velocity are calculated using Equations (1) and (2), respectively, and the data used are mentioned in Table 2. The glacier displacement and velocity for the years 2021–2022 are given in Figure 3(a) and 3(b), respectively.
Figure 3

(a) Glacier displacement of the year 2021–2022, (b) glacier velocity of the year 2021–2022, and (c) glacier displacement vector of the year 2021–2022.

Figure 3

(a) Glacier displacement of the year 2021–2022, (b) glacier velocity of the year 2021–2022, and (c) glacier displacement vector of the year 2021–2022.

Close modal

For the years 2001–2022, the maximum displacement and velocity are located in the upper central trunk with a value of 89.63 m and 98.63 m/year, respectively. The displacement and velocity gradually decrease toward the glacier terminus. The findings of our study on glacier displacement and velocity are consistent with previous research conducted on Himalayan glaciers by Gantayat et al. (2014) and Singh et al. (2021).

The graph (Figure 4) shows the year-wise variation of the glacier velocity pattern from 1976 to 2022 for seven major glaciers of Bhutan. From Figure 4, it is observed that the velocity of all the glaciers progressively increases from lower reach (<20 m/year) to upper central trunk (>30 m/year) for all seven glaciers. From the graph (Figure 4), the glacier velocities are abrupt and irregular between lower reach to the central trunk for 1976–1978. The abrupt and irregular glacier velocity changes observed between the lower reach and central trunk from 1976 to 1978 may be attributed to faster glacier dynamics activities in the 1970s and 1980s or potential errors induced by low-quality data. Further investigation and data quality improvement would be beneficial to better understand this period of irregular velocity changes. Comparing our results with Singh et al. (2021), who studied seven glaciers in the Dhauliganga basin of Uttarakhand, India, we find a similar pattern of decreasing glacier velocity over time. Singh et al. (2021) reported a 40% decrease in glacier velocity between 1993 and 2015, while our study shows a general decrease in high glacier velocity dynamics from 1976 to 2022. This suggests a consistent trend of declining glacier velocities in both the Bhutanese and Indian Himalayan regions.
Figure 4

Graph showing the velocity variation from 1976 to 2022 for the study area.

Figure 4

Graph showing the velocity variation from 1976 to 2022 for the study area.

Close modal

Table 3 shows the detailed comparison of individual glacier velocities. The average velocity for individual glaciers was calculated by averaging different year glacier velocities. The average velocity for individual glaciers was calculated to compare the degree of glacier dynamics between the different glaciers. From Table 3, it is observed that glacier ID 5 undergone the highest glacier average velocity (25.58 m/year), followed by glacier ID 1 (19.62 m/year), glacier ID 3 (18.28 m/year), glacier ID 2 (15.56 m/year), and glacier ID 6 (13.44 m/year) and it indicates these five glaciers are undergoing tremendous glacier change between 1976 and 2022. Glacier ID 7 and glacier ID 4 are undergoing the least alteration with an average velocity of 7.70 and 9.19 m/year, respectively.

Table 3

The average displacement and velocity of the study area

Glacier IDImage acquisition dateTime interval (days)Displacement (m)Velocity (m/year)Average velocity (1976–2022)
17/12/1976–7/12/1978 720 43.66 47.86 19.62 
11/11/1989–14/11/1990 368 20.76 19.73 
9/11/1994–30/11/1996 751 27.94 13.28 
30/9/2000–7/11/2002 768 20.19 9.81 
12/10/2013–31/10/2014 384 20.21 20.05 
19/11/2021–21/10/2022 333 13.81 7.00 
17/12/1976–7/12/1978 720 45.49 23.06 15.56 
11/11/1989–14/11/1990 368 19.98 19.82 
9/11/1994–30/11/1996 751 22.10 10.74 
30/9/2000–7/11/2002 768 31.67 15.05 
12/10/2013–31/10/2014 384 14.34 13.63 
19/11/2021–21/10/2022 333 10.11 11.08 
17/12/1976–7/12/1978 720 39.07 19.81 18.28 
11/11/1989–14/11/1990 368 23.16 22.97 
9/11/1994–30/11/1996 751 25.85 12.56 
30/9/2000–7/11/2002 768 28.42 13.51 
12/10/2013–31/10/2014 384 20.63 19.61 
19/11/2021–21/10/2022 333 19.37 21.23 
17/12/1976–7/12/1978 720 37.65 19.09 9.19 
11/11/1989–14/11/1990 368 11.79 11.69 
9/11/1994–30/11/1996 751 13.26 6.44 
30/9/2000–7/11/2002 768 13.06 6.21 
12/10/2013–31/10/2014 384 7.74 7.36 
19/11/2021–21/10/2022 333 3.96 4.34 
17/12/1976–7/12/1978 720 34.47 17.47 25.58 
11/11/1989–14/11/1990 368 35.21 34.92 
9/11/1994–30/11/1996 751 40.13 19.50 
30/9/2000–7/11/2002 768 42.43 20.17 
12/10/2013–31/10/2014 384 33.52 31.86 
19/11/2021–21/10/2022 333 26.94 29.53 
17/12/1976–7/12/1978 720 26.44 13.40 13.44 
11/11/1989–14/11/1990 368 18.92 18.77 
9/11/1994–30/11/1996 751 27.14 13.19 
30/9/2000–7/11/2002 768 25.71 12.22 
12/10/2013–31/10/2014 384 12.67 12.04 
19/11/2021–21/10/2022 333 10.03 10.99 
17/12/1976–7/12/1978 720 25.19 12.77 7.70 
11/11/1989–14/11/1990 368 11.58 11.49 
9/11/1994–30/11/1996 751 12.60 6.12 
30/9/2000–7/11/2002 768 13.73 6.53 
12/10/2013–31/10/2014 384 6.69 6.36 
19/11/2021–21/10/2022 333 2.68 2.94 
Glacier IDImage acquisition dateTime interval (days)Displacement (m)Velocity (m/year)Average velocity (1976–2022)
17/12/1976–7/12/1978 720 43.66 47.86 19.62 
11/11/1989–14/11/1990 368 20.76 19.73 
9/11/1994–30/11/1996 751 27.94 13.28 
30/9/2000–7/11/2002 768 20.19 9.81 
12/10/2013–31/10/2014 384 20.21 20.05 
19/11/2021–21/10/2022 333 13.81 7.00 
17/12/1976–7/12/1978 720 45.49 23.06 15.56 
11/11/1989–14/11/1990 368 19.98 19.82 
9/11/1994–30/11/1996 751 22.10 10.74 
30/9/2000–7/11/2002 768 31.67 15.05 
12/10/2013–31/10/2014 384 14.34 13.63 
19/11/2021–21/10/2022 333 10.11 11.08 
17/12/1976–7/12/1978 720 39.07 19.81 18.28 
11/11/1989–14/11/1990 368 23.16 22.97 
9/11/1994–30/11/1996 751 25.85 12.56 
30/9/2000–7/11/2002 768 28.42 13.51 
12/10/2013–31/10/2014 384 20.63 19.61 
19/11/2021–21/10/2022 333 19.37 21.23 
17/12/1976–7/12/1978 720 37.65 19.09 9.19 
11/11/1989–14/11/1990 368 11.79 11.69 
9/11/1994–30/11/1996 751 13.26 6.44 
30/9/2000–7/11/2002 768 13.06 6.21 
12/10/2013–31/10/2014 384 7.74 7.36 
19/11/2021–21/10/2022 333 3.96 4.34 
17/12/1976–7/12/1978 720 34.47 17.47 25.58 
11/11/1989–14/11/1990 368 35.21 34.92 
9/11/1994–30/11/1996 751 40.13 19.50 
30/9/2000–7/11/2002 768 42.43 20.17 
12/10/2013–31/10/2014 384 33.52 31.86 
19/11/2021–21/10/2022 333 26.94 29.53 
17/12/1976–7/12/1978 720 26.44 13.40 13.44 
11/11/1989–14/11/1990 368 18.92 18.77 
9/11/1994–30/11/1996 751 27.14 13.19 
30/9/2000–7/11/2002 768 25.71 12.22 
12/10/2013–31/10/2014 384 12.67 12.04 
19/11/2021–21/10/2022 333 10.03 10.99 
17/12/1976–7/12/1978 720 25.19 12.77 7.70 
11/11/1989–14/11/1990 368 11.58 11.49 
9/11/1994–30/11/1996 751 12.60 6.12 
30/9/2000–7/11/2002 768 13.73 6.53 
12/10/2013–31/10/2014 384 6.69 6.36 
19/11/2021–21/10/2022 333 2.68 2.94 

Uncertainty of glacier velocity

The overall uncertainty of glacier velocity was calculated using Equation (6) for all seven glaciers, and the results are presented in Table 4. The velocity uncertainties vary across different periods, reflecting the quality and availability of data. In the year 1976–1978, a higher velocity uncertainty of ±27.10 m/year is observed, followed by 5.78 m/year in 1989–1990, 2.57 m/year in 1994–1996, 2.08 m/year in 2000–2002, 1.24 m/year in 2013–2014, and 2.15 m/year in the year 2021–2022. This can be attributed to several factors, including the limited availability of data during that period and the coarser spatial resolution compared to more recent Landsat data.

Table 4

Velocity uncertainty for different years

Sl NoYearMEDSDNo. of pointsSE (m/year)Error (m/year)
(m/year)(m/year)
1976–1978 25.87 98.98 135 8.08 27.10 
1989–1990 5.76 4.46 150 0.36 5.78 
1994–1996 2.56 2.94 150 0.24 2.57 
2000–2022 2.07 1.20 150 0.10 2.08 
2013–2014 1.24 0.92 150 0.08 1.24 
2021–2022 2.15 2.15 150 0.18 2.15 
Sl NoYearMEDSDNo. of pointsSE (m/year)Error (m/year)
(m/year)(m/year)
1976–1978 25.87 98.98 135 8.08 27.10 
1989–1990 5.76 4.46 150 0.36 5.78 
1994–1996 2.56 2.94 150 0.24 2.57 
2000–2022 2.07 1.20 150 0.10 2.08 
2013–2014 1.24 0.92 150 0.08 1.24 
2021–2022 2.15 2.15 150 0.18 2.15 

The graph (Figure 4) and Table 3 clearly demonstrate the abruptness in estimated velocity and the uncertainty of velocity during the years 1976–1978. It is important to note that the reliability of the data during this period may be compromised due to limited availability and lower quality. The spatial resolution of the data used for that time period might have been coarser, further contributing to the higher uncertainty observed.

To improve the accuracy and reduce uncertainty in glacier velocity estimation, it is crucial to utilize high-quality and high-resolution data. Recent advancements in remote sensing technology, such as the availability of higher-resolution satellite imagery, can provide more precise measurements and reduce uncertainties in glacier velocity calculations.

Estimation of glacier thickness

The estimation of glacier thickness for the year 2022 is depicted in Figure 5 for the seven glaciers under study. Glacier ID 2 exhibits the highest mean glacier thickness of 86.29 m, followed by glacier ID 5 with 81.54 m. Glacier ID 7 has a mean thickness of 70.54 m, glacier ID 4 measures 67.48 m, and glacier ID 6 has a mean thickness of 61.39 m. Glacier ID 1 and glacier ID 3 show the lowest mean thicknesses of 57.75 and 65.16 m, respectively. Figure 5 reveals that the central part of the main trunk of the glaciers exhibits the maximum thickness for all seven glaciers. This observation aligns with the findings of a study conducted by Thakur et al. (2023) on the Gangotri glacier in the Indian Himalayas, which also demonstrated that the middle and upper regions of glaciers have higher thicknesses.
Figure 5

Glacier thickness for the year 2022.

Figure 5

Glacier thickness for the year 2022.

Close modal
To assess the changes in glacier thickness from 2000 to 2022, three cross-sectional lines (A–Aʹ), (B–Bʹ), and (C–Cʹ) were plotted for all seven glaciers, as shown in Figure 5. The thickness curve (Figure 6) was generated from these cross-sectional lines, indicating a considerable reduction in glacier thickness for all seven glaciers between 2000 and 2022. Notably, the central flow line exhibited the maximum changes in thickness compared to the outer periphery of the glaciers. For instance, as depicted in Figure 6(c) for glacier ID 3, the thickness was approximately 500 m in 2000, reduced to around 400 m by 2022. The uncertainty in glacier thickness for the years 2000 and 2022 was calculated using Equation (6), yielding a value of ±1.99 m.
Figure 6

Cross-section of the individual glacier for the years 2000 and 2022.

Figure 6

Cross-section of the individual glacier for the years 2000 and 2022.

Close modal
Figure 7 is constructed by subtracting the glacier thickness of the year 2022 from the thickness of the year 2000. It reveals that both ablation (negative values) and accumulation (positive values) occurred within individual glaciers. However, the ablation area (negative elevation difference) is larger than the accumulation area (positive elevation difference) for all seven glaciers. The central trunk and lower reach of the glaciers predominantly exhibit negative elevation differences, indicating a recession process leading to a reduction in thickness between 2000 and 2022. Among the glaciers, glacier ID 6 experienced the greatest decline in thickness, with a reduction of −192.3 m, while glacier ID 5 had the smallest reduction of −51.1 m. On the other hand, the upper reaches of the glaciers show positive elevation differences, indicating an increase in thickness (accumulation). Notably, glacier ID 5 displayed an increase in thickness of 223.6 m in a certain portion of the upper reach.
Figure 7

Change in glacier thickness in between the years 2000 and 2022.

Figure 7

Change in glacier thickness in between the years 2000 and 2022.

Close modal

Glacier bed topography and prediction of future glacier formation sites

The actual ground elevation below the glacier is the glacier bed topography and it is shown in Figure 8 for seven glaciers of Bhutan. The elevations of the glacier beds are comparably low in the region with large ice thickness. From Figure 8, it is seen that six glacier snouts start their bed topography at around 4,000 m above mean sea level, except glacier ID 3 whose bed topography starts at 3,845 m above mean sea level.
Figure 8

Glacier bed topography and predicted glacier lake formation sites with their depth.

Figure 8

Glacier bed topography and predicted glacier lake formation sites with their depth.

Close modal

The formation of glacier lakes is influenced by factors such as glacier recession and the topography of the glacier's base. Remya et al. (2019) observed that glacier lakes typically arise when glaciers recede and the bottom topography becomes overdeepened. The size and depth of the glacier lake vary depending on the depth of the depression in the bed topography and the size of the glacier (Sergienko 2013).

Glacier lakes that are adjacent to a glacier's toe can accelerate local glacier melt (Nick et al. 2010). However, the presence of glacier lakes also poses a significant risk of glacier lake outburst, which is a major disaster threat downstream (Bhushan et al. 2017). Therefore, it is crucial to monitor the development of glacier lakes to provide early warnings to downstream residents and enable adequate preparation.

In this study, it is anticipated that glacier lakes will form in the lower reaches of all seven glaciers examined. Based on the findings presented in Table 5, glacier ID 6 is expected to give rise to the largest glacier lake, with an area of 2.572 km2 and a depth of 208.5 m. On the other hand, glacier ID 2 is expected to form a smaller lake, with an area of 0.194 km2 and a depth of 41.3 m.

Table 5

Predicted area of lake formation site and its maximum depth

Glacier IDPredicted area of the glacier lake (km2)Predicted maximum depth of the glacier lake (m)
0.863 135.7 
0.194 41.3 
0.563 143.3 
0.413 71.4 
1.348 80.3 
2.572 208.5 
0.550 122.4 
Glacier IDPredicted area of the glacier lake (km2)Predicted maximum depth of the glacier lake (m)
0.863 135.7 
0.194 41.3 
0.563 143.3 
0.413 71.4 
1.348 80.3 
2.572 208.5 
0.550 122.4 

To validate these expectations, a comparison was made between the anticipated glacier lake and panchromatic and Google images. This comparison revealed the presence of previously discovered small lake patches, suggesting that a large lake will indeed form in the future.

Validation of glacier velocity

In the absence of high-resolution data on glacier depth, the validation of glacier thickness was not conducted in this study. However, global data for glacier velocity was available and downloaded from the website https://nsidc.org/apps/itslive/. Table 6 presents the global velocity data alongside the estimated velocity derived from the current study.

Table 6

Global velocity and the estimated velocity derived from the study

Glacier IDEstimated velocityGlobal velocity
16.92 14.58 
15.56 9.99 
18.28 16.68 
9.2 5.97 
25.58 33.22 
13.44 8.5 
7.7 2.98 
Glacier IDEstimated velocityGlobal velocity
16.92 14.58 
15.56 9.99 
18.28 16.68 
9.2 5.97 
25.58 33.22 
13.44 8.5 
7.7 2.98 

To assess the agreement between the estimated and global glacier velocities, a comparison curve was plotted (Figure 9) to examine the trend of the estimated velocity in relation to the global velocity. It was observed that the estimated and global glacier velocities exhibited a synchronous trend, indicating a similarity in their patterns.
Figure 9

Comparison of global glacier velocity and estimated glacier velocity.

Figure 9

Comparison of global glacier velocity and estimated glacier velocity.

Close modal

Furthermore, the Pearson correlation coefficient was calculated using Equation (8) to quantify the correlation between the estimated and global velocities. The obtained correlation coefficient was found to be 0.965, which is close to 1. This high correlation coefficient indicates a strong relationship between the estimated and global glacier velocities, providing additional evidence that the estimated velocity is accurate.

While the study did not validate the glacier thickness directly due to the lack of high-resolution data, the validation of the estimated velocity using global data and correlation analysis lends support to the reliability of the estimated glacier velocity in this research:
(8)
where r is the correlation coefficient, xi is the estimated velocity, is the mean estimated velocity, yi is the global velocity, and is the mean of global velocity.

Glaciers are sources of many rivers in Bhutan that are used for hydropower generation, agriculture, industry, and water for consumption. Apart from Bhutan, some parts of India also benefit from rivers that flow from Bhutan. The glaciers are important for the Bhutanese economy but its undergoing tremendous changes due to climate change. Yet, there is no concrete evidence that climate change is causing the glaciers to retreat. This study offers small evidence of glacier alteration in Bhutan's seven largest glaciers. Since there are no names for these glaciers in any database or papers, we named these glaciers from glacier ID 1 to glacier ID 7. The study presented the glacier changes in terms of displacement, velocity, and thickness. This study also covers predicted future possible glacier lake formation sites with their depth.

From the study, it was found that the glaciers are undergoing major changes in velocity in 1976–1978, and the effect of the change is decelerated in recent times. The speed of glaciers increases from the lower reach to the main central trunk for all the glaciers. Among the seven glaciers, glacier ID 5 is undergoing major changes in terms of speed. It is also noted the glaciers shrunk considerably in thickness between 2000 and 2022. Glacier ID 6 has undergone a maximum decline in its thickness with a reduction of −192.3 m between 2000 and 2022. From the study, it is observed that all the glaciers will form a glacier lake at its base. Glacier ID 6 is expected to form the largest lake with an area of 2.572 km2 and 208.5 m depth.

The next phase of the study is to study the variation of ELA and glacier mass balance. The annual ELA and annual glacier mass balance will give information on the health of the glaciers in the mountains. It is also recommended to do the study on glacier lake outburst flood risk assessment of the potentially dangerous glacier lake in Bhutan.

The authors are thankful to the Indian Science and Research Fellowship for sponsoring 3 months fellowship and the Indian Institute of Remote Sensing for hosting 3 months to do research. Their gratitude also goes to the National Centre for Hydrology and Meteorology, Bhutan, and the International Centre for Integrated Mountain Development (ICIMOD) for sharing glacier inventory, earth explorer for providing free Landsat data, and https://nsidc.org/apps/itslive/ for proving free dataset for global glacier velocity. The authors also thank organizations and individuals who developed and shared COSI-Corr and HIGTHIM software. The principal author also would like to sincerely thank Dr Praveen K. Thakur for his guidance throughout the research period.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.

Ageta
Y.
&
Iwata
S.
1998
The Assessment of Glacier Lake Outburst Flood (GLOF) in Bhutan. Report of Japan-Bhutan Joint Research
.
Bhushan
S.
,
Syed
T. H.
,
Kulkarni
A. V.
,
Gantayat
P.
&
Agarwal
V.
2017
Quantifying changes in the Gangotri Glacier of Central Himalaya: evidence for increasing mass loss and decreasing velocity
.
IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing
10
(
12
),
5295
5306
.
Dyurgerov
M. B.
&
Meier
M. F.
2000
Twentieth century climate change: evidence from small glaciers
.
Proceedings of the National Academy of Sciences
97
(
4
),
1406
1411
.
Gantayat
P.
,
Kulkarni
A. V.
&
Srinivasan
J.
2014
Estimation of ice thickness using surface velocities and slope: case study at Gangotri Glacier, India
.
Journal of Glaciology
60
(
220
),
277
282
.
Garg
V.
,
Thakur
P. K.
,
Rajak
D. R.
,
Aggarwal
S. P.
&
Kumar
P.
2022
Spatio-temporal changes in radar zones and ELA estimation of glaciers in NyÅlesund using Sentinel-1 SAR
.
Polar Science
31
,
100786
.
Gopika
J. S.
,
Kulkarni
A. V.
,
Prasad
V.
,
Srinivasalu
P.
&
Raman
A.
2021
Estimation of glacier stored water in the Bhaga basin using laminar flow and volume-area scaling methods
.
Remote Sensing Applications: Society and Environment
24
,
100656
.
Hasnain
S. I.
1999
Himalayan glaciers: hydrology and hydrochemistry. Allied Publishers, Mumbai
.
Hooker
B. L.
&
Fitzharris
B. B.
1999
The correlation between climatic parameters and the retreat and advance of Franz Josef Glacier, New Zealand
.
Global and Planetary Change
22
(
4
),
39
48
.
Iwata
S.
2010
Glaciers of Asia: Glaciers of Bhutan – An Overview. Satellite Image Atlas of Glaciers of the World, pp. F321–F334
.
Kulkarni
A. V.
,
Goswami
A.
,
Singh
G.
,
Srinivasalu
P.
,
Ramsanakaran
R.
&
Dashora
A.
2018
HIGHTHIM User Manual
.
Divecha Centre for Climate Change, Indian Institute of Science
,
Bangalore
,
India
.
Kumar
R.
,
Manzoor
S.
,
Vishwakarma
D. K.
,
Al-Ansari
N.
,
Kushwaha
N. L.
,
Elbeltagi
A.
&
Kuriqi
A.
2022
Assessment of climate change impact on snowmelt runoff in Himalayan region
.
Sustainability
24
(
3
),
1150
.
Lal
M.
2000
Climatic change-implications for India's water resources
.
Journal of Social and Economic Development
3
,
57
87
.
NCHM
2019
Bhutan Glacier Inventory 2018
.
Nick
F. M.
,
Van der Veen
C. J.
,
Vieli
A.
&
Benn
D. I.
2010
A physically based calving model applied to marine outlet glaciers and implications for the glacier dynamics
.
Journal of Glaciology
56
(
199
),
781
794
.
Paul
F.
,
Kääb
A.
&
Haeberli
W.
2007
Recent glacier changes in the Alps observed by satellite: consequences for future monitoring strategies
.
Global and Planetary Change
56
(
1–2
),
111
122
.
Racoviteanu
A. E.
,
Williams
M. W.
&
Barry
R. G.
2008
Optical remote sensing of glacier characteristics: a review with focus on the Himalaya
.
Sensors
8
(
5
),
3355
3383
.
Remya
S. N.
,
Kulkarni
A. V.
,
Pradeep
S.
&
Shrestha
D. G.
2019
Volume estimation of existing and potential glacier lakes, Sikkim Himalaya, India
.
Current Science
116
(
4
),
620
627
.
Shreve
R. L.
1985
Esker characteristics in terms of glacier physics, Katahdin esker system, Maine
.
Geological Society of America Bulletin
96
(
5
),
639
646
.
Singh
D. K.
,
Thakur
P. K.
,
Naithani
B. P.
&
Dhote
P. R.
2021
Spatio-temporal analysis of glacier surface velocity in Dhauliganga basin using geo-spatial techniques
.
Environmental Earth Sciences
80
,
1
16
.
Thakur
P. K.
,
Ambika
A. K.
,
Bisht
S. M.
,
Stein
A.
,
Mahagaonkar
A.
,
Kumar
U.
,
Garg
V.
,
Khajuria
V.
,
Chouksey
A.
,
Snehmani
,
Chauhan
P.
&
Aggarwal
L. S.
2023
Gangotri glacier dynamics from multi-sensor SAR and optical data
.
Advances in Space Research
72
(
2
),
309
326
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).