Articles | Volume 25, issue 6
https://doi.org/10.5194/nhess-25-1841-2025
https://doi.org/10.5194/nhess-25-1841-2025
Research article
 | 
03 Jun 2025
Research article |  | 03 Jun 2025

Exploring implications of input parameter uncertainties in glacial lake outburst flood (GLOF) modelling results using the modelling code r.avaflow

Sonam Rinzin, Stuart Dunning, Rachel Joanne Carr, Ashim Sattar, and Martin Mergili
Abstract

Modelling complex mass flow processes, such as glacial lake outburst floods (GLOFs), for hazard and risk assessments requires extensive data and computational resources. Researchers often rely on low-resolution, open-access datasets and parameters derived from plausibility due to the difficulty involved in conducting direct measurements. This results in considerable uncertainties in forward modelling, potentially limiting the accuracy and reliability of predictions. To determine the sensitivity of the model outputs stemming from input parameter uncertainties in the forward modelling, we selected 9 parameters relevant to GLOF modelling and performed a total of 84 simulations, each representing a unique GLOF scenario in the physically based r.avaflow model. Our results indicate that mass-movement-triggered moraine-dammed GLOF modelling outputs are notably sensitive to five parameters, which are, in order of importance: (1) volume of mass movement entering the lake, (2) DEM datasets, (3) origin of mass movement, (4) entrainment coefficient, and (5) basal friction angle. The GLOF output parameter resulting from the volume of mass movement entering the lake has the greatest coefficient of variation (CV) (47 %), while the internal friction angle had the lowest CV (0.4 %). For future GLOF modelling, we recommend carefully considering the output uncertainty stemming from the sensitive input parameters identified here, some of which cannot be constrained before a GLOF and which must be addressed using statistical approaches.

Share
1 Introduction

Glacial lakes can store millions of cubic metres of water (Zhang et al., 2024, 2023a; Shugar et al., 2020): as of 2020, it is estimated that glacial lakes (0.002 km2) store about 1280.6 ± 354.1 km3 of water around the world (Zhang et al., 2024). The glacial lakes in High Mountain Asia (HMA) have experienced the greatest expansion (46 %) between 1990 and 2018 (Shugar et al., 2020). Furthermore, over 28 % of glacial lakes in HMA are dammed by loose/destabilizing moraines (Fujita et al., 2013; Zheng et al., 2021b) and the majority of glacial lakes (70 %) are exposed to mass inputs in the form of ice/snow avalanches, rockfalls, and landslides from the surrounding slopes (Dubey et al., 2023). Meltwater resulting from the shrinkage of glaciers leads to the formation of new glacial lakes and the expansion of existing ones (Zhang et al., 2015; Wang et al., 2020). This process sometimes exposes them to mass movement from the slopes above and increases the total volume of stored water (Rounce et al., 2016). Additionally, the degradation of permafrost destabilizes the slopes surrounding the glacial lakes, increasing the likelihood of mass movement entering the lake (Huggel, 2009).

Recent work has documented 3151 glacial lake outburst flood (GLOF) events between 850 and 2022 CE globally (Lützow et al., 2023) and 682 GLOF events in HMA between 1833 and 2022 (Shrestha et al., 2023). In HMA alone, these GLOF events have resulted in 6907 human deaths (Allen et al., 2015) and caused damage to more than 2200 buildings, 71 km2 of agricultural land, 163 MW of hydropower capacity, 2000 heads of livestock, and numerous other structures including bridges and roads (Shrestha et al., 2023). However, the reported deaths and damage are significantly underestimated because of patchy documentation (Carrivick and Tweed, 2016). The risk from GLOFs is expected to rise in the future with the anticipated expansion of glacial lakes (Zheng et al., 2021b; Zhang et al., 2023b) compounded by a growing population and the construction of structures in areas prone to GLOFs (Taylor et al., 2023; Nie et al., 2023).

Most GLOF events in HMA start with mass movements entering the lake from surrounding slopes, leading to the displacement of water and waves overtopping the dam (Shrestha et al., 2023; Lützow et al., 2023; Nie et al., 2018). Mass movements such as rock or ice avalanches and landslides entering the lake constitute 70 % of known causes of historical HMA GLOF events (Shrestha et al., 2023). The overtopping waves cause moraine dam incision and dam failure, resulting in a sudden discharge of lake water. To a lesser extent, GLOF events are also triggered by factors such as excess hydrostatic pressure from runoff snow and ice melt, intense rainfall and cloud outbursts, and dam settling caused by the melting of ice cores or internal piping. However, it is important to note that, in some cases, these triggering factors may not necessarily result in a complete moraine dam failure. As the flood propagates further downstream, it can transform into a debris flow and/or a hyper-concentrated flow depending on the geologic and topographic characteristics of the river channel as well as depending on the availability of erodible sediment and its grain size distribution (GAPHAZ, 2017; Schneider et al., 2014; Westoby et al., 2015; Westoby et al., 2014). These complex GLOF process chains are difficult to accurately capture in numerical models, given the large number of processes and parameters, which limits our ability to model the impacts of the hazard cascade as a whole.

Numerical modelling of GLOFs

Previous studies have used various modelling codes such as HEC-RAS (Sattar et al., 2021b), BASEMENT (Worni et al., 2013, 2012; Byers et al., 2018), FLO-2D (Somos-Valenzuela et al., 2015), RAMMS (Lala et al., 2018), and r.avaflow (Mergili et al., 2020b). Most of these codes, however, cannot model the evolution of the GLOF process chain through interaction at the boundary of different processes involved (e.g. the interaction of mass movements with the lake) and dynamic transformation of flow through sediment entrainment and deposition. To address this limitation, some of the studies modelled each component separately and then fed the results of each modelling component into the next stage (Lala et al., 2018; Schneider et al., 2014; Frey et al., 2018). For example, Lala et al. (2018) have used RAMMS to model mass movement from the surrounding slope into the Imja Lake, Nepal; the Heller–Hager model and BASEMENT to model wave propagation across the lake surface; and BASEMENT to model the subsequent downstream hydrodynamic evolution of GLOFs. In contrast, the r.avaflow model (Mergili et al., 2017; Mergili and Pudasaini, 2024) enables the integration of all components of the GLOF process chain and their interactions and transformation without the need to combine the results of different models. It enables the detailed modelling of the GLOF process chain, covering everything from the initial trigger to the downstream propagation. r.avaflow is an open-source, GIS-based tool for simulating mass flows over arbitrary terrain. Furthermore, r.avaflow allows for modification of all input parameters, which sets it apart from many modelling codes (where most of the parameters remain hidden within a “black box”), making it suitable for conducting GLOF parameter sensitivity analysis (Mergili et al., 2017; Mergili and Pudasaini, 2024).

r.avaflow utilizes the total variation diminishing non-oscillatory central differencing (NOC-TVD) numerical scheme (Wang et al., 2004) to solve an enhanced version of the Pudasaini multiphase flow model (Pudasaini and Mergili, 2019). It also offers added features for entrainment, deposition, dispersion, and phase transformation. Because of these features, r.avaflow can model the full process chain of a GLOF and flow transformation due to erosion of bed material and deposition of entrained material (Mergili et al., 2017; Mergili and Pudasaini, 2024). However, the precision of this model output parameters such as peak flow, depth, and velocity depends on the accuracy of various input parameters and initial conditions, such as the release volume of mass, the resolution and vertical accuracy of the digital elevation model (DEM), density, entrainment, and frictional parameters (Mergili et al., 2017).

Because of the significant logistical challenges associated with collecting field data and the financial costs involved in acquiring high-resolution remote sensing data, many of the parameters in GLOF modelling are derived from open-access data, leading to considerable uncertainties in the resultant discharge, inundation extent, and arrival times. Also, certain input parameters such as the volume of mass movement entering the lake are impossible to measure accurately before a GLOF event. For example, the global-scale DEM SRTMGL1 (Shuttle Radar Topography Mission Global 1 arcsec), with a ground resolution of 30 m, is commonly employed in GLOF modelling without adequately considering the inherent uncertainty due to horizontal and vertical inaccuracies in this DEM (Rinzin et al., 2023). Similarly, the origin of avalanches and other mass movements is determined using low- to medium-resolution remote sensing imagery and DEMs, often supplemented by secondary datasets like permafrost data (Obu et al., 2019), which can introduce notable uncertainties (Sattar et al., 2023; Allen et al., 2016). When estimating the volume of avalanches entering the lake, DEM differencing between pre- and post-event conditions can be advantageous for reconstructing historical events (Baggio et al., 2021; Zheng et al., 2021a), although the accuracy is contingent upon the vertical and horizontal accuracy and resolution of the data and the temporal interval between data acquisition. Likewise, when ice is considered the sole source of the mass movement, ice thickness is employed to calculate the avalanche volume (Allen et al., 2022), for which the accuracy of the computed volume relies on the resolution and availability of ice thickness data from the ice-flow inversion model in the region of interest. When the depths of landslides and avalanches are not known, arbitrary thicknesses of 10, 30, and 50 m based on past events (Dubey et al., 2023) are often utilized for forward modelling, further contributing to significant uncertainties in the modelling results (Rounce et al., 2017, 2016; Dubey and Goyal, 2020).

https://nhess.copernicus.org/articles/25/1841/2025/nhess-25-1841-2025-f01

Figure 1Study area. (a) The location of Thorthormi Tsho and its downstream condition in Bhutan. (b) Elevation (as per SRTMGL1) and the overview of glacial lakes in Lunana and settlements along the Phochu and Punatsangchu basins, downstream of Thorthormi Tsho. The downstream settlement is divided into 17 zones (1–17), each 10 km long. (c) Area of Thorthormi Tsho between 1960 and 2020 and the surrounding slope with topography potential (TPP) for mass movement entering Thorthormi Tsho. The downstream settlements in the (d) Lunana and (e) Punakha and Wangdue Phodrang regions. (f) The change in the area of Thorthormi Tsho between 1960 and 2020 (Rinzin et al., 2023) and (g) the buildings (count) and road (km) within 1 km on either side of the river centreline as per the latest version of OpenStreetMap (as of 30 March 2024). Loc-1 to Loc-6 are the locations of the origin of mass movement entering the lake considered for this study. The background map for the panels is hillshade from SRTMGL1.

Moreover, the flow parameters in r.avaflow are adjusted and optimized based on the fit of the model's results to well-documented past events (Mergili et al., 2017, 2020a; Vilca et al., 2021) and the physically plausible range suggested by Mergili et al. (2017, 2018a, b). Efforts to fine-tune parameters to fit with historical events of varying magnitude, temporality, and spatiality have led to the use of wide-ranging values. For example, Mergili et al. (2020b) used an internal solid friction angle of 28° to reconstruct the 1941 GLOF process chain of Lake Palcacocha in the Cordillera Blanca, Peru. In contrast, Vilca et al. (2021) used 45° to model the 2020 glacial lake outburst process chain of Lake Salkantaycocha located in Cordillera Vilcabamba of Peru. Likewise, the value of the basal friction angle ranges between 6–18° (Baggio et al., 2021; Mergili et al., 2020a) (Supplement Fig. S1). Because each GLOF event is inherently distinct, even when originating from the same glacial lake (Emmer and Cochachine, 2013; Lala et al., 2018), inferring reconstructed values from past events for forward modelling introduces substantial uncertainties (GAPHAZ, 2017; Mergili et al., 2020b). Finally, r.avaflow model outputs are extremely sensitive to parameters like the entrainment coefficient value, basal friction angle, and initial release volume (Mergili et al., 2018b, a; Baggio et al., 2021). While all these indicate that the value of input parameters is highly variable depending on specific events, to our knowledge, how changes in the values of these input parameters affect the model output (for example, peak and total flow, flow depth, flow velocity, and arrival time) is not known.

To determine the relative contribution of uncertainties in different input parameters to variability in GLOF extent, we identified 9 out of 38 input parameters and initial conditions relevant to GLOF flow modelling: digital elevation model, the volume of mass movement entering the lake, the origin of mass movement entering the lake, grain density of mass movement entering the lake, volume of the lake, entrainment coefficient, internal friction angle, basal friction angle, and fluid friction number (Table S1 in the Supplement). Our selection was motivated by the recognition that these input parameters are considered essential and have been frequently adjusted in previous studies to align with values inferred from the observed past events (Baggio et al., 2021; Mergili et al., 2018b; Allen et al., 2022; Zheng et al., 2021a; Vilca et al., 2021). We believe that these parameters are the most likely to influence the results of future modelling efforts, making it critical to evaluate their impacts on model outputs. We assessed the sensitivity of the model output to each of these parameters by conducting up to 10 r.avaflow simulations per parameter and varying their values within the range determined from the literature that employed r.avaflow modelling (Fig. S1). We investigated the impact of variation of these parameter values on the model outputs and used the following diagnostic variables: peak discharge, total discharge, flow arrival time, flow depth, flow velocity, and reach distance. We then calculated the coefficient of variation for each parameter and ranked them based on this metric.

2 Study site

Here, our sensitivity analysis was conducted on Thorthormi Tsho (28.10° N, 90.27° E) located in the Lunana region of the Bhutan Himalaya (Fig. 1). The area of Thorthormi Tsho has expanded by  192 % since 1990, evolving into the largest proglacial lake (area of 4.35 km2) in Bhutan by 2020 (Rinzin et al., 2023) (Fig. 1b and f). Although the lake level was lowered by 5 m by artificially draining out the water between 2008 and 2012 (NCHM, 2019a), Thorthormi Tsho is considered the most dangerous glacial lake as it is determined to be susceptible to triggering factors (like mass movement entering the lake) and potentially unstable moraine damming the lake (NCHM, 2019a; Rinzin et al., 2021) (Fig. 1b). In recent years, Thorthormi Tsho has produced two GLOF events (NCHM, 2023), with the first occurring on 20 June 2019 (NCHM, 2019b) and the latest occurring on 30 October 2023. Also, modelling of future predicted GLOFs from Thorthormi Tsho shows it can produce a flood with flow volume of up to 3 × 108 m3 of water with a peak discharge of up to 75 000 m3 s−1, affecting more than 16 000 people and various infrastructure downstream of this glacial lake (Rinzin et al., 2023). This high level of outburst susceptibility and hazard makes Thorthormi Tsho an ideal candidate for GLOF modelling to improve our modelling output of GLOF uncertainty.

Additionally, the Phochu and Punatsangchu basins, located downstream of Thorthormi Tsho, are the most populated basins in Bhutan. The latest updated version of OpenStreetMap (https://www.openstreetmap.org, last access: 30 March 2024), although it does not have 100 % coverage, shows that there are over 7000 buildings, 50 bridges, 4 schools, 687 km of road, and a large area of agricultural land within a 1 km radius of the Phochu and Punatsangchu rivers (Fig. 1c, d, g). Besides, located downstream are the two biggest hydropower plants (Punatsangchu-I and Punatsangchu-II) in Bhutan, poised to become key contributors to the nation's gross domestic product (GDP). Also, the Punakha Dzong, which has great historical and cultural significance to Bhutan, is located downstream of Thorthormi Tsho. This high downstream exposure to GLOF hazard further highlights the importance of understanding GLOF characteristics from Thorthormi Tsho (Fig. 1).

3 Methods

3.1 r.avaflow model framework

r.avaflow is a comprehensive GIS-based open-source computational framework for modelling mass movement from one or more release areas over the defined basal topography (Mergili et al., 2017; Mergili and Pudasaini, 2024). It can model the entire GLOF process chain from the release of avalanches through the dynamic interaction of the avalanche and lake water and overtopping and retrogressive moraine dam erosion and finally to the downstream evolution of the resulting flow (Mergili et al., 2020b; Vilca et al., 2021; Sattar et al., 2023). It can also robustly consider the interactions between the phases as well as erosion and deposition (Mergili et al., 2017). Furthermore, it is equipped with a built-in function for visualization and validation. Because of this capability, r.avaflow has been widely used to model process chains such as GLOFs in high mountains all over the world, mostly to reconstruct past events (Zheng et al., 2021a; Mergili et al., 2020b; Vilca et al., 2021) and to a lesser extent to predict future hazards (Sattar et al., 2023; Allen et al., 2022).

In r.avaflow, the evolution of the flow in space and time is solved using an enhanced version of the Pudasaini multiphase flow model (Pudasaini and Mergili, 2019; Pudasaini, 2012). The flow is computed through depth-averaged conservation of mass and momentum equations for solid and fluid components. These equations involve six differential equations accounting for solid (Ds) and fluid (Df) flow depths, solid (Msx) and fluid (Mfx) momentum in the x direction (Msx=Ds×vsx,Mfx=Df×vfx), and Msy and Mfy in the y direction (Msy=Ds×vsy,Mfy=Df×vfy, where v is the flow velocity) (Mergili et al., 2017). Mohr–Coulomb plasticity is used to compute solid stress, while fluid is subjected to solid volume-fraction-gradient-enhanced non-Newtonian viscous stress. r.avaflow also considers other factors like virtual mass force, viscous drag, and buoyancy. These factors collectively facilitate momentum transfer between the solid and fluid phases, enabling the simultaneous deformation, separation, and mixing of phases as they propagate across the mountain topography (Pudasaini and Mergili, 2019; Pudasaini and Krautblatter, 2014b; Mergili et al., 2020b; Pudasaini, 2012). To numerically solve these differential equations and propagate flow over time and space, r.avaflow uses a high-resolution total variation diminishing non-oscillatory central differencing (TVD-NOC) scheme, a commonly used numerical scheme to handle the advection of quantities, whilst minimizing numerical artefacts like oscillations (Wang et al., 2004). The internal friction angle and basal friction angle, which are crucial factors governing the frictional forces influencing flow rheology, are adjusted based on the solid fraction of the flow material (Mergili et al., 2018b, 2017; Pudasaini and Mergili, 2019).

r.avaflow has three different models, namely, a single-phase shallow-water model with Voellmy friction relation, an enhanced version of the multiphase flow model of Pudasaini and Mergili (2019), and an equilibrium-of-motion model for slow-flow process (Mergili and Pudasaini, 2024). Here, we chose an enhanced version of the multiphase flow model considering an erodible moraine dam and mass movement entering the lake as the solid component and lake water as the fluid component. The multiphase mass flow model can simulate the propagation of three different elements – solids (coarse material including boulders, cobbles, and gravel), fine solids (including sand and particles larger than clay and silt), and fluids (including water and very fine particle including clay, silt and colloids) – and assign each of them a distinct flow rheology (Pudasaini and Mergili, 2019).

Furthermore, r.avaflow has six specific optional functions including conversion of release height to depth, diffusion control, surface control, entrainment, stopping, and dynamic adaption of frictional parameters. The latest version of r.avaflow has four options to compute erosion and entrainment: (i) a model calculated by multiplying the entrainment coefficient with flow momentum, (ii) a simplified entrainment–deposition numerical model of Pudasaini and Krautblatter (2014a), (iii) a combination of models i and ii, and (iv) an acceleration–deceleration entrainment and deposition model. Since models ii to iv are in the experimental phase, here, we used model i, where the amount of entrainment is computed dynamically by multiplying with the user-defined entrainment coefficient (CE) with the total momentum of the flow at the given raster cell and time step (Mergili et al., 2017) (Eqs. 1 and 2).

(1)qE,s=CE|Ms+Mf|αs,Emax,(2)qE,f=CE|Ms+Mf|(1-αs,Emax),

where qE,s and qE,f are the entrainment rates of solids and fluids, respectively; CE is user-defined entrainment coefficient (kg−1); and αs,Emax uses user-defined the solid entrainable material height (m).

We utilized r.avaflow direct, a web-based tool (Mergili and Pudasaini, 2024), to initially generate the sample model script. We modified the sample script by inputting parameters relevant to each experimental set-up and wrote a Bash shell script for all simulations in each experiment to test various parameter values within our predefined range. We developed one master Bash script for each experiment that allowed us to run all experiments in parallel, leveraging the Rocket high-performance computing (HPC) facilities at Newcastle University. All the GLOF simulations were done for Thorthormi Tsho and were run for 1500 s when the flow reached up to  24 km downstream of the lake, depending on values of various parameters defined here. The flow propagation beyond this point and its interaction with the downstream component are beyond the scope of this study.

https://nhess.copernicus.org/articles/25/1841/2025/nhess-25-1841-2025-f02

Figure 2Schematic view of Thorthormi Tsho, surrounding terrain (study area) and input parameters employed for investigating r.avaflow model output sensitivity used in this study. Points 1–6 show the origin of mass movement entering the lake (considered for this study) with the vertical drop and horizontal distance (separated with a slash) to the lake shoreline.

Table 1Key parameters tested in this study to investigate model output sensitivity. Detailed parameters for r.avaflow modelling are provided in Table S1. The constant value is the common parameter value used across experimental setups.

Download Print Version | Download XLSX

Table 2Characteristics of DEM datasets employed in this study to investigate the impact of DEM dataset variation on GLOF modelling results.

Download Print Version | Download XLSX

3.2 Model input parameterization and experimental setups

In the context of r.avaflow, a parameter is a (often user-defined) variable influencing the physical characteristics of the movement or the numerical behaviour of the flow. Parameters can be based e.g. on physics (such as friction angles) or empirical knowledge. Parameters can be represented by global values, by individual values for each raster cell, or by time-dependent values. r.avaflow has a large choice of parameters and initial conditions, such as a DEM representing initial basal topography, the volume of the solid and liquid phase, entrainment and stopping parameters, and density and frictional parameters (Mergili and Pudasaini, 2024) (Table S1). The values specified for these parameters influence crucial aspects of modelled GLOF flow characteristics, including the impacted area, travel distance, travel time, and volume of sediment deposited at the various downstream locations (Mergili et al., 2017). In this study we selected a total of nine parameters which are identified as important and commonly modified in the previous studies (e.g. Mergili et al., 2020a): (1) DEM dataset, (2) volume of mass movement entering the lake (VS), (3) origin of mass movement entering the lake, (4) volume of the glacial lake (VL), (5) grain density of mass movement entering the lake (ρS), (6) entrainment coefficient (CE), (7) basal friction angle (δ), (8) internal friction angle (ϕ), and (9) fluid friction number (CFF). To investigate the impact of DEM dataset variation, we modelled GLOFs by employing freely available global and regional DEM datasets with differing spatial resolution and vertical accuracy (Table 2). For the origin of mass movement entering the lake, we first computed the topographic potential for mass movement entering the lake (Allen et al., 2019) (Fig. 1b) and selected six different sites by considering the topographic potential values and geometry of the lake (Fig. 2). The VS was varied between 106–107 m3. The ρS value range was considered based on assumed combinations of rock and ice parts in mass movement following the approach used in the earlier studies (Allen et al., 2022; Sattar et al., 2023). For parameters 5–9, we gathered various values employed in previous studies (Allen et al., 2022; Mergili et al., 2020a, b; Vilca et al., 2021) and computed descriptive statistics and established the median, upper quantile, and lower quantile for each parameter using these collated values (Fig. S1). We then varied these parameter values between the calculated upper quartile and lower quartile, yielding 10 equally spaced values in total. This range of 10 values was utilized in a total of 10 experiments for each respective parameter whilst holding other parameter values constant at the median value. For example, for the ϕ experiment, the value of ϕ was varied between the upper and lower quantiles, with 10 increments in total, whilst holding constant the other parameter values (Table 1). An overview of the employed parameters and workflow is shown in Fig. 2 (schematic view) and Table 1, while further details on the parameter range used for each experiment are provided in the following section.

3.2.1 Digital elevation model

Here our goal is to constrain model output uncertainty stemming from the use of freely available global and regional DEM datasets. We performed a series of GLOF simulations using four open-access DEM data sources of various resolutions, vertical accuracy, and elevation derivation methods, namely, the High Mountain Asia DEM (HMA-DEM; 8 m) (Shean, 2017), ALOS (Advanced Land Observing Satellite) global digital surface model (AW3D30; 30 m) (JAXA, 2021), NASADEM ( 30 m) (NASA JPL, 2021), and SRTMGL3 ( 90 m) (SRTM, 2013) (Table 2). The details about the DEM datasets used here are provided in the Supplement.

3.2.2 Volume of the glacial lake and mass movement entering the lake

r.avaflow has the option to define the initial release volume of different phases involved in the GLOF process chain. Here, we assume the GLOF was initiated by rock–ice-mixed mass movement entering the lake followed by a tsunami wave hitting the moraine, damming the lake and causing moraine dam failure. Accordingly, we defined the frontal moraine damming Thorthormi Tsho as phase 1 (rock component with ρ=2700 kg m−3), mass movement entering Thorthormi Tsho as phase 2 (rock–ice component), and Thorthormi Tsho as phase 3 (fluid part).

Field-based bathymetric data are essential for accurately estimating VL, which is a fundamental input for the GLOF modelling. However, the remote and harsh environment of glacial lakes in the mountain regions makes it highly challenging to conduct field-based bathymetric surveys. In the case of Thorthormi Tsho, conducting the bathymetric survey is particularly difficult as the lake is filled with debris and icebergs, which makes it hazardous to conduct a survey. In the absence of field-based bathymetry, the area–volume empirical equations are commonly employed by past studies to approximate VL. However, these area–volume scales are based on sparse field data and may not represent the actual volume, leading to substantial uncertainties in the volume estimates (Zheng et al., 2021a; Gantayat et al., 2024). In this study, we calculated the volume of Thorthormi Tsho considering eight area–volume empirical equations which are commonly used for computing the volume of moraine dam lakes (Table S2). Uncertainty in using these empirical equations where evident as the VL estimate for Thorthormi Tsho ranges between 2.1 × 108 and 3.8 × 108 m3. So, to investigate the impact of this uncertainty in modelled GLOF hydraulics, we modelled 10 scenarios of GLOFs by varying VL from 2.0 × 108 to 3.8 × 108 m3 with increments of 20 × 106 m3. For all other experiments, the median value (2.9×108 m3) of VL estimated from area–volume scaling equations was used uniformly. r.avaflow requires spatially varying lake bathymetry to be used as fluid release height rather than the absolute VL. With Thorthormi being a recently formed lake, it has ice thickness data covering the extent of the lake (Farinotti et al., 2019). Therefore, we computed the bathymetry of Thorthormi Tsho by subtracting ice thickness data from the surface DEM (Linsbauer et al., 2016, 2017). Assuming that the present-day lake has been formed by filling the overdeepening, this ice-thickness-derived bathymetry was adjusted to match the volume we employed for each modelled GLOF. For example, for the experiment with a VL of 3 × 108 m3, the ice thickness bathymetry was adjusted to get the same volume (Table S2).

The volume of mass movement entering the lake serves as a fundamental parameter for defining various scenarios in the forward modelling of a GLOF (for example, Allen et al., 2022, and Sattar et al., 2023). However, it is difficult to predict how big or small the mass movement will be. The historical record of VS is sparse but mostly exceeds 1×106 m3 (for example, Zheng et al., 2021a, and Byers et al., 2018), although a recent event in Sikkim Himalaya was known to have been caused by a volume of up to 14.7×106 m3 (Sattar et al., 2025). Additionally, earlier studies in Nepal based on topographic potential and assumed three depths based on avalanches in Russia and the Swiss Alps (10, 30 and 50 m) (Rounce et al., 2016) has estimated the volume between 2.7 × 106 to 6.7 × 106 m3. Here we selected the most likely zone of mass movement entering the lake based on the topographic potential of the slope surrounding the lake. However, accurately constraining VS remains highly challenging and subject to considerable uncertainty. Considering these uncertainties, to test the effect of mass movement of various volumes, we conducted a series of 10 experiments considering VS ranging from 1×106 to 10×106 m3. For the other experiments, we used a uniform VS of 5 × 106 m3 (Table 1).

3.2.3 Origin of mass movement entering the lake

To account for uncertainties involved in the origin of mass movement entering the lake, we identified a total of six mass movement areas, each characterized by different directions and distances to the lake (Figs. 1 and 2). To do this, we first computed topographic potential for mass movement (rock–ice avalanche and landslide) entering the lake based on slope and runout trajectory criteria (Allen et al., 2019). Based on this first-order estimate, we identified the six potential instances of origin of mass movement each with a unique direction, height, and distance relative to the lake (Fig. 2): Loc-1 (slope  1 km away from the headwall, aligned along the longitudinal axis of the lake), Loc-2 (headwall about 100 m away from the lake along the longitudinal axis), Loc-3 (slope 1 km from the right moraine dam), Loc-4 (right moraine dam), Loc-5 (slope  1.5 km from the left moraine dam), and Loc-6 (left moraine dam) (Figs. 1 and 2). We then ran one scenario of a GLOF from each of these locations (total of 6), keeping all input variables constant for each model run. For the rest of the experiments, we used a mass movement scenario from Loc-1.

3.2.4 Grain density of mass movement entering the lake

Our goal here is to assess the impact of ρS entering the lake, which serves as a proxy for the ice-to-rock ratio. Accordingly, we consistently set the grain density of phase 1 (moraine) at 2700 kg m−3 across all the experiments, whilst the fluid density of phase 3 was also held constant at 1000 kg m−3. In the earlier studies, ρS had been used as a proxy for the portion of a rock–ice component of mass movement entering the lake, which is highly uncertain (Vilca et al., 2021; Allen et al., 2022). The phase separation of rock and ice components of the mass movement with different densities is not well established in r.avaflow (Vilca et al., 2021). Therefore, in this study, following Sattar et al. (2023), a portion of snow and ice in the mass movement is considered fluid by adjusting the ρS of the phase 2 represented by the mass (Table S3). In our experiment set-up, this is executed by varying the ρS value between 2700 kg m−3 (representing 100 % rock) and 1000 kg m−3 (representing 100 % ice) (Table 1).

3.2.5 Entrainment coefficient

Material entrainment due to bed erosion can make the flow more concentrated and thus increase the volume, resulting in spatial and temporal variation in flow. In the r.avaflow model, the user must define entrainment height in the form of a raster covering the entire model domain, which can be identified using either remote sensing imagery or fieldwork (Mergili and Pudasaini, 2024). However, accurately measuring the depth and spatial extent of these erodible materials covering the entire model domain is highly challenging, even with field surveys, and remains infeasible using remote sensing techniques alone. In this study, we focused on the frontal moraine, as it is measurable using remote sensing data. The extent of the moraine was mapped using high-resolution Google Earth imagery, while its height was determined using HMA-DEM (Fig. 1). The amount of entrainment itself is dependent on the user-defined CE. In r.avaflow the base-10 logarithm of CE must be entered (Mergili et al., 2018a, 2017). Here, we modelled 10 scenarios of GLOFs by varying CE between 10−6.95 and 10−5.85 kg−1 (Table 1).

3.2.6 Frictional parameters

The internal friction angle (ϕ), basal friction angle (δ), and fluid friction number (CFF) mechanically control the basal shear stress, internal deformation, anisotropy of the stresses, and hydraulic pressure gradient of the solid constituents (Pudasaini and Krautblatter, 2014b), which are essential attributes influencing flow runout distance and time. Within the r.avaflow model set-up, a user can use either spatially varying values for these frictional parameters using a raster map or one absolute value (Mergili and Pudasaini, 2024; Mergili et al., 2017). In this study, we computed 10 experiments for each of these frictional parameters. Specifically, by varying ϕ between 25 and 35°, δ between 10 and 14°, and CFF between 0.027 and 0.050 (Table 1).

3.3 Sensitivity analysis

Here we use sensitivity analysis to determine how variations in the initial values for key input parameters impact the modelled GLOF hydraulics (Saltelli et al., 2004). Thus, our goal is not to determine the “correct” value for each parameter but to determine the r.avaflow input parameter(s) that cause the most variation in the model output. To constrain this variability, we mainly focused on examining the peak discharge, total discharge, and flow arrival time as the output metrics. The flow for all the experiments was measured from the profile immediately beneath the moraine dam (profile 1 in Fig. 2). We calculated the peak and total discharge based on the flow data obtained from the same profiles (Fig. S2). The flow arrival time was considered the average value across the time recorded from the profiles located 3, 6, and 9 km downstream of Thorthormi Tsho (profiles 2, 3, and 4 in Fig. 2). For the scalable parameters, we also computed simple linear regression considering input parameters as the independent variable and model output as the dependent variable. To ascertain the sensitivity of the three model output metrics defined earlier for variations in values across all parameters, we computed the coefficient of variation (CV) for individual parameters and subsequently ranked them based on this metric. The CV is a statistical measurement of the dispersion of data points around the mean, regardless of the units used to measure it. The CV is deemed suitable here since the r.avaflow output variability is caused by input parameters that are measured in different units. To calculate the CV, we took the standard deviation of the output value range of a particular experiment (e.g. peak discharge) and divided it by the mean of the same output range (Abdi, 2010). We multiplied each CV value by 100 to express them in percentage form.

https://nhess.copernicus.org/articles/25/1841/2025/nhess-25-1841-2025-f03

Figure 3The line graph showing the (a) rate of flow (Q), (b) maximum flow depth (for total and debris part), (c) flow arrival time, and (d) flow velocity profiles along the river centreline generated by conducting a sequence of r.avaflow simulations, employing different types of DEM datasets. The solid and dashed lines (in panel a) show the maximum flow depth and debris part of the flow. The profiles are plotted as rolling averages every 100 m.

Download

Table 3Percentage change in flow velocity, depth, and arrival time resulting from variation in values of different input parameters we employed in this study. The total percentage (%) change represents the output variation resulting from the maximum and minimum input parameter values used in the experiment. The average percentage (%) change is calculated as the mean change across all incremental steps employed in setting up the experiment. The arrival time is the average of the record from three locations (profiles 2, 3, and 4) (Fig. 2). Flow velocity and depth are mean values taken from the river centreline. The negative () value indicates that increasing input parameters decreases output parameters. The detailed flow pattern along the river centreline is provided in Figs. S3, S4, S5, and S6.

Download Print Version | Download XLSX

4 Results

4.1 Effect of the DEM dataset

When the GLOF was modelled employing freely available global and regional DEM datasets (HMA-DEM, AW3D30, NASADEM, SRTMGL3), our results showed a variation in peak and total discharge of a GLOF from Thorthormi Tsho by almost 100 % and 400 %, respectively (Fig. 3). Specifically, HMA-DEM consistently produced the lowest-magnitude GLOF, while SRTMGL3 consistently produced the highest. Although NASADEM and AW3D30 have a similar spatial resolution, notable differences (65 %) in peak discharge and other hydraulic properties emerged between simulations done using these datasets (Fig. 3a).

We observed a significant fluctuation in the mean flow depth (82 %) and velocity (43 %) along the flow path resulting from the change in DEM datasets (Table 3). For instance, the mean flow depth along the river centreline ranged from 39 m (HMA-DEM) to 54 m (SRTMGL3) (Table 3) and the flow reach distance increased from 15.5 km (HMA-DEM) to 24.2 km (SRTMGL3). Once again, NASADEM and AW3D30 yielded significantly different maximum flow depths (8.5 %) and reach distances (72 %) (Fig. 3b). The flow arrival time varied by around 59 %. Flows derived from SRTM data always arrived earlier, while those using HMA-DEM were consistently most delayed (Fig. 3c). There was also a similar variation in the solid concentration of the flow in response to changes in input DEM datasets (Fig. 3b).

https://nhess.copernicus.org/articles/25/1841/2025/nhess-25-1841-2025-f04

Figure 4Flow rate and depth resulting from mass movement entering the lake from different locations: (a–f) Loc-1 to Loc-6.

4.2 Effect of the origin of mass movement entering the lake

We modelled mass movement with a volume of 5 × 106 m3 entering the lake from the surrounding slopes at various locations identified in Sect. 3.2.3. Our study found that the GLOF process chain initiated by mass movements from Loc-1 to Loc-6 results in a significant fluctuation in the GLOF output (Fig. 4). The peak discharge varied by approximately 200 %, with the total discharge varying by 55 % (Fig. 4). Likewise, the mean flow depth, velocity, and arrival time fluctuated by 65 %, 82 %, and 20 %, respectively (Table 3). By comparison, the flow resulting from the GLOF initiated by mass entering from Loc-1 (Fig. 4a) (1 km from the headwall) and Loc-5 ( 1.5 km above the left moraine) (Fig. 4e) produced the highest-magnitude GLOF, while that from Loc-4 (right lateral moraine) (Fig. 4d) was the lowest. For example, the highest peak (180 × 103 m3) and total discharge (11 × 106 m3) occurred from Loc-1, while the lowest peak (60 × 103 m3) and total discharge (7 × 106 m3) were from Loc-4 (Fig. 4a and d). The longest flow reach distance (25 km) was produced by Loc-1 and Loc-5, while the shortest was from Loc-3 (10 km) (Fig. 4c). The solid volumetric portion did not exhibit significant fluctuation, with concentration ranging from 4 % (Loc-4) to 5 % (Loc-2) (Fig. S3).

https://nhess.copernicus.org/articles/25/1841/2025/nhess-25-1841-2025-f05

Figure 5Linear regression between input parameter value variation and percentage change in the resulting peak discharge (Q). The linear regression is computed only for the scalable parameters (a) volume of mass movement entering the lake, (b) mass movement density, (c) volume of the glacial lake, (d) entrainment coefficient, (e) basal friction angle, (f) internal friction angle, and (g) fluid friction number. The shaded envelope around the line shows a 95 % confidence interval.

Download

4.3 Effect of the volume and grain density of mass movement entering the lake

To separate the effect of variation in the volume (VS) and grain density (ρS) of mass movement entering the lake, we simulated all 10 scenarios of the GLOF event for each volume and density variation using the mass movement initiated from Loc-1. Here we observed that only VS variation leads to a very large variation in the resulting peak (1160 %) (Fig. 5a) and total discharge (2500 %) (Fig. 6a). Subsequently, this resulted in maximum variation in flow characteristics such as mean flow depth (923 %) and flow arrival time (143 %) (Table 3, Figs. S3–S6). Conversely, the ρS variation showed the least impact on both peak (5 %) and total discharge (24 %) (Figs. 5b, 6b, and 7b) and subsequent characteristics such as flow depth (31 %) and velocity (2 %) (Table 3 and Figs. S3–S6). Both VS and ρS variation did not result in significant fluctuation in the solid volumetric concentration of the flow (Fig. S3).

https://nhess.copernicus.org/articles/25/1841/2025/nhess-25-1841-2025-f06

Figure 6Linear regression between input parameter value variation and the resulting total discharge. The linear regression is computed only for the scalable parameters (a) volume of mass movement entering the lake, (b) mass movement density, (c) volume of the lake, (d) entrainment coefficient, (e) basal friction angle, (f) internal friction angle, and (g) fluid friction number. The shaded envelope around the line shows a 95 % confidence interval.

Download

4.4 Lake volume

Variation in volume of the glacial lake (VL) from 2 × 108 to 3.8 × 108 m3 results in total peak discharge variations of 25 % with an average of 8.7 % (Fig. 5c). Likewise, the total discharge varied by 18 % with a mean variation of 3 % (Fig. 6c). However, the resulting flow characteristics including flow depth, velocity (Table 3), and arrival time (Fig. 7c) all showed minimal fluctuation, as the total percentage change was well below 10 % (Table 3).

https://nhess.copernicus.org/articles/25/1841/2025/nhess-25-1841-2025-f07

Figure 7Linear regression between input parameter value variation and flow arrival time. The linear regression is computed only for the scalable parameters (a) volume of mass movement entering the lake, (b) mass movement density, (c) volume of the lake, (d) entrainment coefficient, (e) basal friction angle, (f) internal friction angle, and (g) fluid friction number. The shaded envelope around the line shows a 95 % confidence interval.

Download

4.5 Effect of the entrainment coefficient

Variations in the entrainment coefficient (CE) substantially impact the resulting GLOF output, causing fluctuations in total discharge by 123 %, although its impact on peak discharge was minimal (13 %) (Figs. 5d and 6d). These changes also affect the flow characteristics including mean depth (49 %) and reach distance (20 %) (Table 3) but had a minimal effect on arrival time (3 %) (Fig. 7d). Interestingly, CE variation exhibited a threshold effect on peak discharge, total discharge, and flow arrival time. The total discharge increases linearly with increasing CE up to 10−5.97 kg−1 but decreases after this. On the other hand, peak discharge increased when CE increased from 10−6.95 to 10−6.83 kg−1, remained constant between 10−6.83 and 10−6.09 kg−1, and decreased from there on. The arrival time decreased when the CE value was increased from 10−6.46 to 10−6.46 kg−1 but increased upon a further increase in the CE value to 10−5.85 kg−1. Also, unlike other parameters, entrainment variation affected the solid concentration of the flow (Fig. S3). An increase in the entrainment coefficient from 10−6.95 to 10−5.85 kg−1 led to a 30 % increase in the mean solid volumetric concentration of the flow.

4.6 Effect of the frictional parameters

Among the frictional parameters, the basal friction angle (δ) resulted in a significant fluctuation in GLOF magnitude (Figs. 5e, 6e, 7e) and subsequent flow characteristics (Figs. S3–S6). While the variation in the fluid friction number (CFF) had a minimal impact on the resulting peak and total flow (Figs. 5g, 6g), it notably altered other flow characteristics, such as flow velocity (55 %) and depth (70 %) (Table 3). The increase in δ from 10 to 14° resulted in a peak and total discharge decrease of 36 % (Fig. 5d) and 32 % (Fig. 6d), respectively. Likewise, the flow velocity decreased by 23 %, resulting in a delay in flow arrival by 68 % (Table 3). Conversely, the peak flow decreased by 2 % only in response to an increase in the internal friction angle (ϕ) from 25 to 35° (Fig. 5f). The variation in all values of frictional parameters did not result in a significant change in the solid volumetric concentration of the flow (Fig. S5).

https://nhess.copernicus.org/articles/25/1841/2025/nhess-25-1841-2025-f08

Figure 8The coefficient of variation for (a) peak flow, (b) volume, (c) time, (d) average flow depth along the river centreline, (e) flow velocity along the river centreline, and (f) average across all these output parameters.

Download

4.7 Comparison of the effects of all parameters

To compare output sensitivity resulting from all input parameters considered here, we calculated the percentage of the coefficient of variation (CV) for peak flow, total discharge, arrival time, flow depth, and flow velocity. We further computed the average CV (avg CV) across all these output variables and examined the overall impact of each input parameter variation (Fig. 8). Comparing all these output indicators, VS had the greatest impact (avg CV of 47 %), followed by DEM datasets (avg CV of 35 %) and the origin of mass movement (avg CV of 21 %). Other input parameters like δ and CE also caused significant variation in resulting GLOFs. Notably, CFF significantly impacted flow depth with its CV of 16 % despite having minimal impact on other flow characteristics (Fig. 8).

For the seven scalable parameters, we computed linear regression (Figs. 5–7). The linear regression analysis revealed that the five parameters: VS (R2=0.99), ρS (R2=0.96), VL (R2=0.81), δ (R2=0.96), and CFF (R2=0.83) demonstrated strong explanatory power in accounting for the variability in GLOF peak discharge (Fig. 5). Among these sets of parameters, VS (m=1.6), ρS (m=0.085), and VL (m=0.07) showed a positive relationship, while δ (m=-0.347), CE (m=-0.091), and CFF (m=-0.09) exhibited a negative relationship (Fig. 5). All seven parameters (R2>0.86) except for ϕ (R2=0.59) and VL (R2=0.25) indicated a high level of explanatory power regarding the variation in resulting total discharge (Fig. 6). Comparing all seven input parameters, VS exhibited the highest R2 value for both peak and total discharge, signifying strong explanatory power regarding the modelled flow compared to the other parameters. Additionally, the rate of change in peak discharge and total discharge for every unit increase in VS was comparatively higher than other parameters as indicated by the steepness of the slope (Figs. 5a and 6a).

Basal friction angle and CFF demonstrated a high level of explanatory power concerning the variability in flow arrival time, as evidenced by their R2 values of 0.98 and 0.97, respectively. VS variation also exhibited a moderately high explanatory power with a negative relationship, with R2 of 0.81 and a slope (m) of −0.019. However, the linearity became less pronounced within the volume range of 4 × 106 to 10 × 106 m3. In contrast, other parameters, including CE and ϕ, did not exhibit a definitive linear relationship (Fig. 7). However, CE variation showed a threshold effect on arrival time; increasing CE from 10−5.95 to 10−6.42 kg−1 decreased arrival time, while further increases towards 10−5.85 kg−1 led to a linear increase in arrival time.

5 Discussion

Our primary aim was to investigate the sensitivity of the modelled GLOF outputs to a range of values for key model input parameters using the r.avaflow model. Previous studies have underscored the sensitivity of GLOF model outputs to some input parameters, including the basal friction angle, entrainment coefficient, and volume of mass movement entering the lake (Baggio et al., 2021; Mergili et al., 2018b, 2020a). This study advances our understanding of GLOF modelling by conducting a comprehensive sensitivity analysis across nine essential parameters and multiple (84) GLOF simulations. As a result, we have, for the first time, ranked these nine GLOF input parameters based on their contributions to model output variabilities. Our results showed that modelled GLOF output parameters are substantially sensitive to five of the nine parameters we tested here (DEM dataset, volume of mass movement entering the lake, origin of mass movement entering the lake, entrainment coefficient, and basal friction angle), suggesting that GLOF modelling results are subject to uncertainty from multiple sources. The findings offer valuable perspectives on the uncertainty in GLOF modelling results and complexities inherent in modelling the GLOF process chain within rugged mountain terrain such as in the Himalaya.

https://nhess.copernicus.org/articles/25/1841/2025/nhess-25-1841-2025-f09

Figure 9A comparison of the elevation profiles from four DEM datasets and the corresponding flow depths along the river centreline, generated through r.avaflow modelling. Panels (a) and (b) show the flow depths and elevation profiles along the river centreline. Panels (c) and (d) illustrate elevation profiles for two specific sections marked in (a) and (b). DEM and flow depth profiles from the vertical cross-sections at various distances are also provided in the Supplement (Fig. S7). The DEM datasets were co-registered using the algorithm developed by Shean et al. (2016).

Download

5.1 DEM datasets

A DEM is one essential piece of data for modelling GLOFs and other types of floods (Hawker et al., 2018; Saksena and Merwade, 2015; Schumann and Bates, 2018; Westoby et al., 2014). The impact of DEM resolution is even more pronounced in the steep and complex topographic conditions prevalent in high-mountain regions like the Himalaya (Liu et al., 2019). Our study provides the quantification of the effect of DEM in such environments for the first time. Our results suggest that the use of global and regional DEM datasets with a resolution ranging from 8 m (HMA-DEM) to 90 m (SRTMGL3) leads to over 2-fold and 4-fold variations in peak and total discharge, respectively, and causes successive significant fluctuations in flood height, reach distance, and flow arrival time. This likely results from the low-resolution DEMs not fully resolving the river channel compared to higher-resolution DEMs. This limitation results in a reduction in surface roughness and river channel conveyance (carrying capacity of channel). Thus, the flow spreads out more, leading to an increase in the modelled flow extent and reach (Figs. 9 and S7) (Muthusamy et al., 2021). This was supported by comparing the DEM profile and flow depth along the river centreline (Fig. 9) and across the multiple vertical cross-sections along the river channel (Fig. S7). The analysis showed that GLOF output from SRTMGL3, where river channels are poorly resolved, was comparatively higher than that from HMA-DEM with the better-resolved channel. Also, DEM datasets were acquired at different times, meaning the topographic features they captured will also differ depending on natural geomorphological change or human-made alteration of the earth's surface over time (Khosh Bin Ghomash et al., 2019; Bishop et al., 2002; Schumann and Bates, 2018; Watson et al., 2015). Such discrepancies will have a substantial influence on flow characteristics and uncertainty if a DEM with inadequate spatial and temporal resolution is used in GLOF modelling.

Overestimation of flood maps stemming from reductions in DEM resolution has been reported in urban flood modelling (Muthusamy et al., 2021; McClean et al., 2020). However, the impact of DEM data on GLOF modelling, especially in a complex topographic setting such as in the Himalaya, has been rarely documented (Wang et al., 2011). Our results show the substantial variation in GLOF model output stemming from DEM dataset variation, even when employing DEMs with comparable spatial resolutions, underscoring the criticality of high-quality DEM data in GLOF modelling (Fig. 9). DEM datasets covering rugged high-mountain terrain, where GLOFs typically occur, are likely to have more errors due to geometric distortion and data loss due to challenges involved in data acquisition for DEM production (Hugonnet et al., 2021; Liu et al., 2019). Therefore, using global- to regional-scale open-access DEMs, such as SRTM for GLOF modelling, due to the absence of high-resolution alternatives (Wang et al., 2011) may only be suitable for a first-order assessment of GLOFs at large scales (Zhang et al., 2023b). This is important as uncertainty stemming from DEM datasets is often overlooked and/or not well addressed in the previous basin-specific GLOF modelling work (Rinzin et al., 2023; Sattar et al., 2023, 2021b).

5.2 Variation in the origin of mass movement

Our study indicated that different origins of mass movement initiation produced GLOFs with approximately 2-fold variations in their peak discharge, volume, and reach distance (Fig. 4). These variations can be explained based on the lake geometry and the direction, horizontal distance, and vertical distance at which the mass movement enters the lake. The r.avaflow model provides detailed output parameters such as kinetic energy associated with the flow and a flow depth map for each time step, which allowed us to better understand the cause of this variation. For example, the mass movement originating from Loc-1, which is located at the slope above the headwall, directly impacts the head end of the lake with the highest kinetic energy (50 714 GJ) of all sources of mass movement. This head-end impact, coupled with its high energy, facilitates the direct forward propagation of waves toward the frontal outlet, causing the lake water to overtop the frontal moraine and resulting in a higher peak and total discharge (Fig. S8). Thorthormi Tsho is roughly crescentic in shape and curves toward the west, with its maximum curvature facing the origin of mass movement of Loc-5. This shape also allows for the impact wave generated from mass movement from Loc-5 to move almost unimpeded along the flow line, resulting in greater GLOF discharge. Also, Loc-1 and Loc-5 are located at greater distances and higher elevations from the lake compared to other locations (Fig. 2). This spatial location might have also enabled them to generate greater force to impact the lake. In contrast, the direct wave of impact generated by the mass movement from Loc-3, located on the slope above the right moraine dam, is deflected towards the left lateral moraine and only a secondary wave proceeds towards the lake outlet, resulting in a comparatively lower peak and total discharge (Fig. S8) (Emmer et al., 2024). This finding implies that the geometry of the glacial lake and the surrounding source slope plays a vital role in GLOF output. Thus, we underscore the importance of considering catchment shape in GLOF modelling, although we cannot assume that two identical basins will have the same flood properties due to the influence of other factors, such as the involved volume of solid and fluid parts.

Earlier studies (Mergili et al., 2017, 2020b) have explained the interaction between landslides and reservoirs (lakes) and their influence on the resulting hydrograph. However, these studies did not consider the variables such as directions and distances from which the mass impacts the lake. To fill this gap, here we enhanced our understanding of the interplay between the resulting GLOF magnitude and mass movement attributes including the direction and distance from which the mass enters the lake, the amount of kinetic energy the mass possesses, and the geometry of the lake. Our results emphasize the significant impact on the resulting GLOF events caused by the uncertainty in pinpointing the specific location of the origin of mass movement entering the lake. The thawing of permafrost and destabilization of the slope surrounding the lake due to climate warming (Gruber et al., 2017; Kääb et al., 2018; Sattar et al., 2025) combined with the expansion of the glacial lake towards the mountain flank (Rounce et al., 2016) are likely to increase the frequency of mass movement entering the lake, further exacerbating this uncertainty. By effectively identifying and quantifying the uncertainties related to the origin of mass movement entering the lake in addition to the volume, we propose considering the origins of mass movement in designing the scenario-based GLOF modelling in future studies. This approach goes beyond the existing practice, which primarily focuses on the volume of mass movement entering the lake. This is important because, as demonstrated in our study, the impact of mass movement can vary significantly depending on its origin. In some cases, what may be considered a worst-case scenario from one origin might represent only a low-magnitude event from another (GAPHAZ, 2017; Sattar et al., 2021a; Allen et al., 2022).

5.3 Mass movement volume, grain density, entrainment coefficient, and lake volume

Our investigation revealed that variation in GLOF magnitude is most sensitive to the volume of mass movement (VS) entering the lake. It also exhibits a significant level of sensitivity to the entrainment coefficient (CE), whilst the grain density (ρS) and lake volume (VL) exhibit a negligible impact. For example, the variation in VS between 1×106 and 10×106 m3 leads to peak and total discharge fluctuation of 1160 % and 2500 %, respectively, and subsequent variation in maximum flow depth and arrival time (Figs. 5 to 7). The dominant impact of VS and CE on GLOF magnitude could be due to their direct influence on the overall magnitude and intensity of flood events. The total discharge during the GLOF cascade event is a function of VS entering the lake. This is further corroborated by the near-perfect linear relationship of peak discharge (R2=0.99) and total discharge (R2=1) with VS entering the lake observed here. Likewise, the volume of solid content in the flow is solely contributed by the entrainment of frontal moraine material, primarily determined by CE. Additionally, this correlation could be attributed to the amount of energy and associated momentum of the flow, which increases significantly with a corresponding increase in VS. Also, it could be due to the longer timing and duration of the flow as evident in Fig. S2. Most GLOF events in high mountains across HMA and other alpine regions are caused by moraine dam breaches triggered by mass movement entering the lake from the surrounding mountain flank (Shrestha et al., 2023; Lützow et al., 2023; Emmer and Vilímek, 2014). As a result, VS is considered a primary basis for scenario development (Allen et al., 2022). Thus, we believe this finding provides useful insights for improving the development of different scenarios of GLOFs with higher confidence or is a basis for ensemble testing, with the caveat that the range of outputs may be too wide to be of practical use.

The volume of a glacial lake is one of the key parameters in GLOF modelling and subsequent hazard mapping. Without field-measured bathymetry, we applied eight different empirical equations to estimate the volume of Thorthormi Tsho, which ranged between 2.1 × 108 and 3.80 × 108 m3. Interestingly, this variation in VL did not significantly influence the modelled GLOF hydraulics; VL was ranked as the third least sensitive parameter in our analysis. This limited impact may be attributed to the large size of Thorthormi Tsho ( 4.35 km3) and its substantial storage capacity, as estimated by the empirical equations used for calculating lake volume here. Even the lower bound of the calculated VL exceeds 95 % of the measured volume of the glacial lake, and the upper bound is greater than 58 of the 59 measured lake volumes in the greater Himalaya (Zhang et al., 2023a). We believe that Thorthormi Tsho's large volume and expansive area enable it to absorb and dissipate the energy from potential mass movement impacts, effectively acting as a buffer (Emmer et al., 2024). Nevertheless, while our findings suggest that VL uncertainty may have a relatively minor impact on GLOF modelling for large lakes like Thorthormi Tsho, we acknowledge that VL remains a critical parameter in GLOF modelling, particularly for smaller lakes where volume changes could have a more pronounced effect on flood dynamics (Sattar et al., 2023).

5.4 Variations in frictional parameters

Among the frictional parameters, our result showed that GLOF magnitude is most sensitive to the basal friction angle (δ). For example, the variation in total discharge (47.5 %) resulting from fluctuation in δ within the experiment range was 30 times greater than that of the internal friction angle (ϕ) and over 4 times greater than that of the fluid friction angle (CFF). δ plays a dominant role in flow dynamics and the interaction between the flowing material and the channel bed. This direct contact means that even minor changes in δ can have substantial effects on the resistance encountered by the flowing material, thereby influencing the mobility of the flow (Pudasaini and Krautblatter, 2014a; Mergili et al., 2018a, b). ϕ on the other hand primarily affects particle interactions within the flowing material, whilst CFF is a coefficient which quantifies the overall flow resistance within the flow path, mainly depending on surface roughness. Our findings indicate that prioritizing the consideration of δ over the other two frictional parameters is advisable. While the back-calculated values might seem to be reasonable initiation values for δ as measuring it in the field is practically not feasible, we recommend conducting a statistically substantial sensitivity analysis using adequate sampling size and an appropriate statistical model. Despite the relatively low overall impact on GLOF magnitude, CFF notably increased the flow's mobility, especially beyond 12 km downstream, when the flow became fluid-dominant (Fig. S4). This is because CFF controls the mobility of the fluid part (Mergili and Pudasaini, 2024; Mergili et al., 2017). This suggests that CFF could exert a substantial influence, particularly in modelling scenarios encompassing longer flow distances.

5.5 Key points from the comparison of all parameters and the way forward

Identifying the most accurate parameter values or optimal datasets can be achieved through validation with well-constrained historical events (Zheng et al., 2021a; Schneider et al., 2014; Mergili et al., 2020b; Shugar et al., 2021; Sattar et al., 2025), but there are limitations in the transferability of these findings due to unique characteristics and initial conditions of each GLOF, such as varying proportions of solid and liquid parts. Because of this we effectively constrained how variation in input parameter values within commonly used ranges influences the GLOF modelling results, instead of trying to determine the correct value for each parameter. Or, put another way, if we apply our study to a specific example, we may determine that certain factors are more important than others for this specific example, but it would be unclear how applicable our results are to other events. Thus, we see our approach as the least biased towards any event and hence the most generally applicable approach. Our results indicating modelled GLOFs being subjected to uncertainty from multiple input parameters implies that the parameter value determined from one modelled GLOF event might not accurately predict the behaviour of GLOFs in different regions or under different circumstances (Mergili et al., 2018a, 2020b). Therefore, while these back-analysed parameter values can provide valuable insights, they need to be applied with caution and adapted to the specific context of each new GLOF scenario. As a result of these multiple sources of uncertainty in modelled GLOFs, it could pose challenges in effectively communicating risks with communities and other stakeholders (Thompson et al., 2020). We highlight that more sensitive parameters should be treated carefully in future GLOF modelling works by robustly considering associated uncertainties.

Due to the high sensitivity of the model output on DEM resolution, we emphasize the critical importance of a high-resolution and good-quality DEM (Uuemaa et al., 2020; Schumann and Bates, 2018), especially when modelling is aimed at producing hazard maps with higher granularity at the specific basin scale. DEMs should have a high spatial resolution, have a high vertical accuracy, and be recently produced, especially in areas of high relief and rapid landscape change such as in the Himalaya (Schumann and Bates, 2018). Previous studies have indicated that flood modelling accuracy can be improved by addressing the effect of DEM resolution and accuracy (Saksena and Merwade, 2015) or by merging with other high-resolution and accurate DEMs (Muthusamy et al., 2021). These methods appear viable in the context of highly sparse coverage of high-resolution DEMs and the unaffordability of high-resolution commercial DEMs, but the modelling results may be subjected to substantial uncertainty.

The mass movement volume and δ exhibit a strong linear relationship with all output parameters. Whilst the linear relationship does not negate the influence these parameters have on flow characteristics, it suggests that model output errors resulting from uncertainties in these parameters might be predictable. This is essential since predicting the VS movement involved in the forward modelling is highly challenging and determining an accurate value is impossible – the current challenge is rather to establish a likely envelope of volumes. However, such a prediction should be bespoke to the particular events based on the initial parameters like estimated ice thickness, slope, and the presence of permafrost. Furthermore, such predictions must also consider other factors, such as equifinality arising from the interaction of multiple parameters (Mergili et al., 2018a, b, 2020b).

The entrainment coefficient exhibits a threshold effect for peak discharge, arrival time, and total discharge, but these thresholds occur at different CE values for each of these output parameters. The observed threshold effects in flow characteristics can likely be attributed to the initial dominance of the fluid component, wherein the contribution of erosion to the flow dynamics is negligible. As the erosion rate progresses with increasing CE values, the increasing concentration of eroded material reduces the mobility of the flow, leading to a decline in mobility of the flow. These threshold effects suggest that once these thresholds are surpassed, the resulting flow and arrival time demonstrate a heightened sensitivity to variations in entrainment. Consequently, this sensitivity may translate to the flow characteristics such as flow depth and arrival time which are essential for hazard and risk assessments. It is important to note that this threshold value may vary across different GLOF events due to the diverse combinations of other parameters.

The GLOF simulations were conducted using the r.avaflow model due to its capability to model the entire GLOF process chain (Mergili and Pudasaini, 2024; Mergili et al., 2017). While we present the uncertainty involved in the full GLOF process chain from mass movement entering the lake to downstream propagation, we specifically explored the uncertainty in the GLOF input parameters relevant to r.avaflow modelling. Input parameters such as DEM datasets and the volume and density of mass movement involved in a GLOF event might be similar across different models. However, it should be noted that all parameters tested here do not necessarily apply to all models used for GLOF modelling.

The flow arrival time was measured from the profile located 3, 6, and 9 km downstream of the lake since some of our modelled GLOF terminates before proceeding further downstream. This is a reasonable point as human settlement downstream of the lake is mostly concentrated around this area. The variation in flow arrival time might be underestimated if the location is farther downstream from the lake.

Here we focused on nine essential parameters in r.avaflow, which are relevant to GLOF modelling. As an open-source modelling code, the r.avaflow model offers flexibility by allowing the users to manipulate all parameters, a level of transparency that sets it apart from many modelling codes. However, including inbuilt modules, initial conditions, and all flow parameters, r.avaflow has more than 30 tuneable parameters (Mergili and Pudasaini, 2024) (Table S1). Our choice of parameter in this study is essentially motivated by those identified as critical and commonly modified in the previous GLOF modelling studies and which are frequently modified to fit with the parameter from back-calculated events. Because of this, our sensitivity analysis might have potentially overlooked the complexity of r.avaflow stemming from the effect of all these parameters, let alone their interaction effect. Also due to the enormous computational time and resources required for running r.avaflow, we used one-at-a-time sensitivity analysis conducting up to 10 simulations per parameter. While this number of simulations for each parameter produced substantially conclusive results, we do not discount the robustness of the global sensitivity analysis employing an adequate sampling size which can adequately account for the influence of parameter interactions (Saltelli et al., 2004). While we consider for the first time the importance of the uncertainty in these parameters in GLOF modelling, in the light of all these limitations, we provide the following recommendations. (1) We recommend that future studies include a broader range of parameters, further advancing the understanding of their roles in GLOF modelling. (2) This can be achieved by initially focusing on the most sensitive parameters identified in this study and progressively expanding the analysis to include other parameters. Adequate sample sizes and appropriate statistical models should be employed to ensure robust findings. (3) Given that r.avaflow is an open-source tool under continuous development, we underscore the implementation of features that support parallel computation for parameter interaction analyses. This advancement would enable future studies to fully utilize high-performance computing (HPC) resources and address complex modelling challenges more effectively.

6 Conclusions

GLOFs present substantial dangers to communities residing in valleys downstream of glacial lakes. GLOFs involve complex cascading processes and typically occur across rugged mountain terrains. Due to these complexities, modelling GLOFs necessitates extensive input data, parameters, and complex modelling codes for accurate hazard and risk assessments, which is inherently challenging. However, previous studies have mostly relied on open-access data and are grounded in a historical event introducing significant uncertainties to the modelling results. In this study, we have, for the first time, conducted a sensitivity analysis considering multiple GLOF parameters and ranked these inputs based on how their uncertainties in input values apportion to the variation in modelling output by employing cutting-edge modelling code, r.avaflow. Our results suggested GLOF modelling outputs such as peak and total discharge are substantially sensitive to variation in input values of six out of nine parameters we tested here. Specifically, the modelling outputs are the most sensitive to the volume of mass movement entering the lake followed by the variation in DEM datasets and the location of the origin of mass movement entering the lake. Other parameters like the basal friction angle and entrainment coefficient also showed significant sensitivity. Although limited to GLOF modelling with the r.avaflow model, our study emphasizes that GLOF modelling results are influenced by uncertainties stemming from various sources, underscoring the need for careful interpretation of the modelling results. By ranking the model parameters according to their impact on model output, our study prioritizes model input parameters for future modelling efforts, given the challenge of adequately constraining multiple parameters. Additionally, this study lays the groundwork for a thorough investigation into the most sensitive parameters to improve our understanding of GLOF modelling.

Code and data availability

The r.avaflow modelling code we used here for simulating all scenarios of GLOFs can be accessed at https://www.landslidemodels.org/r.avaflow/ (Mergili and Pudasaini, 2024). SRTMGL3 can be downloaded from https://doi.org/10.5069/G9445JDF (SRTM, 2013), NASADEM can be downloaded from https://doi.org/10.5069/G93T9FD9 (NASA JPL, 2021), and AW3D30 can be downloaded from OpenTopogragphy https://doi.org/10.5069/G94M92HB (JAXA, 2021). HMA-DEM can be downloaded from the National Snow and Ice Data Center at https://doi.org/10.5067/KXOVQ9L172S2 (Shean, 2017).

Supplement

The supplement related to this article is available online at https://doi.org/10.5194/nhess-25-1841-2025-supplement.

Author contributions

SR, SD, and RJC conceptualized the study. SR undertook data processing and analysis. SD and RJC secured the funding, supervised, and contributed equally to the work. AS helped with modelling and was involved as an external supervisor. MM revised and provided an expert opinion on the study. All authors wrote and revised the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We thank Sonam Wangchuk of ICIMOD for helping us set up HPC experiments.

Financial support

This work was supported by the IAPETUS Doctoral Training Partnership (grant no. IAP2-21-267) funded by the Natural Environment Research Council (NERC).

Review statement

This paper was edited by Daniele Giordan and reviewed by Takashi Kimura and one anonymous referee.

References

Abdi, H.: Coefficient of variation, Encyclopedia of research design, 1, https://doi.org/10.4135/9781412961288.n56, 2010. 

Allen, S. K., Rastner, P., Arora, M., Huggel, C., and Stoffel, M.: Lake outburst and debris flow disaster at Kedarnath, June 2013: hydrometeorological triggering and topographic predisposition, Landslides, 13, 1479–1491, https://doi.org/10.1007/s10346-015-0584-3, 2015. 

Allen, S. K., Linsbauer, A., Randhawa, S. S., Huggel, C., Rana, P., and Kumari, A.: Glacial lake outburst flood risk in Himachal Pradesh, India: an integrative and anticipatory approach considering current and future threats, Nat. Hazards, 84, 1741–1763, https://doi.org/10.1007/s11069-016-2511-x, 2016. 

Allen, S. K., Zhang, G., Wang, W., Yao, T., and Bolch, T.: Potentially dangerous glacial lakes across the Tibetan Plateau revealed using a large-scale automated assessment approach, Sci. Bull., 64, 435–445, https://doi.org/10.1016/j.scib.2019.03.011, 2019. 

Allen, S. K., Sattar, A., King, O., Zhang, G., Bhattacharya, A., Yao, T., and Bolch, T.: Glacial lake outburst flood hazard under current and future conditions: worst-case scenarios in a transboundary Himalayan basin, Nat. Hazards Earth Syst. Sci., 22, 3765–3785, https://doi.org/10.5194/nhess-22-3765-2022, 2022. 

Baggio, T., Mergili, M., and D'Agostino, V.: Advances in the simulation of debris flow erosion: The case study of the Rio Gere (Italy) event of the 4th August 2017, Geomorphology, 381, 107664, https://doi.org/10.1016/j.geomorph.2021.107664, 2021. 

Bishop, M. P., Shroder, J. F., Bonk, R., and Olsenholler, J.: Geomorphic change in high mountains: a western Himalayan perspective, Global Planet. Change, 32, 311–329, https://doi.org/10.1016/S0921-8181(02)00073-5, 2002. 

Byers, A. C., Rounce, D. R., Shugar, D. H., Lala, J. M., Byers, E. A., and Regmi, D.: A rockfall-induced glacial lake outburst flood, Upper Barun Valley, Nepal, Landslides, 16, 533–549, https://doi.org/10.1007/s10346-018-1079-9, 2018. 

Carrivick, J. L. and Tweed, F. S.: A global assessment of the societal impacts of glacier outburst floods, Global Planet. Change, 144, 1–16, https://doi.org/10.1016/j.gloplacha.2016.07.001, 2016. 

Crippen, R., Buckley, S., Agram, P., Belz, E., Gurrola, E., Hensley, S., Kobrick, M., Lavalle, M., Martin, J., Neumann, M., Nguyen, Q., Rosen, P., Shimada, J., Simard, M., and Tung, W.: NASADEM GLOBAL ELEVATION MODEL: METHODS AND PROGRESS, Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci., XLI-B4, 125–128, https://doi.org/10.5194/isprs-archives-XLI-B4-125-2016, 2016. 

Dubey, S. and Goyal, M. K.: Glacial Lake Outburst Flood Hazard, Downstream Impact, and Risk Over the Indian Himalayas, Water Resour. Res., 56, e2019WR026533, https://doi.org/10.1029/2019wr026533, 2020. 

Dubey, S., Sattar, A., Goyal, M. K., Allen, S., Frey, H., Haritashya, U. K., and Huggel, C.: Mass Movement Hazard and Exposure in the Himalaya, Earth's Future, 11, e2022EF003253, https://doi.org/10.1029/2022ef003253, 2023. 

Emmer, A. and Cochachine, L.: The causes and mechanisms of moraine-dammed lake failures in the Cordillera blanca, North American Cordillera, and Himalayas, AUC Geographica, 48, 5–15, https://doi.org/10.14712/23361980.2014.23, 2013. 

Emmer, A. and Vilímek, V.: New method for assessing the susceptibility of glacial lakes to outburst floods in the Cordillera Blanca, Peru, Hydrol. Earth Syst. Sci., 18, 3461–3479, https://doi.org/10.5194/hess-18-3461-2014, 2014. 

Emmer, A., Vilca, O., Salazar Checa, C., Li, S., Cook, S., Pummer, E., Hrebrina, J., and Haeberli, W.: Causes, consequences and implications of the 2023 landslide-induced Lake Rasac GLOF, Cordillera Huayhuash, Peru, EGUsphere [preprint], https://doi.org/10.5194/egusphere-2024-2316, 2024. 

Farinotti, D., Huss, M., Fürst, J. J., Landmann, J., Machguth, H., Maussion, F., and Pandit, A.: A consensus estimate for the ice thickness distribution of all glaciers on Earth, Nat. Geosci., 12, 168–173, https://doi.org/10.1038/s41561-019-0300-3, 2019. 

Frey, H., Huggel, C., Chisolm, R. E., Baer, P., McArdell, B., Cochachin, A., and Portocarrero, C.: Multi-Source Glacial Lake Outburst Flood Hazard Assessment and Mapping for Huaraz, Cordillera Blanca, Peru, Front. Earth Sci., 6, https://doi.org/10.3389/feart.2018.00210, 2018. 

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

Gantayat, P., Sattar, A., Haritashya, U. K., Watson, C. S., and Kargel, J. S.: Bayesian Approach to Estimate Proglacial Lake Volume (BE-GLAV), Earth Space Sci., 11, e2024EA003542, https://doi.org/10.1029/2024ea003542, 2024. 

GAPHAZ: Assessment of Glacier and Permafrost Hazards in Mountain Regions. Technical Guidance Document, Standing Group on Glacier and Permafrost Hazards in Mountains (GAPHAZ) of the International Association of Cryospheric Sciences (IACS) and the International Permafrost Association (IPA), 72, https://doi.org/10.13140/RG.2.2.26332.90245, 2017. 

Gruber, S., Fleiner, R., Guegan, E., Panday, P., Schmid, M.-O., Stumm, D., Wester, P., Zhang, Y., and Zhao, L.: Review article: Inferring permafrost and permafrost thaw in the mountains of the Hindu Kush Himalaya region, The Cryosphere, 11, 81–99, https://doi.org/10.5194/tc-11-81-2017, 2017. 

Hawker, L., Bates, P., Neal, J., and Rougier, J.: Perspectives on Digital Elevation Model (DEM) Simulation for Flood Modeling in the Absence of a High-Accuracy Open Access Global DEM, Front. Earth Sci., 6, https://doi.org/10.3389/feart.2018.00233, 2018. 

Huggel, C.: Recent extreme slope failures in glacial environments: effects of thermal perturbation, Quaternary Sci. Rev., 28, 1119–1130, https://doi.org/10.1016/j.quascirev.2008.06.007, 2009. 

Hugonnet, R., McNabb, R., Berthier, E., Menounos, B., Nuth, C., Girod, L., Farinotti, D., Huss, M., Dussaillant, I., Brun, F., and Kaab, A.: Accelerated global glacier mass loss in the early twenty-first century, Nature, 592, 726–731, https://doi.org/10.1038/s41586-021-03436-z, 2021. 

Japan Aerospace Exploration Agency (JAXA): ALOS World 3D 30 meter DEM (V3.2), OpenTopography [data set], https://doi.org/10.5069/G94M92HB, 2021. 

Kääb, A., Leinss, S., Gilbert, A., Bühler, Y., Gascoin, S., Evans, S. G., Bartelt, P., Berthier, E., Brun, F., Chao, W.-A., Farinotti, D., Gimbert, F., Guo, W., Huggel, C., Kargel, J. S., Leonard, G. J., Tian, L., Treichler, D., and Yao, T.: Massive collapse of two glaciers in western Tibet in 2016 after surge-like instability, Nat. Geoscie., 11, 114–120, https://doi.org/10.1038/s41561-017-0039-7, 2018. 

Khosh Bin Ghomash, S., Caviedes-Voullieme, D., and Hinz, C.: Effects of erosion-induced changes to topography on runoff dynamics, J. Hydrol., 573, 811–828, https://doi.org/10.1016/j.jhydrol.2019.04.018, 2019. 

Lala, J. M., Rounce, D. R., and McKinney, D. C.: Modeling the glacial lake outburst flood process chain in the Nepal Himalaya: reassessing Imja Tsho's hazard, Hydrol. Earth Syst. Sci., 22, 3721–3737, https://doi.org/10.5194/hess-22-3721-2018, 2018. 

Linsbauer, A., Frey, H., Haeberli, W., Machguth, H., Azam, M. F., and Allen, S.: Modelling glacier-bed overdeepenings and possible future lakes for the glaciers in the Himalaya–Karakoram region, Ann. Glaciol., 57, 119–130, https://doi.org/10.3189/2016AoG71A627, 2016. 

Linsbauer, A., Paul, F., Machguth, H., and Haeberli, W.: Comparing three different methods to model scenarios of future glacier change in the Swiss Alps, Ann. Glaciol., 54, 241–253, https://doi.org/10.3189/2013AoG63A400, 2017. 

Liu, K., Song, C., Ke, L., Jiang, L., Pan, Y., and Ma, R.: Global open-access DEM performances in Earth's most rugged region High Mountain Asia: A multi-level assessment, Geomorphology, 338, 16–26, https://doi.org/10.1016/j.geomorph.2019.04.012, 2019. 

Lützow, N., Veh, G., and Korup, O.: A global database of historic glacier lake outburst floods, Earth Syst. Sci. Data, 15, 2983–3000, https://doi.org/10.5194/essd-15-2983-2023, 2023. 

McClean, F., Dawson, R., and Kilsby, C.: Implications of Using Global Digital Elevation Models for Flood Risk Analysis in Cities, Water Resour. Rese., 56, e2020WR028241, https://doi.org/10.1029/2020wr028241, 2020. 

Mergili, M. and Pudasaini, S. P.: r.avaflow – The mass flow simulation tool, r.avaflow.direct Web interface, https://www.landslidemodels.org/r.avaflow/ (last access: 16 May 2025), 2024. 

Mergili, M., Fischer, J.-T., Krenn, J., and Pudasaini, S. P.: r.avaflow v1, an advanced open-source computational framework for the propagation and interaction of two-phase mass flows, Geosci. Model Dev., 10, 553–569, https://doi.org/10.5194/gmd-10-553-2017, 2017. 

Mergili, M., Frank, B., Fischer, J.-T., Huggel, C., and Pudasaini, S. P.: Computational experiments on the 1962 and 1970 landslide events at Huascarán (Peru) with r.avaflow: Lessons learned for predictive mass flow simulations, Geomorphology, 322, 15–28, https://doi.org/10.1016/j.geomorph.2018.08.032, 2018a. 

Mergili, M., Emmer, A., Juricova, A., Cochachin, A., Fischer, J. T., Huggel, C., and Pudasaini, S. P.: How well can we simulate complex hydro-geomorphic process chains? The 2012 multi-lake outburst flood in the Santa Cruz Valley (Cordillera Blanca, Peru), Earth Surf. Proc. Land., 43, 1373–1389, https://doi.org/10.1002/esp.4318, 2018b. 

Mergili, M., Jaboyedoff, M., Pullarello, J., and Pudasaini, S. P.: Back calculation of the 2017 Piz Cengalo–Bondo landslide cascade with r.avaflow: what we can do and what we can learn, Nat. Hazards Earth Syst. Sci., 20, 505–520, https://doi.org/10.5194/nhess-20-505-2020, 2020a. 

Mergili, M., Pudasaini, S. P., Emmer, A., Fischer, J.-T., Cochachin, A., and Frey, H.: Reconstruction of the 1941 GLOF process chain at Lake Palcacocha (Cordillera Blanca, Peru), Hydrol. Earth Syst. Sci., 24, 93–114, https://doi.org/10.5194/hess-24-93-2020, 2020b. 

Muthusamy, M., Casado, M. R., Butler, D., and Leinster, P.: Understanding the effects of Digital Elevation Model resolution in urban fluvial flood modelling, J. Hydrol., 596, 126088, https://doi.org/10.1016/j.jhydrol.2021.126088, 2021. 

NASA JPL: NASADEM Merged DEM Global 1 arc second V001, OpenTopography [data set], https://doi.org/10.5069/G93T9FD9, 2021. 

NASA Shuttle Radar Topography Mission (SRTM): Shuttle Radar Topography Mission (SRTM) Global [data set], https://doi.org/10.5069/G9445JDF, 2013. 

NCHM: Reassessment of Potentially Dangerous Glacial Lakes in Bhutan, National Centre for Hydrology and Meteorology, Royal Government of Bhutan, National Center for Hydrology and Meteorology, Royal Government of Bhutan, PO Box: 2017, Thimphu, Bhutan, 54, https://www.nchm.gov.bt/attachment/ckfinder/userfiles/files/Re-assessment of Potentailly Dangerous Glacial Lakes.pdf (last access: 16 May 2025), 2019a. 

NCHM: Detailed assessment report on GLOF hazard from Thorthormi glacial lakes and asssociated glaciers, Thimphu, https://www.nchm.gov.bt/attachment/ckfinder/userfiles/files/Report_on_Thorthormi_assessment_2019.pdf (last access: 16 May 2025), 2019b. 

NCHM: Action Taken Report Thorthormi Flood Incident 30th October 2023, National Center for Meteorology and Hydrology, Thimphu, 17, https://www.nchm.gov.bt/attachment/ckfinder/userfiles/files/Action Taken Report Thorthormi Flood Incident 30th October 2023.pdf (last access: 16 May 2025), 2023. 

Nie, Y., Liu, Q., Wang, J., Zhang, Y., Sheng, Y., and Liu, S.: An inventory of historical glacial lake outburst floods in the Himalayas based on remote sensing observations and geomorphological analysis, Geomorphology, 308, 91–106, https://doi.org/10.1016/j.geomorph.2018.02.002, 2018. 

Nie, Y., Deng, Q., Pritchard, H. D., Carrivick, J. L., Ahmed, F., Huggel, C., Liu, L., Wang, W., Lesi, M., Wang, J., Zhang, H., Zhang, B., Lü, Q., and Zhang, Y.: Glacial lake outburst floods threaten Asia's infrastructure, Sci. Bull., 68, 1361–1365, https://doi.org/10.1016/j.scib.2023.05.035, 2023. 

Obu, J., Westermann, S., Bartsch, A., Berdnikov, N., Christiansen, H. H., Dashtseren, A., Delaloye, R., Elberling, B., Etzelmüller, B., Kholodov, A., Khomutov, A., Kääb, A., Leibman, M. O., Lewkowicz, A. G., Panda, S. K., Romanovsky, V., Way, R. G., Westergaard-Nielsen, A., Wu, T., Yamkhin, J., and Zou, D.: Northern Hemisphere permafrost map based on TTOP modelling for 2000–2016 at 1 km2 scale, Earth-Sci. Rev., 193, 299–316, https://doi.org/10.1016/j.earscirev.2019.04.023, 2019. 

Pudasaini, S. P.: A general two-phase debris flow model, J. Geophys. Res.-Earth, 117, https://doi.org/10.1029/2011jf002186, 2012. 

Pudasaini, S. P. and Krautblatter, M.: A two-phase mechanical model for rock-ice avalanches, J.Geophys. Res.-Earth, 119, 2272–2290, https://doi.org/10.1002/2014jf003183, 2014a. 

Pudasaini, S. P. and Krautblatter, M.: A two-phase mechanical model for rock-ice avalanches, J. Geophys. Res.-Earth, 119, 2272–2290, https://doi.org/10.1002/2014jf003183, 2014b. 

Pudasaini, S. P. and Mergili, M.: A Multi-Phase Mass Flow Model, J. Geophys. Res.-Earth, 124, 2920–2942, https://doi.org/10.1029/2019jf005204, 2019. 

Rinzin, S., Zhang, G., and Wangchuk, S.: Glacial Lake Area Change and Potential Outburst Flood Hazard Assessment in the Bhutan Himalaya, Front. Earth Sci., 9, https://doi.org/10.3389/feart.2021.775195, 2021. 

Rinzin, S., Zhang, G., Sattar, A., Wangchuk, S., Allen, S. K., Dunning, S., and Peng, M.: GLOF hazard, exposure, vulnerability, and risk assessment of potentially dangerous glacial lakes in the Bhutan Himalaya, J. Hydrol., 619, 129311, https://doi.org/10.1016/j.jhydrol.2023.129311, 2023. 

Rounce, D., Watson, C., and McKinney, D.: Identification of Hazard and Risk for Glacial Lakes in the Nepal Himalaya Using Satellite Imagery from 2000–2015, Remote. Sens., 9, 654, https://doi.org/10.3390/rs9070654, 2017. 

Rounce, D. R., McKinney, D. C., Lala, J. M., Byers, A. C., and Watson, C. S.: A new remote hazard and risk assessment framework for glacial lakes in the Nepal Himalaya, Hydrol. Earth Syst. Sci., 20, 3455–3475, https://doi.org/10.5194/hess-20-3455-2016, 2016. 

Saksena, S. and Merwade, V.: Incorporating the effect of DEM resolution and accuracy for improved flood inundation mapping, J. Hydrol., 530, 180–194, https://doi.org/10.1016/j.jhydrol.2015.09.069, 2015. 

Saltelli, A., Tarantola, S., Campolongo, F., and Ratto, M.: Sensitivity analysis in practice: a guide to assessing scientific models, Wiley Online Library, https://doi.org/10.1002/0470870958, 2004. 

Sattar, A., Haritashya, U. K., Kargel, J. S., Leonard, G. J., Shugar, D. H., and Chase, D. V.: Modeling lake outburst and downstream hazard assessment of the Lower Barun Glacial Lake, Nepal Himalaya, J. Hydrol., 598, 126208, https://doi.org/10.1016/j.jhydrol.2021.126208, 2021a. 

Sattar, A., Goswami, A., Kulkarni, A. V., Emmer, A., Haritashya, U. K., Allen, S., Frey, H., and Huggel, C.: Future Glacial Lake Outburst Flood (GLOF) hazard of the South Lhonak Lake, Sikkim Himalaya, Geomorphology, 388, 107783, https://doi.org/10.1016/j.geomorph.2021.107783, 2021b. 

Sattar, A., Allen, S., Mergili, M., Haeberli, W., Frey, H., Kulkarni, A. V., Haritashya, U. K., Huggel, C., Goswami, A., and Ramsankaran, R.: Modeling Potential Glacial Lake Outburst Flood Process Chains and Effects From Artificial Lake-Level Lowering at Gepang Gath Lake, Indian Himalaya, J. Geophys. Res.-Earth, 128, e2022JF006826, https://doi.org/10.1029/2022jf006826, 2023. 

Sattar, A., Cook, K. L., Rai, S. K., Berthier, E., Allen, S., Rinzin, S., Van Wyk de Vries, M., Haeberli, W., Kushwaha, P., Shugar, D. H., Emmer, A., Haritashya, U. K., Frey, H., Rao, P., Gurudin, K. S. K., Rai, P., Rajak, R., Hossain, F., Huggel, C., Mergili, M., Azam, M. F., Gascoin, S., Carrivick, J. L., Bell, L. E., Ranjan, R. K., Rashid, I., Kulkarni, A. V., Petley, D., Schwanghart, W., Watson, C. S., Islam, N., Gupta, M. D., Lane, S. N., and Bhat, S. Y.: The Sikkim flood of October 2023: Drivers, causes and impacts of a multihazard cascade, Science, 0, eads2659, https://doi.org/10.1126/science.ads2659, 2025. 

Schneider, D., Huggel, C., Cochachin, A., Guillén, S., and García, J.: Mapping hazards from glacier lake outburst floods based on modelling of process cascades at Lake 513, Carhuaz, Peru, Adv. Geosci., 35, 145–155, https://doi.org/10.5194/adgeo-35-145-2014, 2014. 

Schumann, G. J. P. and Bates, P. D.: The Need for a High-Accuracy, Open-Access Global DEM, Front. Earth Sci., 6, https://doi.org/10.3389/feart.2018.00225, 2018. 

Shean, D.: High Mountain Asia 8-meter DEM mosaics derived from optical imagery (1), NASA National Snow and Ice Data Center Distributed Active Archive Center [data set], https://doi.org/10.5067/KXOVQ9L172S2, 2017. 

Shean, D. E., Alexandrov, O., Moratto, Z. M., Smith, B. E., Joughin, I. R., Porter, C., and Morin, P.: An automated, open-source pipeline for mass production of digital elevation models (DEMs) from very-high-resolution commercial stereo satellite imagery, ISPRS J. Photogramm., 116, 101–117, https://doi.org/10.1016/j.isprsjprs.2016.03.012, 2016. 

Shrestha, F., Steiner, J. F., Shrestha, R., Dhungel, Y., Joshi, S. P., Inglis, S., Ashraf, A., Wali, S., Walizada, K. M., and Zhang, T.: A comprehensive and version-controlled database of glacial lake outburst floods in High Mountain Asia, Earth Syst. Sci. Data, 15, 3941–3961, https://doi.org/10.5194/essd-15-3941-2023, 2023. 

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

Shugar, D. H., Jacquemart, M., Shean, D., Bhushan, S., Upadhyay, K., Sattar, A., Schwanghart, W., McBride, S., de Vries, M. V. W., Mergili, M., Emmer, A., Deschamps-Berger, C., McDonnell, M., Bhambri, R., Allen, S., Berthier, E., Carrivick, J. L., Clague, J. J., Dokukin, M., Dunning, S. A., Frey, H., Gascoin, S., Haritashya, U. K., Huggel, C., Kääb, A., Kargel, J. S., Kavanaugh, J. L., Lacroix, P., Petley, D., Rupper, S., Azam, M. F., Cook, S. J., Dimri, A. P., Eriksson, M., Farinotti, D., Fiddes, J., Gnyawali, K. R., Harrison, S., Jha, M., Koppes, M., Kumar, A., Leinss, S., Majeed, U., Mal, S., Muhuri, A., Noetzli, J., Paul, F., Rashid, I., Sain, K., Steiner, J., Ugalde, F., Watson, C. S., and Westoby, M. J.: A massive rock and ice avalanche caused the 2021 disaster at Chamoli, Indian Himalaya, Science, 373, 300–306, https://doi.org/10.1126/science.abh4455, 2021. 

Somos-Valenzuela, M. A., McKinney, D. C., Byers, A. C., Rounce, D. R., Portocarrero, C., and Lamsal, D.: Assessing downstream flood impacts due to a potential GLOF from Imja Tsho in Nepal, Hydrol. Earth Syst. Sci., 19, 1401–1412, https://doi.org/10.5194/hess-19-1401-2015, 2015. 

Taylor, C., Robinson, T. R., Dunning, S., Rachel Carr, J., and Westoby, M.: Glacial lake outburst floods threaten millions globally, Nat. Commun., 14, 487, https://doi.org/10.1038/s41467-023-36033-x, 2023. 

Thompson, I., Shrestha, M., Chhetri, N., and Agusdinata, D. B.: An institutional analysis of glacial floods and disaster risk management in the Nepal Himalaya, Int. J. Disast. Risk Re., 47, 101567, https://doi.org/10.1016/j.ijdrr.2020.101567, 2020. 

Uuemaa, E., Ahi, S., Montibeller, B., Muru, M., and Kmoch, A.: Vertical Accuracy of Freely Available Global Digital Elevation Models (ASTER, AW3D30, MERIT, TanDEM-X, SRTM, and NASADEM), Remote Sens., 12, 3482, https://doi.org/10.3390/rs12213482, 2020. 

Vilca, O., Mergili, M., Emmer, A., Frey, H., and Huggel, C.: The 2020 glacial lake outburst flood process chain at Lake Salkantaycocha (Cordillera Vilcabamba, Peru), Landslides, 18, 2211–2223, https://doi.org/10.1007/s10346-021-01670-0, 2021. 

Wang, W., Yang, X., and Yao, T.: Evaluation of ASTER GDEM and SRTM and their suitability in hydraulic modelling of a glacial lake outburst flood in southeast Tibet, Hydrol. Process., 26, 213–225, https://doi.org/10.1002/hyp.8127, 2011. 

Wang, X., Guo, X., Yang, C., Liu, Q., Wei, J., Zhang, Y., Liu, S., Zhang, Y., Jiang, Z., and Tang, Z.: Glacial lake inventory of high-mountain Asia in 1990 and 2018 derived from Landsat images, Earth Syst. Sci. Data, 12, 2169–2182, https://doi.org/10.5194/essd-12-2169-2020, 2020. 

Wang, Y., Hutter, K., and Pudasaini, S. P.: The Savage-Hutter theory: A system of partial differential equations for avalanche flows of snow, debris, and mud, ZAMM, Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 84, 507–527, https://doi.org/10.1002/zamm.200310123, 2004. 

Watson, C. S., Carrivick, J., and Quincey, D.: An improved method to represent DEM uncertainty in glacial lake outburst flood propagation using stochastic simulations, J. Hydrol., 529, 1373–1389, https://doi.org/10.1016/j.jhydrol.2015.08.046, 2015. 

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

Westoby, M. J., Brasington, J., Glasser, N. F., Hambrey, M. J., Reynolds, J. M., Hassan, M. A. A. M., and Lowe, A.: Numerical modelling of glacial lake outburst floods using physically based dam-breach models, Earth Surf. Dynam., 3, 171–199, https://doi.org/10.5194/esurf-3-171-2015, 2015. 

Worni, R., Stoffel, M., Huggel, C., Volz, C., Casteller, A., and Luckman, B.: Analysis and dynamic modeling of a moraine failure and glacier lake outburst flood at Ventisquero Negro, Patagonian Andes (Argentina), J. Hydrol., 444–445, 134–145, https://doi.org/10.1016/j.jhydrol.2012.04.013, 2012. 

Worni, R., Huggel, C., and Stoffel, M.: Glacial lakes in the Indian Himalayas–from an area-wide glacial lake inventory to on-site and modeling based risk assessment of critical glacial lakes, Sci. Total Environ., 468–469 S71–S84, https://doi.org/10.1016/j.scitotenv.2012.11.043, 2013. 

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

Zhang, G., Bolch, T., Yao, T., Rounce, D. R., Chen, W., Veh, G., King, O., Allen, S. K., Wang, M., and Wang, W.: Underestimated mass loss from lake-terminating glaciers in the greater Himalaya, Nat. Geosci., 16, 333–338, https://doi.org/10.1038/s41561-023-01150-1, 2023a. 

Zhang, T., Wang, W., An, B., and Wei, L.: Enhanced glacial lake activity threatens numerous communities and infrastructure in the Third Pole, Nat. Commun., 14, https://doi.org/10.1038/s41467-023-44123-z, 2023b. 

Zhang, T., Wang, W., and An, B.: Heterogeneous changes in global glacial lakes under coupled climate warming and glacier thinning, Commun. Earth Environ., 5, 374, https://doi.org/10.1038/s43247-024-01544-y, 2024. 

Zheng, G., Mergili, M., Emmer, A., Allen, S., Bao, A., Guo, H., and Stoffel, M.: The 2020 glacial lake outburst flood at Jinwuco, Tibet: causes, impacts, and implications for hazard and risk assessment, The Cryosphere, 15, 3159–3180, https://doi.org/10.5194/tc-15-3159-2021, 2021a. 

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

Download
Short summary
We modelled multiple glacial lake outburst flood (GLOF) scenarios (84 simulations) and tested the effect of nine key input parameters on the modelling results using r.avaflow. Our results highlight that GLOF modelling results are subject to uncertainty from the multiple input parameters. The variation in the volume of mass movement entering the lake causes the highest uncertainty in the modelled GLOF, followed by the DEM dataset and the origin of mass movement.
Share
Altmetrics
Final-revised paper
Preprint