the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 16 Jun 2026
High-resolution glacier mapping reveals inventory biases and terrain controls on debris-covered glaciers in the Karakoram
Abstract. Accurate glacier inventories are fundamental for quantifying glacier change, estimating ice volume and assessing meltwater resources. However, medium-resolution inventories often fail to resolve critical glacier characteristics in topographically complex and debris-rich mountain environments. Here, we present the 2 m Karakoram Glacier Inventory (2mKGI), developed from high-resolution ZY-3 optical imagery, a co-registered ZY-3 digital elevation model (DEM) and auxiliary optical datasets through an integrated deep-learning and manual-refinement framework. The inventory identifies ~13,900 glaciers covering 21,261.8 ± 278 km², including 2,239.9 ± 82.8 km² of supraglacial debris, with an overall mapping uncertainty of ±4.7 %. Comparison with existing inventories reveals that previous medium-resolution products commonly underestimate glacier numbers while simultaneously overgeneralizing debris-covered glacier margins. These biases substantially influence glacier-count statistics, estimates of debris-covered area, and the interpretations of glacier-change signals. The newly identified glaciers are mostly <0.1 km². Despite their limited area, their thin ice and shorter response times may make them particularly sensitive to warming. Topographic analysis further demonstrates that supraglacial debris is preferentially distributed across low-elevation and low-slope glacier tongues, highlighting the strong controls of valley geometry, ice transport processes, and ablation-zone morphology on debris persistence. The 4 m ZY-3 DEM highlights that high-resolution topographic information improves the delineation of glacier-units, accumulation–ablation zone structure and debris-covered tongues by preserving steep headwalls, slope discontinuities, tributary junctions, local relief and low-gradient terrain. The 2mKGI will provide a high-resolution geometric and topographic benchmark for glacier-change assessment, ice-thickness inversion, glacier-evolution modelling and next-generation automated glacier mapping in the Karakoram.
- Preprint
(4175 KB) - Metadata XML
-
Supplement
(749 KB) - BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2026-2836', Anonymous Referee #1, 08 Jul 2026
-
AC2: 'Response to Reviewer 1', Xin Yang, 05 Sep 2026
Response to Reviewer 1
We sincerely thank the reviewer for the positive assessment of our manuscript and for the constructive comments. In response, we have revised the manuscript to clarify the model configuration, LST processing, manual delineation criteria, treatment of nunataks, and identification of surge-type glaciers. We have also improved the reproducibility and usability of the dataset by updating the code-availability statement and attribute-table documentation and by adding quantitative assessments of imagery-source uncertainty and DEM-source sensitivity. The manuscript and dataset statistics have been updated to match the final revised inventory. Detailed responses to each comment are provided below. For readability, the reviewer’s comments are shown in blue and our responses in red.
Comment 1
Figure 3 provides a clear overview of the U-Net+CBAM workflow. However, the specific activation functions (e.g., ReLU, Sigmoid) used in the convolutional blocks are not explicitly stated. Including these details in the figure caption or the main text would help the reader better understand the model's non-linear fitting capabilities.
Response
We thank the reviewer for this helpful suggestion. We have revised Figure 3 and Section 3.3.1 to specify the activation functions used in the network. All standard convolutional layers use ReLU activation. Within the CBAM, the channel-attention branch uses a shared MLP with a ReLU-activated reduction layer followed by sigmoid gating, and the spatial-attention branch uses a 7×7 convolution followed by sigmoid gating. The final segmentation layer uses softmax to generate pixel-wise probabilities for the three classes.
Comment 2
As the glacier boundary extraction relies heavily on deep learning, I strongly recommend providing the core training and inference code (or pseudocode) via a public repository such as GitHub. This transparency is crucial for the reproducibility of the methodology.
Response
We agree that access to the implementation is important for reproducibility. We have created a publicly accessible GitHub repository containing the core training and inference scripts for the U-Net+CBAM workflow, together with instructions for preparing the input data and reproducing the main processing steps. The repository link has been added to the Code Availability statement (https://github.com/YNUYangX/Debris-covered-glacier-extraction.git).
Comment 3
The manuscript mentions that Land Surface Temperature (LST) was retrieved from Landsat 8 imagery. Given that the ZY-3 imagery was acquired between July and October 2020–2021, could you clarify whether the LST represents a single scene matching the specific ZY-3 acquisition date, or if it is a composite/average product for the ablation season (July–October)?
Response
We appreciate this request for clarification. The LST layer was not derived from one Landsat 8 scene matched exactly to each ZY-3 acquisition. Instead, for each corresponding year we generated an end-of-ablation-season fused LST product from quality-controlled Landsat 8 observations. This product was co-registered with the optical and topographic layers and used as an auxiliary input to reduce the influence of individual-scene noise and residual cloud contamination. We have clarified this processing choice in the revised manuscript.
Comment 4
In the manual refinement phase, the manuscript mentions using features such as lateral moraines and meltwater outlets as geomorphological indicators. However, the distinction between debris-covered ice and lateral moraines is not explicitly detailed. Please clarify the specific visual texture features or criteria used to manually differentiate debris-covered glaciers from lateral moraines.
Response
We thank the reviewer for identifying this missing methodological detail. We have expanded the manual-refinement description. Debris-covered ice was interpreted from its continuity with up-glacier ice, flow-compatible valley morphology, ice cliffs and/or meltwater outlets, and the associated slope breaks and surface-roughness pattern. In contrast, lateral moraines were treated as discontinuous, ridge-like debris accumulations adjacent to the glacier margin that lacked continuity with the ice-flow path and the diagnostic ice-surface features. ZY-3 imagery, DEM-derived slope, LST, and Google Earth imagery were considered jointly for ambiguous margins.
Comment 5
The legend for Figure 6 includes ‘Karakorum’, but there are no corresponding boundaries, nor are any required here. Could this be an error in the legend?
Response
We thank the reviewer for noticing this error. The 'Karakorum' legend item was introduced inadvertently during figure preparation and does not correspond to a required boundary layer. We have removed it from the revised Figure 6 legend.
Comment 6
The manuscript provides a total glacier area of 21,261.8 km². According to standard glacier inventory protocols (e.g., RGI standards), internal bedrock outcrops (Nunataks) should typically be excluded. Please confirm whether the 2mKGI dataset excludes these Nunataks. If they are retained, what is their approximate proportion of the total area?
Response
We thank the reviewer for raising this point. We confirm that internal bedrock outcrops (nunataks) were excluded from the glacier polygons. Exposed bedrock was classified as non-glacier during pixel-level segmentation, and internal holes and exposed-rock areas were additionally checked during manual refinement.
We also rechecked the area reported in the previous manuscript against the final revised vector dataset. The previously reported value of 21,261.8 km² was not consistent with the final revised inventory. The corrected inventory contains 13,272 glacier units with a total mapped glacier area of 19,582.3 km². This value represents clean ice plus supraglacial debris-covered ice and excludes mapped nunataks.
Comment 7
Regarding Section 4.4, which identifies 192 surge-type glaciers, please clarify the data source for these "known" surge-type glaciers. Were they referenced from existing literature, or were they identified independently based on surface morphology observed in the ZY-3 imagery?
Response
We appreciate this request for clarification. The surge-type glacier attribute was compiled from previously published inventories and literature rather than being identified independently from the ZY-3 imagery alone (Bhambri et al., 2017; Guillet et al., 2022; Guo et al., 2023). We rechecked the final revised attribute table and identified 201 glaciers classified as confirmed surge-type glaciers. The previous value of 192 has therefore been corrected in Section 4.4.
Bhambri, R., Hewitt, K., Kawishwar, P., and Pratap, B.: Surge-type and surge-modified glaciers in the Karakoram, Scientific Reports, 7, 15391, 2017.
Guillet, G., King, O., Lv, M., Ghuffar, S., Benn, D., Quincey, D., and Bolch, T.: A regionally resolved inventory of High Mountain Asia surge-type glaciers, derived from a multi-factor remote sensing approach, The Cryosphere, 16, 603-623, 2022.
Guo, L., Li, J., Dehecq, A., Li, Z., Li, X., and Zhu, J.: A new inventory of High Mountain Asia surging glaciers derived from multiple elevation datasets since the 1970s, Earth Syst. Sci. Data, 15, 2841-2861, 2023.
Comment 8
The study notes that 21.5% of the glacier boundaries were supplemented by 10m-resolution Sentinel-2 imagery where ZY-3 coverage was unavailable. Could you provide a breakdown of whether the overall extraction uncertainty of ±4.7% differs significantly between the ZY-3 (2m) and Sentinel-2 (10m) data sources?
Response
We thank the reviewer for raising this important point. The reported +/-4.7% uncertainty is the estimate for the final, manually refined regional inventory; it is not a separately validated uncertainty value for the ZY-3 and Sentinel-2 subsets. We have revised the manuscript to make this distinction explicit. Both data-source subsets underwent the same manual quality-control procedure, including review of clean-ice, debris-covered, and shadow-affected boundaries. The Sentinel-2-only areas are concentrated mainly in the eastern Karakoram, where debris-covered terrain is relatively limited; therefore, the principal source of local uncertainty remains debris-covered and topographically complex margins rather than the data-source label alone. We now record the boundary-imagery source in the attribute table.
Comment 9
The 2mKGI dataset provides high-resolution benchmarks for approximately 13,900 glaciers. As this paper is intended to facilitate data sharing, I recommend stating in the text or supplementary material whether the attribute table of the final vector Shapefile includes standard GLIMS or RGI fields (e.g., standard glacier ID format, specific acquisition dates, area). This would greatly enhance the usability and interoperability of the dataset for the global glaciological community.
Response
We agree that interoperable attributes are essential for reuse. We have updated the data documentation and attribute table. Each glacier now has a persistent 2mKGI identifier generated from its central geographic coordinates, together with area, acquisition-year information, geographic location, elevation, slope, aspect, DEM source, and supraglacial-debris area. Where a spatial correspondence exists, the relevant GLIMS and/or RGI identifier is provided as a cross-reference; newly delineated glaciers without a pre-existing record retain the 2mKGI identifier. The revised data-release documentation explains these fields and the matching procedure.
Comment 10
Glacier unit subdivision uses ZY-3 4 m DEM where available and ASTER GDEM V3 elsewhere. The two DEMs differ in resolution, vertical accuracy, and ridge-line definition, which will introduce systematic inconsistencies in glacier unit boundaries across the mosaic edge. Please assess ridge-line positional differences between the two DEMs in overlapping areas, and report how many glacier units straddle the DEM mosaic boundary. Also discuss whether the use of ASTER GDEM introduces detectable biases in glacier count, elevation statistics, or slope metrics for the affected subregions.
Response
We thank the reviewer for highlighting the potential inconsistencies introduced by the use of two DEM sources. We have now assessed this issue explicitly from the perspectives of DEM-boundary effects, ridge-line position, glacier-unit subdivision, and derived topographic attributes.First, we delineated the valid-data footprint of the original 4 m ZY-3 DEM and intersected its boundary with the revised glacier inventory. Of the 13,272 glacier units, 337 (2.54%) straddle the ZY-3–ASTER mosaic boundary. These boundary-crossing units were manually re-examined using the optical imagery together with RGI v7.0 and KGI-2020s, and the initial DEM-derived drainage divides were corrected where necessary to avoid artificial splitting or merging caused by the change in DEM source. A local sensitivity analysis showed no significant discontinuity in glacier-unit size across the DEM boundary within either 5 km (p = 0.821) or 10 km (p = 0.139), providing no evidence of a systematic boundary-related effect on the final glacier-unit subdivision or count. Second, we quantified the differences between the two DEMs in their overlapping area. A paired comparison at 11,128 glacier-unit centres showed a median ZY-3-minus-ASTER elevation difference of −26.7 m, indicating a detectable vertical offset between the two elevation products. After harmonizing both DEMs to the same spatial scale, the median slope difference was only +0.28°, although larger local differences remained in steep terrain. Ridge-position analysis further showed that 91.1% of corresponding ridge features were within 120 m and 95.9% were within 300 m. These results indicate that ASTER GDEM introduces a measurable elevation offset and stronger local smoothing of slopes and secondary ridge features, but these differences do not translate into a detectable systematic change in the final glacier-unit count after manual correction.We also evaluated the broader effect of DEM resolution on glacier-scale morphometric metrics using 1,781 glaciers ≥1 km² with at least 95% valid coverage in both DEMs. Using identical glacier outlines, we compared differences in mean slope, 90th-percentile slope, 90 m surface roughness, and 90 m local relief. The results confirm that the 30 m DEM produces a systematic smoothing effect in local terrain representation, but this effect is widespread across glacier sizes rather than concentrated in a particular size class. We have added these quantitative comparisons and the corresponding discussion to the revised manuscript.Comment 11
The inventory includes glaciers down to 0.005 km² (~1250 pixels at 2 m resolution). No justification is given for this low threshold, and there is no discussion of whether very small features may be perennial snow patches.
Response
Thank you for raising this important point. The 0.005 km² threshold was not adopted as a universal physical definition of a glacier, but was determined empirically from the smallest glacier bodies that could still be identified with reasonable confidence in the revised inventory. During our glacier-by-glacier re-examination, we found several very small residual ice bodies that correspond to glacier units already represented in RGI v7.0 but had undergone substantial shrinkage by around 2020. The smallest of these confidently identified residual glaciers had an area of approximately 0.005 km². We therefore used 0.005 km² as the operational minimum mapping unit so that such genuine remnant glaciers would not be removed simply because their present-day area had fallen below 0.01 km².
At 2 m spatial resolution, 0.005 km² corresponds to approximately 1,250 pixels, which also allows these features to be inspected reliably at the image scale. However, area alone was not used to determine whether a feature was a glacier. All very small polygons were re-examined using multi-temporal Sentinel-2 imagery from 2020, 2022, and 2024 together with the available ZY-3 imagery. Features that appeared only as transient snow, or for which persistent glacier ice could not be supported by the multi-temporal observations, were removed. For retained features, persistence through time and spatial correspondence with previously mapped glacier units provided additional evidence for their interpretation as residual glacier ice.
We nevertheless acknowledge that some very small, persistently snow-covered features remain difficult to distinguish unequivocally from perennial or long-lasting snow patches using remote sensing alone. We have therefore clarified in the revised manuscript that 0.005 km² is a dataset-specific operational threshold rather than a general minimum glacier size, and that the identification of very small glaciers relies on multi-temporal and contextual evidence rather than area alone.
Citation: https://doi.org/10.5194/egusphere-2026-2836-AC2
-
AC2: 'Response to Reviewer 1', Xin Yang, 05 Sep 2026
-
RC2: 'Comment on egusphere-2026-2836', Frank Paul, 07 Aug 2026
Please find my review in the attached document.
-
AC1: 'Response to Professor Frank Paul', Xin Yang, 05 Sep 2026
We sincerely thank you for reviewing our manuscript and for the considerable time and effort you devoted to evaluating both the manuscript and the 2mKGI dataset. We greatly appreciate the detailed mapping examples, critical observations, and constructive recommendations provided throughout your review. These comments were highly valuable in identifying important issues in the original submission and in guiding a thorough revision of both the dataset and the manuscript. We have carefully considered all of your comments and made substantial revisions accordingly. Following completion of the revision, we will also upload the updated version of the glacier inventory to the data repository so that the publicly available dataset is consistent with the revised manuscript. Our detailed point-by-point responses and the corresponding changes are provided in the attached response document.
-
AC1: 'Response to Professor Frank Paul', Xin Yang, 05 Sep 2026
-
RC3: 'Comment on egusphere-2026-2836', Anonymous Referee #3, 11 Aug 2026
The authors present 2mKGI, a Karakoram glacier inventory mapped from 2 m ZY-3 imagery using a U-Net+CBAM deep learning model. They identify more than 13,900 glaciers covering 21,261.8 km², including 2,239.9 km² of supraglacial debris, with an overall mapping uncertainty of ±4.7%. The paper is well written, and a meter-scale glacier inventory for the Karakoram will be useful to the glaciological community. However, I have several concerns below.
Major comments
- At L305-307 the authors state that approximately 21.5% of the glacier-outline coverage was derived from Sentinel-2 imagery in areas without ZY-3 coverage, with the remainder mapped from ZY-3. Please provide a map showing the spatial coverage of each source. This would allow the reader to assess how outline quality varies between the two, and I would encourage the authors to report glacier number, size distribution, and debris fraction separately for the ZY-3- and Sentinel-2-mapped areas.
- At L445-448 the authors argue that RGI v7.0 and KGI-2020s overgeneralise debris-covered termini. However, 2mKGI excludes stagnant debris and non-glacial debris-covered surfaces, whereas RGI and KGI-2020s follow the convention of including debris-covered ice that remains connected to the glacier. Under these two criteria, 2mKGI will be systematically smaller than the other inventories, so the difference at debris-covered termini may reflect the mapping criterion rather than an error in the earlier products. Please state explicitly how debris-covered ice was identified and removed in 2mKGI, and what evidence was used to decide that a given debris surface is no longer ice-covered
- The comparison between inventories spans a wide range of data acquisition periods. The authors identify 192 surge-type glaciers in 2mKGI, and rapid terminus advance during surge events would produce area differences of similar magnitude to the "credible differences" described at L474. Comparisons with other inventories should therefore exclude surge-type glaciers, or at least report the results with and without them.
- In Section 5.2, it is not surprising that the 4 m ZY-3 DEM resolves higher roughness and more short-wavelength variability than the 30 m ASTER DEM. This is largely caused by the grid spacing and is not by itself evidence of a data-quality difference. I recommend aggregating the ZY-3 DEM to 30 m and repeating the comparison at matched resolution. Any residual difference can then be attributed to ASTER error.
- L215-217. The text implies that the training, validation, and test subsets were split after data augmentation, which would place augmented copies of the same source tile on both sides of the split and lead to data leakage and an optimistic evaluation. Please clarify the order of these steps and describe the split protocol in detail. I would also encourage a spatially blocked split rather than a random one, since adjacent tiles from the same scene are strongly correlated. In addition, please confirm that the 4,595 km² test set is independent of the four scenes used for label generation, and clarify its relationship to the 10% test subset described here.
Specific comments
- L120-125. Please give the exact number of Sentinel-2 images used rather than "approximately 10", and provide a table listing all ZY-3, Sentinel-2, and Landsat scene IDs with their acquisition dates.
- L196-199. Please state the ATL06 date range used for validation, confirm whether the DEM was co-registered before the assessment, and define the slope thresholds separating the flat, hilly, mountainous, and alpine terrain classes.
- L397. Please give the source for the known surge-type glaciers.
- Table 1. Please give the resolution of KGI-2020s and RGI v7.0. The date of "2000" for RGI v7.0 is also inaccurate, as the source imagery for this region spans several years.
- Figure 1. Please give the source of the Karakoram boundary polygon.
- Figure 5. The red ridgelines derived from the ZY-3 DEM are difficult to distinguish from the background DEM. I recommend a different colour.
- Please release the training and inference code in a public repository and include the code link in the data availability section, together with the training labels, so that the study can be reproduced.
Citation: https://doi.org/10.5194/egusphere-2026-2836-RC3 -
AC3: 'Response to Reviewer 3', Xin Yang, 06 Sep 2026
We sincerely thank the reviewer 3 for the careful assessment of our manuscript and for the constructive comments. We have carefully considered all of the reviewer’s suggestions and revised both the manuscript and the 2mKGI dataset accordingly. Our detailed point-by-point responses and the corresponding revisions are provided in the attached response document.
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 143 | 44 | 25 | 212 | 27 | 15 | 15 |
- HTML: 143
- PDF: 44
- XML: 25
- Total: 212
- Supplement: 27
- BibTeX: 15
- EndNote: 15
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
The authors present a comprehensive and high-resolution glacier inventory (2mKGI) for the Karakoram range, integrating multi-source satellite imagery and a deep learning approach (U-Net+CBAM). The methodology is innovative, and the resulting dataset holds significant value for cryospheric research. The manuscript is well-structured; however, there are several areas where additional methodological clarification, transparency regarding data processing, and adherence to established glaciological standards would significantly improve the clarity, reproducibility, and utility of the dataset. I have compiled the following specific comments to assist in strengthening the manuscript for publication.
1. Figure 3 provides a clear overview of the U-Net+CBAM workflow. However, the specific activation functions (e.g., ReLU, Sigmoid) used in the convolutional blocks are not explicitly stated. Including these details in the figure caption or the main text would help the reader better understand the model's non-linear fitting capabilities.
2. As the glacier boundary extraction relies heavily on deep learning, I strongly recommend providing the core training and inference code (or pseudocode) via a public repository such as GitHub. This transparency is crucial for the reproducibility of the methodology.
3. The manuscript mentions that Land Surface Temperature (LST) was retrieved from Landsat 8 imagery. Given that the ZY-3 imagery was acquired between July and October 2020–2021, could you clarify whether the LST represents a single scene matching the specific ZY-3 acquisition date, or if it is a composite/average product for the ablation season (July–October)?
4. In the manual refinement phase, the manuscript mentions using features such as lateral moraines and meltwater outlets as geomorphological indicators. However, the distinction between debris-covered ice and lateral moraines is not explicitly detailed. Please clarify the specific visual texture features or criteria used to manually differentiate debris-covered glaciers from lateral moraines.
5. The legend for Figure 6 includes ‘Karakorum’, but there are no corresponding boundaries, nor are any required here. Could this be an error in the legend?
6. The manuscript provides a total glacier area of 21,261.8 km². According to standard glacier inventory protocols (e.g., RGI standards), internal bedrock outcrops (Nunataks) should typically be excluded. Please confirm whether the 2mKGI dataset excludes these Nunataks. If they are retained, what is their approximate proportion of the total area?
7. Regarding Section 4.4, which identifies 192 surge-type glaciers, please clarify the data source for these "known" surge-type glaciers. Were they referenced from existing literature, or were they identified independently based on surface morphology observed in the ZY-3 imagery?
8. The study notes that 21.5% of the glacier boundaries were supplemented by 10m-resolution Sentinel-2 imagery where ZY-3 coverage was unavailable. Could you provide a breakdown of whether the overall extraction uncertainty of ±4.7% differs significantly between the ZY-3 (2m) and Sentinel-2 (10m) data sources?
9. The 2mKGI dataset provides high-resolution benchmarks for approximately 13,900 glaciers. As this paper is intended to facilitate data sharing, I recommend stating in the text or supplementary material whether the attribute table of the final vector Shapefile includes standard GLIMS or RGI fields (e.g., standard glacier ID format, specific acquisition dates, area). This would greatly enhance the usability and interoperability of the dataset for the global glaciological community.
10. Glacier unit subdivision uses ZY-3 4 m DEM where available and ASTER GDEM V3 elsewhere. The two DEMs differ in resolution, vertical accuracy, and ridge-line definition, which will introduce systematic inconsistencies in glacier unit boundaries across the mosaic edge. Please assess ridge-line positional differences between the two DEMs in overlapping areas, and report how many glacier units straddle the DEM mosaic boundary. Also discuss whether the use of ASTER GDEM introduces detectable biases in glacier count, elevation statistics, or slope metrics for the affected subregions.
11. The inventory includes glaciers down to 0.005 km² (~1250 pixels at 2 m resolution). No justification is given for this low threshold, and there is no discussion of whether very small features may be perennial snow patches.