Multitemporal glacier inventory revealing four decades of glacier changes in the Ladakh region

Multi-temporal inventories of glacierised regions provide an improved understanding of water resource availability. In this study, we present a Landsat-based multi-temporal inventory of glaciers in four Upper Indus sub-basins and three internal drainage basins in the Ladakh region for the years 1977, 1994, 2009 and 2019. The study records data on 2257 glaciers (of individual size >0.5 km²) covering an area of $\sim \mathrm{7923}\pm \mathrm{106}$ km² which is equivalent to ∼30 % of the total glacier population and ∼89 % of the total glacierised area of the region. Glacier area ranged between 0.5±0.02 and 862±16 km², while glacier length ranged between 0.4±0.02 and 73±0.54 km. Shayok Basin has the largest glacierised area and glacier population, while Tsokar has the least. Results show that the highest concentration of glaciers is found in the higher elevation zones, between 5000 and 6000 m a.s.l., with most of the glaciers facing towards the NW–NE quadrant. The error assessment shows that the uncertainty, based on the buffer-based approach, ranges between 2.6 % and 5.1 % for glacier area, and 1.5 % and 2.6 % for glacier length with a mean uncertainty of 3.2 % and 1.8 %, respectively. This multitemporal inventory is in good agreement with previous studies undertaken in parts of the Ladakh region. The new glacier database for the Ladakh region will be valuable for policy-making bodies, and future glaciological and hydrological studies. The data can be viewed and downloaded from PANGAEA, https://doi.org/10.1594/PANGAEA.940994 (Soheb et al., 2022).

1 Introduction

The Himalaya is the largest storehouse of snow and ice outside the Polar Regions. This large reserve of water plays a crucial role in the hydro-economy of the region (Bolch, 2019; Frey et al., 2014; Maurer et al., 2019; Pritchard, 2019). Any change to the Himalayan cryosphere would have a direct impact on the hydrology, further influencing the communities downstream whose livelihood and economy relies on, and are supported by, the major river systems e.g. the Brahmaputra, Ganges and Indus, among others. In high altitude arid regions like Ladakh, where the majority of glaciers are small and restricted to higher altitudes, meltwater serves as an important driver of the economy, especially in years with low winter precipitation when glacier melt becomes the major (or only) source of water (Schmidt and Nüsser, 2012, 2017). Recent studies have reported that Himalayan glaciers are retreating at an alarming rate (Azam et al., 2021; Bolch, 2019; Kääb et al., 2015; Maurer et al., 2019; Pritchard, 2019; Shean et al., 2020, among others) with glaciers of the Western Himalayas showing less shrinkage than the glaciers of the central and eastern parts (Azam et al., 2021; Shukla et al., 2020; Singh et al., 2016). Glaciers in the nearby Karakoram region display long-term irregular behaviour with frequent glacier advances/surges and minimal shrinkage, which is yet to be fully understood (Azam et al., 2021; Bhambri et al., 2013; Bolch et al., 2012; Kulkarni, 2010; Liu et al., 2006; Minora et al., 2013; Negi et al., 2021). Glaciers of the Karakoram region experienced an increase in area post-2000, due to surge-type glaciers. In just the upper Shayok Valley, as many as 18 glaciers, occupying more than one-third of the glacierised area, showed surge-type behaviour (Bhambri et al., 2011, 2013; Negi et al., 2021). However, not all regions of Ladakh have been analysed at the same level of spatio-temporal detail. In particular, our knowledge of glacier dynamics and their response to climate change is still incomplete in the cold–arid, high-altitude Ladakh region (∼105 476 km²) comprising both the Himalayan and Karakoram ranges. Few studies have focused on the glaciers of this region (e.g. Bhambri et al., 2011, 2013; Chudley et al., 2017; Negi et al., 2021; Nüsser et al., 2012; Schmidt and Nüsser, 2012, 2017; Shukla et al., 2020).

The advent of remote sensing technologies has permitted the mapping and measuring of various glacier attributes even in the absence of sufficient in situ observations (Bhardwaj et al., 2015). Glacierised area estimations have often relied on global and regional glacier inventories such as the Randolph Glacier Inventory (RGI), Global Land Ice Measurements from Space (GLIMS), Geological Survey of India (GSI) inventory and Space Application Centre India (SAC) inventory, among others (Chinese Glacier Inventory (CGI), Glacier Area Mapping for Discharge from the Asian Mountains (GAMDAM), International Centre for Integrated Mountain Development (ICIMOD)). However, given the large scale of these inventories, automated techniques are employed, in most of the cases, to map and calculate glacier extent with differing levels of success. Additionally, the varying quality of satellite imagery acquired from different time periods are sometimes necessitated in high mountain areas, such as Ladakh. Together, these two factors can lead to over- or under-estimation of glacier areas leading to erroneous information on temporal change. Moreover, there is no multi-temporal glacier inventory available for the entire Ladakh region, which can inform us on the changes in the natural frozen water reserves which have put the water security of this entire cold–arid region under significant stress during recent years. The residents of Ladakh have witnessed a decrease in agricultural yields, the main driver of economic development of the region, due to a decrease in water resources (Barrett and Bosak, 2018). The water scarcity together with an increase in tourism footprint (four times more tourists (327 366) in 2018 than 2010, a number that is more than the entire population of Ladakh) has led to a shift in livelihood from agriculture to other commercial activities (Müller et al., 2020), though even the latter relies heavily on water resources. In order to cope with water scarcity, some people of Ladakh have developed new water management techniques, commonly known as “ice reservoirs” or “ice stupas”, to supplement agricultural activities (Nüsser et al., 2019a, b).

This study presents a new multi-temporal glacier inventory for the Union Territory of Ladakh, India, covering 42 years of change between 1977 and 2019. This new dataset and analyses of glacier distribution will help to improve understanding of the glacier dynamics and the impact of ongoing climate change on water resources in the Ladakh region, where glaciers are the only source of water in the dry season. The inventories are entirely based on Landsat images acquired mostly during late summer with additional quality control provided through high-resolution PlanetScope and Google Earth imagery. We further establish a comparison with the existing inventories and data available in recent studies from the region. The dataset produced in this study can be viewed and downloaded from PANGAEA, https://doi.org/10.1594/PANGAEA.940994 (Soheb et al., 2022).

2 Study area

This study focuses on glaciers in the Upper Indus Basin (UIB) upstream of Skardu and three internal drainage/endorheic basins (IDBs) within Ladakh, namely Tsokar, Tsomoriri and Pangong basins. The geographic extent of the study area lies within a latitude of 31.1 to 35.6^∘ N and a longitude of 75.1 to 81.8^∘ E and covers a vast region of the Karakoram and Western Himalayan ranges. UIB has an area of ∼105 476 km², of which ∼8302 km² (8 %) is glacierised by ∼6300 glaciers spanning elevations between ∼3400 and ∼7500 m a.s.l. (as per RGI 6.0). IDBs of Tsokar (1036 km²), Tsomoriri (5462 km²) and Pangong (21 206 km²) house ∼30, 345 and 812 glaciers, comprising a glacierised area of ∼7 (0.6 %), 185 (3.4 %) and 437 (2.1 %) km², respectively (as per RGI 6.0). The glaciers of IDBs are at a comparatively higher elevation, spanning from ∼4800 to ∼6800 m a.s.l. Meltwater from these glaciers drains into the lakes within each basin. Pangong Lake (a saline lake), situated at an elevation of ∼4241 m a.s.l., is the largest with an area of ∼703 km². Both Tsomoriri (freshwater lake at ∼4522 m a.s.l.) and Tsokar (saline lake at ∼4531 m a.s.l.) lakes are designated Ramsar sites which occupy areas of ∼140 and ∼15 km², respectively. Since the majority of the investigated area (UIB and IDBs combined) falls within Ladakh, the combined area of UIB and IDBs will be referred to as “Ladakh region” hereafter.

Figure 1 Location map of the study area: the boundaries of studied Upper Indus Basin and internal drainage basins are outlined in black and red on the digital elevation model (DEM) and in the inset map. Inset map shows the study area with respect to the Himalayan and Karakoram region. Black dots and stars represent the respective basins' major settlements and field-investigated glaciers. The background image (ASTER GDEM) courtesy: NASA/METI/AIST/Japan Space Systems and U.S./Japan ASTER Science Team, https://earthdata.nasa.gov/ (last access: 17 June 2020).

The Ladakh region has a cold–arid climate due to the rain shadow and elevation effects of the Himalaya and Karakoram mountains (Schmidt and Nüsser, 2017). Mean annual air temperature and annual precipitation range between 0 to 10 ^∘C and 20 to 145 mm, respectively (Hersbach et al., 2020; Fig. 2). This region is inhabited by ∼700 000 people (as per Census of India, 2011; Census of China, 2020), most of which are directly, or indirectly, dependent on snow and glacier meltwater to support hydropower generation, irrigation and domestic needs.

Figure 2 Mean annual (a, c) and monthly (b, d) temperature and precipitation at Shiquanhe and Leh stations.

3 Data and methods

3.1 Data

This study utilises multiple Landsat level-1 precision and terrain (L1TP) corrected scenes (63 scenes in total) from four different periods: 1977±5 (hereafter 1977), 1994±1 (hereafter 1994), 2009 and 2019±1 (hereafter 2019). Scenes from the 1970s are majorly (12 out of 17) from the year 1976 and 1977, however due to higher cloud cover and less availability of imagery during the earlier Landsat period, five scenes from 1972, 1979 and 1980 were also included to aid the digitisation of glaciers (Table S1 in the Supplement). Images from late in the ablation season (July–October), having least snow and cloud cover (<30 % overall, and not over the glacierised parts), were selected and used for glacier identification and boundary delineation. Advanced Spaceborne Thermal Emission and Reflection Radiometer global digital elevation model (ASTER GDEM) scenes were also used for basin delineation and calculating slope, aspect and elevation metrics of the glaciers. Glacier digitisation, basin delineation and calculation of area were all performed in ArcGIS 10.4. Details of the imagery used in this study are presented in (Tables 1 and S1).

Table 1 Information on the satellite imagery used in this study (detailed information in Table S1).

3.2 Basin delineation

Basin delineation was carried using ASTER GDEM V003 and the hydrology tool in ArcGIS. The input DEM was first analysed to fill in all sinks with careful consideration of the potential for basin area over-estimation (Khan et al., 2014). UIB was delineated using a pour point selected at the Indus River in Skardu as we aimed to assess all the tributary basins of the Ladakh region. The UIB obtained by this approach was further divided into second-order tributary basins, i.e. Shayok, Suru, Zanskar and Leh basins. A small portion of the leftover area from UIB after second-order tributary basin delineation was merged into the Leh Basin in order to investigate the UIB upstream of Skardu. Delineation of the three endorheic basins (IDBs) that lie partially or completely in the Ladakh region, i.e. Tsokar, Tsomoriri and Pangong basins, was also carried out using the same method with the help of respective lakes as a pour point. The digitisation of the three lakes (Tsokar, Tsomoriri and Pangong lakes) was carried out manually for the years 1977, 1994, 2009 and 2019 using Landsat imagery.

3.3 Glacier mapping

Glaciers were mapped using a two-way approach, closely following the Global Land Ice Measurements from Space (GLIMS) guidelines (Paul et al., 2009): (1) automatic mapping of the clean glacier and (2) manually correcting the glacier outlines and digitisation of debris cover. First, a band ratio approach between NIR (near infrared) and SWIR (shortwave infrared) (as suggested by Paul et al., 2002, 2015; Racoviteanu et al., 2009; Bhardwaj et al., 2015; Schmidt and Nüsser, 2017; Smith et al., 2015; Winsvold et al., 2014, 2016) with a threshold of 2.0 (NIR/SWIR $>\mathrm{2}=$ ice/snow) was used on 2019 Landsat OLI images to delineate the clean part of glaciers. A median filter of kernel size 3×3 was applied to remove the isolated and small pixels outside the glacier area. The NIR and SWIR band ratio approach is good at distinguishing glacier pixels from water features with similar spectral reflectance values (Racoviteanu et al., 2009; Zhang et al., 2019). This approach failed in areas with high snow/cloud cover, shadows, frozen channels/lakes and debris cover. The snow/cloud cover and frozen lakes/stream problems were addressed by selecting Landsat scenes from the ablation period (July–October) with the cloud cover <30 %. The issue with the snow-covered regions in accumulation zones, where the delineation was the most challenging, was resolved using the best available imagery of any time between 1977 and 2019 because glaciers are not expected to change their shape significantly in the higher accumulation zones. One of the major issues was the debris-covered glaciers, which had to be manually digitised, with the support of high-resolution Google Earth and PlanetScope imagery from 2019±2. The result was then used as a basis for manual digitisation of debris-covered glaciers in other years where high-resolution images are not available. In most cases, identification of the glacier terminus was made with certain contextual characteristics at the snout, e.g. the emergence of meltwater streams, proglacial lakes, ice walls, end moraines (Fig. S1 in the Supplement).

The glacier outlines from 2019 were used as a starting point for the subsequent digitisation of glacier areas in 2009, 1994 and 1977. Glacier length was measured using a semi-automatic approach, by employing the DEM to identify a central flow line for each mapped glacier (Ji et al., 2017; Le Bris and Paul, 2013). Further manual corrections were undertaken to account for the flow lines of glaciers that have multiple tributaries and multiple highest/lowest points. Furthermore, some mapping errors are still expected to be present in this inventory due to a possible misinterpretation of glacier features, and the quantification of such errors are difficult owing to the lack of reliable reference in situ data in the Ladakh region. Such errors were minimised by keeping a fixed map–scale of 1:10 000 in most cases, and undertaking a quality check on glacier outlines using high-resolution images. In case of MSS images and smaller glaciers, a map–scale of 1:25 000 was also used whenever required.

Other specific glacier attributes were also extracted including new glacier IDs, Global Land Ice Measurements from Space (GLIMS)-0IDs, Randolph Glacier Inventory (RGI 6.0)–IDs, coordinates (latitude and longitude), elevation (maximum, mean and minimum), aspect (mean), slope (mean), area, length (maximum), area uncertainty and length uncertainty.

3.4 Uncertainty

This study involves the use of satellite imagery to extract various glacier parameters. It is therefore subject to uncertainties which may arise mainly from four different sources: (1) the quality of the image (with potential issues due to seasonal snow, shadows and cloud cover), (2) sensor characteristics (spatial/spectral resolution), (3) interpretation of glacial features and methodology used and (4) post-processing techniques (Le Bris and Paul, 2013; Paul et al., 2013, 2017; Racoviteanu et al., 2009, 2019). Error due to sources 1, 3 and 4 are generally minor and can be visually identified and corrected (Sect. 3.3), but an exact quantification is difficult due to the lack of reference data available from the region (Racoviteanu et al., 2009; Shukla et al., 2020). Type 4 errors are significant and have an impact on both glacier area and length estimation. Therefore, we applied a buffer-based assessment to glacier areas with the buffer width set to one pixel for debris-covered and a half pixel for clean ice (Bolch et al., 2010; Granshaw and Fountain, 2006; Mölg et al., 2018; Paul et al., 2017; Racoviteanu et al., 2009; Shukla et al., 2020; Tielidze and Wheate, 2018), given that the level 1TP Landsat images were corrected to sub-pixel geometric accuracy (Bhambri et al., 2013). A buffer-based method provides the maximum and minimum estimates of uncertainty with respect to glacier size, where the values vary with size of the glacier and spatial resolution of the imagery used. Thus, it is more specific to the dataset and most recommended when there are no reliable reference data available (Paul et al., 2017; Racoviteanu et al., 2009; Shukla et al., 2020). The same approach was also followed to estimate the uncertainties in lake areas with one pixel as the buffer width.

The associated uncertainty for smaller glaciers (< 0.5 km²) amounts to ∼12 %–25 %. Therefore, all the glaciers with an area of less than 0.5 km², which comprise ∼70 % and ∼10 % of the total glacier count and glacierised area, respectively, are not included in this study. For the remaining glaciers, the uncertainty in glacier area ranged between ±2.1 % and ±7.2 % depending on the spatial resolution of the satellite imagery and the individual glacier size. The highest uncertainty was for the year 1977 due to the coarser spatial resolution of Landsat MSS data when applied to the smallest glaciers (0.5–1 km²). For most of the glaciers, lengths are assumed to be accurate to ±1 pixel at the terminus (Le Bris and Paul, 2013). Therefore, a buffer of one pixel was set to determine the uncertainty in glacier length. The length uncertainty ranged between ±1.5 % and ±2.6 % with maximum uncertainty observed for the smallest glacier category (0.5–1 km²). The methods yielded an overall uncertainty of 4.2 %, 1.8 % and 1.5 % for glacier area, glacier length and lake area, respectively (Table S2).

Uncertainties related to other attributes (mean elevation, mean slope and mean aspect) of the inventory are difficult to estimate due to the use of the ASTER GDEM product in this study, which was developed using a collage of archived scenes acquired between 2000 and 2013. In addition, the local undulations and surface change over time will have only marginal effects on parameters (elevation, slope and aspect) that are averaged over the entire glacier as averaging compensates for most of the changes (Frey and Paul, 2012). However, for parameters like maximum and minimum elevations, where one cell is used and no averaging is applied, the uncertainty is $\sim \pm \mathrm{9}$ m, as the vertical accuracy of ASTER GDEM is ±8.55 m for glacierised areas of high Asia (Yao et al., 2020) and ±8.86 m elsewhere (Mukherjee et al., 2013).

4 Results

4.1 General statistics

In total, 2257 glaciers (>0.5 km²) were compiled in the current inventory (Table 2), with a total glacierised area of $\sim \mathrm{8511}\pm \mathrm{430}$, 8173±215, 8096±214 and 7923±106 km² for the years 1977, 1994, 2009 and 2019, respectively. The glacierised area corresponds to ∼6 % of the Ladakh region with individual areas ranging between 0.5±0.02 and 862±16 km². Glacier length in the Ladakh region varies between 0.4±0.02 and 73±0.54 km with a mean length of 2.9±0.05 km. About 90 % of the glaciers are shorter than 5 km in length while only 6 % of glaciers have a length of <1 km. Larger glaciers are mainly located in the Shayok and Zanskar basins with the Siachen Glacier being the largest (862±16 km²), longest (73±0.54 km) and covering the greatest elevation range of ∼3616 m (3702–7318 m a.s.l.). The major lakes in each endorheic basins of Pangong, Tsokar and Tsomoriri occupy an area of 3 %, 2 % and 2.5 %, respectively. The lake areas for the year 1977, 1994, 2009 and 2019 were 610±14, 619±8, 669±8 and 705±8 km² for Pangong, 13.5±0.9, 17±0.7, 18.3±0.7 and 18.8±0.6 km² for Tsokar, and 140±2.6, 141±1.3, 141±1.3 and 141±1.1 km², respectively.

4.2 Glacier distribution in the Ladakh region

Glacierised areas and population in the Ladakh region vary across basins. Shayok Basin has the largest distribution of glacierised area and population (74 % and 56 %), whereas the Tsokar Basin has the least (0.04 % and 0.1 %), respectively (Table 2). Based on size distribution, the glacier area category of 1–5 km² comprises the highest area (28 % of the total), while the category of 50–100 km² occupies the least glacierised area (9 %) of the region. Most glaciers (∼90 % of the total) in the Ladakh region have an area of <5 km² but occupy only 37 % of the total glacierised area. The population and area of glaciers in each area class are different in each basin but the proportion of glaciers, smaller than 5 km², is greater than 87 % in all basins. Glaciers larger than 100 km² (n=7, <1 % of the total) are only present in the Shayok Basin and occupy ∼24 % and 32 % of the total glacierised area of Ladakh and Shayok Basin, respectively.

4.3 Glacier hypsometry, slope and aspect

Figure 3iii and iv show the glacier elevations and hypsometry with 100 m elevation intervals of seven basins of the Ladakh region. The highest and lowest glacier elevation are 7740 and 3249 m a.s.l., both in the Shayok Basin. Whereas mean elevation of the glacier ranges between 4345–6355 m a.s.l. (Fig. 3iii). Small glaciers mainly occupy the higher elevations above 5500, and vice versa. The majority (73 %, 5810 km²) of the glacierised area is distributed in the 5000–6000 m a.s.l. elevation range, while only 14 % is located below 5000 m and 13 % above 6000 m a.s.l. (Fig. 3iv). The mean slope of these glaciers ranges between 8 and 46^∘, and is found to decrease with increasing glacier area. Glaciers with an area greater than 100 km² (n=7, <1 % of the total) have the lowest mean slope of 13^∘ wherea higher mean slopes (23^∘) are found for smaller glaciers (43 % of the total). Overall, the mean glacier slope is ∼21^∘ (Fig. 3v). Around 74 % (1665) of the glaciers face the northern quadrant (NW–NE) amounting to ∼50 % (3940 km²) of the glacierised area. While 9 %, 5 %, 3 %, 3 % and 4 % of the glaciers face East, South-East, South, South-West and West which constitute 24 %, 6 %, 8 %, 6 % and 6 % of the glacierised area, respectively. However, the orientation and respective area coverage of glaciers vary within individual basins (Fig. 3i and ii).

Table 2 Basin-wide glacier information of Ladakh region based on present study for the year 2019.

Figure 3 General statistics of the glaciers in the Ladakh region: orientation of glaciers (i) and associated area distribution (ii), maximum, minimum and mean elevation of glaciers (iii), hypsometry of glacierised area (iv) and slope against glacier area (v) and elevation (vi).

5 Discussion

5.1 The produced dataset and limitations

The multitemporal inventory of glaciers (>0.5 km²) in the Ladakh region for the years 1977, 1994, 2009 and 2019 is available at PANGAEA portal (https://doi.org/10.1594/PANGAEA.940994; Soheb et al., 2022). The dataset is provided in two different GIS-ready file formats, i.e. GeoPackages (*.gpkg) and Shapefiles (*.dbf, *.prj, *.sbn, *.sbx, *.shp, *.shx) to support wider end users. GeoPackage is a relatively new and open-source file format which is now being widely used and supported, whereas Shapefile format is one of the most widely used proprietary but open file formats for vector datasets, supported by open-source GIS tools such as QGIS. The outlines of glaciers, basins and lakes are all referenced to the WGS 84/UTM zone 43N datum. For each region, there is one file for basin outlines, and four files for glacier and lake (if present) outlines for 1977, 1994, 2009 and 2019. Each glacier outline file contains glacier IDs (New glacier IDs, Randolph Glacier Inventory 6.0 IDs, and Global Land Ice Measurements from Space initiative IDs), coordinates (latitude and longitude), elevation (maximum, mean and minimum), aspect (mean), slope (mean), area, length (maximum), area uncertainty and length uncertainty. Whereas the lake outline file contains coordinates, area, elevation and area uncertainty.

When using this dataset it is important to understand the key limitations of such regional-scale glacier inventories. Some of the key user limitations of the dataset are: (1) glaciers smaller than 0.5 km² (which comprise ∼70 % and ∼10 % of the total glacier population and glacierised area, respectively) were not included in this inventory due to the higher uncertainty (∼12 %–25 %) associated with these glacier outlines; (2) inventories produced in this study are entirely based on the medium-resolution Landsat imagery, in the same way as other global or regional-scale glacier inventories. Although the uncertainty associated with these inventories do not considerably impact regional-scale analyses, care should be taken while using these data for a small subset of glaciers. It should also be noted that it is not feasible to produce multitemporal inventories regionally using high-resolution datasets due to the paucity and high costs of such high-resolution datasets; (3) the inventories of 1977±5, 1994±1 and 2019±1 are produced using images with a range of acquisition dates due to the lack of data continuity within a particular year (more details in Sect. 3.1); and (4) the time periods chosen in this study are based on the availability of datasets and sufficient temporal gaps between the datasets to allow multitemporal glacio-hydrological analyses.

Table 3 Basin and class-wise comparison of the glacierised area between the present study and other inventories (RGI 6.0, ICIMOD and GAMDAM).

5.2 Significance of the present inventory

The glacier inventory presented here has several improvements compared to the existing regional and global inventories. Firstly, it covers the glaciers (>0.5 km²; n=2257; $\sim \mathrm{7923}\pm \mathrm{106}$ km²) for the entire Ladakh region with manual correction and quality control undertaken using freely available high-resolution images. The analyses were further extended to estimate the distribution of ice masses at the sub-basin scale. Secondly, the temporal aspect of the glacierised area will aid hydrological and glaciological modelling aimed at understanding past and future system evolution. Finally, the new inventory will aid both the scientific community studying the glaciers and water resources of the Ladakh region, and the administration of the Union Territory of Ladakh, Government of India in developing efficient mitigation and adaptation strategies by improving the projections of change on timescales relevant to policy makers.

5.3 Comparison of inventories in the Ladakh region

Differences in estimates of the glacierised areas are meaningful as they can lead to an over or underestimation of the available water resources. Therefore, correctly estimating glacier area over time is necessary for understanding glacier dynamics, future response to climate forcing and the water resources they provide. Table 3 presents a comparison between the present inventory and the Randolph Glacier Inventory (RGI) 6.0 (Pfeffer et al., 2014), the International Centre for Integrated Mountain Development (ICIMOD) inventory (Bajracharya et al., 2011, 2019; Williams, 2013) and the Glacier Area Mapping for Discharge in Asian Mountains (GAMDAM) inventory (Guo et al., 2015; Nuimura et al., 2015; Sakai, 2019) for the Ladakh region. The comparison involves glacier outlines for 2009 from the present study and excludes glaciers smaller than 0.5 km² from regional inventories to achieve the closest match temporally and for glacier size categories. This should be taken as a first-order comparison, given the fact that the uncertainties have been estimated with different approaches for the different inventories. Specifically, the uncertainty estimated for the GAMDAM and ICIMOD inventories differs only slightly to the one applied here, given that they used a normalised standard deviation approach on the datasets produced by several operators on the same subset of glaciers (Bajracharya et al., 2011, 2019; Guo et al., 2015; Nuimura et al., 2015; Sakai, 2019). Whereas, in case of RGI 6.0 inventory, the uncertainty estimation approach differs significantly from the one presented here, because their errors were calculated on a collection of glaciers due to the vast quantity of data acquired from multiple sources and approaches used to produced them (Pfeffer et al., 2014). Figure 4 presents a comparison of the only three inventories (present, RGI 6.0 and ICIMOD) for the five field-investigated glaciers of Ladakh region because RGI and GAMDAM inventories share the same outlines for these glaciers.

The comparison showed a higher glacierised area in the RGI/GAMDAM inventories and lower in the ICIMOD inventory (Table 3) than the present inventory, with most of the differences contributed by the basins having the higher glacierised areas (Shayok and Zanskar) and from the larger glaciers (>10 km²). Such inconsistencies among the inventories are a product of several factors, e.g. (1) absence of change in glaciers over time due to the use of imagery with a wide range of acquisition years (Fig. 4a, c, d); (2) misinterpretation of the glacier terminus due to icing, debris, snow and cloud cover (Nagai et al., 2016); and (3) the methodology used. The smaller difference between the present and the ICIMOD inventory is mainly due to the adoption of a similar technique (i.e. a semi-automated approach) and the shorter time frame of the analysis that generated the ICIMOD inventory (i.e. 2002–2009).

Figure 4 Comparison of inventories on the field-investigated glaciers of the Ladakh region: (a) Parkachik Glacier, Suru Basin; (b) Pensila Glacier, Suru Basin; (c) Lato Glacier, Leh Basin; (d) Khardung Glacier, Shayok Basin; (e) Stok Glacier, Leh Basin. The background image (Landsat) courtesy: U.S. Geological Survey/NASA, https://earthexplorer.usgs.gov/ (last access: 17 June 2020).

5.4 Comparison with recent studies

The data from the recent spatio-temporal change studies from different sub-regions of Ladakh (Fig. 5) are not in the public domain, except from Shukla et al. (2020). Hence, it is not possible to use these to validate our results. Therefore, our comparison mostly focuses on the rate of change for some of the individual glaciers (n=21, Fig. 5) from the literature and the bulk properties of a set of glaciers in different regions (Table S3, Fig. 6). Our results agree well with the studies conducted by others (Bhambri et al., 2013; Chudley et al., 2017; Garg et al., 2022a, b, 2021; Negi et al., 2021; Schmidt and Nüsser, 2012, 2017; Shukla et al., 2020) on individual glaciers of various sizes and on a set of glaciers, respectively (Fig. 6, Table S3). However, the results differ significantly only on some glaciers and especially in a part of the Shayok Basin (e.g. Kumdan (D), Aktash (E) and Thusa glaciers(I)). In the Shayok Basin, surge-type glaciers are common (Bhambri et al., 2013, 2017); the difference in analysis period between the present and other studies is the likely cause of the difference in glacier area statistics. Figure S2 presents the dynamics of the Kumdan and Aktash glaciers as an example of surge-type glacier of this region.

No significant difference was observed in rate of change of glacierised areas between the present study and other studies in the Leh, Tsomoriri, Zaskar and Suru basins. In contrast, the number of glaciers and glacierised area vary among these studies (present and others) but paint a similar picture of relatively lower retreat in the Shayok Basin (Bhambri et al., 2013; Negi et al., 2021), higher in Leh, Tsokar, Tsomoriri (Chudley et al., 2017; Schmidt and Nüsser, 2012, 2017), and moderate in Zanskar and Suru basins (Garg et al., 2022a, b; Shukla et al., 2020).

Figure 5 Presenting the spatial extent of different studies undertaken in Ladakh region. Black stars represent the individual glaciers. The background image (ASTER GDEM) courtesy: NASA/METI/AIST/Japan Space Systems and U.S./Japan ASTER Science Team, https://earthdata.nasa.gov/ (last access: 17 June 2020).

Figure 6 Comparison between the present study and other studies undertaken in different basins of Ladakh region over different time periods.

6 Data availability

The entire dataset of the Landsat-based multitemporal inventory of glaciers, larger than 0.5 km², in Ladakh region for the years 1977, 1994, 2009 and 2019 is available at PANGAEA, https://doi.org/10.1594/PANGAEA.940994 (Soheb et al., 2022).

7 Conclusions

We compiled a new glacier inventory of the Ladakh region for 1977, 1994, 2009 and 2019 based on 63 Landsat (MSS, TM and OLI) images, with least cloud/snow cover, acquired during the summertime (July–October). The inventory includes 2257 glaciers, larger than 0.5 km², covering an area of $\sim \mathrm{7923}\pm \mathrm{106}$ km² which is ∼14 % and ∼11 % less than the RGI 6.0 and the GAMDAM, and 7 % more than the ICIMOD inventories. The glacierised area accounts for ∼6 % of the Ladakh region with individual glacier areas ranging between 0.5±0.02 and 862±16 km². About 90 % of the glacier population are smaller than 5 km² but combined they occupy only 37 % of the glacierised area. The seven largest glaciers, larger than 100 km², account for ∼1879 km² or 24 % of the total. The Shayok Basin and glacier area category 1–5 km² hosts the highest number of glacier population and glacierised area; whereas, Tsokar Basin accounts for the least. More than 70 % of the glaciers are in the north-facing quadrant (NW–NE) and are concentrated in the higher elevation zones, between 5000 and 6000 m a.s.l. The error assessment shows that the uncertainty, based on a buffer-based approach, ranges between 2.6 % and 5.1 % for glacier area, and 1.5 % and 2.6 % for glacier length with a mean uncertainty of 3.2 % and 1.8 %, respectively. The uncertainty varies depending on the quality of the images and size of the glaciers. Our results also show a good agreement with other studies undertaken in parts of the Ladakh region for individual glaciers (n=21) and bulk properties of a set of glaciers.

The new multi-temporal inventory presented here will assist in planning the management of water resources, and for guiding scientific research focusing on glacier mass balance, hydrology and glacier change within the region. The detailed information and multi-temporal nature of this inventory will also aid in improving the existing global and regional glacier inventories especially in the cold–arid Ladakh region where the majority of the population is highly dependent on glacier-derived melt water resources for domestic, irrigation and hydropower generation needs.