1. Introduction
Glaciers in High Mountain Asia (HMA) contain ~98 000 km2 of ice (Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022), making them the largest ice mass outside the polar regions These glaciers have been losing mass in recent decades at a rate of −19.0 ± 2.5 Gt a−1 (−0.19 ± 0.03 m w.e. a−1) between 2000 and 2018 (Shean and others, Reference Shean2020), and this ice loss is forecast to continue during the 21st century (Rounce and others, Reference Rounce, Hock and Shean2020). Glacier mass balance has been particularly negative across central and eastern regions of the Himalaya and has coincided with rising air temperatures (e.g. Bolch and others, Reference Bolch2012; Brun and others, Reference Brun, Berthier, Wagnon, Kääb and Treichler2017; Maurer and others, Reference Maurer, Schaefer, Rupper and Corley2019b; Shean and others, Reference Shean2020; Hugonnet and others, Reference Hugonnet2021). Quantifying Himalayan glacier loss is crucial, as these glaciers are the source of Asia's largest river systems and ~220 million people depend on glacier melt to maintain water supplies throughout the dry season (e.g. Immerzeel and others, Reference Immerzeel, van Beek and Bierkens2010; Pritchard, Reference Pritchard2019; Immerzeel and others, Reference Immerzeel2020; Nie and others, Reference Nie2021). Consequently, downstream populations are highly vulnerable to climate change-induced changes in glacial meltwater provision. Furthermore, populations in HMA are the most vulnerable globally to Glacial Lake Outburst Floods (GLOFs) (Carrivick and Tweed, Reference Carrivick and Tweed2016; Taylor and others, Reference Taylor, Robinson, Dunning, Carr and Westoby2023), which occur when water impounded behind a natural dam (e.g. a moraine) is suddenly released (e.g. Fujita and others, Reference Fujita2013; Westoby and others, Reference Westoby2015). GLOFs can rapidly generate large volumes of water and sediment, which can cause major damage to downstream infrastructure, property and agricultural land, and loss of life (e.g. Carrivick and Tweed, Reference Carrivick and Tweed2016; Huggel and others, Reference Huggel, Dimri, Bookhagen, Stoffel and Yasunari2020; Zheng and others, Reference Zheng2021b; Emmer and others, Reference Emmer2022). Rapid ice loss has led to a dramatic increase in the number and volume of glacial lakes globally in recent years (Shugar and others, Reference Shugar2020), and this has been particularly marked in the eastern Himalaya (e.g. Komori, Reference Komori2008; Wang and others, Reference Wang, Xiang, Gao, Lu and Yao2015; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016; Zheng and others, Reference Zheng2021a). Given the potential impacts on water supplies and natural hazards, there is therefore an urgent need to quantify recent changes in glacial water storage in the eastern Himalaya.
GLOFs can occur from water impounded at the front/margins of glaciers or on their surface, and can be triggered by a variety of mechanisms, including: (i) failure of the dam impounding the lake; (ii) drainage of supraglacial ponds supra- or englacially; (iii) inputs from surrounding hillslopes and/or glaciers; (v) ice calving and/or rapid meltwater inputs; (v) earthquakes; (vi) intense rainfall and/or high temperatures (e.g. Richardson and Reynolds, Reference Richardson and Reynolds2000; Westoby and others, Reference Westoby2014; Westoby and others, Reference Westoby2015; Rounce and others, Reference Rounce, McKinney, Lala, Byers and Watson2016; Emmer and others, Reference Emmer2022; Shrestha and others, Reference Shrestha2023). Supraglacial ponds also enhance ablation, by absorbing solar radiation and transmitting it to exposed ice (e.g. Quincey and others, Reference Quincey2007; Komori, Reference Komori2008; Carrivick and Tweed, Reference Carrivick and Tweed2013; Narama and others, Reference Narama2017). Once formed, ponds expand via enhanced ablation and calving (Sakai and others, Reference Sakai, Takeuchi, Fujita and Nakawo2000, Reference Sakai, Nishimura, Kadota and Takeuchi2009; Benn and others, Reference Benn, Wiseman and Hands2001; Miles and others, Reference Miles2016, Reference Miles2017b) and positive feedbacks develop: increased ablation enhances surface lowering, which reduces the ice surface gradient and promotes further ponding (Benn and others, Reference Benn2012; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016; Narama and others, Reference Narama2017). These positive feedbacks can increase the size and number of ponds, which may coalesce to form large proglacial lakes (Richardson and Reynolds, Reference Richardson and Reynolds2000; Miles and others, Reference Miles2016; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016). Supraglacial ponds can also drain on seasonal to interannual timescales, when they intersect with the englacial hydrological system and/or overtop (Benn and others, Reference Benn, Wiseman and Hands2001; Wessels and others, Reference Wessels, Kargel and Kieffer2002; Miles and others, Reference Miles, Willis, Arnold, Steiner and Pellicciotti2017a; Taylor and others, Reference Taylor, Carr and Rounce2022). Ice surface slope is the primary control on supraglacial pond distribution and characteristics: shallower slopes promote widespread ponding, coalescence and lake formation, whereas steeper slopes lead to more isolated ponds (Reynolds, Reference Reynolds2000). Ice velocities also influence pond size and distribution, as faster flow encourages fractures to open, leading to pond drainage and hence limiting pond expansion and persistence (Salerno and others, Reference Salerno2012; Miles and others, Reference Miles2016, Reference Miles2017b).
Previous work attributes 11% (72 total) of GLOFs recorded in the HMA between 1833 and 2022 with an identifiable cause to supraglacial lake drainage events, compared to 49% (330 total) from moraine damned lakes and 32% (234 total) from ice dammed (Shrestha and others, Reference Shrestha2023). However, recurring GLOFs were often supraglacial in origin and almost all casualties were associated with failure of supraglacial or moraine dammed lakes, with the former accounting for 101 deaths from 72 recorded incidents (Shrestha and others, Reference Shrestha2023). Furthermore, at least four of the 24 historical GLOFs observed in Bhutan originated from supraglacial lakes (Komori and others, Reference Komori, Koike, Yamanokuchi and Tshering2012). Jointly with Nepal, Bhutan is the most vulnerable country globally to GLOFs (Liu and others, Reference Liu, Mayer and Liu2015; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016, Reference Watson2017; Miles and others, Reference Miles, Willis, Arnold, Steiner and Pellicciotti2017a), as its glaciers are shrinking rapidly (RGoB, 2011; Carrivick and Tweed, Reference Carrivick and Tweed2016) and its proglacial lakes are expanding (Brun and others, Reference Brun, Berthier, Wagnon, Kääb and Treichler2017; Maurer and others, Reference Maurer, Schaefer, Rupper and Corley2019b; Shean and others, Reference Shean2020). Furthermore, its main population centres, agricultural land, cultural and administrative centres lay downstream of large, rapidly growing lakes (Dorji and others, Reference Dorji, Olesen, Bøcher and Seidenkrantz2016). Bhutan also relies on Hydroelectric Power (HEP) to generate 100% of its power and ~45% of its GDP (UNDPBhutan, 2011), while subsistence agriculture employs ~79% of its population (Dorji and others, Reference Dorji, Olesen, Bøcher and Seidenkrantz2016).
To date, most work on GLOFs in Bhutan has focused on the Lunana region, which has numerous large, growing glacial lakes that have generated multiple historical GLOFs (e.g. Richardson and Reynolds, Reference Richardson and Reynolds2000; Fujita and others, Reference Fujita, Suzuki, Nuimura and Sakai2008; Narama and others, Reference Narama2017; Maurer and others, Reference Maurer2020; Rinzin and others, Reference Rinzin2023). However, prior work has highlighted the nearby Tshojo Glacier as potentially dangerous (Fig. 1): in 2009, a previously undetected ice surface pond burst and caused the downstream city of Punakha to be evacuated, although no damage was recorded (Komori and Tshering, Reference Komori and Tshering2010; Komori and others, Reference Komori, Koike, Yamanokuchi and Tshering2012; Yamanokuchi and others, Reference Yamanokuchi, Tadono, Komori, Kawamoto and Tomiyama2011). Tshojo Glacier is located on the southern side of the main Himalayan ridge in the Lunana region and flows southwards (28°6′21″N, 90°9′52″E; Fig. 1). Tshojo Glacier's tongue has extensive debris cover and is situated between ~4000 and 4500 m a.s.l (Fig. 1), while the glacier as a whole is ~15 km long, making it one of the largest glaciers in Bhutan (Komori and others, Reference Komori, Koike, Yamanokuchi and Tshering2012). Meltwater from Tshojo Glacier is main source of the Gotey chhu, one of main tributaries of Phochu River, and the glacier is located ~2.5 km upstream of the Lunana Gewog centre, which is the administrative centre for the Lunana region. There are two major hydropower plants located downstream of Tshojo Glacier: the Basochhu power plants I (capacity = 24 MW) and II (capacity = 40 MW) are located approximately 95 km downstream and the Punatsangchhu-I (capacity = 1200 MW) and Punatsangchhu-II (capacity = 1020 MW) plants, located ~85 km downstream. Thus, there is a pressing need to quantify the variability in water volumes stored in Tshojo Glacier's ice surface ponds, given its location directly upstream major settlements, cultural heritage sites and HEP plants. In this study, we aim to: (i) quantify the multi-decadal evolution and interannual variability of supraglacial water storage on Tshojo Glacier between 1987 and 2020; (ii) to determine controls on observed supraglacial water storage patters, including ice velocities, longitudinal gradient of the glacier surface and air temperatures; and (iii) determine the area and infrastructure downstream that would be inundated for different pond drainage scenarios.
(a) Overview map, showing the location of Tshojo Glacier within Bhutan and the surrounding region of the Himalaya. Inset imagery is from the SRTM DEM and is colour coded by elevation. Country boundaries are indicated in yellow and the red box indicates the extent of ‘b’. (b) location of Tshojo Glacier in relation to the Lunana area and other major glaciers and proglacial lakes. The glacier outline is indicated in light blue and the glacier tongue in dark blue. Both outlines were manually digitised from the background image, which is from Landsat 8, acquired on 16th Dec 2019, from USGS Earth Explorer.
2. Methods
2.1. Imagery sources
We used multiple image sources to quantify multidecadal patterns in supraglacial pond volumes and locations on Tshojo Glacier. For each time period, we used the highest resolution data that were consistently available and ensured overlap between datasets to evaluate the impacts of image resolution on pond statistics. Specifically, for the period 2017–2020, we used 4-band PlanetScope scenes at 3 m resolution and for 2012–2017 we used RapidEye Ortho tiles (corrected for surface reflectance) at 5 m resolution (Table S1). Between 1987 and 2020, we utilised imagery from the Landsat Collection 2 Level-1 data, provided by the USGS Earth Explorer, specifically Landsat 4–5 TM, Landsat 7 ETM+, and Landsat 8 OLI (Table S1). To extend the data further back in time (1962–1976) we used declassified spy satellite data supplied by USGS Earth Explorer, specifically the Corona KH-4 imagery and Hexagon KH-9 Lower Resolution Mapping and 70 mm Panoramic film datasets (Table S1). Unfortunately, georeferencing could not be conducted accurately enough to reliably quantify pond areas, due to both the lack of camera information for the Corona imagery and the changes on the steep slopes surrounding Tshojo Glacier in recent decades that would otherwise be used for tie points (e.g. ridge edges). Consequently, the spy satellite imagery was used only to determine the longevity of Tshojo Glacier's large ponds and not to calculate pond area.
To improve comparability between the datasets of different spatial resolution, and because errors in pond area scale inversely with pond area i.e. errors are proportionally higher for smaller ponds (Salerno and others, Reference Salerno2012; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016; Taylor and others, Reference Taylor, Carr and Rounce2022), we exclude ponds less than 1000 m2 in area. This threshold was chosen to minimise differences associated with image resolution, while including as much of the ponded area as possible, specifically 73% of the total ponded area for RapidEye, 81% for Planet Scope and 90% for Landsat. We evaluated differences in pond statistics (Table S2) and the spatial distribution of ponds (Fig. S2) for Planet Scope (3 m resolution) and RapidEye (5 m resolution) data using images from the same date (13 December 2017), to determine comparability between these commonly used sensors.
To minimise fluctuations in pond area due to seasonal cycles, we aimed to select images in early to mid-December (Table S1 & Fig. S1). This date was chosen to avoid monsoon cloud cover (May – August) and heavy winter snowfall (January – March), which can obscure frozen ponds. Previous work from nearby glaciers (Taylor and others, Reference Taylor, Carr and Rounce2022) has demonstrated that maximum ponded area occurs in either winter or during the monsoon. Furthermore, Taylor and others, Reference Taylor, Carr and Rounce2022 showed that pond drainage events occurred more frequently during the monsoon (average: 8 times/month) than the winter (average: 3 times/month), although cloud cover during the monsoon adds uncertainty to these seasonal patterns (Taylor and others, Reference Taylor, Carr and Rounce2022). Thus, the limited data available on seasonal variations in pond area in Bhutan suggest that ponds are likely to be comparatively full and close to their maximum seasonal extent in December. For the Planet and RapidEye imagery, all dates were in December. For Landsat 8, three image dates were in November (from 14 November onwards), for Landsat 7, one image was acquired on 26 November and three images were acquired in Nov for Landsat 4–5 (Table S1 & Fig. S1). Scenes were only used if Tshojo Glacier was fully visible and unobscured by cloud. Where required, individual scenes were mosaicked and clipped to the study area extent to generate a single image of Tshojo Glacier and its tributaries for each image date.
2.2. Quantification of pond area
To identify supraglacial ponds on Tshojo Glacier from PlanetScope imagery (2017–2020) and Landsat data (1987–2020) we used Normalised Difference Water Index (NDWI) (McFeeters, Reference McFeeters1996). The workflow is summarised in Fig. S3. Water surfaces have high relectance in the green band and low reflectance in the near-infrared band (Fig. S3a). NDWI uses a ratio to magnify this difference in reflectance to identify water from other surface types (i.e. ponds from supraglacial debris cover) by calculating ${{\Delta \;\;Reflectance} \over {\sum \;\;Reflectance}}$
(Lillesand, Reference Lillesand2008; Bolch and others, Reference Bolch2011; Jha and Khare, Reference Jha and Khare2017; Miles and others, Reference Miles, Willis, Arnold, Steiner and Pellicciotti2017a), as follows:
For the PlanetScope images we used Band 2 (Green: 0.50–0.59 μm) and Band 4 (Near-Infrared [NIR]: 0.78–0.86 μm; Table S1). The same bands were used for Landsat TM 4–5 (Band 2: 0.52–0.60 μm; Band 4: 0.76–0.90 μm) and Landsat ETM + 7 (Band 2: 0.52–0.60 μm; Band 4: 0.77–0.90 μm; Table S1). For the Landsat 8 data, we utilised Bands 3 and 5, which equate to the Green and Near-Infrared bands, at wavelengths of 0.525–0.600 μm and 0.845–0.885 μm respectively (Table S1).
After calculating the NDWI (Fig. S3b), a threshold of water bodies vs nonwater needs to be set. Here, we made an initial estimate of the NDWI threshold by determining the NDWI value of ponds that were clearly identifiable in each image. Following previous studies (i.e. Bolch and others, Reference Bolch2011; Miles and others, Reference Miles, Willis, Arnold, Steiner and Pellicciotti2017a), the initial threshold was then adjusted manually for each image, to find the optimum NDWI threshold that most accurately identified ponds while minimising misclassification errors (e.g. exposed ice cliffs and crevasses incorrectly classified as water). A final threshold of 0.07 was applied to all PlanetScope NDWI raster outputs to produce classified maps of ponds vs no ponds (Fig. S3c), and areas with NDWI values > 0.07 were exported to create polygons for each individual pond (Fig. S3e). For the Landsat data, the most appropriate NDWI threshold varied between Landsat satellites and between images, and so the threshold used for each image is provided in Table S1. Initially, we used NDWI to classify the RapidEye imagery, but the results were poor when manually compared with the imagery, and so the ponds were manually digitised. To estimate uncertainties due to the NDWI classification, we compared our PlanetScope and Landsat derived NDWI results to the manually digitised Rapid Eye data. Specifically, we compared the total ponded area obtained for an individual image date (13 December 2017), for the PlanetScope NDWI results (205 605 m2) and the manually digitised RapidEye data (206 160 m2), which resulted in a percentage difference of 0.27%. We then compared total ponded area from the Landsat NDWI classification from 10th December 2017 (155 583 m2) with all ponds from the manually digitised RapidEye data that intersected with the Landsat ponds (171 900 m2): this was necessary to facilitate direct comparison because the coarser resolution of Landsat means that only larger ponds are detected. The percentage difference in total ponded area between the manually-digitised data and the Landsat NDWI was 10.5%.
2.3. Ice velocities, longitudinal gradient of the glacier surface and thinning rates
Ice surface velocities were determined from the Inter-Mission Time Series of Land Ice Velocity and Elevation (ITS_LIVE) project (Gardner and others, Reference Gardner2018), provided publicly as part of NASA's MEaSUREs 2017 program. The data provide regional-scale maps of annual ice velocities for the world's major glaciated regions, and we use the HMA data product. Velocity maps were produced annually, from 1985 to 2021, and we used data from 1988 to 2020: the first three years of data (1985–1987) had large coverage gaps over Tshojo Glacier and 2020 corresponded to our most recent pond data. In our study, we utilised individual data years, which had a spatial resolution of 240 m, and the 120 m resolution composite, which provides the average velocity across all years.
The longitudinal gradient of the glacier surface was calculated from the HMA 8-metre DEM Mosaics Derived from Optical Imagery, Version 1 (Shean, Reference Shean2017). The dataset was generated from very high-resolution satellite imagery (GEO-EYE, QUICKBIRD and the WORLDVIEW satellites) and has a final spatial resolution of 8 m. The dataset was clipped to Tshojo Glacier's outline, which was determined manually from PlanetScope imagery. Following approaches developed in (Quincey and others, Reference Quincey2007; Miles and others, Reference Miles, Willis, Arnold, Steiner and Pellicciotti2017a), we calculated the longitudinal gradient of the glacier surface, as opposed to surface slope, to remove the impacts of the hummocky topography often found on debris covered glaciers (i.e. undulations at the scale of tens – hundreds m). Specifically, we calculated the longitudinal surface gradient as the angle from horizontal of each point on Tshojo Glacier's surface, relative to a base contour (Fig. S4). The base contour (4060 m) was selected as the first ice surface contour up glacier from the terminus that traversed the entire glacier tongue (Fig. S4). We first calculated the vertical distance, or difference in elevation, between the base contour and all points on the glacier surface, by subtracting 4060 m from the HMA DEM data (Fig. S4A). We then determined the horizontal distance of all points on the glacier from the base contour by calculating the Euclidean distance. Next, we calculated the longitudinal surface gradient in degrees using:
The longitudinal surface gradient data was divided into four categories, which are associated with differing pond characteristics on Bhutanese glaciers, as defined by Reynolds (Reference Reynolds2000).
Ice surface thinning rates were determined from the HMA Gridded Glacier Thickness Change from Multi-Sensor DEMs, Version 1 (Maurer and others, Reference Maurer, Rupper and Schaefer2019a). The dataset provides thickness change mosaics for selected Himalayan glaciers for the periods 1975–2000 and 2000–2016 and has a spatial resolution of 30 m. The thickness change data were created by combining digital elevation models (DEMs) derived from HEXAGON KH-9 and ASTER DEMs. The data were clipped to the manually-digitised outline of Tshojo Glacier. Changes in long-profile glacier thickness were calculated from Pléiades DSMs, generated using the Ames Stereo Pipeline (Beyer and others, Reference Beyer, Alexandrov and McMichael2018) and made available via the Pléiades Glacier Observatory (Berthier and others, Reference Berthier2024). DSMs were available for Tshojo Glacier from 7 November 2017 and 19 October 2022. We utilised the 20 m resolution product to extract long profiles, to remove small-scale undulations on the glacier surface, although visual comparison of the 2 and 20 m products showed a very similar overall profile. We used the product generated using the Semi-Global Matching (SGM) correlator (Hirschmuller, Reference Hirschmuller2007), as it is able to highly-resolve topography and minimises data gaps; no areas of noise were apparent upon visual inspection of the data (Berthier and others, Reference Berthier2024). Ice surface elevations were extracted along the glacier centreline from the 2017 and 2022 DSMs, at an interval of 20 m, to correspond with the data resolution.
2.4. Pond outburst flood scenarios
To determine the downstream impact from potential future scenarios of outbursts from Tshojo Glacier's supraglacial ponds, HEC-RAS 2D model (v6.4.1) with 2D shallow water equation was used to perform unsteady flow routing of 11 pond outburst scenarios from Tshojo Glacier (CEIWR-HEC, 2021). The scenarios were based on the range of individual pond areas observed during the study period, from 10 000 to 100 000 m2, with an increment of 10 000 m2, which corresponds to our smallest (scen1) to largest magnitude scenarios (scen10; Table 1). Detailed data from the 2009 event demonstrated full pond drainage occurred suddenly (Komori and Tshering, Reference Komori and Tshering2010; Komori and others, Reference Komori, Koike, Yamanokuchi and Tshering2012; Yamanokuchi and others, Reference Yamanokuchi, Tadono, Komori, Kawamoto and Tomiyama2011), meaning that Scenarios 1 to 10 could occur approximately annually, depending on the exact cycle of the largest ponds. We also include a worst-case scenario: this uses the total area of all ponds visible on the 2008 image, which was the largest total ponded area observed on Tshojo Glacier during the study period. It represents the ‘worst-case’ scenario, whereby ponded area is at its observed maximum and all ponds drain at once, which represents an unlikely, but plausible, end-member scenario potentially triggered by a rare event, such as a catastrophic mass movement event, as observed in Chamoli (Shugar and others, Reference Shugar2021) and/or sudden release of water following blocking of a subglacial conduit, as observed on the Changri Shar and Khumbu glaciers (Miles and others, Reference Miles2018). For each scenario, we used a different flow hydrograph, which is the function of pond volume for that scenario as an upstream boundary condition, while keeping other parameters constant (Table 2). To estimate pond volumes from observed pond areas, we applied an empirical area-volume relationship that was developed based on data from supraglacial ponds on debris-covered glaciers in the Himalaya (Watson and others, Reference Watson, Quincey, Carrivick and Smith2016). For the worst-case scenario, we converted the area of each pond detect to volume using the same empirical equation and then summed the estimated volumes to get a total volume for 2008. Specifically:
where:
Area of supraglacial pond and corresponding empirically derived flood volume and peak discharge (Qp) calculated from Eqn (4) for each scenario considered in the study and then inputted into the HEC-RAS 2-D hydrodynamic model.
Amount of agricultural land, footpaths, buildings and bridges exposed to our range of scenarios of supraglacial pond outburst flood from Tshojo Glacier. Pond area and volume for each scenario are given in Table 1.
Based on observations of previous supraglacial pond drainage on Tshojo Glacier (Yamanokuchi and others, Reference Yamanokuchi, Tadono, Komori, Kawamoto and Tomiyama2011), and in the neighbouring central Himaya (Kropáček and others, Reference Kropáček2015), we assumed that each pond drains fully in each scenario. Thus, the estimates presented are upper bounds for potential flood volume and resultant downstream extent.
2.5. Pond outburst model parameters
The model domain was established by creating a 500 m buffer zone on both sides of the centreline of the Phochu River, downstream of Tshojo Glacier's terminus (Fig. 1). Within this model domain, a computational mesh with a grid resolution of 30 × 30 m was generated. ALOS-PALSAR DEM, which is a SRTM 30 GL1 resampled to a ground resolution of 12.5 m, was used to derive the terrain information for the model setup. The DEM data was hydrologically corrected by filling the sinks and burning the artificial channel in the areas where deep gorges are detected as sinks, as per the approach used by Rinzin and others (Reference Rinzin2023). The HMA 8-metre DEM used to determine Tshojo Glacier's longitudinal surface gradient could not be used for the GLOF modelling across the entire domain, as its coverage is limited and patchy. However, HMA DEM data were available within ~10 km of Tshojo Glacier's snout and we therefore used this smaller domain to evaluate the sensitivity of our results to DEM resolution. Specifically, we ran Scenario 10 (i.e. pond area = 100 000 m2) with both the ALOS-PALSAR DEM and HMA 8 m DEM and compared their resulting flow hydrodynamics.
The other essential model input parameters are boundary conditions and the Manning's roughness coefficient n. The inflow hydrograph was imposed as upstream boundary condition, while the downstream boundary condition at the outlet was assumed normal depth (i.e. average slope at the downstream boundary condition) with the discreet energy slope of 0.01/The latter was determined by taking the bed slope across a profile length of 500 m, perpendicular to the downstream boundary condition. The supraglacial pond volume calculated using Eqn (3) was used as the total flood volume for the inflow hydrograph construction (Table 2). The initial condition for the 2D flow area was established by setting 2D initial condition ramp-up time of 10 hours out of which 10% was allocated for ramping up the 2D boundary conditions up from zero to their first value (CEIWR-HEC, 2021).
The Manning's n value represents the resistance to flow from channel and flood plain and plays an influential role in the propagation of flow downstream and resulting flow dynamics (CEIWR-HEC, 2021). A spatially distributed pixel-based Manning's value was assigned based on land cover data from the International Centre for Integrated Mountain Development (ICIMOD) (Uddin, Reference Uddin2021). In certain areas, the river channel was not well resolved by the landcover map, so we manually mapped the river channel using high resolution, publicly available Google Earth imagery, which was from Pleiades Neo 4, with a date of 27 December 2023 and a horizontal resolution of 30 cm, and assigned a Manning's value of 0.035 to this manually mapped channel (Arcement and Schneider, Reference Arcement and Schneider1989). The Manning's n values defined from the landcover dataset and Google Earth ranged between 0.02 and 0.12. Due to the potentially importance influence of this value on our results, we conducted further sensitivity analysis. As for our DEM sensitivity assessments, we utilised scenario 10 and varied the value for Manning's n between our calculated range of 0.02 and 0.12, in increments of 0.01. We then compared flood hydrodynamic data (flow arrival time and peak flow) at three sites, located at distances of 5, 10 and 20 km downstream of Tshojo Glacier respectively.
The water released from supraglacial lakes can drain either into existing sub- or englacial conduits, and eventually emerges at the glacier snout (Richardson and Reynolds, Reference Richardson and Reynolds2000), or can flow supraglacially, if it overwhelms the sub- and englacial drainage networks (Miles and others, Reference Miles2018). The exact routing depends on a variety of factors, including the glacier geometry near the terminus and the configuration of the drainage network, which are very difficult to quantify in the absence of directly measured data. Thus, we assumed the worst case scenario (i.e. the most rapid water delivery to the snout) and calculated peak discharge using the empirical equation of tunnel drainage events proposed by Clague and Mathews (Reference Clague and Mathews2017) (Eqn (4)). Due to the difficulties involved in collecting essential parameters such as drainage tunnel size, temperature or filling level of the lake, using a numerical approach to generate the flow hydrograph was not possible (Kropáček and others, Reference Kropáček2015). Instead, we estimated the hydrograph by fitting the peak flow and volume of each scenario to a log-normal distribution curve with a sigma value of 4.5 and mean of 1 following the convention used by (Kropáček and others, Reference Kropáček2015).
where:
2.6. Pond outburst impact assessment
To estimate the extent of flood inundation downstream caused by an outburst from Tshojo Glacier's supraglacial ponds, we collated flow depth from all scenarios and assigned a consistent value of 1. The final flood inundation map was produced by summing the resulting inundation extents from all 11 scenarios. This map depicted the range of pixel values from 1 to 11, corresponding to the number of scenarios that affected a specific area. For instance, a pixel with a value of 1 indicated that the area was impacted solely from one scenario, while a pixel with a value of 11 indicated impact from all 11 scenarios.
To assess the exposure of downstream infrastructure and settlements to our modelled potential outburst floods from Tshojo Glacier, we extracted the buildings, roads, and bridges that intersect with our flood inundation maps from the OpenStreetMap. We also overlaid these data on Google Earth imagery from 27 January 2021 (i.e. the closest date to the Open Street map data) and manually checked its accuracy. Fortunately, OpenStreetMap had good coverage of buildings and infrastructure within the inundation map, except for the two suspension bridges in the Lunana region, which we added manually. However, the existing OpenStreetMap and land-use and land cover maps at the global and Hindu Kush Himalaya (HKH) scales (Uddin, Reference Uddin2021; Brown and others, Reference Brown2022) do not adequately resolve the agricultural land, so we manually mapped the agricultural land within our model domain using Google Earth imagery.
3. Results
3.1. Comparison of pond statistics between PlanetScope, RapidEye and Landsat imagery
We observed limited difference between for total ponded area (0.27% difference), area of the largest pond (0.76%), and percentage of the glacier tongue covered by ponds (0.54%; Table S2 & Fig. S2), and so we treat these statistics as directly comparable between the two datasets. Differences were somewhat higher for the number of ponds (7.23%) and the mean area across all ponds (6.69%; Table S2) but remained small enough to be broadly comparable. Landsat data were not available on the same dates as PlanetScope or RapidEye and so we compared images within three days of each other, to minimise any temporal variability. For PlanetScope and Landsat imagery, the total ponded area (0.94%), area of the largest pond (4.63%) and percentage of the glacier tongue covered by ponds (0.94%) were thus directly comparable (Table S3, Fig. S2). However, differences in the number of ponds (41.13%) and mean pond area (42.02%) were large, so are not directly compared in this study (Table S3). For RapidEye and Landsat data, the total ponded area (14.95%), area of the largest pond (10.63%) and percentage of the glacier tongue covered by ponds (14.95%) were broadly comparable (Table S4, Fig. S2), but the number of ponds (29.65%) and the mean pond area (40.09%) were not (Table S4). We attribute these differences in mean pond area and number of ponds for Landsat vs PlanetScope/RapidEye to the difference in image resolution: Landsat does not detect the substantial number of smaller ponds (Fig. S2), leading to a lower count and a higher mean pond size, but total ponded area is largely unaffected by resolution, as it is dominated by the larger ponds.
3.2. Temporal variability in bulk pond statistics
We used Landsat data to investigate temporal variability in total pond area, percentage of the glacier tongue covered by ponds and the area of the largest pond, because values for these variables were comparable across the different image resolutions (Tables S3 & S4) and Landsat enables us to construct the longest temporal record. Total ponded area on Tshojo Glacier showed two distinct periods during our study period: it doubled between the earlier portion (1987–2003) and the latter portion (2007–2020), increasing from 104 529 m2 to 213 943 m2 (Fig. 2a). Interannual variability in total ponded area also increased markedly between these two time periods: the largest change between consecutive years in the earlier period was 85 500 m2, between 2001 and 2002, compared to 165 600 m2 between 2012 and 2013 (Fig. 2a). Over the entire study period, total ponded area varied by approximately an order of magnitude, from 34 200 m2 in 1997 to 331 200 m2 in 2008 (Fig. 2a). The percentage of the glacier tongue covered by ponds was derived from total pond area, so follows the same temporal pattern, and ranges between 0.23 and 2.2% of Tshojo Glacier's area (Fig. 2b). As for total ponded area, the area of the largest pond increased markedly between 1987 and 2003 (average area = 16 457 m2) and 2007–2020 (average area = 60 879 m2; Fig. 2c). Likewise, the interannual variability increased by an order of magnitude between these two time intervals, with the largest year-on-year change during the earlier period being 8100 m2 (1997–1998), compared to 83 700 m2 during the later period (2014–2015; Fig. 2c). Finally, we observed no apparent up-glacier migration in pond locations of any size at multi-decadal time scales (Fig. S5).
Variability in (a) Total area of all ponds on Tshojo Glacier; (b) the percentage of Tshojo Glacier's tongue covered by ponds; and c) the area of the largest pond detected in a given year on Tshojo Glacier, for the period 1987 to 2020, derived from Landsat imagery (30 m resolution). Error bars (lighter colours) show the 10.5% error calculated for the Landsat data (Section 2.2.). Values are given for ponds over 1000 m2 only. Data are from Landsat 4–5 TM (1987–2003), Landsat 7 ETM + (2007–2012) and Landsat 8 (2013–2020) and specific image dates and details are given in Table S1.
To assess variability in the number of ponds and mean pond area for ponds above the 1000 m2 threshold, we use data obtained from PlanetScope and RapidEye, as these variables are sensitive to image resolution (Tables S3 & S4) and so the highest image resolution available was used. Between 2012 and 2020, the number of ponds showed a net reduction over time (Fig. 3a). The number of ponds showed considerable interannual variability between 2012 and 2017, ranging between 51 in 2013 and 35 in 2017 (Fig. 3a). Pond number was lowest in 2020, at 30 ponds (Fig. 3a). Mean pond area increased overall between 2012 and 2020 (Fig. 3b) but was temporarily variable: mean pond area increased from 2837 m2 in 2012 to 5256 m2 in 2015, before decreasing to 3727 m2 in 2016. Thereafter, mean pond area increased again to 5230 m2 in 2018 and remained at similar values until the end of the dataset in 2020 (Fig. 3b).
Variability in (a) the number of ponds and (b) mean area of the ponds on Tshojo Glacier between 2012 and 2020. Values are given for ponds over 1000 m2 only. Magenta indicates data derived from RapidEye imagery (5 m spatial resolution) and dark blue incites PlanetScope (3 m resolution). Specific image dates and details are given in Table S1. Error in the PlanetScope data, relative to our manually digitised reference ponds from RapidEye, was 0.27%.
3.3. Temporal variability in large ponds
Our data demonstrate that both the number and area of large ponds increased dramatically between 1987–2003 and 2007–2020 (Fig. 4). Specifically, from 2007 onwards, we observed seven ponds over 50 000 m2, two between 500 000 and 40 000 m2 and nine between 40 000 and 30 000 m2, whereas no ponds were found in these categories between 1987 and 2003 (Fig. 4). Furthermore, ponds with an area between 20 000 and 30 000 m2 increased from four in 1987–2003 to 16 in 2007–2020 and those between 10 000 and 20 000 m2 increased from 21 to 29 for the same time intervals (Fig. 4). Finally, the number of ponds between 1000 and 10 000 m2 increased from 248 to 343 (Fig. 4).
Location and area of large ponds detected on Tshojo Glacier between 1987 and 2020. (a) Location of ponds on Tshojo Glacier. Ponds larger than 10 000 m2 are symbolised by size and colour, according to pond area category. Ponds smaller than 10 000 m2 are symbolised by a cyan dot. Ponds larger than 20 000 m2 are named. (b) Individual pond area by year. For ponds larger than 20 000 m2, symbol type and colour indicate the pond letter, located in (a). Note that B and D indicate separate ponds and BD denotes incidences when B and D merge. Ponds smaller than 20 000 m2, are symbolised with a blue dot (10 000–20 000 m2) or a black dot (less than 10 000 m2). (c) violin plots of pond area for 1987–2003 and 2007–2020, with the mean indicated by the black line and the median by the do-dashed red line.
Focusing on individual ponds, the largest (117 000 m2 in 2008) and the second largest (116 100 m2 in 2015) individual pond area observed during the time series occurred when Pond B and D merged and the third largest was Pond A, reaching 100 800 m2 in 2013 (Fig. 4; Table S5). The area of an individual pond exceeded 50 000 m2 at three different ponds during period: three times at BD, twice at A and once at V (Fig. 4; Table S5). Pond A reached over 40 000 m2 twice during the study period, while the ponds exceeding 30 000 m2 were: B three times, A and V twice, and C and H once (Fig. 4; Table S5). Ponds that exceeded over 20 000 m2 during the study period were: A five times, Y four times, J twice and Ponds B-F, H, M, S and V once (Fig. 4; Table S5). Thus, our data indicate that the largest ponds consistently form at the same locations during our study period. This is supported by data from the declassified satellite imagery, in which we identified most of our large ponds (A, B, D, F, H, M, V and Y) in images from 1967, 1974 and 1976 (Fig. S6). This demonstrates that some of the large ponds have formed in the same location on Tshojo Glacier for at least 53 years. Of the ponds over 20 000 m2 (Fig. 1), we did not observe ponds C, E, J or S in the declassified imagery (Fig. S6). We could not definitively identify any ponds in declassified imagery from 1966 or 1962 (Fig. S6).
3.4. Ice velocities, longitudinal gradient of the glacier surface and thickness change
Overall, Tshojo Glacier's ice velocities increased with distance from the terminus and the glacier margins (Fig. 6a). Velocities were low (≤5 m a−1) over the lowest 4.5 km of the glacier tongue and increased to over 50 m a−1 at ~12 km from the terminus (Fig. 6a & c). Tshojo Glacier slowed down between 1998 and 1993, particularly in its middle section (~4–10 km from the terminus), where ice velocities reduced by ~10 ma−1 (Fig. 6c). Smaller magnitude reductions in ice velocities were also observed over the lower 4 km and the upper ~1 km of the glacier tongue (Fig. 6c). Ice velocities remained low in 2003, before increasing slightly in the middle section of the glacier tongue (~4–10 km) for the remainder of the study period (Fig. 5c).
Ice velocities, longitudinal gradient of the glacier surface, and supraglacial pond distribution on Tshojo Glacier. (a) Map of ice velocities from ITS LIVE 120 m composite, compiled from velocity data from 1988–2020. (b) Map of the longitudinal gradient of the glacier surface, derived from High Mountain Asia 8 m resolution DEM (data date: 16 July 2017). Panels a & b show the location and area of all ponds on Tshojo Glacier detected during the study period (1987–2020). Ponds larger than 10 000 m2 are categorised and symbolised by size and colour. Ponds smaller than 10 000 m2 are symbolised by a cyan dot. Ponds larger than 20 000 m2 are named. In panels a & b the thin black line represents the glacier centreline, with distance markers (black circle) numbered by distance in kilometres from the glacier terminus, along the centre line. In panel b, the thick black line represents the 4060 m contour, which was the baseline from which slope was calculated (See Section 2.4). (c) Centreline ice velocities from selected years. Velocities were sampled every 120 m along the centreline (black line) shown in a. Lines are colour coded by year and displayed every 4 years from 1988. The black line shows the velocity profile from the 120 m composite. (d) Centreline surface elevation profiles, sampled every 20 m from Pléiades DSMs from 7 November 2017 (blue line) and 19th October 2022 (magenta line). The DSMs were provided by the Pléiades Glacier Observatory (Berthier and others, Reference Berthier2024).
Map showing ice thickness change and its correspondence with supraglacial pond distribution on Tshojo Glacier, for (a) 1975–2000 and (b) 2000–2016. Ice thickness change is derived from the High Mountain Asia Gridded Glacier Thickness Change from Multi-Sensor DEMs, Version 1 (Maurer and others, Reference Maurer, Rupper and Schaefer2019a). Gaps in the ice thickness change data are symbolised in mid-grey. (c) Location and area of all ponds on Tshojo Glacier detected during the study period (1987–2020). Ponds larger than 10 000 m2 are categorised and symbolised by size and colour. Ponds smaller than 10 000 m2 are symbolised by a cyan dot. Ponds larger than 20 000 m2 are named.
The longitudinal surface gradient of Tshojo Glacier's tongue ranged between 0 and 6° in 2017 (Fig. 5b). The majority of the glacier was classified in the 3–6° category and about 20% of the glacier was classified as >2° (Fig. 5b). Longitudinal surface gradient did not show a clear trend with distance up glacier. However, patches of lower longitudinal surface gradient occurred along the glacier centreline line from the terminus to ~10 km up glacier (Fig. 5b). Areas of lower longitudinal surface gradient were also apparent along the eastern margin of Tshojo Glacier's tongue, between ~5 and 11 km up glacier (Fig. 5b). Surface elevation at Tshojo Glacier showed a broadly linear increase along its centreline, from ~4000 m.a.s.l. at the terminus to ~4450 m.a.s.l. at the top of the glacier tongue, approximately 12 km up glacier, and was characterised by undulations of tens to hundreds meters in amplitude (Fig. 5d).
Rates of ice thinning on Tshojo Glacier increased markedly between 1975–2000 and 2000–2016, particularly over the middle of the ice tongue (~3–10 km from the terminus; Figs 6a, b). In 1975–2000, rates of thinning on the lower 8 km of Tshojo Glacier's tongue were generally between −0.25 and −0.75 m a−1, with isolated patches of higher thinning rates (between −0.75 and −1.25 m a−1) and net thickening (Fig. 6a). Thinning rates were somewhat higher up glacier of 8 km, at between −0.75 and 1.25 m a−1 (Fig. 6a). For 2000–2016, ice thinning rates were generally >-1 m a−1 over the middle portion of the tongue (~3–10 km), with substantial patches reaching over −2 m a−1, and maximum rates of thinning occurred between 3 and 8 km from the terminus (Fig. 6b). Over the lower 3 km of Tshojo Glacier's tongue, ice thinning rates for 2000–2016 were similar to 1975–2000, at −0.25 to −0.75 m a−1 and rates of thinning were also comparable between the two time periods above ~10 km (Fig. 6b). For both time intervals and all areas of the glacier, the spatial pattern of ice thinning was heterogenous and patchy (Figs 6a, b). Limited changes in glacier long profile were observed between 2017 and 2022 (Fig. 5d).
3.5. Pond outburst hydrodynamics and downstream impacts
Applying Eqn (3) to our drainage scenarios gave pond volumes of between 0.06 × 106 m3 (10 000 m2) and 1.37 × 106 m3 (100 000 m2), with a worst-case scenario of 3.03 × 106 m3 i.e. the total pond volume recorded in 2008, which was the maximum observed during our study period (Table 1). In our potential breach scenarios for these ponds, the peak discharge ranged between 10.84–92.52 m3 s−1 for scen1 to scen10 and was 119 m3 s−1 for the ‘worst-case’ scenario (Table 1). Our results suggest that the downstream propagation of the flow wave varied from ~6 km (scen1) to ~35 km (scen10) and reached a maximum of 47 km in the worst-case scenario (Figs 7–9). All scenarios reached the Lunana area and, although the worst-case scenario travelled further downstream, it attenuated before reaching the next settlement at Wolathang village (Figs 7–9). In all scenarios, the flood waves arrived at Lunana within an hour of initiation at the upstream boundary, with the worst-case scenario arriving ~5 minutes earlier than the more moderate scenarios (Fig. 8). For scenarios 1–10, the flow magnitude at Lunana was similar in magnitude to the initial discharge, ranging from 5–86 m3 s−1 (Figs 8a, b). The mean flow velocity ranged between 0.1–2 m s−1 while the maximum velocity was up to 9 m s−1 (Figs 7, 10d). Meanwhile, the mean flow depth varied between 0.1–1 m, with a maximum depth of up to 10 m (Figs 7, 10e). Our inundation map shows that the outflow from all modelling scenarios was largely confined within the pre-existing flood plain (Fig. 9). An exception to this was observed in Lunana village (Fig. 9d), where approximately 90 m2 of agricultural land was inundated in scenario 1 to 10 and 1410 m2 was exposed in the worst-case scenario. All the scenarios intersected with the network of footpaths in Lunana, impacting between 0.3 (scenario 1) to 1.2 km (worst-case scenario) of the paths and two suspension bridges.
Flow depth and velocity profile every 200 m along the horizontal river centreline under the various scenarios of supraglacial pond outburst floods from Tshojo Glacier. Scenarios are based on pond area, in increments from 10 000–100 000 m2 (sc1–sc10) and the worst-case (wc) scenario, which is derived from the maximum total ponded area recorded during the study period. The vertical grey dashed lines show the start and end of the Lunana, i.e. the downstream community inundated under some of the model scenarios.
(a) Hydrographs used as the upstream boundary condition, located at the terminus of Tshojo Glacier, and (b) flow hydrograph at Lunana resulted from our multiple scenario pond outburst flood modelling. The locations of Tshojo Glacier and Lunana are shown in Figure 9.
Flood inundation maps for different scenarios of supraglacial pond outburst floods from Tshojo Glacier, detailed in Table 1. (a) Overview of inundation over the entire model domain for each scenario, where a value of 11 means the area is flooded in all scenarios and 1 means the area is only flooded in one scenario. Inundation maps as per (a), focused on: (b) the Lunana region (5–10 km downstream of Tshojo Glacier). (c–d) Flow depth and velocity under the worst-case scenario at the specific locations in Lunana.
Focusing on our sensitivity analysis, varying the manning n value had limited effect on peak and total flow close to Tshojo Glacier but its impact increased with distance downstream. Specifically, at 5 km downstream of Tshojo Glacier, the peak flow decreased by 2% (between 94 and 93.4 m3 s−1) when Manning's number was increased from 0.02 to 0.12 (Fig. S7). This peak flow decreased by 3% (between 94 and 92 m3 s−1) at 10 km downstream and vy 52% (between 78 and 38 m3 s−1) at 20 km downstream (Fig. S7). The changes in total flow and duration exhibited a similar pattern in response to changes in the Manning's n value (Fig. S7). The flow arrival time was sensitive to the choice of Manning's n value and this sensitivity increases with distance from Tshojo Glacier (Fig. S7). For example, at a distance of 5 km downstream, every 0.1 increase in Manning's n value results in to flow delay of 3.6 minutes, totalling almost 40 minutes when Manning's n increased from 0.02 to 0.12. This delay increased to 46 minutes per 0.1 increase in Manning's n and reached a total delay of 8 hour at 20 km downstream. Furthermore, we observed a linear relationship (R 2 = 0.9 to 0.99) between flow arrival time and Manning's n variations, with linearity strengthening as the distance increased from 5 km (R 2 = 0.9) to 20 km (R 2 = 0.99) downstream of Tshojo Glacier.
There was minimal variation in simulated peak and total flow from our model runs using the SRTM DEM vs the high-resolution HMA 8 m DEM (Fig. S8). Likewise, there was no substantial discrepancy in flow arrival at 5 km and 10 km downstream of Tshojo Glacier, but at 20 km downstream, flow modelled with HMA 8 m DEM arrived slightly earlier than the simulations using the SRTM DEM (Fig. S8). The inundation extent from both the DEMs remained within the flood plain in both cases (Fig. S8), although the inundation pattern from the SRTM DEM was patchier at the fine scale.
4. Discussion
4.1. Temporal patterns in pond characteristics
Overall, our data show an increase in total ponded area on Tshojo Glacier between 1987 and 2020 (Fig. 2a). This is consistent with supraglacial lake expansion observed elsewhere in the eastern Himalaya in recent decades, which has been attributed to the combined impacts of climate warming and glacier thinning (e.g. Gardelle and others, Reference Gardelle, Arnaud and Berthier2011; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016; Nie and others, Reference Nie2017; Khadka and others, Reference Khadka, Zhang and Thakuri2018). We also observe large interannual variation in total ponded area on Tshojo Glacier, particularly from 2007 onwards (Fig. 2a). Similar temporal variability has been previously reported elsewhere in Bhutan (Taylor and others, Reference Taylor, Carr and Rounce2022) and the eastern and central Himalaya (e.g. Gardelle and others, Reference Gardelle, Arnaud and Berthier2011; Wang and others, Reference Wang, Xiang, Gao, Lu and Yao2015; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016; Nie and others, Reference Nie2017; Steiner and others, Reference Steiner, Buri, Miles, Ragettli and Pellicciotti2019; Kneib and others, Reference Kneib2021), suggesting that the temporal evolution of ponded area on Tshojo mirrors the regional pattern of large interannual variability, within an overall trend of decadal-scale pond expansion.
Our data show that interannual variations in ponded area and the percentage of the glacier tongue covered by ponds were predominantly driven by changes in the area of the largest ponds (Figs 2c, 4) and that these large ponds occupied similar locations throughout the study period (Fig. 4). This is supported by declassified spy satellite imagery (Fig. S6), which showed that most of Tshojo Glacier's major ponds occurred at the same locations for at least 53 years. This duration of consistent pond locations is among the longest reported for the Himalaya, with previous work demonstrating pond recurrence over ~10-year time intervals in the Everest and Langtang region (Watson and others, Reference Watson, Quincey, Carrivick and Smith2016; Steiner and others, Reference Steiner, Buri, Miles, Ragettli and Pellicciotti2019), reaching up to 28 years at the termini of large Everest glaciers (Chand and Watanabe, Reference Chand and Watanabe2019). We suggest that ponds repeatedly formed at the same locations on Tshojo Glacier initially due to local ice characteristics (e.g. basal topography, surface geometry and uneven melt) and are then maintained via a series of melt-related feedbacks, as observed elsewhere in the Himalaya (e.g. Benn and others, Reference Benn2012; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016). Furthermore, the regular drainage events observed on Tshojo Glacier may also aid pond persistence and the growth of the largest ponds: data from Ngozumpa Glacier, Nepal, suggest that englacial pond drainage may remove debris, causing localised enhancement of melt rates, and prevent ponds from overflowing, thus encouraging higher melt rates and incision at the pond margins (Strickland and others, Reference Strickland, Covington, Gulley, Kayastha and Blackstock2023).
4.2. Controls on pond location and characteristics
On debris-covered glaciers, surface slope is the primary control on the location and characteristics of supraglacial ponds (e.g. Reynolds, Reference Reynolds2000; Quincey and others, Reference Quincey2007; Sakai and Fujita, Reference Sakai and Fujita2010; Salerno and others, Reference Salerno2012; Miles and others, Reference Miles, Willis, Arnold, Steiner and Pellicciotti2017a; Miles and others, Reference Miles2020). Specifically, longitudinal surface gradients of less than 2° promote the formation of large ponds and encourage individual ponds to merge, while those between 2 and 6° are usually associated with the formation of individual ponds, which are transient and cover substantial areas of the ice surface (Reynolds, Reference Reynolds2000). The majority Tshojo Glacier has a longitudinal surface gradient of between 2 and 6° (Fig. 5b), which we suggest is sufficiently flat to encourage widespread pond formation and growth, but too steep for large, linked ponds to form and persist. Furthermore, approximately 20% of Tshojo Glacier's surface has a longitudinal surface gradient of >2°, and these areas tend to be associated with the larger ponds (Fig. 5b), supporting past inferences that localised areas of lower longitudinal surface gradient provide preferential conditions for pond formation and growth.
Ice velocities also influence pond distribution, via their impact on the likelihood of drainage (Miles and others, Reference Miles2020): faster ice flow encourages crevasses to open and ponds to drain, whereas slower flow causes compression, leading to closure of fractures and thus enables pond persistence and growth (Salerno and others, Reference Salerno2012; Miles and others, Reference Miles2016, Reference Miles2017b). Ice velocities on Tshojo Glacier (Figs 5, c) were higher than those observed on the majority of Everest glaciers (Quincey and others, Reference Quincey2007) and elsewhere in Bhutan (Taylor and others, Reference Taylor, Carr and Rounce2022), and Tshojo Glacier's large ponds drained regularly during the study period (Fig. 4). At the same time, the majority of the large, recurrent ponds occur in areas with relatively low surface velocities (<20 m a−1), meaning that ice inflow and down-glacier pond advection are likely to be limited, and hence ponds occur at comparable locations over time (Fig. 4, Fig. S6). Thus, we suggest that ice velocities on Tshojo Glacier (Fig. 5a) were slow enough to enable large, perched ponds to form, grow and persist in similar locations with minimal down-glacier advection, but were sufficiently fast to regularly open fractures and enable perched ponds to connect to englacial drainage networks (Benn and others, Reference Benn, Wiseman and Hands2001; Liu and others, Reference Liu, Mayer and Liu2015; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016, Reference Watson2017; Miles and others, Reference Miles, Willis, Arnold, Steiner and Pellicciotti2017a). This periodic drainage of Tshojo Glacier's ponds may put a limit on pond growth and coalescence (Benn and others, Reference Benn, Wiseman and Hands2001; Thompson and others, Reference Thompson, Benn, Dennis and Luckman2012), resulting instead in numerous large, individual, ponds (Fig. 4). Thus, the life cycle of Tshojo Glacier's large ponds is likely influenced by the configuration of the englacial hydrological system and the frequency with which drainage conduits are connected/blocked (Benn and others, Reference Benn2017). Observations from across the Himalaya document similarly rapid (days to weeks) pond drainage events to those seen on Tshojo Glacier, but pond refilling and drainage recurrence is generally slower, occurring over multiple years, compared to approximately annually at Tshojo Glacier (Benn and others, Reference Benn, Wiseman and Hands2001; Wessels and others, Reference Wessels, Kargel and Kieffer2002; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016; Miles and others, Reference Miles, Willis, Arnold, Steiner and Pellicciotti2017a). We suggest that the more regular connection of perched ponds to Tshojo Glacier's englacial hydrology may result from its faster flow, which would produce more fractures and may encourage more rapid evolution, but mapping of its hydrological system would be required to confirm this.
Overall, our data suggest that Tshojo Glacier's ice velocity and longitudinal surface gradient patterns (Fig. 5) enabled large perched ponds to expand over time in the same locations (Fig. 4), with relatively limited down-glacier advection, and regularly connect with the englacial system (Fig. 4), meaning that their drainage and refilling patterns dominated the total ponded area and percentage area coverage on Tshojo Glacier from 2007 onwards (Fig. 2). This is supported by spatial patterns of ponding on the glacier: most of the large ponds occurred between ~3.5 and 8.5 km from the terminus and the small ponds between ~3 and 10 km, which we attribute to the generally higher longitudinal surface gradients closer to the terminus and higher ice velocities up glacier of 10 km (Figs 5a, b). Furthermore, the limit on expansion resulting from the drainage events meant that the ponds did not coalesce to form proglacial lakes, as seen in the Everest region (Watson and others, Reference Watson, Quincey, Carrivick and Smith2016, Reference Watson2017) and elsewhere in Bhutan (Komori, Reference Komori2008). Thus, Tshojo Glacier has remained at Stage 1 of the Komori (Reference Komori2008) classification, i.e. the appearance and growth of supraglacial lakes on the lower ablation area, for at least 53 years (Fig. S6). This is much longer than its neighbouring glaciers, both in the Lunana area and in the Bhutan-China boarder area more generally, despite the first supraglacial lakes emerging at a similar time across the region (Komori, Reference Komori2008). Given their proximity and similar aspects, we suggest that differences in climate are unlikely to be sufficient to produce these divergent lake growth trajectories. Instead, we speculate that this may result from Tshojo Glacier's considerably larger total area and high-elevation accumulation area, with multiple tributaries delivering ice to the glacier tongue (Fig. 1), which could offset thinning and pond growth.
4.3. Relationship between ice thickness change and supraglacial ponds
We observed a marked increase in thinning rates between 1975–2000 and 2000–2016 (Figs 6a, b) and the areas of most rapid thinning between 2000 and 2016 coincided with Tshojo Glacier's largest ponds, i.e. between ~3.5 and 8 km up glacier of the terminus (Figs 6b, c). Thus, we suggest that Tshojo Glacier's large ponds contributed to glacier thinning and ice loss, as observed elsewhere in the Himalaya (Liu and others, Reference Liu, Mayer and Liu2015; Pellicciotti and others, Reference Pellicciotti2015; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016, Reference Watson2017; Salerno and others, Reference Salerno2017; Miles and others, Reference Miles, Willis, Arnold, Steiner and Pellicciotti2017a). Furthermore, Tshojo Glacier shows a marked reduction in ice velocities with distance down-glacier (Fig. 5) and maximum thinning rates between ~3.5 and 8 km up glacier from the terminus (Fig. 6), which is indicative of an inverted mass balance gradient (Benn and Lehmkuhl, Reference Benn and Lehmkuhl2000; Bisset and others, Reference Bisset2020). Inverted mass balance gradients can form on debris-covered glaciers because debris cover is generally thickest nearer the terminus, which supresses melt, and thins with distance up glacier, which enhances melt at the mid-elevations and leads to a flattening of the surface elevation profile (Quincey and others, Reference Quincey, Luckman and Benn2009; Jouvet and others, Reference Jouvet, Huss, Funk and Blatter2011; Rowan and others, Reference Rowan, Egholm, Quincey and Glasser2015). However, we also note that Tshojo Glacier has a relatively uniform surface profile with distance up glacier (Fig. 5d) (which has remained similar between 2017 and 2022), whereas as many debris-covered glacier tongues experiencing inverted mass balance gradients have a flat longitudinal surface gradient in their mid-elevations, which becomes increasingly flat over time (Quincey and others, Reference Quincey2007). Furthermore, the area of Tshojo Glacier's tongue in which the ponds are situated has not visibly expanded in the +30 years of Landsat coverage (Fig. S5). Thus, we suggest that Tshojo Glacier has an inverted pattern of thinning with respect to altitude (i.e. peak thinning occurs at mid-elevations), but we did not observe a concomitant upstream expansion of the zone of supraglacial ponding.
4.4. Dominance of larger ponds after 2007
After 2007, there was a marked increase in the area of the largest ponds (Fig. 4): ponds over 30 000 m2 were only observed from 2007 onwards and the number of ponds between 20 000 and 30 000 m2 quadrupled, compared to 1987–2003 (Fig. 4). The number of ponds of between 1000 and 10 000 m2 also increased from 284 (1987–2003) to 343 (2007–2018; Fig. 4), although we did see reduction in pond number from 2017–2020 (Fig. 3b). We suggest a number of potential explanations for this switch to the dominance of larger pond areas after 2007:
(i) Differences in satellite imagery: The change in largest pond area coincides with changing from using Landsat 4–5 TM imagery up to 2003 and Landsat 7 and 8 imagery thereafter (Table S1). We suggest that this is an unlikely cause, as we visually inspected each NDWI classification output and both the spatial resolution and band wavelengths for Landsat 4–5 and Landsat 7 are identical for bands utilised in the NDWI classification (Table S1). We thus discount this explanation.
(ii) Changes in longitudinal gradient: Surface slope is a key control on pond distribution and characteristics (e.g. Reynolds, Reference Reynolds2000; Quincey and others, Reference Quincey2007; Watson and others, Reference Watson, Quincey, Carrivick and Smith2016), meaning that surface thinning observed on Tshojo Glacier (Fig. 6) could have reduced longitudinal gradients and promoted ponding. However, comparison of longitudinal profiles from Pléiades DSM suggest this is unlikely to be a factor, as the overall surface slope changed little between 2017 and 2022 and does not have a flat profile in the mid elevations.
(iii) Changes in ice velocity: We observed a slow-down in ice velocities between 1988 and 2003, with a subsequent slight speed up, particularly between ~4 and 10 km from the glacier terminus (Fig. 5c). This earlier slowdown may have encouraged pond expansion, by promoting compressional flow, closing fractures, and thus reducing the frequency of drainage events (Salerno and others, Reference Salerno2012; Miles and others, Reference Miles2016; Reference Miles2017b). These larger depressions may then have been retained for the rest of the study period, providing the accommodation space for large ponds to form, even once drainage events began occurring regularly (Fig. 4).
(iv) Increased air temperature and melt rates: meteorological conditions are likely to be an important control on pond characteristics at interannual and seasonal timescales (e.g. Miles and others, Reference Miles2016; Narama and others, Reference Narama2017; Watson and others, Reference Watson2017; Miles and others, Reference Miles2017b, Reference Miles, Willis, Arnold, Steiner and Pellicciotti2017a; Taylor and others, Reference Taylor, Carr and Rounce2022), but the closest meteorological station with data extending back to 2003 is Gyetsa, which is located 80 km from Tshojo Glacier and is therefore not utilised, as it is unlikely to be representative.
Based on the available evidence, we suggest that the emergence of much larger ponds on Tshojo Glacier after 2003 was at least partly facilitated by a slowdown in ice velocities over the preceding years causing conduits to close, reducing the number of drainage events and facilitating pond growth. Once these hollows were formed, then were then reoccupied in subsequent years, and could then be enlarged by feedbacks between pond presence, glacier ablation and pond drainage. However, we underscore the need for meteorological observations proximal to Tshojo Glacier, to comprehensively assess controls on its subglacial pond characteristics, their potential future evolution under a warming climate, and the implications for future GLOF risk.
4.5. Pond outburst flood impacts
In HMA, ~11% of the cumulative historical instances of GLOFs are from supraglacial lakes and almost all GLOF casualties are due to floods from moraine-dammed lakes and supraglacial ponds (Shrestha and others, Reference Shrestha2023), making it vital to quantify the potential impact of supraglacial pond outbursts. This is particularly important for Tshojo Glacier, as a recent outburst event, on 29 April 2009, necessitated the evacuation of residents within the Punakha region, although no damage was recorded (Komori and others, Reference Komori, Koike, Yamanokuchi and Tshering2012). Our modelling scenarios suggest that, aside from the worst-case scenario, floods generated by supraglacial ponds on Tshojo Glacier are small in magnitude (11–93 m3 s−1) and floodwater typically remains within the existing flood plains, resulting to minimal inundation of structures and agricultural land (Fig. 9, Table 1). While this result is consistent with the minimal damage caused by previous supraglacial pond outburst floods in nearby areas, such as from Halji glacier in north western Nepal (Kropáček and others, Reference Kropáček2015), it is important to note that some of the most damaging floods in HMA are small in magnitude. For example, a small, inconspicuous moraine-dammed lake failure event in the transboundary Shakhimardan catchment of Kyrgyzstan and Uzbekistan in the Central Asia had and initial discharge of only 2–3 m3 s−1 but increased to 150–200 m3 s−1 as it propagated downstream into the valley, resulting in deaths of more than 100 people, while 500–600 went missing and ~ 14 000 people were evacuated (Petrakov and others, Reference Petrakov2020). Similarly, the rapid coalescence and subsequent drainage of a series of ponds through the Changri Shar and Khumbu glaciers caused substantial downstream geomorphic change and destroyed major walking routes (Miles and others, Reference Miles2018). These seemingly low magnitude GLOF possess remarkable mobility and are able to carry a substantial amount of debris (Cook and others, Reference Cook, Andermann, Gimbert, Adhikari and Hovius2018), and, on occasion, can combine with other flood events, such as cloud bursts (Allen and others, Reference Allen, Rastner, Arora, Huggel and Stoffel2015). Thus, it is important to recognise the potential for damage resulting from outbursts from Tshojo Glacier's ponds, particularly if similar cascading events occur. This is particularly important given the large proglacial lakes that currently sit upstream of Tshojo Glacier in the Lunana area (Fig. 1) and their potential to generate substantial flood volumes. Furthermore, despite their comparatively small magnitude, our data show that these pond drainage events occur almost annually, meaning they have the potential to impact downstream communities far more regularly than more irregular outbursts from Lunana's large proglacial lakes. Finally, all scenarios impacted footpaths and bridges in the Lunana area, with the worst-case scenario effecting 1.2 km of footpath and 2 bridges (Table 1). Although this damage is minor in the broader context it could have major impacts on local communities: the Lunana area is the remotest part of Bhutan, and so loss of these routes could lead to further isolation from the nearest road access and create challenges in delivering essential aid and responding to a potential flood event.
4.6. GLOF modelling assumptions, sensitivity and recommendations
The shape of the input hydrograph is an essential initial condition in our pond outburst modelling and can significantly affect the resulting flood characteristics and impacts. Despite 11% of reported GLOF events originating from supraglacial ponds (Shrestha and others, Reference Shrestha2023), there is a paucity of recorded flood characteristics and empirical relationships to construct hydrographs from such outbursts. Consequently, we used a log-normal distribution curve to construct the hydrograph, assuming a positively skewed flow rate that gradually decreases after peaking. While this assumption is associated with notable uncertainties, we considered multiple (10) outburst flood modelling scenarios from Tshojo Glacier to account for the uncertainty associated with hydrograph. The relatively small predicted magnitudes of the pond outburst events can be at least partly attributed to our assumption of a tunnelled-like drainage system, based on lack of available data to the contrary. However, peak discharge and resulting flood magnitude could greatly surpass the flow we have calculated here if the glacial conduit became obstructed and then rapidly released large amounts of water (Rounce and others, Reference Rounce, Byers, Byers and McKinney2017; Miles and others, Reference Miles2018). Furthermore, we cannot discount the worst-case scenario for flooding from Tshojo Glacier, i.e. that all ponds drain at once, due to a rare, but catastrophic, event such as a large mass movement and/or blockage of a glacier conduit and sudden water release. A recent example of such an event is the Chamoli rock-ice-avalanche on 7 February 2021 in Uttarakhand, India, which had a magnitude ~27 × 106 m3 (Shugar and others, Reference Shugar2021). We include this extreme, but unlikely scenario, as catastrophic events are difficult to predict but the frequency of these extreme events is expected to increase in the future, due to climate warming and the resultant deglaciation and destabilization of surrounding slopes.
Our area-volume scaling is based on data that includes comparatively few supraglacial ponds (Cook and Quincey, Reference Cook and Quincey2015; Watson and others, Reference Watson2017) and the equation used to estimate peak discharge (Clague and Mathews, Reference Clague and Mathews2017) is predominantly focused on large, ice dammed lakes. As such, we highlight the need to collect data on area-volume relationship for supraglacial ponds, both to refine our results at Tshojo Glacier and to enable more accurate modelling of supraglacial pond outburst more broadly. We could not incorporate baseflow in our pond outburst modelling scenarios due to the absence of gauged data within our modelling domain. However, base flow plays a critical role in flood propagation, contributing significantly to channel conveyance and overall flow behaviour such as peak flow, flood duration, and floodplain inundation (Sharma and Mujumdar, Reference Sharma and Mujumdar2024). Thus, we recommend installation of water level sensors to capture base flow variations, to inform GLOF modelling both at Tshojo Glacier and the lakes upstream at Lunana. Future GLOF modelling work at Tshojo Glacier should also be conducted in conjunction scenarios for Lunana's largest lakes, to identify the possibility for cascading GLOF events and their potential impacts.
The flood volume resulting from Tshojo Glacier's supraglacial pond outburst flood on 29 April 2009 was estimated to be ~0.5 × 106 m3, which is approximately equivalent to our Scenario 5 (Table 2). However, the 2009 estimation was made ~70 km downstream, in Punakha, (Komori and others, Reference Komori, Koike, Yamanokuchi and Tshering2012), so the actual flood volume and discharge from the event at Tshojo Glacier itself is unknown. Available data suggest that the pond drained rapidly: it was present in satellite imagery from 24 April 2009 image, fully drained by the next available image from 19 May 2009, and had started to refill in December 2009 (Yamanokuchi and others, Reference Yamanokuchi, Tadono, Komori, Kawamoto and Tomiyama2011). This supports our modelling approach of sudden drainage of individual ponds, informed by the range of pond sizes we observed on Tshojo Glacier during the study period. Future work should use high temporal resolution satellite imagery (e.g. PlanetScope) to further constrain pond drainage timeframes, and hence refine our modelled outburst hydrographs. Furthermore, we lack comprehensive data on vital model input parameters, such as flow rate, depth, and velocity, even though the 2009 event is considered a contemporary example of a GLOF in Bhutan. This data deficiency makes it challenging to calibrate our modelled outburst floods and so our scenario-based approach should be considered as a guide to potential inundation extent and location, based on the data that are currently available.
In this study, we assume clear water flow in HEC-RAS 2-D model, as our aim is to provide a first-pass assessment of potential pond outburst flood volumes and extents and HEC-RAS provides a computationally affordable tool to conduct this assessment, which has been used in many comparable studies (e.g. Klimeš and others, Reference Klimeš, Benešová, Vilímek, Bouška and Cochachin Rapre2014; Maskey and others, Reference Maskey, Kayastha and Kayastha2020). However, outburst floods may contain significant amount of debris (Miles and others, Reference Miles2018), which may mean that our clear water results underestimate both flow volume and flow rate, and hence inundation extent. Models like r.avaflow (Mergili and others, Reference Mergili, Fischer, Krenn and Pudasaini2017) and RAMMS (Christen and others, Reference Christen, Kowalski and Bartelt2010) have the capability to model bed erosion and the resulting debris/mud flow, but we lack the necessary input information for Tshojo Glacier. Thus, future work could evaluate the impact of clear vs debris-laden flow, in combination with field data to inform parameter choices.
Due to the absence of open access, high-resolution DEMs that adequately cover the Lunana region, our outburst simulations are based on medium resolution ALOS-PALSAR DEM data. Our sensitivity analysis for part of our model domain suggest that key modelled GLOF output parameters (total and peak flow, and inundation extent) showed limited variation between the ALOS-PALSAR DEM (30 m resolution) and the high-resolution HMA 8 m DEM (Fig. S8). However, our comparison is confined to the Lunana region, which is a relatively flat part of our model domain and so discrepancies could be higher further downstream, where the topography is steeper (Wang and others, Reference Wang, Liu, Liu, Wei and Jiang2016; Liu and others, Reference Liu2019; Rinzin and others, Reference Rinzin2023). Furthermore, the detailed pattern of inundation and parameters such as flood depth does vary between the different DEMs. Thus, our modelling results provide an overview of potential flood routing from Tshojo Glacier, using the datasets currently available, but could be further improved and refined by acquiring and incorporating high-resolution DEMs, e.g. from site-specific UAV surveys. Our sensitivity analysis shows that flow arrival time is highly sensitive to the Manning's n value and that this impact increases with distance downstream (Fig. S7), which has important implications for developing and locating early warning systems (Carey and others, Reference Carey, Huggel, Bury, Portocarrero and Haeberli2011). Thus, we caution against only using coarse-scale (i.e. global to regional scale) land-use products to derive Manning's n values and instead recommend generating basin-specific data from high-resolution satellite imagery sources such as PlanetScope.
Despite the limitations noted above, our inundation map, based on first-order hydraulic modelling, provides a basis for further modelling of outburst scenarios from Tshojo Glacier's supraglacial ponds. Future modelling work should be informed by monitoring of changes in ponded area and, in particular the size of Tshojo Glacier's largest ponds, as climate warms and the glacier evolves. This is particularly important given the frequent drainage of Tshojo Glacier's largest ponds (Figs 4, 5), their impact on total ponded area (Fig. 2), their persistence (Fig. S6), and their potential impacts on downstream populations and HEP stations (Fig. 9 and Table 1). Overall, our work underscores the threat posed by GLOFs to communities in the Lunana region (Rinzin and others, Reference Rinzin2023) and Bhutan more broadly (Taylor and others, Reference Taylor, Robinson, Dunning, Carr and Westoby2023) and highlights the urgent need to quantify their evolving threat in a warming world.
5. Conclusions
Overall, our data show that ponded area and its variability increased substantially between 1987 and 2020 on Tshojo Glacier, with a marked shift occurring between 2003 and 2007. These temporal patterns were predominantly driven by the area of the largest ponds: from 2007 onwards, we observed much larger ponds, which drained regularly. Declassified satellite imagery indicates that these large ponds have occupied the same locations since at least 1967 and that Tshojo Glacier has remained in the first stage of proglacial lake development (i.e. isolated ponds on its surface) during this time, in contrast to the neighbouring glaciers in the Lunana region. We attribute the limited pond coalescence on Tshojo Glacier to its moderate longitudinal surface gradient and ice velocities: these factors are low enough to promote the formation of large ponds, but are insufficient for ponds to coalesce, while velocities are high enough to facilitate regular drainage events, thus limiting pond growth. We suggest that the step-change in the area of Tshojo Glacier's largest ponds between 2003 and 2007 may have been linked to reduced ice velocities during the 1990s, which reduced drainage events and thus promoted growth. However, changes in climate are a potentially important driving factor that could not be investigated, as the closest available meteorological station with data pre-2003 was located ~80 km from Tshojo Glacier. Numerical modelling shows that the flood wave from our pond outburst scenarios could reach between ~6 km (scenario1) and 47 km (worst-case scenario) downstream of Tshojo Glacier's terminus. In all cases, the flood waves reached the settlement of Lunana within an hour of initiation and damage occurred to bridges, agricultural land and footpaths. Although these impacts may appear minor, they could isolate already very remote communities and prevent the delivery of emergency aid. Furthermore, pond outbursts from Tshojo could potentially be magnified by combining with other flood events, such as cloud bursts or drainage of Lunana's large glacier lakes upstream. Our work highlights a number of key areas where we could improve GLOF simulations, both at Tshojo Glacier and from supraglacial outburst floods in general, including: incorporating base flow data; improved knowledge of water routing from burst ponds; specific volume-area scaling relationships for supraglacial ponds; evaluation of the impact of clear vs debris-laden flow; and high resolution data to determining how rapidly ponds drain. Sensitivity analysis suggested that DEM resolution did not markedly impact our results in the Lunana area, but could not be evaluated across the entire domain due to data gaps, and that the Manning's n value is an important control on arrival time and hence for the development of early warning systems. Overall, our results highlight the need to continue monitoring the evolution of Tshojo Glacier and its large ponds, to better understand GLOF threat in Bhutan and how it will evolve as climate warms.
Data availability
Data on ice velocities, the longitudinal gradient of the glacier surface and glacier thickness change are freely available online and the sources are specified in the methods section. Furthermore information on ITS LIVE data processing and resultant dataset are available in Gardner and others, Reference Gardner2018 and via this web link: http://its-live-data.jpl.nasa.gov.s3.amazonaws.com/documentation/ITS_LIVE-Regional-Glacier-and-Ice-Sheet-Surface-Velocities.pdf. Further details on the glacier thickness change dataset are available at: https://nsidc.org/sites/default/files/hma_glacier_dh_mosaics-v001-userguide_1.pdf. Landsat images are freely available online via the USGS Earth Explorer. Shapefiles of supraglacial pond area are available on request from the corresponding author (Rachel.carr@newcastle.ac.uk).
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2024.62.
Open access
