Glacial lake outburst flood hazard assessment by satellite Earth observation in the Himalayas (Chomolhari area, Bhutan)

A case study of glacial lakes outburst flood (GLOF) hazard assessment by satellite Earth observation (EO) and numerical modelling is presented for the supraglacial and ice-contact lakes on Thangothang Chhu glacier, Chomolhari area (Bhutan). Detailed geomorphological mapping, including landslide and rock glacier inventories, as well as surface displacement determination using an interferometric SAR (InSAR) satellite, allowed a GLOF hazard assessment for lake Wa-007 to be performed. Outburst scenario modelling was achieved by combining both empirical and numerical modelling approaches, revealing that only a flood wave can have an impact on the two human settlements located downslope of Wa-007 lake. The worst-case scenario, modelled thanks to r.damflood, allowed the wave-front arrival time, the maximum water depth and the arrival time of maximum water height for the two human settlements to be quantified. A long-term monitoring strategy based entirely on EO data, with an update cycle of 5 years, is proposed to assess the future evolution of the area.

1 Introduction

Thanks to satellite Earth observation (EO), hazard assessments of glacial lake outburst floods (GLOFs) were performed in several mountain chains of the world, including the Alps (e.g. Haeberli, 1983; Huggel et al., 2004; Haeberli et al., 2001), the Andes (e.g. Schneider et al., 2014; Schaub et al., 2015; Frey et al., 2016) and the Himalayas (e.g. Budhathoki et al., 2010; Haemmig et al., 2014) in the last years. Due to its extension, importance in the south-Asiatic geopolitical context and natural sensitivity to climate change, the Himalayas were recently of great interest in research projects related to landslide- (e.g. Ambrosi et al., 2018) and glacier-related (e.g. Allen et al., 2016) hazard assessments, including morphological studies trying to quantify future lake formation by glacier melting (e.g. Schaub et al., 2013; Linsbauer et al., 2016) or quantify permafrost degradation (e.g. Gruber et al., 2017). In this framework, the GLOF hazard assessment is to be considered to be a part of the general topic concerning the connection between high-mountain hazards and climate change (e.g. GAPHAZ, 2017), which integrates the impact of permafrost degradation on the stability of rock slopes above glaciers (e.g. Gruber and Haeberli, 2007; Krautblatter et al., 2013) or the increasing hazard related to ice-rock avalanches falling into new lakes from icy peaks with degrading permafrost (e.g. Salzmann et al., 2004; Haeberli et al., 2017).

This paper presents a GLOF hazard assessment for the supraglacial lake Wa-007 on the Thangothang Chhu glacier, Chomolhari area in the Bhutan Himalayas. This assessment was carried out in the framework of a European Space Agency–World Bank (ESA–WB) collaboration project on satellite EO in the Himalayas, where landslides, floods and GLOF hazard assessments were performed in Nepal and Bhutan (Sauerbier et al., 2015a; Scapozza, 2015; Ambrosi et al., 2018). Within this framework, the GLOF hazard assessment has made it possible to combine observations from flood and landslide hazard mapping, as well as to valorise the products derived from land use and land cover mapping and from the generation of high-resolution digital-elevation models (DEMs).

2 Study area

2.1 Site description

The Chomolhari area is located in north-western Bhutan, in the Thimphu District (Thimphu is also the capital of Bhutan and the largest city of the state) and a few kilometres north of the village of Jangothang (Fig. 1). Chomolhari (7326 m a.s.l., also known as Jomolhari or Chomo Lhari) is the second highest peak in Bhutan, lying on the western border with Tibet.

Figure 1 Google Earth image with the Glacial Lake Inventory of Bhutan using ALOS (Daichi) data (light blue) (GLIB, 2000). The supraglacial lake object of this assessment report is Wa-007.

According to the climate diagram based on 30 years of hourly weather model simulations on the Chomolhari area and calculated for the 4156 m a.s.l. of Jantogang Camp (Meteoblue, 2018), the sum of annual precipitation is 317 mm. Minimum temperatures during cold nights can reach −21 ^∘C in January, whereas maximum temperatures during hot days can reach 17 ^∘C in summer. Mean daily minimum temperature is between −13 ^∘C in January and 3 ^∘C in July (mean annual value of −4.5 ^∘C) and mean daily maximum temperatures is between 0 ^∘C in January and February and 14 ^∘C in summer (mean annual value of 7.75 ^∘C). Mean annual air temperature (MAAT) is then 1.6 ^∘C. Considering a mean rate lapse of 0.6 ^∘C 100 m⁻¹, the 0 ^∘C MAAT is located at about 4420 m a.s.l. As a rough estimation from Google Earth images, the equilibrium line altitude (ELA) on glaciers is at about 5500–5700 m a.s.l. for south and south-east expositions; this corresponds to a MAAT at the ELA of about −6.5 to −7.7 ^∘C. Concerning the permafrost distribution, Gruber (2012) tried to derivate a high-resolution estimation of permafrost zonation, whereas Gruber et al. (2017) reviewed the importance of permafrost degradation in the Hindu Kush region of Himalayas.

The glacial lake object of this assessment report, Wa-007 according to the glacial and glacial lake inventory of Bhutan (GLIB, 2000), is located in the frontal part of the 8.5 km long debris-covered glacier that originated from the glacial basin located on the south-western part of Jichu Drake peak (6790 m a.s.l.). This glacier will be named Thangothang Chhu glacier, according to the river generated at this front (RGI Consortium, 2017). The Wa-007 lake lies directly on the Thangothang Chuu debris-covered glacier and can be considered as a supraglacial lake. Several other small supraglacial lakes are present on the Thangothang Chhu glacier (the second biggest lake is Wa-006), as shown in Figs. 1 and 2. Two other important ice-contact lakes are present on the northern side of the glacier, outside of its historical moraine ridges (lakes Wa-008 and Wa-009).

Figure 2 TerraSAR-X backscattering intensity image of 2 June 2014 showing many supraglacial ponds in black. The names of the lakes are according to GLIB (2000).

The main settlements present on the valley starting from the Chomolhari area is the town of Paro (ca. 20 000 inhabitants), located about 50 km from the Thangothang Chhu glacier tongue. Closest to the glacier lake area, only two settlements are present. Jantogang Camp (4150 m a.s.l.) is located at 2.5 km and the village of Jangothang (4000 m a.s.l.) at 4.2 km downslope of the Wa-007 lake. No other human settlements are present downslope of Jangothang for more than 20 km. Jantogang Camp is the base camp for hiking and climbing in the Bhutan part of Chomolhari area and is composed of three sparse permanent buildings (visible in Google Earth) and by a variable number of tents according to the season. The village of Jangothang presents 16 buildings (also visible in Google Earth) located on two fluvial terraces of the Thangothang Chhu River and several livestock fences. Several other very small hamlets and isolated buildings are located on the valley floor up to 4 km downslope of the main settlement.

2.2 Historical events

On the glaciological mapping in 1962 by Reynolds (2000), the two main ice-contact lakes (Wa-008 and Wa-009) were already present, whereas the supraglacial lake Wa-007 on the terminal part of the Thangothang Chhu glacier was absent. In contrast to other Bhutan glacial lakes (see Richardson and Reynolds, 2000; Nayar, 2009), there is no historical information about past events on the Thangothang Chhu glacier. In the scientific literature, the surface gradients of this glacier were analysed by Reynolds (2000) to explain the formation of supraglacial lakes on debris-covered glaciers, resulting from the stagnation of very slow-moving ice (negative mass balance) with surface gradients lower than 2^∘ and highlighting the supraglacial lake formation since 1966.

3 Material and methods

3.1 Satellite EO mapping

Detailed geomorphological mapping took place by combining three EO base maps and considering the pros and cons of every product. Direct mapping of geomorphological units was done on a VHRO KOMPSAT-2 orthophoto from 27 November 2010, with a ground resolution of 2 m. Shape reconnaissance was improved thanks to hillshades of TanDEM-X DEM, with absolute and relative vertical accuracies of 5 and 2 m (Ambrosi et al., 2018). The joint use of orthophotos and hillshaded DEMs for geomorphological mapping, as well as the differentiation of slope, slope foot, glacial, glaciofluvial and periglacial landforms and deposits, was done according to the recommendations of Ambrosi and Scapozza (2015). Google Earth images were used as ancillary information for supporting the direct mapping and for the interpretation of the lake level fluctuations.

Satellite interferometric SAR (InSAR) – using both differential interferometric SAR (DInSAR, e.g. Bamler and Hartl, 1998) and persistent scatterer interferometry (PSI, e.g. Ferretti et al., 2001) – was used to assess the surface displacements. DInSAR has been applied to stacks of TerraSAR-X (Fig. 2), Sentinel-1, ERS-1/2 SAR, ENVISAT ASAR, ALOS-1 PALSAR-1 and Radarsat-2 images acquired between 1998 and 2015 from ascending and descending orbits. PSI analysis was conducted only with the stacks of ENVISAT ASAR images of the descending orbit and of ALOS-1 PALSAR-1 images of the ascending orbit, because for the other image stacks the number of acquisitions was not sufficient for this kind of processing. For details about satellite SAR acquisitions, processing and interpretation please refer to Ambrosi et al. (2018). The complete set of images and data available for GLOF hazard assessment in the Chomolhari area is listed in Table 1.

Table 1 Data available for this hazard assessment of the Chomolhari area processed within this study and additional sources.

Landslide and rock glacier inventories, comprising degrees of activity, were created from the joint analysis of detailed geomorphological mapping and of surface displacement rates quantified from InSAR, as described in detail by Ambrosi et al. (2018). Ice surface velocities were computed using offset tracking (Strozzi et al., 2002; Paul et al., 2015).

3.2 Hazard assessment

The GLOF hazard assessment was performed based on the checklist for glacier lake hazard assessments provided by the ESA's S:GLA:MO project (Frey et al., 2015). The S:GLA:MO project (Slope Stability and Glacial Lake Monitoring, http://sglamo.gamma-rs.ch, last access: 7 January 2019), funded by ESA and conducted between 2014 and 2015 by Gamma Remote Sensing (Switzerland), the Department of Geography of the University of Zürich (Switzerland), the Department of Geosciences of the University of Oslo (Norway) and Asiaq (Greenland) aimed to assess the hazard potential of glacial lakes based on EO data and products, in situ data and flow modelling. An integrated assessment of hazards related to glacial lakes was developed that address the widely acknowledged dependencies, integrated detection, monitoring and modelling of glacial lakes together with detection, monitoring and modelling of slope instabilities and glacier conditions and behaviour that potentially affect the glacial lakes. In order to provide a reproducible and transparent hazard assessment, the generic checklist-like structure was followed, according to which the analyses were conducted for the specific study case.

3.3 Potential outburst scenarios modelling

Two kinds of models were applied to the potential outburst scenario defined on the basis of the hazard assessment. The first one is based on empirical relationships proposed by Huggel et al. (2002, 2004) and Huggel (2004). Because the equation for lake volume quantification as proposed by Huggel et al. (2002: 31) and Huggel et al. (2004: 1069) is based on self-correlation relating the mathematical products between lake area and lake depth and the factor (the lake area) from which this product had been calculated, here we made a lake volume estimation [V] based on the mean depth [D], which was determined as follows starting from the lake surface [A] (Huggel et al., 2002: 319):

$$\begin{array}{}\text{(1)}& D=\mathrm{0.104}\times A^{\mathrm{0.42}}.\end{array}$$

In the case of ice-dammed lakes, the probable maximum discharge [Q_max] is a function of the volume [V]. Regarding subglacial drainage of ice-dammed lakes, Q_max can be calculated as follows (Huggel et al., 2004: 1071):

$$\begin{array}{}\text{(2)}& Q_{\max }=\mathrm{46}\times {\left({\displaystyle \frac{V}{{\mathrm{10}}^{\mathrm{6}}}}\right)}^{\mathrm{0.66}}.\end{array}$$

The maximum discharge [Q_max] for a worst-case scenario, related to the ice dam suddenly breaking, is calculated by assuming mean drainage duration [t] of 1000 s for maximum estimates as proposed in Haeberli's (1983) empirical relationship (Huggel et al., 2004: 1071).

$$\begin{array}{}\text{(3)}& Q_{\max }={\displaystyle \frac{V}{t}}\end{array}$$

According to Huggel et al. (2002: 322), the lake volume can be doubled in the case of outbursts from moraine-dammed lakes:

$$\begin{array}{}\text{(4)}& Q_{\max }={\displaystyle \frac{\mathrm{2}V}{t}}.\end{array}$$

Finally, the minimum average slope for debris flows from lake outbursts [α], which conditioned the probable maximum travel distance of the debris flow, is a function of the maximum discharge [Q_max] (Huggel, 2004: 35):

$$\begin{array}{}\text{(5)}& \mathit{\alpha }=\mathrm{18}\times Q_{\max }^{-\mathrm{0.07}}.\end{array}$$

The second model adopted here is based on numerical modelling using a two-dimensional dam break flooding simulation and applying a GIS-embedded approach. Here we applied the numerical model r.damflood, a GRASS GIS-integrated two-dimensional dam break model that solves the conservative form of the 2-D shallow water equations using a finite volume method (Cannata and Marzocchi, 2012). The intercell flux is computed by a one-sided upwind conservative scheme extended to a two-dimensional problem. Data used for the simulation include the DEM, reclassified land use to assign Manning's friction coefficient values and a dam-breach map indicating the depth of the dam after collapse (Table 2).

Table 2 Input and output parameters of r.damflood model for estimating the area potentially inundated in case of a dam break. For more detail, see https://grass.osgeo.org/grass74/manuals/addons/r.damflood.html (last access: 19 December 2018).

4 Results

4.1 Satellite EO mapping

Detailed geomorphological mapping is reported in Fig. 3. The areal extent data of the four main glacial lakes on and around the Thangothang Chuu glacier, derived from this mapping, are reported in Table 3. The slightly different values reported on the GLIB (2000) and our own mapping indicate a growth or shrinkage of the lakes between 2000 and 2010. The surface increase of lake Wa-007 is limited (+4.0 %) and rather negligible for Wa-006 (+1.2 %), whereas it was very important for Wa-008 (+16.3 %). In only one case are we in the presence of lake shrinkage, with a decrease of 8.4 % for the Wa-009 lake surface.

Table 3 Surface of the four main glacial lakes of the Thangothang Chhu glacier.

¹ From TanDEM-X DEM produced within this project.

Figure 3 Glacial and periglacial landforms of the lower part of the basin of the Thangothang Chhu glacier. An active rock glacier is also present on the slope above lake Wa-009, but it is outside the frame of the image. Base map: VHRO KOMPSAT-2 orthophoto from 27 November 2010.

The Wa-007 lake does not present evidence of important fluctuations of the water level, for example terraces and erosion scars in the freeboard. A visual comparison between five satellite images taken between 2006 and 2014 and visualized in Google Earth (Fig. 4) allowed the lake levels to be estimated on the basis of the width of the small island present in the south-western part of the lake. The island's extent was the largest on the oldest image (14 February 2006), indicating a low lake level rising until 2009 and 2010 (the two dates with the smallest island widths). In recent times, the lake level was the lowest in the observation period on January 2014 (largest measurement of the island width), reaching a width comparable to 2009 and 2010 in December 2014. It is nevertheless too simplistic to interpret such data as an increase in water level between 2006 and 2010, a decrease until the beginning of 2014 and a new increase. Considering the season of each analysed image, it is probable that the lower lake level in January and February compared to November and December is due to the lower specific runoff in consideration of the typical discharge regime of Himalayan glacier catchments (e.g. Thayyen and Gergan, 2010). We can thus conclude that the lake level does not present important fluctuations so far.

Figure 4 Google Earth images of the small island on the south-western part of the Wa-007 supraglacial lake (location of the white rectangle), used for a qualitative assessment of the lake level fluctuations.

As Wa-007 is a supraglacial lake, it is directly influenced by glacier dynamics. According to the TerraSAR-X images of 2 and 13 June 2014 (Fig. 5) and the ALOS PSI analysis (Fig. 6a), glacier velocities close to the lake are very low. This indicates a stagnant glacier, also explaining the very important cover in debris, which is the consequence of a positive balance between the accumulation of debris on the glacier surface and their evacuation by ice flow. Only the very steep upper part of the glacier moves with velocities larger than 100 m a⁻¹. Calving processes on the lake are possible, as shown by floating ice blocks on the 20 December 2014 image reported in Fig. 4. The detachment of small pieces of ice is also visible by the presence of several gullies on the ice-debris scarp dominating the northern side of the lake.

Figure 5 Horizontal ice surface velocity determined from TerraSAR-X images of 2 to 13 June 2014. Image background is a shaded relief of the DEM. Glacier boundaries (yellow) from RGI Consortium (2017). Lake boundaries (blue) from GLIB (2000).

Figure 6 SAR interferometry of the debris-covered glacier tongue. (a) ALOS PSI between 18 February 2007 and 1 February 2011. Within the red circle the lack of points indicates that velocities might be higher than a few cm a⁻¹. (b) Radarsat-2 InSAR between 4 and 28 November 2014 (24 days). The signal is coherent, indicating a quite uniform movement with a rate between 10 and 50 cm a⁻¹.

The dam of the supraglacial lake Wa-007 is composed of the debris-covered glacier itself: we can then consider it to be an ice dam. Just downstream of the glacier front, several moraine ridges are present (Fig. 3). In case of failure of the ice dam, these ridges can act as moraine dams because of their minimum altitude above lake level. Considering the kinematics of the glacier front, analysed by DInSAR from ALOS PSI and Radarsat-2 images, the ice dam of lake Wa-007 is not stable and presents movements of several decimetres per year on the right side of the glacier (Fig. 6). The freeboard is between 2 and 5 m depending on the location. As the lake is not located directly close to the front of the debris-covered glacier, the width-to-height ratio of the dam is not critical. Considering that no morphological signs of surficial runoff from the lake are present (e.g. gullies), the freeboard is probably still sufficient to contain impact waves generated from ice falls. This is exclusively a morphological interpretation, because a numerical modelling of the possible run-up height of an impact wave (e.g. Schaub et al., 2015) has not been carried out. The lack of morphological signs also indicates that a surface outlet of the lake does not exist. Even in the absence of indications of percolation through the dam, it is probable that a percolation of water is present within the debris-covered glacier or, at least, at the contact between ice and the debris cover. As a consequence, the possibility of a sudden drainage of the lake through the probably polythermal glacier should be considered.

Table 4 Topographical characteristics and magnitude of velocity of the rock glaciers located on the lower part of the Thangothang Chhu glacier basin. Name refers to the landforms mapped in Fig. 3.

In the lower part of the Thangothang Chhu glacier basin, eight rock glaciers were identified and mapped (Fig. 3). An assessment of their state of activity was performed based on ALOS and Radarsat-2 SAR interferometric analyses, which revealed significant surface displacements on the slope dominating the right side of the Thangothang Chhu glacier tongue (Fig. 7). The joint interpretation of the morphological characteristics and of the magnitude of movement allowed the differentiations between active (i.e. containing ice and moving downslope), inactive (i.e. containing ice but not moving) and relict (without ice) rock glaciers (Table 4), according to the conventional classification proposed by Barsch (1996). Only the rock glacier TC-05 can be considered to be a relict. This is located lower in altitude (from 180 to 280 m) with respect to TC-2, TC-3 and TC-4, which present the same slope exposition. Rock glacier TC-6 is polymorphic sensu stricto (Frauenfelder and Kääb, 2000). The presence of two superimposed lobes is also visible in the kinematics of this complex landform, with the upper lobe that presents a magnitude of velocity higher than the lower one (Table 4). Estimated MAAT at the front of active or inactive rock glaciers is between −2.6 and −0.8 ^∘C, whereas for the only relict rock glacier it is +0.2 ^∘C (Table 3). Considering that the front altitude of the lowest active rock glaciers in a region allows for an estimation of the lower limit of discontinuous permafrost in loose debris (e.g. Barsch, 1996; Frauenfelder and Kääb, 2000; Boeckli et al., 2012), the presence of permafrost is probable above an altitude of 4600–4700 m a.s.l., corresponding to a MAAT of about −1.1 to −1.7 ^∘C. The occurrence of active rock glaciers on the right side of the Wa-007 lake indicates a certain presence of permafrost in the catchment of the lake. Considering that the mapped rock glaciers are mainly exposed to the E and SE, and in accordance with the high-resolution permafrost zonation carried out by Gruber (2012), it is probable that the permafrost lower limit on north-exposed slopes may reach lower altitudes.

Except for Wa-006, all the lakes are located below 4600 m a.s.l.; therefore, according to the climate parameters presented above, MAAT at their elevation is between −1.1 and +0.5 ^∘C (Table 3). As a consequence, for the two ice-contact lakes partially dammed by moraine ridges (Wa-008 and Wa-009), the presence of permafrost on their dams is improbable. Despite this, it is nevertheless important to consider that the presence of dead ice, as well as the presence of very sporadic local permafrost conditions, is also possible several hundreds of metres below the lower limit of discontinuous permafrost.

Regarding unstable terrains identified upstream of the glacier lake, several landslides and rockfall deposits were mapped (Fig. 8). The state of activity of the mapped phenomena was evaluated from SAR interferometric data. A certain number of active soil slides affect in particular both lateral moraine ridges. Small rock slides are also present on the slopes above the glacier tongue. The movements of these slides are in all the cases of the order of magnitude of 1 cm a⁻¹ and cannot influence the stability of the Wa-007 lake. Small volumes of fractured rock mass were mapped in the area upstream of the glacier. These mapped volumes could evolve as rockfall but, in any case, they cannot influence the stability of the lake. In addition, areas of active erosion and areas affected by diffuse shallow landslides were recognized on the left slope above the glacier. Considering their limited volume and the phenomena restricted to the shallow surface, they do not represent a hazard factor for the stability of the lake.

The region of the lake itself is not directly exposed to eventual ice-rock avalanches related to hanging glaciers. These kinds of glaciers are, however, present in the upper part of the Thangothang Chhu glacier, with very steep slopes between 4800 and 6000 m a.s.l. and the fastest glacier velocities (Fig. 5). They can have an indirect effect on the Wa-007 lake in the case of major ice-rock avalanches related to the collapse of the main exposed part of these hanging glaciers (sectors with a pervasive presence of crevasses mapped in Fig. 8 as potential detachment zones of an ice avalanche), by causing a potential surge of the glacier related to a wave effect produced by the increase in load on the accumulation zone caused by the collapse of a hanging glacier (e.g. Gardner and Hewitt, 1990). However, as the lake is located at the middle of the glacier tongue, separated from the valley slopes by significant morainic ridges (Fig. 3), it is not directly subject to potential avalanches or rockfall trajectories. The presence of important moraine ridges provides a kind of natural protection from mass movement impacts from the lateral slopes for lakes Wa-006 and Wa-007 (but not for Wa-008 and Wa-009). The impact of ice-rock avalanches can, however, be important for lakes Wa-008 and Wa-009, located outside of the left lateral moraine ridge. Due to their position, these two lakes have a relatively high chance of being impacted by ice-rock avalanches. Deposits probably related to ice-rock avalanches are present in contact with these two lakes at the base of the slope.

Figure 7 SAR interferometry of the slope dominating the Wa-007 supraglacial lake. (a) ALOS PSI between 18 February 2007 and 1 February 2011, with red points indicating velocities higher than 1 cm a⁻¹. (b) Radarsat-2 InSAR between 4 and 28 November 2014 (24 days). Unfilled ovals represent the magnitude of velocity for selected rock glaciers.

Figure 8 Google Earth image reporting a detailed inventory of landslides and other potential unstable areas present in the catchment of the Wa-007 lake excluding rock glaciers.

4.2 Hazard consideration

Three other supraglacial lakes and two ice-contact lakes, not directly linked with the Wa-007 supraglacial lake, are present in the lower part of the Thangothang Chhu glacier catchment. The other supraglacial lakes are not sufficiently large to cause a GLOF and a subsequent debris flow impacting directly on the Wa-007 lake (the larger one is Wa-006, with a surface of 11 100 m²). However, lake Wa-006 in particular can be a potential source of hazard for the Wa-007 supraglacial lake because, if it drains within the glacier, this could lead to a very rapid increase in the lake volume at Wa-007 and cause drainage of this lake. The two ice-contact lakes present a larger surface (17 700 m² for Wa-008 and 25 800 m² for Wa-009), but an eventual debris flow related to a GLOF starting from these two lakes will be channelized between the left side of the valley and the left lateral moraine of the Thangothang Chhu glacier and not attend the Wa-007 supraglacial lake.

Considering the results of geomorphological mapping, the probability of (i) an ice-rock avalanche directly impacting the lake, (ii) large debris flows caused by GLOF related to the other glacial lakes located upstream, (iii) the presence of rapid landslides directly impacting the lake or (iv) very large lake fluctuations, appears to be low. As the flow velocities of the glacier around the lake are close to zero, a blockage of the subglacial drainage causing a rapid increase in the water level by intense precipitation or ice and snow melt is improbable. The main probable lake outburst trigger is related to the thinning of the ice dam related to important glacier retreat in the next years or decennia (significant ablation at the glacier front).

Another aspect that must be considered to be an outburst mechanism is a subglacial drainage by progressive enlargement of the subglacial channels within the glacier. This is the most common drainage process for supraglacial lakes, as highlighted for the large supraglacial lake in 2002–2003 on Ghiacciaio del Belvedere above Macugnaga in Italy (Kääb et al., 2004) and in 2004 on a lake on the Gorner Glacier in the Valais Alps (Sugiyama et al., 2008). The connection between supraglacial lakes should not be forgotten. A potential drainage of Wa-006 lake within the glacier, for example, can lead to a very rapid increase in Wa-007 lake volume and therefore cause drainage of this lake. As a comparable situation, in June 2015, only about 45 km north-east of the Chomolhari area, a GLOF happened: the drainage of a smaller lake caused an outburst of the lower lake (Palden, 2015).

Downstream of the Thangothang Chhu glacier tongue, abundant loose material is available, first by the presence of several frontal moraines and second by the presence of a glaciofluvial plain just outside of the historical moraines. However, considering that the main valley bottom presents an average slope of ca. 4^∘ and that the critical channel slope for erosion is about 8^∘ (Huggel et al., 2004), it is improbable that an eventual flood wave can transform into a debris flow. In other words, material deposition will occur on the first hundreds of metres downslope of the glacier front. Other lakes are not present downstream of Wa-007 and the valley bottom is probably sufficiently large (100 to 300 m) to prevent damming of the main river.

4.3 Potential outburst scenarios modelling

Considering the situation during the time of this assessment procedure and the consequent current hazard evaluation described above, critical conditions triggering a GLOF from the Wa-007 supraglacial lake are probably not present. In future, therefore, radical modifications in the Thangothang Chhu glacier itself (in particular by glacial melting) or in the surrounding terrain (related to permafrost warming and ground ice thawing and activation and acceleration of shallow landslides and erosion zones) must to be considered as realistic. As a consequence, a worst-case scenario considering a complete rupture of the front of the glacier (ice dam failure) and a subsequent lake shrinkage was adopted for the modelling.

The results of empirical modelling are calculated, starting with a Wa-007 lake surface of 92 308 m² (Table 3). The mean lake depth calculated following Eq. (1) is 12.7 m. Considering the scatter in regression presented by Huggel et al. (2002), however, the lake volume was quantified by adopting a mean lake depth of 15 m. The resulting lake volume calculated by multiplying the lake area by the mean depth is 1 384 620 m³. For the ice dam it is necessary to use the hydraulic relationship proposed in Eq. (2). The probable maximum discharge in the case of subglacial draining of the ice-dammed lake is 57 m³ s⁻¹. The maximum discharge increase significantly considering a sudden break of the ice dam; following Eq. (3), it is 1384 m³ s⁻¹. As a consequence, the minimum average slope for debris flows from lake outbursts calculated following Eq. (5) is 13.6 and 10.9^∘ considering the maximum discharge calculated following Eqs. (2) and (3).

The results presented above include an extreme uncertainty because based on regressions providing mean-value statistics (Huggel, 2004). If we assume that a worst-case scenario must work with extreme values, the volume calculations are then based on a maximal mean value of lake depth. A very large uncertainty range is also present for the travel distance determined by the minimum average slope, which is based on a very large range of maximum discharge considering the diverse mechanics of breach formation. In the case of an outburst from moraine-dammed lakes, which is not the case for the Thangothang Chhu glacier, considering the estimated volume of lake Wa-007, the maximum discharge might even reach 2769 m³ s⁻¹, as calculated according to Eq. (4). The subsequent minimum average slope is 10.3^∘. Even considering the scatter between these values, as the main slope of the valley floor is ca. 4^∘, the formation of a debris flow is highly improbable. As the channel slope is lower than 8^∘, considered to be a starting value for erosion (Huggel et al., 2004), it is assumed to only be the effect of a flood wave from the lake outburst.

The wave-front arrival time (Figs. 9 and 10), the maximum water height (Fig. 9) and the arrival time of maximum water height (Fig. 10) were simulated with the numerical model r.damflood. The results show that the wave front reaches the Jantogang Camp in less than 1.5 min and the village of Jangothang ca. 3 min after the dam breach, corresponding to a maximal flow velocity between 23 and 28 m s⁻¹. The maximum water height is reached at about 2.5 km distance from the lake (Fig. 9) in ca. 3 min (Fig. 10), corresponding to a mean flow velocity of 13.9 m s⁻¹. Considering the lateral narrowing of the valley downslope of Jangothang, an estimation of about 10 m of water was simulated in that area. The only two human settlements reached by the simulated flood wave are Jantogang Camp and the village of Jangothang. In the first case, considering that Jantogang Camp is located about 5–6 m above the Thangothang Chhu River bed, the maximum water height modelled is between 0.1 and 1.5 m, which will reach this settlement less than 3 min after the ice-dam breaching and with a mean flow velocity higher than 11.7 m s⁻¹. For Jangothang village, located in correspondence of a widening of the valley bottom, the maximum modelled water height is between 0.1 m on the upper part of the settlement and 4.5 m on its lower part, which is located on a fluvial terrace located only about 2–3 m above the river bed. In this case, the arrival time of the maximum water height is estimated between 4 and 8 min, corresponding to a mean flow velocity between 8.8 and 17.5 m s⁻¹.

Figure 9 Wave front arrival time and maximum water height simulated with the r.damflood model for a complete collapse of the ice-dam of supraglacial lake Wa-007. In the box, a detail of the land cover map. Coordinates are expressed in WGS 1984 UTM zone 45N.

Figure 10 Wave front and maximum water height arrival times simulated with the r.damflood model for a complete collapse of the ice-dam of supraglacial lake Wa-007. In the box is a detail of the land cover map. Coordinates are expressed in WGS 1984 UTM zone 45N.

5 Discussion

5.1 Future evolution

Considering a surface gradient of less than 2^∘ and further lowering of glacial surfaces indicating a thickness loss observed since the beginning of the 1960s on the glaciers of the Chomolhari area (Reynolds, 2000), we can expect a potential increase in area and volume in the Waa-007 lake due to glacier melting. This will also cause a probable change in the dam, by first modifying the dam width, and then the ratio between lake depth and ice thickness (Westoby et al., 2015). This normally determines the floating equilibrium: the subglacial channels typically start opening when the floating equilibrium is reached, also modifying the opening and close of part of the drainage network for meltwater due to modifications in the flow activity of the glacier (Bælum and Benn, 2011; Westoby et al., 2014). With an increase in glacier melting, we can also expect an increase in calving activity on the lake (Sakai et al., 2009).

In the accumulation zone of the Thangothang Chhu glacier, several hanging glaciers may warm in the next years or decennia and then pass from cold-based conditions to polythermal or even temperate-based conditions (Reynolds, 2000). Thus, their dynamics will accelerate and we can also expect an increase in the meltwater production (Thayyen and Gergan, 2010). This will then increase the probability of ice-rock avalanches in the upper part of the Thangothang Chhu glacier.

With a general increase in temperature and a modification in precipitation regime, we can also expect an evolution of the mass movements affecting the slopes around the lake. Warming and partial thawing of permafrost ice can have an effect both on shallow landslides and on the flow velocity and stability of rock glaciers of the area (Daanen et al., 2012). With increasing water availability, we can expect a reactivation or an acceleration of shallow landslides and areas of erosion present on the left valley slope just above the Wa-007 lake (Ambrosi et al., 2018), with an increasing amount of sediment transfer toward the glacier tongue (although not directly in the lake itself). For rock glaciers, we can also expect a reactivation of inactive rock glaciers or an acceleration of active ones with the increasing temperatures of permafrost ice, as evidenced, for example, in the European Alps (e.g. Kääb et al., 2007; Delaloye et al., 2010; Scapozza et al., 2014).

As for the Wa-007 supraglacial lake, the GLOF hazard under current conditions is low, a long-term monitoring strategy based entirely on EO data is proposed (Table 5). With an update cycle of 5 years, this monitoring strategy has the objective of measuring the future developments of the lake to survey the evolution of glacier morphology and dynamics and to monitor the kinematics of the mass movements (including rock glaciers) surrounding the lower part of the Thangothang Chhu glacier catchment. Satellite EO alone will be normally sufficient for monitoring the evolution of this lake. However, as drainage variations are very difficult to detect by EO data, measuring the discharge at the glacier front might be very helpful to detect a GLOF event related to an increase in subglacial drainage by progressive enlargement of the subglacial channels within the glacier.

Table 5 Suggested future monitoring strategy for Thangothang Chhu glacier and supraglacial lake Wa-007, based on an update cycle of 5 years.

5.2 Outburst scenario modelling

The study was structured methodologically following the workflow proposed by Huggel (2004) and based on the integration of remote sensing and GIS modelling for the detection, assessment and detailed analysis of glacial hazards. Hazard assessment was based on EO mapping, terrain surface modelling (both for the production of DEMs and land cover maps) and determination of surface displacements from InSAR, according to the classical methodology applied for landslide hazard assessment based on EO data (e.g. Ambrosi et al., 2018 and references therein). Outburst scenarios were modelled by combining both empirical and numerical modelling, integrating detailed information from EO mapping, such as the dam characteristics, the loose deposits availability and the main slope (derived from the TanDEM-X DEM). Empirical model was applied to quantify the lake volume, the probable maximum discharge and the minimum average slope for debris flows from lake outbursts. This last parameter, compared with the effective topographical and morphological context of the Thangothang Chhu valley, allowed the most plausible scenario for numerical modelling to be defined. Based on these considerations, a debris-flow formation was thought to be very improbable in this hazard assessment. As a consequence, the applied numerical model was chosen only for the simulation of the flood wave from the lake outburst. This approach is particularly robust because, even if an empirical model can be perceived as relatively simplistic model, the geomorphological analysis and the subsequent hazard assessment is based on the detailed knowledge of the terrain characteristics (even if derived from EO). As a consequence, the GLOF numerical modelling, even if it is physically based and integrates the basic hydraulic principles that describe the mechanics of breach formation, represents only additional information about the entire process of hazard assessment and not a unique and unquestionable answer.

6 Conclusions

For Wa-007, the largest lake in the Thangothang Chhu glacier catchment, the likelihood of an outburst flood under the current morphoclimatical conditions appears to be low. The ice dam is, for the moment, relatively stable and sufficiently thick to prevent a dam breaching. The lack of morphological indicators at the surface of the glacier tongue also exclude lake overflowing because of a significant rise in the lake level, mass movements or ice collapse impacts on the lake.

Both empirical and numerical models indicate that a GLOF from Wa-007 lake is potentially dangerous, only when considering a worst-case scenario related to the complete collapse of the ice dam and then a subsequent almost complete emptying of the lake. This event is possible only under radical modifications of the glaciological conditions in the area: important ablation of the Thangothang Chhu glacier tongue, collapses of hanging glaciers located in the upper part of the catchment, which can originate large ice or ice-rock avalanches, or important acceleration and subsequent potential collapse of the rock glaciers and shallow landslides located directly upslope of the lake due to permafrost warming and thawing.

Model estimates revealed that a flood wave rather than a debris flow can impact the human settlements located downslope of Wa-007 lake. The formation of a large debris flow is very improbable because of the low inclination of the valley bottom. Under the worst-case scenario, the Jantogang Camp will be attained after less than 1.5 min from the ice dam collapse by a flood with a maximal water height between 0.1 and 1.5 m, whereas the lower part of the village of Jangothang will be inundated by a maximum water height of 4.5 m about 4 to 8 min after the glacial lake outburst. Glacial lakes also exist in the next valley of the same catchment, but they have not been analysed. Consequently, this hazard analysis is incomplete for the involved village and is related exclusively to lake Wa-007.

The exclusive application of EO techniques allowed a detailed hazard assessment to be performed. Considering the remoteness of the Thangothang Chhu valley, a future monitoring strategy is also proposed based entirely on EO data, with the objectives of (1) detecting potential critical evolution of the glacial lake and of its surroundings, (2) surveying the Thangothang Chhu glacier and Wa-007 lake surface evolution, (3) assessing the hanging glaciers and slope instability evolution, (4) monitoring the kinematics of permafrost-related landforms (rock glaciers) and of shallow landslides.