Future glacial lakes in High Mountain Asia: an inventory and assessment of hazard potential from surrounding slopes

Abstract

Bedrock overdeepenings exposed by continued glacial retreat can store precipitation and meltwater, potentially leading to the formation of new proglacial lakes. These lakes may pose threats of glacial lake outburst floods (GLOFs) in high mountain areas, particularly if new lakes form in geomorphological setups prone to triggering events such as landslides or moraine collapses. We present the first complete inventory for future glacial lakes in High Mountain Asia by computing the subglacial bedrock for ~100 000 glaciers and estimating overdeepening area, volume and impact hazard for the larger potential lakes. We detect 25 285 overdeepenings larger than 10⁴ m² with a volume of 99.1 ± 28.6 km³ covering an area of 2683 ± 773.8 km². For the 2700 overdeepenings larger than 10⁵ m², we assess the lake predisposition for mass-movement impacts that could trigger a GLOF by estimating the hazard of material detaching from surrounding slopes. Our findings indicate a shift in lake area, volume and GLOF hazard from the southwestern Himalayan region toward the Karakoram. The results of this study can be used for anticipating emerging threats and potentials connected to glacial lakes and as a basis for further studies at suspected GLOF hazard hotspots.

1. Introduction

Globally rising temperatures cause increased glacial melting in nearly all glaciated regions of the world (Hock and others, 2019b). This holds especially true for the mountain ranges of High Mountain Asia (HMA) (Shean and others, 2020). One of the most dynamic and apparent impacts of a warming climate is the expansion of glacial lakes (Haritashya and others, 2018; Farinotti and others, 2019b). However, glacial melting processes are not only connected to the expansion of existing glacial lakes but also to the formation of new ones. Paraglacial bedrock troughs and laterofrontal moraine complexes are abundant in many deglaciating mountain ranges, often providing morphological configurations capable of storing precipitation and meltwater (Cuellar and McKinney, 2017). Newly formed moraine dams consist of unconsolidated sediments that were typically deposited along a steep ice front that recently melted away and are thus highly unstable. Dam failures can be triggered by thawing permafrost or by displacement waves from rockfall, ice avalanches or calving (Kershaw and others, 2005; Harrison and others, 2018). Both types of events can lead to reduced material cohesion and/or subsequent dam erosion (Eriksson and others, 2009). The following sudden and often catastrophic discharge of water characteristically transforms into a debris flow, gaining momentum and range, sometimes even exceeding 100 km, and involving valley damming and secondary outbursts (Lliboutry and others, 1977; Allen and others, 2016). Hence, glacial lake outburst floods (GLOFs) are among the most dangerous natural hazards in high mountain areas and pose an immense threat to people and infrastructures regionally (Clague and O'Connor, 2015; Carrivick and Tweed, 2016). Glaciological and hydrological information about glacial lakes is crucial for early hazard detection and effective risk reduction (Haeberli and others, 2016). To date, most research focuses on evaluating the threat from existing glacial lakes by means of remote sensing (Veh and others, 2018; Scapozza and others, 2019; Treichler and others, 2019; Shugar and others, 2020; Zhang and others, 2020). Goals were, for example, the development of early warning systems (Huggel and others, 2020), decision making strategies (Cuellar and McKinney, 2017) or inventories of lakes and previous flood events (Nagai and others, 2017; Nie and others, 2018).

Over the last few years, another research perspective emerged focusing on anticipating the development of future lakes as increasing temperatures lead to permafrost degradation and increasing slope instability (Deline and others, 2015). These future glacial lakes will most likely differ from existing ones regarding several aspects: first, they will rather form in bedrock overdeepenings in ice-free glacier beds than behind large moraine dams which would require a stationary glacier to build-up (Frey and others, 2010; Farinotti and others, 2019b). Second, they are likely to be located in the vicinity of increasingly unstable slopes (Schaub and others, 2013). Third, they will probably have a higher potential for disasters because of the continued expansion of people and infrastructure into higher elevations (GAPHAZ, 2017). Linsbauer and others (2012) developed the widely used GlabTop model to identify potential sites of lake formation. To date, variations of this method have been applied to individual mountain ranges, such as the European Alps (Haeberli and others, 2016), the Himalayan and Karakoram ranges (Allen and others, 2016; Linsbauer and others, 2016), or the Peruvian Andes (Colonia and others, 2017; Drenkhan and others, 2019). However, we are not aware of any inventory of future glacial lakes covering all mountain ranges of HMA, neither do we know about a future glacial lake hazard analysis of the same regional scope.

With the current study we provide the first inventory of future glacial lakes with associated hazard analyses on a regional-to-continental scale. Our study area covers the entire HMA, including the mountain ranges of the Himalaya, Karakoram and Kunlun, the Hindu Kush, Pamir and Tien Shan. We use global ice thickness data provided by Farinotti and others (2019a) to create a DEM of the study area ‘without glaciers’. Based on this, we calculate subglacial bedrock morphology which is then used to identify potential locations for lake development. As a result, we provide information on location, area and volume of potential future lakes for ~100 000 glaciers in HMA. We categorize the slopes surrounding the overdeepenings according to their geomorphological and hydrological attributes and obtain an approximate assessment of the hazard of rock or ice avalanches impacting the lake. Subsequently, we classify each lake accordingly, from highest to lowest impact hazard. We combine these results in a synoptic analysis to provide a first assessment of mass-movement impact hazards anticipated at future glacial lakes in the entire HMA. Due to the complex interactions between surface processes and landforms, modeling future lakes in high mountain areas is inherently subject to large uncertainties regarding, e.g. lake location and morphology, moraine dam height and water level (GAPHAZ, 2017). Therefore, the goal of this study is not a final assessment of glacial lake hazards but a comprehensive first overview. We aim at narrowing down potentially dangerous future lakes as a basis for more extensive investigations at the respective locations. It is with such early anticipation that planning and adaptation responses can be improved and risk-reducing measures can be defined.

2. Study area

This study focuses on the mountain ranges of HMA – the most glaciated region in the world aside from the polar ice caps (Vaughan and others (2013)). The Randolph Glacier Inventory, version 6.0 (RGI v6), lists 95 536 glaciers in HMA that cover an area of 97 605 km² (RGI Consortium, 2017). Despite the fact that HMA comprises extensive glaciers and glacier systems (e.g. the Fedchenko Glacier in the Pamir or the Siachen, the Baltoro and the Biafo Glaciers in the Karakoram), only ~1.5% of the inventoried glaciers are larger than 10 km² and ~85% are smaller than 1 km² (RGI Consortium, 2017). The RGI separates the glaciers of HMA into three zones: Central Asia (zone 13), South West Asia (zone 14) and South East Asia (zone 15) (Pfeffer and others, 2014). However, to be able to conduct comparative analyses of regional hydro- and topoclimatic regimes, a more detailed dataset is required. Here, we use glacier region outlines by Loibl (2020) which are based on second-order glacier regions provided by the Global Terrestrial Network for Glaciers (GTN-G, 2017). The aim of this new dataset is to refine GTN-G regions to facilitate regional glaciological and climatological analyses, e.g. for patterns in forcing and changes. This is achieved by delineating individual glaciated (sub)mountain ranges using major valleys for separation. The resulting dataset thus allows for analyses of regional geomorphological and hydrological properties relevant in this study at a much higher level of detail. Figure 1 shows the regions of this dataset together with the three RGI zones and the respective glaciers.

Fig. 1. Overview of the study area in HMA with glaciers in blue and orographic subregions in orange. Capital letters after some names stand for the compass directions while C stands for ‘central’. Glacier data are taken from the Randolph Glacier Inventory, version 6.0, orographic regions from Loibl (2020). The country boundaries are part of the TM World Borders dataset (Sandvik, 2008).

In the Tien Shan, the Karakoram and the western Himalaya, mid-latitude cyclones of the Westerlies provide 60–70% of the annual precipitation, while areas in the eastern and central Himalaya receive ~80% from the East Asian and the Indian summer monsoons (Bolch and others, 2012; Maussion and others, 2014). Brun and others (2017) calculate a total annual glacier mass change for HMA of −16.3 ± 3.5 Gt for the period between 2000 and 2016. Shean and others (2020) estimate the annual glacier mass change to be even higher, at ~−19.0 ± 2.5 Gt for the period between 2000 and 2018. Continued glacier wastage of this order of magnitude would substantially accelerate the formation of new glacial lakes during coming decades. HMA is expected to have lost between 36 and 87% of its glacier mass by 2100, depending on the chosen representative concentration pathway scenario, climate forcing and model (Kraaijenbrink and others, 2017; Hock and others, 2019a; Marzeion and others, 2020). Estimates of the number of current glacial lakes range from ~15 000 covering an area of 1600 km² (Chen and others, in press) to ~25 000 lakes with an area of 1800 km² (Wang and others, 2020). For all estimates, it is undebated that glacial lake area has expanded significantly over the last few decades, especially at proglacial lakes, and that new lakes tend to develop in higher elevations (Chen and others, in press; Wang and others, 2020). According to King and others (2018, 2019), this will further accelerate glacier mass loss due to a positive feedback between lowered effective pressure, rising ice velocity, longitudinal strain and increased thinning and ablation at the glacier front.

3. Data and methods

3.1 Input data

To detect subglacial overdeepenings and assess potential future lake area and volume, this study relies on three datasets: DEM data, glacier ice thickness data and glacier outlines. Regarding the DEM, previous studies of glacial lake development in certain parts of HMA mainly relied on SRTM (Linsbauer and others, 2016) or ASTER data (Allen and others, 2016). However, these DEMs are subject to considerable artifacts in high mountain regions (Bolch and Loibl, 2018). In this study, we therefore employ the ALOS World 3D-30m (AW3D30) DEM (version 3.1), which is based on the AW3D DEM with a resolution of 5 m by the PRISM sensor aboard the Advanced Land Observing Satellite (ALOS) operated by the Japan Aerospace Exploration Agency (JAXA) (Takaku and others, 2018). Despite void-filling anomalies still being present in a few mountainous areas, the AW3D30 offers the highest accuracy among publicly available global DEMs (Mo and others, 2018). Especially for studies in HMA, it is currently regarded the most suitable dataset surpassing both ASTER and SRTM elevation datasets regarding the mean absolute error (MAE) and the root mean square error (Liu and others, 2019). Naturally, both error measures show increasing trends with an increase in slope. However, even in very rugged areas with slopes >50°, the MAE for the AW3D30 DEM (>10 m) is relatively low when compared to the SRTM data (>30 m) or ASTER data (>26 m). Also, the AW3D30 DEM outperforms other DEMs particularly well in areas with elevations over 1500 m (Liu and others, 2019) which is especially required in this study.

For the glacier ice thickness, we chose to rely on a global model ensemble instead of one individual model as those tend to suffer from substantial uncertainties (Farinotti and others, 2017). We use data published by Farinotti and others (2019a) who employed a combination of up to five models to increase the accuracy and robustness of the results. Their model ensemble uses RGI v6 data for glacier outlines, SRTM v4 data for the glacier surface topography, glacier mass turnover estimates and principles of ice flow dynamics to produce a glacier ice thickness dataset with a spatial resolution of 100 m. Cross-validation is employed to assess model performance while inverse variance and bias weighting are used for the final composite result (Farinotti and others, 2019a). The difference between single model outputs and measured ice thickness (−17 ± 36%) reduces to +10 ± 24% when considering the average composite solution (Farinotti and others, 2017). This value is only surpassed by manually selecting the most suitable solution for every glacier which is not feasible in the scope of our study.

Glacier outlines for HMA are available through the RGI v6. This globally nearly complete inventory of glaciers provides information on, e.g. a glacier's area, type, elevation, slope, length and aspect. In this dataset, ice caps with multiple outlet glaciers are divided up into separate units (RGI Consortium, 2017).

3.2 Detection of overdeepenings

To locate potential future lakes, we create a DEM ‘without glaciers’, facilitating an assessment of the subglacial bedrock morphology (Paul and Linsbauer, 2012). We do so by subtracting the ice thickness data (Fig. 2a) from the AW3D30 elevation data (Fig. 2b and c). We refrain from coregistering the AW3D30 and the SRTM v4 DEM due to their very similar horizontal offset (Hu and others, 2017), the higher sampling rate in the AW3D30 data (~30 vs 90 m in the SRTM) and the potential to introduce errors during the transformation. Additionally, any horizontal uncertainties can be considered negligible compared to the much higher uncertainties of the ice thickness dataset. To reduce computing time, we decided to exclude all glaciers <1 km² which, according to preliminary tests as well as previous studies (Farinotti and others, 2019a), does not significantly affect the results. Since accumulation areas shared by multiple glaciers are separated along flow divides in the RGI, calculations near the common borders of glaciers can lead to artifacts. To avoid those and other inconsistencies in the data, we group together all adjacent glaciers within a buffer of 100 m. Finally, the bedrock topography for the whole glacier complex is calculated. Subsequently, overdeepenings in the subglacial bedrock are filled using hydrological GIS tools (ESRI, 2019; Hijmans, 2019). The difference between the DEM without glaciers and the filled DEM yields a bathymetry raster which is then used to assess the area and volume of the overdeepenings (Fig. 2d). We exclude small overdeepenings that most likely would fill with sediments rather than with water or not form at all (red areas in Fig. 2d). In choosing this threshold, we are guided by previous studies. Linsbauer and others (2016) and Colonia and others (2017) chose a threshold of 10⁴ m² of surface area for the initial calculation. For HMA, Linsbauer and others (2016) considered only overdeepenings with surface areas larger than 10⁶ m² for most of their evaluations. We also employ 10⁴ m² as a threshold for the computation of overdeepenings (yellow areas in Fig. 2d). However, to investigate future lake development in more detail, we set our analysis threshold one order of magnitude smaller than previous studies in HMA – at 10⁵ m² of surface area (blue/violet areas in Fig. 2d). Additionally, we use a minimum depth threshold to rule out very shallow sinks. Following our assessments, lakes with a depth >10 m can potentially be the source of a significant GLOF event. Therefore, we employ 5 m as a second threshold leaving only those sinks that realistically can form and fill up with water.

Fig. 2. Illustration of the detection of overdeepenings. (a) Ice thickness data for Biafo Glacier, Pakistan, as provided by Farinotti and others (2019a). (b) DEM of the northern part of the glacier and the surrounding mountains. (c) DEM ‘without ice’ depicting the calculated subglacial bedrock. (d) Bathymetry of bedrock overdeepenings >10⁵ m² (blue) as well as overdeepenings >10⁴ m² (yellow) and overdeepenings <10⁴ m² (red).

Uncertainties regarding overdeepening area and volume estimates are determined by two factors: (1) the accuracy of the ice thickness data and, (2) the vertical precision of the AW3D30 DEM. With uncertainties of ±25.99% (RGI zone 13), ±25.78% (RGI zone 14) and ±26.14% (RGI zone 15), the accuracy of the five-model-ensemble ice thickness estimates is nearly consistent for the entire HMA and significantly higher than with previous single model runs (Farinotti and others, 2019a). In this study, we use the mean uncertainty of ±25.97% for the modeled ice thickness in all three RGI zones in HMA. Following NELAK (2013), glacial lakes deeper than 30 m are described as ‘deep lakes’. Few glacial lake inventories include depth estimates but Muñoz and others (2020) describe a maximum glacial lake depth of ~130 m for the Cordillera Blanca. However, in our study we found that some computed overdeepenings exceed this depth by far. As this is partly the effect of erroneous void fillings in the AW3D30 data, we disregard overdeepenings where depth values are strongly affected by such void fillings from our calculation using a primary maximum depth threshold of 1000 m. We use a secondary, still fairly high threshold of 300 m to mark overdeepenings for manual validation, as the overdeepenings mapped in our study are about one order of magnitude larger than the lakes described by Muñoz and others (2020). Each marked overdeepening was checked in GIS for visible erroneous void fillings in the AW3D30 DEM as well as for incongruencies between Landsat surface reflectance images and the shape of the overdeepening in the DEM. Overdeepenings with questionable characteristics were excluded from further analyses.

As we avoid grave outliers in our data with the previous step, we rely on the MAE for describing the DEM accuracy as it is less sensitive to such outliers and, therefore, better fits our adjusted data. For the AW3D30, Liu and others (2019) measured the MAE as 2.79 m for HMA. The resulting uncertainty depends on the respective glacier ice thickness and ranges from ±0.59 to ±69.75%. For all overdeepenings larger than 10⁵ m², the AW3D30 uncertainty averages ±2.87%. Combined with the uncertainty from the ice thickness model (±25.97%) this results in a cumulative uncertainty of ±28.8%. We employ this cumulative value for both lake area and volume, as Cook and Quincey (2015) found a strong correlation (R ² of 0.91) between both parameters that we were able to confirm for our findings (R ² of 0.86, significance level: <0.001).

3.3 Slope hazard assessment

To assess the hazard of rock or ice avalanches hitting the lake we identify all relevant slopes of the DEM without glaciers in a 2 km buffer. Relevant slopes are initially defined as slopes belonging to the catchment area of the lake. However, the flow behaviors of avalanches and surface runoff are substantially different from one another, since avalanche tracks are, although guided by the underlying topography, far more linear (Norem and others, 1989; Fischer and others, 2012a). Ice and rock avalanches in lateral valleys may well originate from a location within the lake's catchment area but are unlikely to directly hit the lake and initiate a GLOF. If avalanches reach the floor of the lateral valley, their mobility will be heavily decreased by a run up on the opposing valley side (Nicoletti and Sorriso-Valvo, 1991). Therefore, if lateral valleys exceed a certain drainage area, we exclude their catchment area from the hazard assessment. This is determined by the number of pixels that are drained through a valley's inlet into the lake. We set the threshold at 200 pixels (~0.18 km²) which, depending on the geomorphology of the catchment, generally excludes lateral valleys with a length of more than 1.5–2 km. With this approach, we provide a catchment area containing exclusively slopes from where mass movements would directly impact the lake.

The relationship between a region's terrain and the morphological characteristics of its future glacial lakes can be investigated using the terrain ruggedness index (TRI) by Riley and others (1999). This index corresponds to the mean elevation difference between adjacent pixels in a DEM and is used to evaluate the average elevation change in a region (Riley and others, 1999). We calculate the TRI for each raster cell and use the mean TRI of all raster cells in an orographic region as a proxy for this area's ruggedness (Evans and others, 2014). In this study, the TRI is used for regional comparisons of the spatial distribution of future glacial lakes in order to find correlations between a region's terrain and the development and morphological characteristics of future glacial lakes. Allen and others (2016) introduced the term ‘lake impact predisposition area’ (LIPA) for the combined area around a glacial lake from where an impact of ice or rock avalanches is possible. We expand this concept by dividing DEM cells in the respective catchment area into single slope units (SU) based on slope angles. Although this may split morphologically connected slopes in the landscape into smaller SU, this approach results in a more detailed approximation of the hazard of mass movements. For the SU delineation, we define four groups (Table 1), following the results of previous studies (Allen and others, 2011; Hermanns and others, 2012; Fischer and others, 2012b).

Table 1. Grouping criteria for the angle-dependent slope division

Areas with a slope ≤20° are considered safe and excluded from further calculations.

The first group contains areas from where mass movements typically are unlikely (slope angle <20°). Although this group is excluded from further analyses, the remaining three groups comprise potentially dangerous areas that are analyzed in the slope hazard classification. Additionally, we consider further topographical data for each SU: area, distance from and height difference to the lake. Thus, we can assign more detailed hazard scores to each SU, yielding a more comprehensive estimation of the overall impact hazard for each lake. The resulting values should, however, not be interpreted as absolute – in the sense of a final hazard assessment – but instead be treated as a means for recognizing similarities and differences between the different regions of HMA. First and foremost, they are meant as a comparative measure for organizing and prioritizing further research and field studies.

The hazard score (H _SU) for each SU is approximated as a weighted mean of normalized relevant topographical parameters, as shown in Eqn (1):

(1)$$H_{S{\rm U}} = \displaystyle{{2A + 2S + D + 0.5E} \over {5.5}}$$

where A is the SU area, S is the mean SU slope, D is the shortest distance from any part of the SU border to any part of the lake border and E is the elevation difference between mean SU elevation and the lake surface. Higher values for A indicate a higher hazard, as more potentially unstable material is to be expected in a larger area. Higher values for S also contribute to a higher hazard classification, as material is more likely to detach from steeper slopes. As we invert the value D, high values for this parameter indicate areas in close proximity to the lake which raises the hazard. Higher values for E indicate a higher hazard as detached material can potentially gather more speed and hit the lake as well as other unstable material in its path with increased force. In accordance with the recommendations of the ‘Standing Group on Glacier and Permafrost Hazards’ of the IACS/IPA (GAPHAZ, 2017), A and S are considered to be the most important parameters and thus have the highest weight factors. D and E were introduced to improve the accuracy of the classification by adding factors that may exert further control on the hazard regimes in the heterogeneous topographic settings of HMA. All four parameters are rescaled. We normalize the slope S between 20 and 90° and use the 2 km buffer as a maximum and 0 as a minimum for the distance D. The elevation E is rescaled between 0 and an arbitrary maximum of 1000 m, while for the area A, we use the ratio between SU and lake area. Figure 3 shows the topographic parameters on the basis of three selected SU in a schematic depiction. The colors attributed to the slopes reflect the potential for rockfalls or icefalls impacting the lake (red = very high potential; orange = moderate potential; yellow = very low potential). Due to the weighting of the parameters, the large, steep area (left in Fig. 3) very close to the lake is attributed a very high hazard despite its smaller elevation difference. The smaller area (in the middle of Fig. 3) receives a moderate hazard score despite its high values for S and E value because of its small values for A and D. The gently sloped area (Fig. 3, right) is attributed a very low hazard score as its slope angle is <20° and would be excluded from further calculations.

Fig. 3. Schematic depiction of the topographic parameters used for the hazard classification of selected SU surrounding a future lake in a 2 km buffer. E = elevation difference between mean SU elevation and lake surface, A = SU area, D = distance to the lake, S = mean SU slope. The colors indicate different hazard scores attributed to the slopes (red = very high potential; orange = moderate potential; yellow = very low potential).

The arithmetic mean of all relevant slopes in a lake's watershed is considered to be a proxy for the hazard of rock or ice avalanches impacting the lake. Resulting values range from 0 to 1, with higher values implying a higher probability of detached material to impact the lake and produce a considerable displacement wave that could lead to a GLOF event. However, as an avalanche or a rockfall from one very hazardous SU would be enough to trigger an event even for a lake with a very low mean hazard score, we additionally calculate the maximum H _SU for each lake. Furthermore, we consider the lake volume to account for the potential magnitude of the GLOF (Eqn (2)). We employ a linear normalization for this parameter between 0 and 1 with 10⁻⁴ km³ as the lower boundary and 10⁻¹ km³ as the upper. The lake hazard level is defined as the product of either the mean slope score (LHL _mean) or the maximum slope score (LHL _max) and the normalized lake volume. These variables are considered a first-order evaluation of the hazard of GLOF events triggered by mass impacts into future glacial lakes:

(2)$$LHL_{{\rm mean}} = \displaystyle{{\mathop \sum \nolimits_{i = 0}^n H_{{\rm S}{\rm U}_i}} \over n} \times V\semicolon \;\;\quad LHL_{{\rm max}} = H_{{\rm S}{\rm U}_{{\rm max}}} \times V$$

where V is the normalized lake volume.

4. Results

4.1 Bedrock overdeepenings

In the computed subglacial bedrock, we find 25 285 overdeepenings with a surface area >10⁴ m² and a total volume of 99.1 ± 28.6 km³ covering 2683 ± 773.8 km². Applying the 10⁵ m² threshold leaves 2700 overdeepenings in the bedrock of 669 glaciers, with 2257 being larger than 10⁶ m². Overdeepenings larger than 10⁵ m² cover an area of 1623 ± 468 km² which is 1.18–2.14% of today's glacierized area. Their volume of 72.6 ± 20.9 km³ equals 0.74–1.33% of the current glacier volume estimated by Farinotti and others (2019a). A total of 1535 future lakes (57%) can be classified as ‘deep lakes’ with a depth of more than 30 m (NELAK, 2013). Altogether, 16 overdeepenings >10⁵ m² (0.6%) were marked for manual validation because of a depth >300 m. With depth values ranging from 1028 to 2223 m, five overdeepenings are automatically excluded for crossing the threshold (depth >1000 m). Supplementary Figure S1 shows some of those overdeepenings that are clearly affected by artifacts in the DEM. Another eight overdeepenings, with depths ranging from 308 to 719 m, were excluded because of clearly visible void fill errors in the AW3D30 DEM data. In case of doubt, comparisons with Landsat surface reflectance images provided additional indications, as some erroneous overdeepenings were placed across mountain ridges. Only three overdeepenings marked for manual validation (300, 308 and 472 m depth) were kept after this step to verify values of extreme overdeepening. Figure 4 shows the spatial distribution of computed overdeepenings >10⁵ m² over HMA. Clearly visible are the regional differences between highly glacierized rugged mountain ranges (e.g. the Central Karakoram, Central Tien Shan and Pamir) and regions with less extreme relief (e.g. the Gangdise, Qilian or Hengduan). Supplementary Table S1 provides a TRI-based terrain comparison for selected regions. Statistics for the most important morphological features of the computed overdeepenings are summarized in Table 2.

Fig. 4. Identified subglacial bedrock overdeepenings (>10⁵ m²) per region and cumulative overdeepening volume.

Table 2. Summary statistics for subglacial bedrock overdeepenings (surface area >10⁵ m²) in HMA

By far the most (653) and biggest (29.5 ± 8.5 km³) overdeepenings are located in the Central Karakoram, accounting for about one-quarter of the total number of overdeepenings. However, owing to the abundance of glaciated surfaces, overdeepenings cover only 5.5% of the currently glacierized area in the Central Karakoram. In comparison, the 51 overdeepenings in eastern Kunlun cover 16% of the local glacierized area. Nevertheless, overdeepenings in Central Karakoram account for 18.7% of the area and 29.6% of the total volume of all overdeepenings in HMA, owing to a number of very large lakes.

Besides the dominant Central Karakoram, there are other regional differences to be found regarding the morphometric characteristics of overdeepenings. Table 3 lists selected regions (see Supplementary Table S2 for all regions) that stand out regarding their relationships between overdeepening area/volume and their current glacierized area/total ice volume. Particularly for the Nun Kun Range, Banderpunch Gangotri and the Dhaulagiri region, the volume ratio is smaller than the area ratio. In general, this indicates shallower overdeepenings for the western part of the Himalaya. The same applies for the northern Hindu Kush and the western Kunlun. In contrast, for the Central Karakoram the volume ratio is larger than the area ratio, suggesting relatively deep overdeepenings. For the Central Tien Shan, eastern Kunlun and Bhutan the same is true in attenuated form. Furthermore, the eastern Kunlun and Bhutan stand out due to a relatively low number of overdeepenings in combination with a comparably high volume and area ratio. In both regions, this can be explained by a relatively small glacierized area in combination with few very large overdeepenings beneath glacier tongues.

Table 3. Statistics for glaciers and overdeepenings (surface area >10⁵ m²) in selected regions of HMA

4.2 Hazard classifications

Over the entire HMA, we classified 8516 km² as hazardous slopes around future glacial lakes. Approximately half of those slopes (48.4%) are mapped within exposed glacier beds. The HMA-wide mean slope hazard forms a mesokurtic distribution with a mean of 0.49, a std dev. of 0.1, and a skewness of −0.4 (Fig. 5). As such, the data relatively closely resembles a normal distribution indicating equilibrated results produced by the algorithm. Expectably, the leptokurtic distribution of the maximum hazard is skewed even further to the left with a mean of 0.7 and a std dev. of 0.17. The secondary peak at zero indicates 46 (1.7%) of the total of 2700 future lakes that were assigned an impact hazard of zero. The average angles of slopes adjacent to those lakes were <20° and, therefore, classified as not relevant. Occurring predominantly in the Kunlun Shan and Inner Tibet, these lakes are characterized by small average area and volume. Its specimens are often located in large cirques or in the central parts of large, flat glacial troughs. With a cumulative area of 17.5 km² and a cumulative volume of 0.3 km³, they only make up 1.1% of the area and 0.4% of the volume of all overdeepenings in HMA and were thus excluded from further analyses.

Fig. 5. Distribution of the mean and maximum lake hazard for all overdeepenings in HMA.

Figure 6 depicts two overdeepenings under Siachen Glacier, Karakoram and the hazard classification of their surrounding slopes to demonstrate the results of our algorithm. Notice that inlets of surface runoff streams from larger tributary valleys and the upper glacier bed (small black circles in Fig. 6) are used to exclude the associated catchments. Slopes belonging to the remaining catchment area of each overdeepening (black outline) are then classified according to the hazard of icefalls or rockfalls impacting the lake. With all areas with a slope <20° excluded, the yellow-to-red color ramp indicates the hazard from moderate (yellow) to very high (red). The division of the mountainsides into SU allows an in-depth analysis of the area surrounding a potential lake as the overall slope morphology can be reproduced following the different hazard scores. Larger, steeper slopes in close distance to the lake with high values for A, S and D receive higher scores. More gently sloping areas with low S values or areas in greater distance (low values for D) receive a lower hazard score.

Fig. 6. Bathymetry of two overdeepenings under Siachen Glacier, Karakoram, and the results of the classification of all relevant SU inside the catchment area (dark outline). SU can be delineated by their gray outline and the different colors indicating their hazard score. Small black circles indicate the inlets of large tributary valleys that were excluded from the calculation.

Figure 7 depicts a very large overdeepening, located beneath the glacier tongues of central and south Rimo Glacier, Karakoram, as another, more extreme example (see also Section 5.1). Noticeable are the tributary valleys of considerable size whose catchment areas are excluded from the hazard assessment. Those valleys are indicated by their inlets (black points in Fig. 7). Smaller, less incised lateral valleys, however, are included in the assessment as their slopes more directly point toward the lake and detached material is more likely to directly impact the lake. In Figure 7, these areas are indicated by black circles. Dashed black rectangles indicate the location of large, steep side moraines discernable by their higher score. They partly block the catchments behind from draining into the lake and would probably also be an obstacle in the way of a potential rockfalls or avalanches.

Fig. 7. Bathymetry (blue colored) for the largest potential lake (11.6 km³) at Rimo Glacier, Karakoram, with hazard scores (yellow-to-red color ramp) for adjacent slopes. Smaller lateral valleys included in the slope assessment are indicated by black circles, excluded larger valleys by their inlets (black points). Dashed black rectangles show steep side moraines distinguishable by their higher hazard classification. Surface runoff channels are indicated by blue lines.

The analysis of the relevant slopes for all lakes on a regional scale, i.e. the cumulative LIPA (Fig. 8), shows a pattern similar to the distribution of the overdeepenings with a contrast between high mountain glaciers and glaciers in less rugged regions. In consequence, most areas with a large number of overdeepenings, such as the Tien Shan, the northern Pamir and the Nun Kun Range, also have comparably large LIPA. Conversely, only very few slopes are classified as LIPA in the eastern Kunlun, despite the region's relatively high number of overdeepenings. Again, the central Karakoram stands out with nearly 9% of its total area being a potential origin for mass movements that could impact future lakes.

Fig. 8. Regional differences between absolute and relative LIPA for subglacial overdeepenings in HMA. The relative values reflect the percentage of LIPA for the whole subregion area.

Spatial patterns of impact hazards are assessed based on a grid of 1° cells (see also Fig. 10). We calculated the mean gridcell TRI as well as the mean gridcell impact hazard. The latter serves as a proxy for the average impact hazard of all lakes within, i.e. their predisposition for mass-movement impacts. Through the analysis of all 112 gridcells, a weak yet statistically significant correlation can be found between the TRI of a gridcell and its mean LIPA (R ² = 0.37; Fig. 9a). Hardly any correlations exist between the TRI and a gridcell's mean impact hazard (Fig. 9b) or its number of overdeepenings (Fig. 9c). However, the relationship between both variables and the TRI is significant. In contrast, overdeepening area, volume and depth show no significant correlation to the TRI, neither does the lake hazard level (LHL _mean and LHL _max), which combines a lake's impact hazard and its volume (see Supplementary Fig. S2 for further correlations).

Fig. 9. Correlation between a gridcell's ruggedness and (a) its mean LIPA, (b) its mean hazard and (c) the number of its future lakes.

Clearly visible is the aforementioned contrast between the high-mountain glaciers of the Karakoram, Tien Shan Pamir and parts of the Himalaya on the one hand and the gently sloped glaciers and ice caps of Inner Tibet, the eastern Kunlun and the Qilian Shan on the other hand (Fig. 10). Lakes with higher impact hazard scores are mostly found in the high mountain ranges. Conversely, low impact hazard scores prevail on the less rugged topography of the Tibetan Plateau.

Fig. 10. Gridcell-based spatial distribution of the mean future lake impact hazard in HMA and the associated cumulative potential future lake volume. Gridcells with fewer than five overdeepenings are symbolized by dashes.

Parameter values of the impact hazard score for both high relief and low relief areas show noticeable similarities within their respective group, yet are very different from each other (Fig. 11). The same holds true for the region-wide mean impact hazard. In regions of less rugged topography, mountain slopes adjacent to overdeepenings tend to be relatively small and exhibit lower overall slope. The overall very high values for D (distance) indicate slopes in close proximity to the overdeepening whereas the very low values for E (elevation difference) reflect only minor differences between mean slope elevation and the calculated lake surface altitude. Therefore, lakes in these regions are expected to have a comparably low predisposition to mass-movement impacts. Moreover, the low maximum impact hazard indicates that most lakes are indeed surrounded exclusively by slopes with only a moderate probability of producing rockfalls or icefalls. In comparison, overdeepenings in higher altitudes are surrounded by slopes of varying areas A with higher overall slope S. In these settings, slopes are classified as hazardous up to greater distances from the overdeepening (high D) and display higher elevation differences to the lake surface altitude (high E). This results in a higher mean hazard for these regions. Also, the very high maximum hazard indicates that most lakes in higher altitudes are located close to at least one slope with a very high potential of generating mass movements in their direction.

Fig. 11. Comparison of the four slope hazard classification parameters (area, slope, distance and elevation) in high-relief (upper half) and less extreme relief (lower half) regions of HMA as well as region-wide mean and maximum slope hazard classifications.

Considering the lake hazard level (both LHL _mean and LHL _max) as the combination of slope hazard score and potential lake volume, smaller lakes in high altitudes with low potential volume are contrasted by large overdeepenings under gently sloped glacier tongues. Although the former are assigned a low impact hazard due to their comparably low volume, the large proglacial lakes receive the highest impact hazard scores, always depending also on the properties of their surrounding slopes. The largest of these lakes are located in the Central Karakoram (20 lakes with a combined volume of 22.58 ± 6.52 km³), the Mahalangur region of the Himalaya (12 lakes with combined volume of 2.93 ± 0.85 km³) and the Northern Pamir (eight lakes with a total volume of 2.57 ± 0.74 km³).

5. Discussion

5.1 Data quality

This study benefits from recently published data that allows for a more robust computing of overdeepenings as well as for an impact hazard classification at a larger scale than previously possible. The quality of the glacier outlines in RGI v6 substantially improved since v4, which was used by Linsbauer and others (2016). Still, there are some glaciers missing in the RGI or other inconsistencies in the data, however, they have only minor influences on the accuracy of the regionally computed glacial lake volumes. Frey and others (2010) proposed a four-level approach to assess the validity of computed subglacial overdeepenings. According to their findings, our results are categorized as a level-3 result, outperformed only by on-site measurements using drills and/or ground-penetrating radar. Despite significant uncertainties regarding the absolute values, level-3 data are expected to be a realistic estimate of lake area and volume as well as the lake bottom topography (Frey and others, 2010).

Figure 12 shows a comparison between SRTM and AW3D30 data in the Karakoram area (see Fig. 14 for an overview of the larger region). Black areas in Figure 12a indicate numerous voids in the SRTM data, which would lead to incomplete hazard classifications of nearby overdeepenings. Only small artifacts remain in the AW3D30 DEM as indicated by the black circles in Figure 12b. In general, this DEM tends to overestimate the terrain elevation over void-filled areas (Liu and others, 2019), as can be seen by the maximum height of 9035 m. Even so, local underestimations can lead to unrealistically high depths of overdeepenings that have to be addressed (see also Supplementary Fig. S1). The difference raster (Fig. 12c) shows a good agreement of both DEMs (±10 m) for some of the wide flat glacier areas with smaller overdeepenings. However, for the large future proglacial lakes in the northeast, there are substantial differences between both DEM. In general, both datasets have a very high correlation, although it depends on the elevation (Fig. 12d). Figures 12c and d indicate that, especially in higher and more rugged areas, the SRTM DEM differs substantially from the AW3D30. The AW3D30 dataset is expected to be a more realistic portrayal of the actual terrain (Liu and others, 2019; Guan and others, 2020).

Fig. 12. Data quality comparison between the SRTM and the AW3D30 DEM in the Karakoram area. (a) Numerous voids in the SRTM DEM are indicated by black areas. (b) Small artifacts in void-filled regions in the AW3D30 DEM are indicated by black circles. (c) Difference raster between both DEM. (d) Raster value scatterplot for both DEM showing growing differences with increasing elevation. SRTM 1 arc-second Global DEM data provided by EROS (2017).

In general, the AW3D30 DEM improves the quality of our computed overdeepenings due to its so far unmatched accuracy in HMA, especially in higher areas. As the lowest future lakes investigated in this study are located ~2200 m, AW3D30 is best suited for our purpose. Due to the higher quality of the DEM, we also eliminate some of the erroneous void filling that posed challenges in previous analyses of future lakes in HMA conducted with SRTM data (Linsbauer and others, 2016) or ASTER data (Allen and others, 2016). Thus, and following Frey and others (2010), we assume the results of our overdeepening mapping to be based on realistic estimations of glacier bed morphology. The vertical accuracy of the AW3D30 is only available HMA-wide as an MAE of 2.79 m (Liu and others, 2019). Thus, the relative uncertainty for each single glacier solely depends on its ice thickness. With glaciers large enough for overdeepenings >10⁵ m² of surface area to form in their bedrock, the mean additional uncertainty of 2.9% is only a minor factor. For glaciers with overdeepenings >10⁴ m², the mean uncertainty rises by 6.5% and for all other, even smaller glaciers in HMA it rises by 10.9%.

Uncertainties in the calculated DEM of subglacial topography mainly originate from the employed ice thickness data. Due to the global scope of their study, Farinotti and others (2019a) do not provide uncertainties for each glacier individually but composite values per RGI region. Therefore, it remains unclear how the given uncertainty values for each RGI region relate to, e.g. glacier type, size or location. Whether the cumulative uncertainty for HMA (28.8%) is an over- or underestimation for larger glaciers as investigated in this study cannot be determined. However, as the accuracy of composite ice thickness data is substantially improved compared to the results of single models (Farinotti and others, 2017), we expect the accuracy of our overdeepening and slope mapping to be comparably high.

Due to the fact that the ice thickness estimates of Farinotti and others (2019a) are based on SRTM v4 data, void-filling errors still factor into the uncertainty of our results, the exclusion of especially affected overdeepenings notwithstanding. Moreover, half of the mapped slopes (48.4%) are located inside the exposed glacier bed. Despite our realistic estimate of lake bottom topography, the morphology of those slopes is subject to higher uncertainties as it is not directly based on the AW3D30 DEM but on SRTM v4 data and the ice thickness modeled by Farinotti and others (2019a). In most cases, such slopes within the exposed glacier bed are adjacent to small lakes in heavily glaciated areas in the high mountains. Conversely, large proglacial lakes in most cases fill up the whole glacier bed and are primarily surrounded by slopes represented by the AW3D30 data. As a result, SRTM uncertainties mainly increase the uncertainty regarding the hazard classifications of smaller lakes in the high mountains and do hardly factor into the uncertainty for larger proglacial lakes. In general, our results confirm morphological characteristics and the surface area for larger overdeepenings in the Himalaya–Karakoram range (RGI zone 14) presented by Linsbauer and others (2016). However, we are not able to reproduce the number of overdeepenings and the total volume calculated in this previous study. Linsbauer and others (2016) mapped ~16 000 overdeepenings with a combined volume of ~120 km³ whereas our study found 8744 (55%) overdeepenings with a total volume of 44.8 km³ (37%) in the Himalaya–Karakoram range. As fully reproducing the results of Linsbauer and others (2016) would exceed the scope of this study, we assume the differences are predominantly caused by two factors: (1) the lower resolution and abundant void-filling artifacts in the SRTM v4 data used by Linsbauer and others (2016) and (2) the different ice thickness estimates of both studies. The ice thickness dataset used in this study was calculated by Farinotti and others (2019a) based on an ensemble of five models, including the GlabTop2 model used for the ice thickness calculations by Linsbauer and others (2016). According to Farinotti and others (2017, 2019a), ice thickness model ensembles significantly improve robustness and accuracy of the findings over individual models that are more prone to suffer from substantial uncertainties. Still, larger errors for individual glaciers are possible despite overall robust and accurate regional thickness estimates. Spatially, higher overdeepening volumes in the dataset by Linsbauer and others (2016) prevail in the western Himalaya and the Karakoram and thus coincide with the region of most abundant void-filling errors in the SRTM dataset. The most pronounced example is the large (50 km²) overdeepening under the central Rimo Glacier (Fig. 7): with a volume of 11 km³ it contains 9.2% of the total overdeepening volume calculated by Linsbauer and others (2016). We estimate this overdeepening's area to be 86.6 km², assuming a confluence of the two overdeepenings under central and south Rimo Glacier. With only 57.7% of the area according to our study, Linsbauer and others (2016) still estimated nearly the same overdeepening volume (11 km³) as we found in our study (11.6 km³). Presumably, this is due to a calculated depth of 588 m in their study compared to 329 m in our study. Due to the smaller surface area, this indicates a much larger overdeepening and excess overdeepening volume by ~4.3 km³ by Linsbauer and others (2016).

5.2 Overdeepenings

The 25 285 overdeepenings with an individual area of >10⁴ m² cover 2683 ± 773.8 km² and indicate a potential increase in the number of glacial lakes in HMA of ~100 and 170% according to current lake estimates by Wang and others (2020) and in press, respectively. This corresponds to an expected increase in lake area of 150% to 170%. The mean glacial lake area, currently estimated to be 0.06–0.1 km², is projected to increase up to 0.6 km². Besides the expected growth of current glacial lakes, this corroborates a trend found by previous studies, indicating future glacial lakes are going to be larger than already existing ones (Linsbauer and others, 2016; Wang and others, 2020). Exact timescales for the development of future lakes are highly uncertain. However, it is possible to provide an approximate time span for lakes located in the glacier ablation area, which in this context can be considered as the area below the mean glacier elevation, following Braithwaite and Raper (2009). With constant rates of glacier retreat, the Karakoram is expected to lose 50–75% (relative to 1985) of its glacier mass by 2035 (Cogley, 2011) which would uncover the glacier bed of most glaciers in this region. Following this estimate, many large proglacial lakes will develop over the next two decades in this region of HMA. This is especially important as these lakes account for a significant fraction of both the total lake area and the total lake volume found by this study and could pose a significant threat to local infrastructure (see also Section 5.4). When calculating the overdeepening volume, we regard every overdeepening to be fully filled which, at times, may overestimate the actual volume. Similar to other studies (Allen and others, 2016; Colonia and others, 2017), we regard the results as a reliable average that accounts for lakes that are not fully filled as well as lakes surrounded by moraine dams that possibly fill up to a higher level. In HMA, the latter case is likely especially close to existing glacier tongues due to high debris availability (Kirkbride, 2011; Kääb and others, 2012).

Currently, the central and eastern Himalaya, particularly Bhutan and Nepal (Fig. 13b) are considered the GLOF hotspots of HMA with glacial lakes being about two orders of magnitude larger than in the less affected western parts, e.g. the Karakoram, Hindu Kush or Pamir (Gardelle and others, 2011; Veh and others, 2019) (Fig. 13a). In the Karakoram, smaller supraglacial lakes prevail (>90%) whereas in Bhutan and Nepal, the majority (~75 and 85%, respectively) are fast-growing proglacial lakes (Gardelle and others, 2011). Our assessment suggests that the hotspot of large glacial lakes (>10⁵ m²) will shift toward the northwest in the future. In part, this is just due to the fact that there are more and larger glaciers yet to melt in the Karakoram region as many of them are stagnant or retreating at a much slower pace than in the Himalaya. Therefore, this region offers a higher potential for the development of future glacial lakes. Another factor could be the comparably high percentage of surging glaciers in this region. However, data on surge-type glaciers in the RGI v6 are still fragmentary. Of the 669 glaciers with overdeepenings, only 84 (13%) are given a data entry regarding their surge classification. Unfortunately, surging was only observed at 36 of those glaciers and could be possible at another 28 glaciers. With data this sparse, a mathematical correlation between overdeepenings and the occurrence of surge-type glaciers cannot be carried out while visual interpretation remains highly speculative. Thus, this issue needs to be investigated further.

Fig. 13. Comparison of (a) the potential future and (b) the current glacial lake hotspots in HMA. The lake distribution is expected to shift from the southeastern Himalaya (b) toward the Karakoram and Pamir (a). The mapped current glacial lakes are a composite of the Hi-MAG database by Chen and others (in press) and the lake inventory of  Wang and others (2020).

With HMA expected to lose up to 87% of its glacier ice mass in the next 70 years (Shean and others, 2020) and the Karakoram transitioning from a positive to a negative mass balance (Brun and others, 2017), our projections indicate that in a scenario with considerably less ice mass the largest glacial lakes will be found in the Karakoram and Pamir (Fig. 13a). Moreover, future glacial lakes in this region will be far more susceptible to mass-movement impacts (see also Fig. 10). The described northwestern shift not only in the number of glacial lakes but also in lake size and impact predisposition could, in turn, further enhance glacier melt and the potential GLOF hazards in that region.

5.3 Impact hazards

For the assessment of the GLOF potential for already existing lakes, parameters such as permafrost conditions (Bolch and others, 2011), lake expansion (Dubey and Goyal, 2020), moraine dam morphology (Wang and others, 2012), vegetation cover (Veh, 2019) and the likelihood of extreme meteorological conditions (Singh and others, 2014) can be included in addition to a slope hazard classification. When investigating future lakes, however, most of these parameters are unsatisfactory due to enormous uncertainties or a general lack of data. Therefore, and since more than 50% of all GLOF events in HMA occur due to displacement waves (Allen and others, 2016), the lakes' predisposition to mass-movement impacts in combination with its volume can serve as a first impression of the potential GLOF hazard in a region. We provide large-scale and detailed information about potential LIPA as well as detailed hazard classifications of related slopes in the catchment. Over the whole of HMA, we classified 8516 km² as potential LIPA for glacial lakes with a surface area >10⁵ m². Notably, this could be considered a conservative estimate due to the exclusion of large lateral valleys. Almost every future lake (98.3%) is surrounded by potentially dangerous slopes with slope angles >20°. Very few small lakes in large cirque basins or flat glacial troughs exhibit adjacent slopes entirely <20° from where mass movements generally are unlikely (Hermanns and others, 2012). Our results show that comparably small future lakes in the higher parts of mountain ranges within HMA tend to have an especially high predisposition for mass-movement impacts. Typically, these high hazard scores are due to steep slopes >40° in close proximity to the respective lake and/or large elevation differences between the mean slope altitude and the calculated lake surface altitude. In regions with less extreme relief, such as Inner Tibet, or parts of the Gangdise and Kunlun Ranges, many of the smaller future lakes show low hazard scores. In these situations, glaciers have mostly retreated to the mountain peaks and few steep slopes remain above possible future lakes. Therefore, these regions seem to be less threatened by GLOF events triggered by mass-movement impacts in general. However, for the few larger lakes in this area a more in-depth analysis of their potential hazard is required, as a significant destructive potential can still be attributed to large lakes at lower altitudes due to their sheer size (GAPHAZ, 2017).

Our results indicate that most of the future lakes in HMA will form in the still heavily glaciated regions of the central Karakoram, Tien Shan and Pamir. In the higher altitudes of these mountain ranges, above the lower permafrost limit, mass movements are dominated by low-magnitude, high-frequency rockfalls and avalanches (Blöthe and others, 2015). With transported material of <10⁶ m³, these events may be small compared to some of the catastrophic landslides of the past. However, they still are capable of generating large displacement waves. The continued warming trend and glacial retreat in nearly every part of HMA increases the slope instability at higher altitudes due to retreating permafrost and the loss of stabilizing pressure at the glacier flanks (Schaub and others, 2013) and, thus, further increases the GLOF hazard in these regions. Currently, landslides of larger magnitudes detach with low frequency and below the mean elevation of a mountain range. They transport much more material (up to 20 km³) over a vertical drop of mostly more than 700–900 m (Blöthe and others, 2015). Therefore, they still must be considered as potential hazards for the few future glacial lakes at lower altitudes. Due to the fact that lower lakes are in general closer to more densely populated areas, the few larger overdeepenings at lower altitudes may still pose a significant threat to local communities. Although the maximum slope hazard is highest in the regions with the most rugged topographies, i.e. the Karakoram, Pamir and Tien Shan (Figs 14a, c), the consideration of lake volume changes the spatial pattern. Although both LHL _mean and LHL _max display this change, it is most visible when considering LHL _max. With this, larger proglacial lakes at lower altitudes can often be expected to pose the highest threat (Figs 14b, d). Due to their size, they are very susceptible to mass-movement impacts and, at the same time, capable of generating large outburst floods relatively close to settlements and infrastructure. Although infrastructure information is fragmentary for the Karakoram region, several roads can be determined that run close to potentially dangerous future lakes (Fig. 14).

Fig. 14. Impact hazard for future lakes in the Central Karakoram mountains based on (a) the maximum slope hazard score and (b) the maximum lake hazard level (LHL _max). Clearly visible is the difference between small lakes at high altitudes and larger proglacial lakes at the glacier tongues. The subsets (c) and (d) illustrate this phenomenon at greater detail. The extent of the subsets is indicated by black rectangles in the overview maps (a) and (b). Roads (black lines) are a composite of data provided by Meijer and others (2018) and the CIESIN/ITOS (2013) dataset.

Additionally, and this has not yet been considered in the scope of this study, some overdeepenings situated below glacier tongues may develop behind current proglacial lakes, eventually coalescing with them and forming even larger lakes. Also, ice avalanches are an alternative GLOF trigger mechanism (Schaub and others, 2016). To date, an assessment of predisposition to ice avalanching is not feasible for large regions, so that this also does not factor into this study. However, the major glacier collapses in the Aru Range in July 2016 (Tian and others, 2017; Kääb and others, 2018) demonstrate the hazard potential such processes may pose in HMA under continuously warming climate. In most glaciated regions of HMA, settlements and infrastructure are still sparse. Nevertheless, vulnerability is often high because the destruction of a single road or bridge poses a significant threat as it can leave entire upstream valleys and settlements cut off (Carrivick and Tweed, 2016). The settlement density is expected to increase by the ongoing expansion of infrastructure into higher elevations for tourism, hydropower or agriculture (GAPHAZ, 2017). With regard to the growing GLOF threat, the long-term lake hazard assessment of this study can serve as the groundwork for future research into regional or glacier specific hazards in order to provide further information for planning and mitigation measures. Himalayan communities often struggle with effectively planning, managing and funding mitigation projects (Thompson and others, 2020). Following the continuing retreat of glaciers in most parts of HMA, the larger overdeepenings under glacier tongues at lower altitudes will probably become ice free and fill up due to glacial melt in the near future. As our study indicates the development of large future lakes will focus on the northwestern part of HMA, adaptation measures taken in the southwestern Himalayan region (Bajracharya, 2009; Mahagaonkar and others, 2017) will need to be considered and adopted by communities in the Karakoram and Pamir as well.

5.4 Proglacial lakes

Currently, there are more than 700 proglacial lakes in HMA whose total area increased by 50% over the last three decades (King and others, 2019). Proglacial lakes pose a significant threat for two reasons: first, the debris cover on glaciers in HMA is higher at lower elevations leading to higher moraine dams with larger potential for failure and increased GLOF risk (Linsbauer and others, 2016). Second, proglacial lakes aggravate glacier mass loss (King and others, 2019) – a factor not included in current ice loss projections. Despite the fact that lake-terminating glaciers currently only account for a small fraction of all glaciers in HMA, our calculations indicate the formation of proglacial lakes at nearly every glacier with overdeepenings >10⁵ m². Some of these lakes will newly form while others will increase the area of already existing ones leading to even larger lakes with a higher potential for catastrophic GLOF events. Lake Merzbacher in the Tien Shan is an example of such a situation. With a maximum depth of 100 m (Narama and others, 2017) and an estimated volume of ~0.17 km³ (Mayer and others, 2008), it is the largest proglacial lake in the region and outbursts almost every year, causing damage to downstream infrastructures and communities (Wortmann and others, 2014). Xie and others (2013) report that the outbursts at Lake Merzbacher are occurring earlier every year due to warming trends. Figure 15 depicts the current lake as well as a computed overdeepening under Engilchek Glacier with a maximum volume of 1.5 km³. Following the ongoing retreat of the glacier, the overdeepening can gradually fill and eventually connect with Lake Merzbacher. The abundance of supraglacial lakes over the lower 7 km of the glacier (indicated by red dots in Fig. 15) could well be an indicator for the development of a massive proglacial lake in the near future (Benn and others, 2012). In case the whole overdeepening will coalesce with Lake Merzbacher, the projected increase in area and volume (by ~400 and 900%, respectively) could lead to a substantial increase in both discharge rate and volume. Whether the lake would still be outbursting regularly at this point or only after rockfall or avalanche events cannot be determined. In both cases, however, a lake of this size would lead to a significantly exacerbated threat downstream.

Fig. 15. Comparison of the current Lake Merzbacher (light blue) and the potential lake extension (shaded dark blue) that would almost quadruple its size. Red dots indicate an area with a multitude of supraglacial lakes which could develop into a proglacial lake in the near future.

6. Conclusions

This study provides the first complete inventory of future glacial lakes >10⁴ m² in area in HMA by computing the bedrock topography of the ~100 000 glaciers in the region using ice thickness data of a five-model-ensemble by Farinotti and others (2019a) and the AW3D30 DEM. A total of 25 285 overdeepenings with a volume of 99.1 ± 29.5 km³ were computed covering 2683 ± 812 km². The location and size of the overdeepenings is in good agreement with previous studies. Particularly for the western Himalaya and the Karakoram range, however, our results suggest a smaller number of overdeepenings (54%) as well as less overdeepening volume (37%) than previously projected. Our analysis of the hazard potential for slopes adjacent to the 2700 lakes larger than 10⁵ m², generally, shows a very high impact predisposition for mainly the smaller lakes in high mountain regions such as the Tien Shan, Pamir, Karakoram and parts of the Himalaya. In regions with less pronounced relief, such as Inner Tibet and the Kunlun range, overdeepenings are less prone to be impacted by mass movements. When additionally considering lake volume, the most hazardous future proglacial lakes are those that are very vulnerable to mass-movement impacts and, at the same time, capable of releasing large outburst floods relatively close to infrastructure and settlements. According to our findings, the number of proglacial lakes is expected to increase substantially, particularly in the northwestern part of HMA. Considering the profound impact of proglacial lakes on glacier ablation at the ice front, we emphasize the need for improving the integration of the impact of proglacial lakes on glacier melt and the projection of glacier-mass loss. At this stage, the results of this study provide the first comprehensive overview on lake development and emerging GLOF hazard in HMA and can serve as a basis for further assessments of risks and opportunities related to future glacial lakes. Investigating the temporal perspective of lake formation and development as well as performing GLOF risk assessments including infrastructure and land-use would be significant future steps toward adapting to the emerging challenges from future glacial lake formation in HMA.