Hazard Assessment for a Glacier Lake Outburst Flood in the Mo Chu River Basin, Bhutan

: The frequency of glacier lake outbursts ﬂoods (GLOFs) is likely to increase with the ongoing glacier retreat, which produces new glacial lakes and enlarges existing ones. Here, we simulate the outburst of a potentially dangerous glacial lake in Bhutan by applying hydrodynamic modelling. Although the lake volume is known, several parameters connected to the dam breach and the routing of the ﬂood are rough estimates or assumptions, which introduce uncertainties in the results. For this reason, we create an ensemble of nine outburst scenarios. The simulation of magnitude and timing of possible inundation depths is an important asset to prepare emergency action plans. For our case study in the Mo Chu River Basin, the results show that, even under the worst case scenario, little damage to residential buildings can be expected. However, such an outburst ﬂood would probably destroy infrastructure and farmland and might even affect the operation of a hydroelectric powerplant more than 120 km downstream the lake. Our simulation efforts revealed that, by using a 30-m elevation model instead of a 5-m raster, ﬂood magnitude and inundation areas are overestimated signiﬁcantly, which highly suggests the use of high-resolution terrain data. These results may be a valuable input for risk mitigation efforts.


Introduction
The formation and extension of proglacial lakes as a consequence of the current glacier recession is a global phenomenon. From 1980 to 2018, the worldwide number and area of proglacial lakes has increased by around 50%, respectively [1]. Where glacier tongues retreat up-valley, they commonly expose depressions which are filled with melt water. These lakes are dammed by moraine walls from the little ice age or younger deposits, which bear the danger of failure, causing so-called glacier lake outburst floods (GLOFs) or 'jökulhlaups' [2]. In the five major river basins of the Hindu Kush Himalaya, 25,614 glacial lakes covering an area of 1444 km 2 were identified [3]. Both the number and the area of glacial lakes increased in recent years, most substantially in the Eastern Himalaya [4]. In the entire Hindu Kush Himalaya region, the number of glacial lakes, normalized by the total glacier area, is highest in Bhutan (1.8 × 10 −2 ha/ha) [4]. This increase in the number and area certainly enhances both the possibility of GLOF occurrence and the risk potential. In the Bhutan Himalayas, a total of 21 past GLOF events has been identified by field evidence [5]. While the youngest incident happened in 2015 [6], the most devastating outburst took place in October 1994, when the collapse of the Lugge Tsho dam caused severe flooding in the Pho Chu Valley and a total of 21 fatalities [7].
GLOFs may appear after rainfall or excessive melt, followed by incision into the moraine. A lowering of the outflow level increases runoff through the outflow, which In the Mo Chu catchment, the glacier area decrease since 1980 was 25.5%, leaving 195 glaciers with an extent of 114 km² in 2010 [20]. Currently, a total of 66 glacial lakes exists in the catchment [21]. Two PDGLs were identified by the NCHM [13], "Mo_gl235" and Sintaphu Tsho (tsho means "lake"), both are located northeast of Gasa village. Since "Mo_gl235" is the smaller lake and has not been surveyed bathymetrically up to now (Phuntsho Tshering, pers. comm., 12. March 2021), we focus on Sintaphu Tsho. This moraine-dammed glacial lake is located at an altitude of 4480 m a.s.l., has an area of 238,314 m², and a volume of 6.2 mio. m³ [13]. A second lake is located at a distance of 100 m up the valley. While the NCHM [13] mentions a growing supraglacial lake in close proximity, the high-resolution In the Mo Chu catchment, the glacier area decrease since 1980 was 25.5%, leaving 195 glaciers with an extent of 114 km 2 in 2010 [20]. Currently, a total of 66 glacial lakes exists in the catchment [21]. Two PDGLs were identified by the NCHM [13], "Mo_gl235" and Sintaphu Tsho (tsho means "lake"), both are located northeast of Gasa village. Since "Mo_gl235" is the smaller lake and has not been surveyed bathymetrically up to now (Phuntsho Tshering, pers. comm., 12. March 2021), we focus on Sintaphu Tsho. This moraine-dammed glacial lake is located at an altitude of 4480 m a.s.l., has an area of 238,314 m 2 , and a volume of 6.2 mio. m 3 [13]. A second lake is located at a distance of 100 m up the valley. While the NCHM [13] mentions a growing supraglacial lake in close proximity, the high-resolution WorldView-2 imagery suggests that the active glacier is located in a horizontal distance of 250 m from the second lake. A supraglacial lake cannot be identified on the recent image. The upper lake is smaller (ca. 70,000 m 2 ) and seems to be shallower, because the areal Appl. Sci. 2021, 11, 9463 4 of 13 extent changes significantly with the varying water level on different Google Earth images. For this reason, we exclude it from our GLOF study.
For Sintaphu Tsho, the lake size and the short distance of the glacier justify the classification as a PDGL. The lake is located ca. 85 km upstream of the confluence with Pho Chu near Punakha, which is also the first and the largest settlement (population: 21,500 according to worldpopulationreview.com) along the route of a potential GLOF. Two more villages follow downstream 22 km and 36 km after the confluence, whereby two large hydropower projects (Punatsangchhu I and II) are under construction, which involve a 137-m-high and 279-m-wide concrete diversion dam across the river, respectively. Although the power houses themselves are underground, the projects might be vulnerable to large GLOFs.

Digital Elevation Models
In order to test the robustness and computational requirements during the model set-up and first model runs, we used NASA Shuttle Radar Topography Mission (SRTM) v. 3.0. This dataset has a ground resolution of 30 m and was created from Synthetic Aperture Radar observations on the Space Shuttle Endeavour in February 2000. In version 3.0, voids are filled by using Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), Global Digital Elevation Model 2 (GDEM2), USGS Global Multiresolution Terrain Elevation Data (GMTED) 2010, and USGS National Elevation Dataset (NED). SRTM elevation products are very common, freely available used in many GLOF studies [22,23].
The resolution of SRTM faces its limits for hydrodynamic routing, especially in steep or narrow valley sections. We therefore used Euro-Maps 3D DEM (EM3D) Level-A Data with a spatial resolution of 5 m for our final GLOF simulations. EM3D is derived from 36 stereo images acquired by the Indian IRS-P5 Cartosat-1 satellite. One advantage of this DEM is that a rough hydrological correction has been applied to the data by default [24]. The dataset is free of charge for scientific use through the ESA Third Party Mission program, which provided the acceptance of a project proposal. In the study area, EM3D is the most detailed of available DEMs. Although we are aware that 5 m resolution is not yet the ne plus ultra when it comes to the representation of terrain details, we refer to EM3D as "high-resolution" data. Compared to 30-m products, which are often used in GLOF studies, this seems to be justified. The best possible DEM resolution at the moment is 30 cm, but this is beyond our financial and computational resources.
In order to get a constant flow, the DEM was edited by erasing bridges and bumps in the river valley. EM3D is delivered as 16-Bit integer pixel type. It was smoothed by converting it to 32-bit floating point and applying a 5 × 5 box filter. In the last step, we inserted the elevation data of the lake floor into the DEM.
Compared to SRTM, the higher spatial resolution of EM3D provides a better representation of the terrain (see Figure 2). It is expected to result in a more realistic simulation of local hydraulic phenomena in a complex topography. We used the software package, Hydrologic Engineering Center's River Analysis System (HEC-RAS), to simulate the flooding scenarios. In addition to other features, this model can simulate two-dimensional unsteady river hydraulics. HEC-RAS was chosen, because it is freely available and has already delivered realistic results in other studies [18,22,23,25,26]. A full model description can be found in [27].
HEC-RAS is capable of simulating multidirectional flows and super-elevation of flow around channel bends, which are typical features of sudden-onset floods such as GLOFs [28]. We used the Saint-Venant equations, a system of partial differential equations implemented in HEC-RAS, which describe incompressible flow in an open channel. We made use of an implicit finite volume algorithm to simulate the flood as an unsteady clear-water flow. In reality, sediment might be entrained into the flow, changing its rheology and shortening the travel distance. For this reason, assuming a clear water flow limits the validity of our results, which is a common problem in GLOF studies [26]. However, this simplification had to be accepted, because the simulation of erosion, transport, and deposition requires extensive field data which are not available.

Hydrodynamic Modelling Using HEC-RAS 5.0.7
We used the software package, Hydrologic Engineering Center's River Analysis System (HEC-RAS), to simulate the flooding scenarios. In addition to other features, this model can simulate two-dimensional unsteady river hydraulics. HEC-RAS was chosen, because it is freely available and has already delivered realistic results in other studies [18,22,23,25,26]. A full model description can be found in [27].
HEC-RAS is capable of simulating multidirectional flows and super-elevation of flow around channel bends, which are typical features of sudden-onset floods such as GLOFs [28]. We used the Saint-Venant equations, a system of partial differential equations implemented in HEC-RAS, which describe incompressible flow in an open channel. We made use of an implicit finite volume algorithm to simulate the flood as an unsteady clearwater flow. In reality, sediment might be entrained into the flow, changing its rheology and shortening the travel distance. For this reason, assuming a clear water flow limits the validity of our results, which is a common problem in GLOF studies [26]. However, this simplification had to be accepted, because the simulation of erosion, transport, and deposition requires extensive field data which are not available.
The model set-up starts with the definition of the spatial domain. The overall size of our simulated flow area is 260 km². HEC-RAS has the capability to capture sub-grid dynamics. For this reason, we have chosen a 20-m resolution for the computational grid. Combining the area with the mesh width results in 644,892 computational cells that have been modelled with a temporal resolution of one second.
The next steps are the preparation of the elevation data (DEM), the selection of upstream and downstream boundary conditions, and the selection of the channel roughness (Manning coefficient or n-value). A dam breach hydrograph served as an upper boundary condition and input for the hydrodynamic model routine. This hydrograph was simulated by the equations of Fröhlich [29], which are derived empirically from dam failures.
The volume of the lake only equals the outburst volume in case of a 100% outburst. The bathymetric survey of Sintaphu Tsho [13] clearly shows that the lake is overdeepened by glacial erosion, which makes a complete drainage highly unlikely. For this reason and considering the topographic situation, we define the maximum possible outburst as 90% of the lake volume. Medium and small outburst scenarios are set to 70% and 40%, respectively. The absolute volumes and other dam breach parameters are presented in Table 1.
Dam breaches are usually modelled as incisions which form in a certain time, releasing the flood wave known as breach hydrograph. The cross-section of the breach can be simplified as a trapezoid, which is defined by its depth, mean width, and side slopes ( Figure 3). The model set-up starts with the definition of the spatial domain. The overall size of our simulated flow area is 260 km 2 . HEC-RAS has the capability to capture sub-grid dynamics. For this reason, we have chosen a 20-m resolution for the computational grid. Combining the area with the mesh width results in 644,892 computational cells that have been modelled with a temporal resolution of one second.
The next steps are the preparation of the elevation data (DEM), the selection of upstream and downstream boundary conditions, and the selection of the channel roughness (Manning coefficient or n-value). A dam breach hydrograph served as an upper boundary condition and input for the hydrodynamic model routine. This hydrograph was simulated by the equations of Fröhlich [29], which are derived empirically from dam failures.
The volume of the lake only equals the outburst volume in case of a 100% outburst. The bathymetric survey of Sintaphu Tsho [13] clearly shows that the lake is overdeepened by glacial erosion, which makes a complete drainage highly unlikely. For this reason and considering the topographic situation, we define the maximum possible outburst as 90% of the lake volume. Medium and small outburst scenarios are set to 70% and 40%, respectively. The absolute volumes and other dam breach parameters are presented in Table 1. Dam breaches are usually modelled as incisions which form in a certain time, releasing the flood wave known as breach hydrograph. The cross-section of the breach can be simplified as a trapezoid, which is defined by its depth, mean width, and side slopes ( Figure 3). The anticipated outburst volumes define the breach depths Hb for the given topography, which in our case range between 17 and 43 m ( Table 1). The side slopes were set to a ratio of z = 0.7 (or 35°), which is the value proposed by [29] for failure modes other than overtopping. Since Sintaphu Tsho is not connected to its mother glacier and not threatened by ice avalanches from hanging glaciers, this choice seemed to be justified. Mean breach width (B;¯)(B)and formation time (tf) were derived by Fröhlich's empiric equations [29], which are based on data from 74 embankment dam failures: where Vw is water volume (m³), and k_0 is a factor that accounts for the effect of failure mode and is set to "1" for modes other than overtopping.
The resulting values for the dam breach are shown in Table 1. With this information, HEC-RAS generates breach hydrographs. A 90% outburst causes a flood peak of 5700 m³/s already 8 min after dam failure begins. For the 70% and 40% outbursts, the corresponding values are 3250 m³/s after 12 min and 816 m³/s after 26 min.
These breach hydrographs serve as input for the two-dimensional hydrodynamic routing. In this step, the channel roughness is the most important parameter, because it controls flow velocity [23]. According to the HEC-RAS manual [27], n-values between 0.025 and 0.075 are typical. In high mountain environments, channel roughness is high and the values for GLOF modelling found in the literature typically range from 0.03 to 0.15. Supported by empirical data [30] and by field observations on a rafting tour in 2018, we decided to test the sensitivity of this parameter by using three different values: 0.04 (smooth), 0.06 (medium) and 0.08 (rough). When no past GLOF events can be used for calibration, the selection of different n-values within a reasonable range is a practical solution, which was also chosen by most other authors [22,23,25,26]. Channel roughness might change along the course of the river and flood plains or forested areas might have n-values different from those of the river bed. However, we decided not to make a differentiation here, but to use mean n-values. In our opinion, defining different areas would pretend a level of detail which is not justified by any data or observation. The anticipated outburst volumes define the breach depths H b for the given topography, which in our case range between 17 and 43 m ( Table 1). The side slopes were set to a ratio of z = 0.7 (or 35 • ), which is the value proposed by [29] for failure modes other than overtopping. Since Sintaphu Tsho is not connected to its mother glacier and not threatened by ice avalanches from hanging glaciers, this choice seemed to be justified. Mean breach width (B) and formation time (tf) were derived by Fröhlich's empiric equations [29], which are based on data from 74 embankment dam failures: where V w is water volume (m 3 ), and k_0 is a factor that accounts for the effect of failure mode and is set to "1" for modes other than overtopping.
The resulting values for the dam breach are shown in Table 1.
With this information, HEC-RAS generates breach hydrographs. A 90% outburst causes a flood peak of 5700 m 3 /s already 8 min after dam failure begins. For the 70% and 40% outbursts, the corresponding values are 3250 m 3 /s after 12 min and 816 m 3 /s after 26 min.
These breach hydrographs serve as input for the two-dimensional hydrodynamic routing. In this step, the channel roughness is the most important parameter, because it controls flow velocity [23]. According to the HEC-RAS manual [27], n-values between 0.025 and 0.075 are typical. In high mountain environments, channel roughness is high and the values for GLOF modelling found in the literature typically range from 0.03 to 0.15. Supported by empirical data [30] and by field observations on a rafting tour in 2018, we decided to test the sensitivity of this parameter by using three different values: 0.04 (smooth), 0.06 (medium) and 0.08 (rough). When no past GLOF events can be used for calibration, the selection of different n-values within a reasonable range is a practical solution, which was also chosen by most other authors [22,23,25,26]. Channel roughness might change along the course of the river and flood plains or forested areas might have n-values different from those of the river bed. However, we decided not to make a differentiation here, but to use mean n-values. In our opinion, defining different areas would pretend a level of detail which is not justified by any data or observation.
In combination with the three possible outburst volumes, we create an ensemble of 9 model runs or scenarios. From here on, they will be denominated in the following style: "outburst volume (%)_Manning value". As an example, 70_0.06 represents an outburst volume of 70% and a Manning value of 0.06.

GLOF Hydrographs
In Figure 4, the GLOF hydrographs for Punakha and Punstsangchu I are shown.
In combination with the three possible outburst volumes, we create an ensemble of 9 model runs or scenarios. From here on, they will be denominated in the following style: "outburst volume (%)_Manning value". As an example, 70_0.06 represents an outburst volume of 70% and a Manning value of 0.06.

GLOF Hydrographs
In Figure 4, the GLOF hydrographs for Punakha and Punstsangchu I are shown. The peak discharges in Punakha show a variation by a factor of more than four: the worst case scenario has a flood peak of almost 950 m³/s, for the average GLOF it is 520 m³/s and for the best case only 200 m³/s. At Punatsangchu I, HEC-RAS predicts a maximum runoff of 365 m³/s.

Timing and Inundation
There are significant differences in the arrival times (Table 2). In the "worst case" scenario (90_0.04), the flood arrives in Punakha 3:45 h after the breach. Both decreasing outburst volume and increasing channel roughness result in a delay of the flood by up to almost 3 h in the "best case" scenario (40_0.08). The peak discharges in Punakha show a variation by a factor of more than four: the worst case scenario has a flood peak of almost 950 m 3 /s, for the average GLOF it is 520 m 3 /s and for the best case only 200 m 3 /s. At Punatsangchu I, HEC-RAS predicts a maximum runoff of 365 m 3 /s.

Timing and Inundation
There are significant differences in the arrival times (Table 2). In the "worst case" scenario (90_0.04), the flood arrives in Punakha 3:45 h after the breach. Both decreasing outburst volume and increasing channel roughness result in a delay of the flood by up to almost 3 h in the "best case" scenario (40_0.08). Table 2 also shows that total inundation areas are sensitive to the outburst volume, rather than to channel roughness. Within the same outburst volume, areas differ within a range of only 2-4%, while between the 40% and 90% outburst volumes, scenarios with corresponding n-values have a variation of inundation areas of around 25%.
The flood depths also depend on the selected outburst volume rather than on channel roughness. Here, increasing n-values cause higher water stages. If the flood is retarded in a rougher channel, water levels rise. This effect is in the order of few decimeters, while there is a difference in the mean inundation depth of more than 2 m between the 40% and 90% outburst scenario. HEC-RAS calculates the maximum velocity for each raster cell. Table 2 shows the average of these maximum velocities in the flow section from the dam to Punakha. Along the flow path, the values generally decrease downstream due to the decreasing slope, but they also have local variations. For the worst case scenario (90_0.04), they increase in narrow valley sections to a maximum of more than 20 m/s and slowing to less than 1 m/s in wide floodplains ( Figure 5).  Table 2 also shows that total inundation areas are sensitive to the outburst volume, rather than to channel roughness. Within the same outburst volume, areas differ within a range of only 2-4%, while between the 40% and 90% outburst volumes, scenarios with corresponding n-values have a variation of inundation areas of around 25%.
The flood depths also depend on the selected outburst volume rather than on channel roughness. Here, increasing n-values cause higher water stages. If the flood is retarded in a rougher channel, water levels rise. This effect is in the order of few decimeters, while there is a difference in the mean inundation depth of more than 2 m between the 40% and 90% outburst scenario.
HEC-RAS calculates the maximum velocity for each raster cell. Table 2 shows the average of these maximum velocities in the flow section from the dam to Punakha. Along the flow path, the values generally decrease downstream due to the decreasing slope, but they also have local variations. For the worst case scenario (90_0.04), they increase in narrow valley sections to a maximum of more than 20 m/s and slowing to less than 1 m/s in wide floodplains ( Figure 5). Some 45 km below the glacial lake, the flood reaches the five bath houses of Gasa Hot Springs (Figure 6), which is a popular site for tourists. These buildings are affected under every scenario, with flood depths between 1 m (best case) and 1.60 m (worst case). Some 45 km below the glacial lake, the flood reaches the five bath houses of Gasa Hot Springs (Figure 6), which is a popular site for tourists. These buildings are affected under every scenario, with flood depths between 1 m (best case) and 1.60 m (worst case). Figure 7 shows an inundation map for the most populated Punakha region. Here, more buildings are affected under the worst case and under the average GLOF scenario, but not under the best case scenario, inclding a small commercial or industrial complex of buildings on the outskirts of Punakha (a in Figure 7), two buildings at the Punakha riverfront (b in Figure 7), and one religious place (Thangzona Choelkhang) located 700 m upstream of Punakha Dzong (c in Figure 7).   Figure 7), two buildings at the Punakha riverfront (b in Figure 7), and one religious place (Thangzona Choelkhang) located 700 m upstream of Punakha Dzong (c in Figure 7).    Figure 7 shows an inundation map for the most populated Punakha region. Here, more buildings are affected under the worst case and under the average GLOF scenario, but not under the best case scenario, inclding a small commercial or industrial complex of buildings on the outskirts of Punakha (a in Figure 7), two buildings at the Punakha riverfront (b in Figure 7), and one religious place (Thangzona Choelkhang) located 700 m upstream of Punakha Dzong (c in Figure 7).

Comparison with SRTM
From the dam to Punakha, the worst case inundation area on the SRTM DEM is 16.86 km 2 ; this is 26% larger than on the 5-m DEM. It is obvious that the coarser DEM delivers a poorer representation of the riverbed topography (see Figure 2), especially steep banks and terraces that constrain lateral water flow and the widening of the river. On the smoother SRTM DEM, the water can run more easily onto floodplains and is less restricted to its minor bed. This effect and thus the difference between the two DEMs is largest in river sections where a mild slope follows a steeper one. In the very first river section, the mean slope is 4.7%. After 23 km, where the influent enters the main stream, the mean slope is suddenly reduced to 2.5% along the following 7 km. In this section, the inundation area on the SRTM DEM is three times larger than on EM3D.
Inundation depths for the Punakha region are displayed in Figure 8.

Comparison with SRTM
From the dam to Punakha, the worst case inundation area on the SRTM DEM is 16.86 km²; this is 26% larger than on the 5-m DEM. It is obvious that the coarser DEM delivers a poorer representation of the riverbed topography (see Figure 2), especially steep banks and terraces that constrain lateral water flow and the widening of the river. On the smoother SRTM DEM, the water can run more easily onto floodplains and is less restricted to its minor bed. This effect and thus the difference between the two DEMs is largest in river sections where a mild slope follows a steeper one. In the very first river section, the mean slope is 4.7%. After 23 km, where the influent enters the main stream, the mean slope is suddenly reduced to 2.5% along the following 7 km. In this section, the inundation area on the SRTM DEM is three times larger than on EM3D.
Inundation depths for the Punakha region are displayed in Figure 8. Although the velocity tends to be higher on a smoother terrain model, the flood is slowed down at the same time, because larger areas are flooded along the route. In other words, the water flows more rapidly, but also more onto floodplains, which again acts decelerating. In our case study, the latter effects overcompensates the first and, as a result, the flood arrives 15 min later in Punakha when the coarser DEM is used. Although the velocity tends to be higher on a smoother terrain model, the flood is slowed down at the same time, because larger areas are flooded along the route. In other words, the water flows more rapidly, but also more onto floodplains, which again acts decelerating. In our case study, the latter effects overcompensates the first and, as a result, the flood arrives 15 min later in Punakha when the coarser DEM is used.

Hazard Potential
Overall, the hazard potential for inhabited areas is comparably small. This is mainly due to the fact that the valley is sparsely populated and that most buildings are built on an elevated location above the river. The river itself is well incised into the valley floor over large parts of its course; this also helps to limit inundation. When EM3D is used, a total of around ten buildings is affected in the worst case scenario, only two of which seem to be built for residential purposes.
Aside from this, infrastructure, such as roads and bridges, are most probably endangered along the flood route. Moreover, economic harm can be expected from the damage of farm and forestry land as well as from the loss of livestock. Our study, however, did not lead to definitive findings regarding affect to the Punatsangchu I hydropower plant. The worst case flood peak of 365 m 3 /s would most likely not threaten the site, but it may be required to implement some measures to ensure the continued safety of operation. If a lot of sediment load is involved, the flood wave could be a challenge for the desilting chamber.
For the worst case scenario based on SRTM, the hazard potential is significantly higher. In the Gasa Hot Springs, all buildings of the bath, more than a few hundred meters upstream, would be affected too. In Punakha, the flood would hit ca. 20 additional buildings, including a small temple and the Dzong Downstream of Punakha, the damage would be even higher with up to 40 buildings in Khurutang, including a temple, a technical institute, a holiday resort, and a dairy farm that lies within the inundation area. Around 30 buildings in Jagathang and Wangdue Phodrang, respectively, would suffer from a worst case GLOF. In total, more than 130 buildings could be damaged or otherwise affected by flooding. However, we are convinced that SRTM overestimates inundation areas and depths and that the GLOF simulates much more realistic with the 5-m DEM.

Discussion and Conclusions
We applied the HEC-RAS hydrodymamic model to simulate GLOF scenarios with a reasonable effort in data preparation and computational time. The duration of one model run ranged between 37.5 h for the best case scenario and almost 42 h for the worst case. This variation is caused by the computational method of HEC-RAS, which is more complex for wet cells. The most challenging steps were the creation of outburst hydrographs and the determination of the channel roughness. The first depends on the outburst volume and on dam parameters, which are generally unknown. The second is difficult to estimate by means of remote sensing as well as by field techniques. For these reasons, several outburst volumes and n-values were selected to create an ensemble of possible combinations. Future research should aim at the quantification of uncertainties, for example by probabilistic dam breach modelling using a Monte Carlo approach. This could help to keep input parameters within realistic ranges. Another limitation of our results is the fact that sediment transport could not be considered. Here, more case studies in regions with better data availability could deliver valuable information on reasonable model input.
Our modelling results in the Mo Chu basin in Bhutan revealed that channel roughness is not essential for the extent of the inundation, but plays an important role for the velocity and arrival time of the flood. These findings are in agreement with the results from other studies [23].
The hazard potential from a jökulhlaup triggered by the outburst of Sintaphu Tsho is comparably small. Even under a worst scenario, which assumes that 90% of the lake is drained and the channel roughness is rather low (n = 0.04), only few buildings will be affected. However, along its route through sections with narrow valley, such a flood would surely damage bridges and roads, causing logistic problems in this mountain region which is already difficult to access. Attempts to reduce the risk potential could aim at a lowering of the lake level by artificial outlets, pumping or tunneling [31]. Since the remoteness of the lake might hamper extensive construction works, a feasible way can also be to widen or deepen the river channel along the GLOF pathway. A different strategy is to reduce the vulnerability of the valley to flood hazards. This can be achieved by technical measures, such as the strengthening of bridges or the construction of dams. Another important adaption method is the prevention of new buildings in threatened areas through appropriate legislation. Here, our results can help to identify prohibitive zones. The vulnerability can also be reduced by early warning components, enabling residents to save valuable goods or at least their lives. Compared to rapid mass movements, the response time for GLOF events is longer and allows warning systems or the design of an evacuation plan. In the Mo Chu Valley, a cost-friendly version could be the nomination of a local observer in Gasa Hot Springs and to equip him, if necessary, with means of communication in order to alert authorities in the case of a flood event. More security can be provided by an automatic warning system, especially in the case of events which happen during nighttime. In the Punakha-Wangdue Valley, where larger glacial lakes exist, an automatic early warning system is in operation since 2011. There, automatic weather and water level stations are connected to 17 sirens which warn the valley population. The NCHM is planning to extend such systems to other valleys (www.hydromet.gov.bt, accessed on 15 August 2021). Hydrodynamic modelling could be a helpful tool in designing such a system or in preparing flood hazard maps.
Integrated disaster management should also undertake awareness and capacity building programs to train the population in identifying hazards and knowing safe evacuation zones. Due to the catastrophic outburst of Lugge Tsho in 1994 and smaller events afterwards, the awareness of GLOF hazards is comparably high in Bhutan.
Our worst case scenario might mark a limit for a magnitude, which already causes some destruction, but cannot yet be considered catastrophic. One should keep in mind that Sintaphu Tsho has a moderate size and volume compared to other lakes in Bhutan. In the Mo Chu basin, a collapse of the moraine between Thorthormi Tsho and Raphstreng Tsho could release a flood volume of 53 mio. m 3 [13]; this is 9.5 times higher than the worst case at Sintaohu Tsho.
The comparison of the results with a model run based on a coarser DEM impressively proves the importance of high-resolution elevation data. In our case, the GLOF magnitude is significantly overestimated if the 30-m DEM is used. However, at present, this resolution seems to be standard for GLOF modelling. Progress in earth observation will make highresolution products freely available in the future. Making use of this in flood routing applications as soon as possible is highly recommendable.