Mapping Glacial Lakes in the Upper Indus Basin (UIB) Using Synthetic Aperture Radar (SAR) Data

1. Introduction

Climate change is dramatically transforming glaciated regions around the world, resulting in widespread glacier mass loss and retreat [1,2,3]. This retreat is driving the accelerated formation and expansion of glacial lakes [4,5,6]. Between 1990 and 2018, the number of global glacial lakes (>0.05 km²) increased by 53%, reaching a total volume of 156.5 km³, which is a sea level equivalent of 0.43 mm [5]. A recent global inventory identified 117,352 glacial lakes (≥0.01 km²), with a cumulative area of approximately 25,000 km² [7]. As ongoing glacier retreat and lake expansion accelerate under anthropogenic climate change, glacial lake outburst flood (GLOF) activity is also expected to increase significantly over the 21st century [8].

The distribution of glacial lakes at regional scales (e.g., High Mountain Asia (HMA), the Alps, the Andes, Greenland, and Alaska) has been extensively studied using remote sensing datasets, specifically optical imagery [4,9,10,11,12,13,14]. Particularly for the HMA region, which includes the Hindu Kush, Himalayas, Karakoram, Tienshan Mountains, and Tibetan Plateau, several glacial lake inventories have been created [6,15,16,17,18]. The first comprehensive glacial lake inventory for HMA manually delineated 27,205 and 30,121 glacial lakes (>0.0054 km²), covering a combined area of 1,806 and 2,080 km² in 1990 and 2018, respectively [6]. Other important findings for the HMA include a total lake area increase of approximately 7% between 2008 and 2017 [18], and an annual expansion rate of 1.38% between 2008 and 2016 [17]. A summary of existing regional and global glacial lake inventories covering the HMA region is provided in the supplementary material (Table S1).

In addition to optical remote sensing approaches, several recent studies have explored the use of synthetic aperture radar (SAR) data for glacial lake mapping. SAR offers significant advantages over optical data, particularly its ability to acquire images regardless of cloud cover or illumination conditions [19]. A variety of methods have been applied to SAR imagery for glacial lake detection, including simple backscatter thresholding, random forest classifiers, and deep learning algorithms [20,21,22,23,24]. To further improve detection accuracy, some studies have combined SAR and optical imagery to leverage the strengths of both sensors [25,26,27].

Despite the abundance of published studies and datasets, several key research gaps persist. One major limitation concerns the temporal scope of the imagery used. To minimize cloud and snow cover, studies often rely on satellite imagery spanning extended seasonal windows, for example, June to November [6], July to November [18], August to December [28], and even multiple years [15]. Such wide timeframes, which are even broader in global inventories, reduce the temporal specificity of lake conditions, making it difficult to capture glacial lake dynamics during a consistent seasonal period.

A second limitation involves the underrepresentation of highly dynamic supraglacial lakes (SGLs), which are typically short-lived and fluctuate at sub-seasonal to seasonal timescales [22,23,29]. These lakes are particularly susceptible to omission due to both their small size and transient nature. The use of extended windows for image selection contributes to their exclusion, as SGLs may appear and disappear between acquisition dates. Additionally, the application of minimum lake area thresholds introduces further bias against smaller water bodies. Consequently, regional-scale inventories that use coarser thresholds and wider acquisition windows systematically underrepresent small and short-lived lakes (Table S1).

While SAR-based techniques address several limitations inherent to optical data, significant challenges remain. Processing SAR data in glaciated mountainous terrain is complicated by the complex nature of radar backscatter [19,30,31]. Moreover, achieving broader spatiotemporal transferability of SAR-based mapping methods remains difficult and often requires extensive training datasets [21,22]. As a result, no regional-scale SAR-based glacial lake inventory currently exists for the HMA region.

This study addresses the limitations of existing glacial lake inventories by utilizing Sentinel-1 SAR imagery to develop a comprehensive inventory for the Upper Indus Basin (UIB), one of the most glacierized regions in the world and a critical water supply for Pakistan. Leveraging the all-weather data acquisition capabilities of SAR [26], the analysis is confined to a narrow observation window in August 2023, to reduce interference from seasonal snow. This narrow window also enhances the temporal consistency of lake detection. We also analyze the spatial distribution and characteristics of the mapped lakes and systematically compare them to existing inventories to assess discrepancies and highlight improvements.

2. Materials and Methods

2.1. Site Description

The UIB (Figure 1) is a transboundary river basin located in the high mountain ranges of the Hindu Kush, Karakoram, Himalayas, and the Tibetan Plateau, spanning parts of Afghanistan, Pakistan, India, and China, with a total area of approximately 425,000 km² [32]. With over 17,000 individual glaciers and a total glacierized area exceeding 25,000 km², the UIB is among the most glaciated areas on the planet [33,34,35]. Glacier melt contributes up to 40% of annual streamflow, with peak discharge occurring during the summer months when melt rates are highest [36,37]. River flows exhibit pronounced seasonality, primarily driven by snow and glacier melt [38]. This seasonal meltwater supply is critical for downstream communities and much of Pakistan, where it supports agriculture, hydropower generation, and industrial activities [32,38,39].

The UIB has a complex climatology due to the interaction of multiple atmospheric circulation systems, including the westerlies and the Indian summer monsoon, and its highly variable topography [32,40,41]. According to the Köppen-Geiger classification, the region falls under a cold desert climate, characterized by warm, dry summers and cold winters [42]. The westerly flows, dominant during winters, are responsible for most of the total annual precipitation, while the Indian summer monsoon dominates during the summer months [43,44].

There is high variability in elevation across the UIB, ranging from 569 m in the valleys to over 8,000 m at mountain peaks, with a mean elevation of 3,750 m [45]. This topographic variation significantly influences spatial patterns of precipitation. Annual precipitation ranges from 100 to 200 mm in lower valleys, increases to approximately 600 mm around 4,400 m, and can exceed 1,500 to 2,000 mm above 5,500 m [43,46].

2.2. Data

For this study, we utilized multiple satellite and geospatial datasets to map glacial lakes. Our primary dataset was Sentinel-1A C-band SAR imagery (5.405 GHz, 5.6 cm), specifically the Ground Range Detected (GRD) collection available in Google Earth Engine (GEE). Sentinel-1A has a repeat cycle of 12 days, and the GRD product provides calibrated and geometrically corrected backscatter data at spatial resolutions of 10, 25, or 40 m and across various polarizations and modes [47,48]. We used the 10 m vertical transmit/vertical receive (VV) polarization and the Interferometric Wide Swath (IW) instrument mode data to delineate glacial lakes. VV polarization reduces interference from surface reflections and enhances the contrast between water and surrounding terrain, making it well-suited for glacial lake mapping [49].

To create normalized difference water index (NDWI), normalized difference snow index (NDSI) maps, and manually delineated validation datasets, we used Sentinel-2 MSI imagery with cloud cover below 10%. Slope masking was performed using digital elevation data from the Shuttle Radar Topography Mission (SRTM) [50]. The Randolph Glacier Inventory (RGI) glacier outlines were utilized to define the glacierized area [34].

2.3. Glacial Lake Mapping

The glacial lake mapping process involved three primary steps: (1) data preprocessing and preparation, (2) water body detection and lake delineation, and (3) accuracy assessment. Each step is described in the following subsections.

2.3.1. Data Preparation

Sentinel-1 GRD imagery in IW mode with VV polarization for August 2023 was filtered by date and study area and imported into GEE. While the GRD collection is already preprocessed to generate a calibrated and ortho-corrected product through thermal noise removal, radiometric calibration, and geometric correction [51], additional preprocessing was required to correct radiometric distortions and mitigate speckle noise. SAR datasets suffer from terrain-induced distortions and speckle noise due to radar backscatter interference [47,51,52]. The terrain flattening involved computing the local incidence angle by combining the radar incidence angle from the SAR image with terrain slope and aspect from the SRTM DEM, followed by gamma-nought normalization [53] using Equation 1,

$${\gamma }^{°}={\displaystyle \frac{{\sigma }^{°}}{\mathrm{c}\mathrm{o}\mathrm{s}\left({\theta }_{LIA}\right)}} \mathrm{E}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} 1$$

where ${\gamma }^{°}$ is the normalized backscatter, ${\sigma }^{°}$ represents the backscatter on a linear scale, and ${\theta }_{LIA}$ is the local incidence angle. The normalized backscatter was converted back to the dB scale and renamed as the new VV band.

To reduce noise while preserving features, we applied a 3x3 median filter, which replaces the center pixel intensity with the median of its eight neighboring pixels. This approach effectively suppresses extreme outliers while maintaining spatial details [23]. Where available, Sentinel-1 images from ascending and descending orbits with similar acquisition dates were combined to generate a median image. This step further reduced speckle noise and mitigated shadow effects, enhancing the quality of the data for glacial lake detection [27,54].

In addition to SAR data, auxiliary datasets were prepared to mask irrelevant features. An NDSI mask was generated from cloud-free Sentinel-2 MSI imagery, when available, to filter snow-covered areas. A slope mask was used to remove regions with slopes greater than 20 degrees, which commonly introduce classification errors due to terrain-induced radar scattering. These preprocessing steps significantly improved computational efficiency and reduced misclassifications caused by wet snow and mountain shadows, which exhibit backscatter values similar to glacial lakes.

2.3.2. Glacial Lake Delineation

The final preprocessed SAR dataset was used for lake delineation. Potential water bodies were initially identified using a lenient SAR backscatter threshold of -13 dB, classifying pixels with backscatter values below this threshold as water. Previous studies have employed a -14 dB threshold for glacial lake delineation in the Himalayas [26]. However, to minimize the exclusion of smaller lakes, a less restrictive threshold of -13 dB was used in this study. This simple thresholding approach led to the inclusion of some non-lake features, such as glacial ice, seasonal snow, and rugged terrain, resulting in classification errors. The water mask was vectorized into individual polygons, and key geometric and topographic attributes were extracted, including surface area, perimeter, and shape index (perimeter-to-area ratio). We applied a minimum shape index filter (0.2) to eliminate highly irregular-shaped non-lake features.

Since a global backscatter threshold is ineffective for accurately delineating glacial lakes of varying sizes, compositions, and surrounding features, we applied an adaptive thresholding method based on the Otus algorithm to reclassify the filtered polygons [55]. The Otsu algorithm analyzes the image histogram to obtain the variance for every possible threshold value and calculates the threshold that minimizes the inter-class variance. First, a buffer equal to the area was applied to each lake. For each buffered lake polygon, the Otsu algorithm was applied to identify the optimal backscatter threshold. We then used the optimized backscatter thresholds for each buffered polygon to reclassify them and obtain the final lake polygons. The associated attributes, including surface area, mean elevation, and centroid coordinates, were recalculated, and the lakes were filtered by a minimum lake area threshold of 1,000 m². This threshold was chosen to balance mapping accuracy with the retention of essential lake information. The lake polygons were overlaid onto false-color Sentinel-2 MSI composites for quality control. Manual inspection was conducted to remove false positives. Additionally, a few overestimated lake boundaries caused by surrounding ice were corrected manually using an NDWI mask. The final lake polygons, along with associated attributes including surface area, mean elevation, and centroid coordinates, were exported as shapefiles for further analysis.

The lakes were categorized into three major groups:

• Supraglacial lakes (SGLs) form in topographical depressions on the surface of glaciers;
• Moraine-dammed lakes (MDLs) form behind terminal or lateral moraines;
• Glacial erosion lakes (GELs) form in topographical depressions in the bedrock.

SGLs were identified using glacier outlines from the RGI inventory, while MDLs and GELs were identified by analyzing their location and distance from current glaciers, comparing with previous inventories, and through visual inspection using a false color composite of Sentinel-2 optical imagery.

2.3.3. Accuracy Assessment

To evaluate the accuracy of the SAR-based glacial lake mapping, we manually delineated 65 lakes using Sentinel-2 MSI false-color composites. This validation dataset included glacial lakes of various sizes and types across the UIB. This dataset served as the reference for validating the lake mapping approach. We assessed the accuracy using mean absolute error (MAE), mean absolute percentage error (MAPE), root mean square error (RMSE), and the coefficient of determination (R²). These metrics quantify both the magnitude of error and the agreement between SAR-derived and reference lake areas. The MAE and MAPE were obtained using Equations 2 and 3, respectively, as follows

$$MAE={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n\left|{Area}_{S2}-{Area}_{S1}\right| \mathrm{E}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} 2$$

$$MAPE={\displaystyle \frac{100}{n}}{\sum }_{i=1}^n\left|{\displaystyle \frac{{Area}_{S2}-{Area}_{S1}}{{Area}_{S2}}}\right| \mathrm{E}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} 3$$

where n is the number of lakes in the validation dataset, Area_S1 is the area of each lake obtained from the Sentinel-1 imagery, and Area_S2 is the area obtained using manual delineation from Sentinel-2 optical imagery.

2.3.4. Glacial Lake Size Distribution

We assessed the statistical distribution of glacial lake sizes in the UIB using the powerlaw Python package, which evaluates the plausibility of power-law models relative to alternative distributions [56]. A power-law model was fitted to the observed lake size data, and the power-law exponent (α) and was estimated for lakes exceeding a minimum area threshold (x_min). This follows the general form of the power-law model:

P(x) ∝ x-α, for x ≥ x_min

where P(x) is the probability of observing a glacial lake of area x.

The Python package automatically determined the minimum lake area threshold through an iterative procedure that identifies the cutoff yielding the best fit, based on minimizing the Kolmogorov-Smirnov (KS) statistic and maximizing the p-value (Clauset et al. 2009). The KS statistic represents the maximum difference between the empirical and fitted cumulative distribution functions. This analysis was applied across the entire UIB and repeated for each of its four major subregions and lake types to assess variations in scaling behavior.

3. Results

3.1. SAR-Based Lake Mapping and Accuracy Assessment

The automated SAR-based lake outlines demonstrate a high level of agreement with the manually digitized lakes from Sentinel-2 optical imagery, as illustrated by the 1:1 plot (Figure 2) and the overlapping lake outlines (Figure 3). The MAE and MAPE are 0.013 km² and 8.18%, respectively, with minimal bias (0.002 km²) (full Willmott statistics in Table S2). The SAR-based estimated lake areas (mean = 0.274 km²) are equivalent to manually delineated lake areas (mean = 0.272 km²) on average. For a few moraine-dammed lakes in contact with glacier ice, areas are overestimated by the SAR-based method, likely due to the low backscatter contrast between water and adjacent ice or wet snow, leading to boundary misclassifications, as illustrated in Figure 3e and Figure 3f.

Otsu-derived backscatter thresholds range from –13 to –28 dB, with a mean of –20 dB. The variability in threshold values, and their relationship with lake characteristics such as area and shape, underscores the importance of adaptive thresholding (Figure 4). A strong negative relationship was observed between lake area and threshold values (Spearman’s ρ = –0.74, p < 0.001), indicating that larger lakes tend to require lower thresholds. This reflects the fact that expansive water surfaces are easier to classify due to their more homogeneous backscatter. In contrast, lake shape shows a significant positive relationship with threshold values (ρ = 0.77, p < 0.001). Lakes with higher shape indices, reflecting more irregular outlines, require higher thresholds to achieve accurate delineation.

3.2. Glacial Lake Area and Elevation Distribution

A total of 6,019 glacial lakes larger than 0.001 km² were identified across the UIB (Figure 5a), covering a cumulative area of 266 km², equivalent to 0.06% of the basin area (Table 1). Lake sizes range from 0.001 to 5.80 km², with a mean area of 0.04 km². Of these, 39.16% of the lakes are classified as small (< 0.01 km²), 53.08% as medium (0.01 to 0.1 km²), and 7.76% as large (> 0.1 km²). Despite their low frequency, large lakes account for 57% of the total lake area, while small lakes contribute only 4%.

The lake elevations range from 2,557 to 6,005 m, with peak lake concentrations between 4,000 and 4,600 m (Figure 5b). The median lake area remains relatively consistent across higher elevation bands. At the lowest elevations (2,400 – 2,600 m), the largest median areas (approximately 0.16 km²) occur, driven primarily by a few large lakes in the Himalayas and the Hindu Kush. Between 2,800 and 3,000 m, the median area declines sharply to its minimum (0.002 – 0.004 km²), reflecting the dominance of numerous small SGLs in the Karakoram. Above 3,400 m, the median area remains relatively small but stable at around 0.015 km² up to 5,800 m.

Figure 5. Glacial lake distribution and characteristics across the study area and sub-regions. (a) Map of the UIB study area, showing sub-region boundaries and delineated lakes (blue). (b) Number of lakes and lake area by elevation for the entire UIB basin, with lake counts shown as bars and lake area as boxplots within 200 m elevation bins. (c-f) Similar elevation-based distributions for the sub-regions: Hindu Kush, Karakoram, Himalayas, and Tibetan Plateau, respectively.

Figure 5. Glacial lake distribution and characteristics across the study area and sub-regions. (a) Map of the UIB study area, showing sub-region boundaries and delineated lakes (blue). (b) Number of lakes and lake area by elevation for the entire UIB basin, with lake counts shown as bars and lake area as boxplots within 200 m elevation bins. (c-f) Similar elevation-based distributions for the sub-regions: Hindu Kush, Karakoram, Himalayas, and Tibetan Plateau, respectively.

Table 1. Lake statistics, including lake number, size distributions, and mean elevation in the UIB and its sub-regions.

Region	Basin Area (km2)	Lake Count	Total Lake Area (km2)	% of Basin Area	Mean Area (km2)	Mean Lake Elevation (m)
UIB	425,000	6,019	266.0	0.06	0.044	4,566
Hindu Kush	134,000	2,171	89.4	0.06	0.041	4,323
Karakoram	40,000	866	11.4	0.03	0.013	4,249
Himalayas	201,000	2,615	129.6	0.06	0.050	4,747
Tibetan Plateau	50,000	367	35.6	0.07	0.100	5,469

Table 1. Lake statistics, including lake number, size distributions, and mean elevation in the UIB and its sub-regions.

Region	Basin Area (km2)	Lake Count	Total Lake Area (km2)	% of Basin Area	Mean Area (km2)	Mean Lake Elevation (m)
UIB	425,000	6,019	266.0	0.06	0.044	4,566
Hindu Kush	134,000	2,171	89.4	0.06	0.041	4,323
Karakoram	40,000	866	11.4	0.03	0.013	4,249
Himalayas	201,000	2,615	129.6	0.06	0.050	4,747
Tibetan Plateau	50,000	367	35.6	0.07	0.100	5,469

The Himalayas and Hindu Kush, which collectively occupy a significant portion of the UIB, have the highest number of lakes with 2,615 and 2,171, and cover total areas of 129.6 km² and 89.4 km², respectively. In both regions, medium-sized lakes dominate by count, comprising 59% of the total, and contribute 38-47% to the cumulative area. Large lakes, although fewer in number (7-9%), account for a substantial share of the lake area (48-58%). Similar to the overall UIB, most lakes are concentrated between 4,000 to 4,600 m elevations in these two sub-regions, with a slight decrease in median area with elevation (Figure 5c and Figure 5e).

The Karakoram has 866 lakes and a comparatively smaller total lake area of 11.4 km². The region is dominated by small lakes (81% by count). However, the few large lakes (3%) contribute 44% to the total lake area. Most lakes are found between 3,600 and 4,200 m, and the median lake area, which is the lowest among the sub-regions, shows slight increase with elevation (Figure 5d).

The Tibetan Plateau’s 367 lakes cover 35.6 km², over three times the lake area of the 866 Karakoram lakes. Large lakes, accounting for 12% of the total number, contribute a substantial 79% of the total lake area, highlighting their dominance. Medium-sized lakes make up 58% of the count but contribute only 18% of the area. This region is characterized by higher elevations, and its lakes are distributed between 4,800 and 6,000 m (Figure 5f).

3.2. Glacial Lake Type Distribution

Glacial lake types across the UIB and its sub-basins show distinct regional patterns (Figure 6). At the basin scale, GELs and MDLs dominate, comprising 54% and 34% of all lakes, and contributing 61% and 37% to the total lake area, respectively. Although SGLs represent 12% of the lake count, they account for only 2% of the total area. In terms of typical size, GELs and MDLs have similar median areas (0.018 and 0.017 km², respectively), whereas SGLs are distinctly smaller with a median area of 0.003 km².

The Himalayas, Hindu Kush, and Tibetan Plateau show broadly similar distributions, all characterized by the dominance of GELs and MDLs, which together account for >96% of the total count and >98% of the cumulative area. SGLs are rare in the Himalayas and Hindu Kush, contributing less than 1% to the total lake area, and are entirely absent in the Tibetan Plateau.

The Karakoram sub-region has a markedly different distribution of glacial lake types. Here, SGLs are the most prevalent, comprising 77% of the lakes and 40% of the total lake area. Although MDLs represent only 20% of the lakes, they account for the largest share of the area (54%), reflecting their larger typical sizes relative to SGLs. GELs are rare in the Karakoram, contributing just 2% of the lake count and about 6% of the cumulative area.

Figure 6. Distribution of glacial lake types across the UIB study area and sub-regions. (a) Map of the study area, showing the distribution of glaciers and glacial lakes in the UIB. Stacked bar plots showing (b) the number of lakes and (c the cumulative lake area by lake type (GEL, MDL, SGL) across the UIB and the four sub-regions (Hindu Kush, Karakoram, Himalayas, and Tibetan Plateau).

Figure 6. Distribution of glacial lake types across the UIB study area and sub-regions. (a) Map of the study area, showing the distribution of glaciers and glacial lakes in the UIB. Stacked bar plots showing (b) the number of lakes and (c the cumulative lake area by lake type (GEL, MDL, SGL) across the UIB and the four sub-regions (Hindu Kush, Karakoram, Himalayas, and Tibetan Plateau).

3.2. Power Law Distribution of Glacial Lake Area

The size distributions of glacial lakes in the UIB follow power-law behavior, indicating that small lakes are far more abundant than large ones (Figure 7). Across the entire basin, the power-law exponent (α) is 2.26 with a minimum lake area of 0.059 km².

The lake size distributions for the subregions have varying scaling behavior and minimum lake area thresholds (Table S3). The Himalayas exhibit a similar power-law distribution (α = 2.20) as the UIB, with a slightly lower area threshold of 0.041 km². In the Hindu Kush, the power-law fit was limited to the largest lakes (x_min = 0.118 km²) with a steep slope (α = 2.67), reflecting the strong dominance of small lakes and scarcity of large ones.

In contrast, the Karakoram displays a much flatter slope (α = 1.91, x_min = 0.003 km²), implying a higher relative frequency of large lakes. The scaling holds for very small lakes compared to other regions. The Tibetan Plateau has the flattest slope overall (α = 1.79, x_min = 0.016 km²), consistent with the prevalence of large lakes.

Scaling also varies by lake type (Figure 7f, g, h). GELs follow a steep slope (α = 2.32, x_min = 0.07 km²), reflecting the dominance of smaller lakes. MDLs show a slightly flatter slope (α = 2.19, x_min = 0.04 km²), indicating a relatively greater contribution from larger lakes, although both types are overall dominated by small lakes. SGLs also show steep scaling (α = 2.30), but with a very low threshold (x_min = 0.005 km²), reflecting their smaller sizes and high abundance.

The scaling analysis highlights two distinct regimes of glacial lake size distributions across UIB. The Himalayas and Hindu Kush exhibit steeper exponents, indicating that lake size distributions decay rapidly and that large lakes are relatively scarce. In contrast, the Karakoram and Tibetan Plateau have flatter slopes, reflecting a more gradual decay of the size distribution and a greater relative contribution from large lakes. In all cases, the lake size distributions deviate from power-law behavior below the estimated x_min, reflecting a relative scarcity of the smallest lakes compared to what would be expected under strict scaling.

4. Discussion

4.1. SAR-Based Glacial Lake Inventory and Comparison with Other Datasets

SAR offers a viable alternative to optical remote sensing-based glacial lake delineation, as it enables delineation under cloud cover and variable illumination conditions. In this study, we combined Sentinel-1 SAR data with adaptive backscatter thresholding to generate a new inventory of glacial lakes in the UIB. The SAR-based outlines show strong agreement with manually digitized reference lakes, confirming the reliability of the approach. Minor overestimation occurred where lakes bordered wet snow or glacier ice, likely due to the similarly low backscatter characteristics of these adjacent surfaces, a challenge also noted in previous SAR-based studies [20,57]. In addition, manual corrections were required to exclude misclassifications in such areas, highlighting a key limitation of SAR-based delineation where surface types have overlapping backscatter signatures.

To evaluate consistency with existing inventories, we compared our results with the dataset of Wang et al. (2020) [6]. After filtering their data to match our study area, Wang et al. (2020) [6] mapped 4,035 glacial lakes (≥ 0.0054 km²) covering a total area of 220.8 km² in 2018. In comparison, our 2023 SAR-based inventory identified 6,019 lakes (≥ 0.001 km²) in 2023, covering a total area of 266 km². When we applied the same minimum area threshold (≥ 0.0054 km²), our dataset contained 4,617 lakes covering 261.8 km², showing 87% agreement in lake number and 84% agreement in cumulative area relative to Wang et al. (2020) [6]. When both inventories were further aggregated to a 0.1° × 0.1° grid, 1:1 comparison plots for lake count and area showed very strong agreement between the two inventories (Figure 8).

The higher lake count and total area in our inventory compared to Wang et al. (2020) [6], even after applying the same minimum area threshold, can be attributed to several factors. First, our inventory is based on 2023 data, whereas Wang et al. (2020) [6] mapped lakes in 2018. Given the continued retreat of glaciers and the observed global trend of glacial lake expansion [3,5], some increase in lake number and area is expected over the five-year period. However, consistent and comparable mapping techniques are needed to identify trends over time in the number and size of lakes. Second, the use of Sentinel-1 SAR imagery in our study enabled consistent lake detection under cloudy or shadowed conditions, where optical sensors often face limitations [26]. However, while post-processing steps minimized major delineation errors, some overestimation of lake extents may persist in our results, potentially contributing to the slightly higher cumulative area compared to the optical-based inventory.

Finally, our analysis was restricted to a single month (August), to coincide with peak meltwater availability and increase the likelihood of capturing seasonally transient SGLs. In contrast, Wang et al. (2020) [6] used a broader observation window (June to November), which may have led to the omission of short-lived lakes. This difference, combined with the all-weather capabilities of SAR, likely contributed to the 223 SGLs identified in our inventory compared to 91 in Wang et al. (2020) [6]. This comparison highlights the strengths of SAR-based mapping, particularly in high-mountain regions with persistent cloud cover and dynamic lake evolution. Given its ability to penetrate cloud cover and capture dynamic lakes, such as the SGLs, our SAR-based approach likely offers a more comprehensive inventory of glacial lakes in the UIB. Additionally, the temporal update to 2023 provides a more current representation of lake distribution compared to earlier inventories.

3.2. Spatial Characteristics of Glacial Lake Distribution

The spatial distribution and characteristics of glacial lakes are influenced by glacier dynamics, topographic controls, and climatic influences [9,58]. The lake size, spatial, and altitudinal distribution in UIB are consistent with observations from other mountain regions. Across the UIB, lake size is skewed toward small and medium lakes by number, while a relatively small number of large lakes (> 0.1 km²) account for more than 50% of the total lake area. An abundance of smaller lakes is common across the glaciated regions. For example, in the Alps, lakes larger than 0.1 km² made up only 14% of the total lakes but accounted for 52% of the total lake area [11]. In the UIB, most glacial lakes are concentrated between 4,000 and 4,600 m in elevation, coinciding with the elevation range where the majority of glaciers are located [59,60]. Large lakes are typically found at lower elevations compared to small and medium lakes, a pattern attributed to the flatter terrain in these areas, which offers more favorable conditions for lake formation and expansion [11,61]. This relationship between terrain and lake development was particularly evident on the flatter Tibetan Plateau, where some of the largest lakes in the study area are located [62]. Additionally, increasing precipitation and warming temperatures since the late 1990s, along with permafrost degradation, have contributed to ongoing lake growth in this region [63,64].

GELs and MDLs are the most prevalent lake types in the UIB. Their dominance has been observed in other regions, including the greater Himalayas [58,65,66,67,68], the Tibetan Plateau [69], the Alps [12,70], the Andes [4,71], and Canada [72]. The abundance of GELs and MDLs is indicative of an advanced stage of glacier retreat, leading to the construction of substantial moraine dams and exposure of overdeepened bedrock basins, where MDLs and GELs form [9,14,73,74]. Most GELs are detached from active glaciers and respond more to changes in regional water balance than to direct glacier melt [14,75].

Across the subregions, glacial lake abundance is highest in the Himalayas and Hindu Kush, dominated by GELs and MDLs. This pattern reflects widespread glacier retreat, the presence of substantial moraines, and suitable topography for lake impoundment. In the Himalayas, approximately 40% of glacier area has been lost since the Little Ice Age (LIA), consistent with the widespread occurrence of GELs and MDLs [76]. Based on the RGI glacier outlines, currently both the Himalayas and Hindu Kush subregions have a 4% glacierized area [34]. In contrast, the Karakoram shows limited glacier area loss, with only 6.6% decline since the LIA, while some glaciers are advancing or showing stable mass balance [77,78]. Glacier cover remains extensive, at 41% of the total area. Under these conditions, the region is predominantly characterized by SGLs. The smaller number of MDLs and GELs in the Karakoram region, as opposed to their dominance in the Himalayas and Hindu Kush, shows the influence of the Karakoram Anomaly on glacial lake formation dynamics. With glacier stability and surge behavior dominating in this region, bedrock depressions and moraine dams are not available for extensive GEL and MDL formation, while the debris-covered flatter ablation zones of the Karakoram glaciers support extensive SGL formation.

The formation of SGLs is closely associated with debris-covered glacier tongues [18,79]. Extensive debris cover in this region leads to differential thinning of glaciers, enabling meltwater to accumulate in depressions on glacier surfaces [80,81]. SGL area is positively correlated with debris-covered glacier area and the glacier length, because longer glaciers provide ideal conditions for the coalescence of supraglacial ponds [18,79,82]. SGLs are typically smaller, more variable, and transient compared to GELs and MDLs [14,54]. However, they can coalesce and evolve into larger moraine-dammed lakes, influencing glacier dynamics and stability [83,84]. For example, lake-terminating glaciers experience enhanced surface lowering, retreat, and increased ice velocities, as observed in parts of the Himalayas [85,86,87]. The varying composition and dynamics of glacial lake types across the UIB thus provide valuable insights into regional glacier behavior and the broader impacts of climate change on high mountain hydrology.

3.2. Power-Law Distribution of Glacial Lake Area

The power-law behavior observed in the UIB and its subregions highlights the scale-invariant nature of glacial lake size distributions. However, variations in the power-law exponent and minimum area threshold across the subregions reflect the influence of regional geomorphology and dominant lake types. For example, the prevalence of smaller-sized SGLs on the debris-covered glaciers of Karakoram led to a power-law adherence that extended to lakes an order of magnitude smaller than other regions. In the Hindu Kush and Himalayas, where GELs and MDLs dominate, scaling applies only above larger minimum sizes compared to the SGL-rich Karakoram.

The power-law distribution of glacial lake sizes in the UIB is consistent with scaling behavior observed in other mountainous and glaciated regions for lake inventories [88,89] and in lake volume-area scaling [90,91], supporting the universality of such patterns in hydrological and geomorphological processes [92,93]. Most of the studies’ minimum areas are an order of magnitude larger than this study. For example, Cael and Seekell (2016) [88] reported α = 2.14 for global lakes ≥ 0.46 km²; by comparison, our study yields α = 2.26 for lakes ≥ 0.06 km². Levenson et al. (2025) [89] found α = 2.23 for Alaskan glacial lakes ≥ 1.01 km².

Modeling lake size distributions using a power-law framework provides both a conceptual and quantitative tool to describe and compare glacial lake systems. These models enable statistical extrapolation beyond observed data, making them particularly valuable in glaciated basins where small lakes are numerous but often underrepresented due to resolution or mapping constraints. However, deviations below minimum area thresholds emphasize the need for careful parameterization when applying scaling models to lake inventories.