Ensemble modelling of ice volume dynamics of Chhota Shigri Glacier in Himachal Pradesh from 2017 to 2022 Open Access

Among the 215,000 glaciers around the world, excluding the Greenland and Antarctic ice sheets (Farinotti et al. 2019), High Mountain Asia (HMA) consists of 97,965 glaciers (RGI Consortium 2017). The HMA hosts the most significant volume of water (7.0 ± 1.8 × 10³ km³ comprising 4.4%) available outside the polar regions, following Alaska (19.0 ± 4.9 × 10³ km³ comprising 12.0%) (Farinotti et al. 2019). Among the glaciers in HMA, the glaciers in the Himalaya–Karakoram region cover approximately 40,000 km² of the area (Bolch et al. 2012) and the Indian Himalayas constitute about 20,000 km² (Kulkarni & Karyakarte 2014; Nie et al. 2021). Hence, the Himalayan region is known as the third pole of the Earth owing to its extensive coverage of snow and glaciers (Gopika et al. 2021). The Himalayas also exhibit the largest repository of glacier systems globally (Bandyopadhyay et al. 2019). However, the Himalayas have become one of the most fragile mountain ranges owing to climate change in recent decades (Gopika et al. 2021). This glacier system is not only known for its socioeconomic importance such as freshwater supply, reservoir projects, or climate barriers but is also under the limelight for geohazards resulting from actively involved deglaciation (Thayyen et al. 2022; Rawat et al. 2024). For instance, Maurer et al. (2019) found that mass loss in the Himalayas for the years 2000–2016 has doubled to −0.43 ± 0.14 m w.e. year⁻¹ (metres of water equivalent per year) as compared with −0.22 ± 0.12 m w.e. year⁻¹ in the years 1975–2000. Such consequences due to the various glaciological processes demonstrate the ongoing climate-change scenarios over the Himalayas and are some of the most significant indicators of climate change (Kääb et al. 2012). Furthermore, the projected change in temperature under various CMIP5 scenarios is expected to be 2.6−4.6 °C by the end of the 21st century. This increase in temperature can significantly alter cryospheric processes, glacier area, and glacier hypsometry and have adverse socioeconomic impacts (Sabin et al. 2020). The Himalayan glacier system also contributes significantly to the major river systems of Asia such as the Ganges, Indus, and Brahmaputra, where often it is the primary source of freshwater for millions of people (Bandyopadhyay et al. 2019). The Himalayas accommodate approximately 54,000 glaciers covering an area of 60,000 km², of which 23,314 km² lie in India (Kulkarni & Karyakarte 2014). Because of its many glaciers, the Himalayan region is also known as the Water Tower of Asia (Bolch et al. 2012), as it stores approximately 7.0 ± 1.8 × 10³ km³ of water in solid form (Farinotti et al. 2009).

In recent decades, the Himalayan region has experienced a warming of temperature with a mean rate of increase of 0.21 ± 0.08 °C/decade (Gautam et al. 2010) resulting in glacier recession and shrinkage that directly affects the magnitude of glacial meltwater. These changes in glaciers lead to instabilities in the river systems (Bolch et al. 2012) and often trigger catastrophic events directly or indirectly, such as glacial lake outburst floods and landslides (Singh et al. 2022; Singh et al. 2023). Subsequently, a comprehensive investigation and monitoring of these glacier systems is needed for better understanding of the anomalies in Himalayan glacier dynamics as a consequence of climate change (Immerzeel et al. 2014).

Several studies have reported that a high rate of glacier thinning has been observed owing to climate change with a mean mass balance of −0.42 m w.e. year⁻¹ globally, where the glacier thickness in the Himalayan regions was observed to decline with a mean of −0.14 m w.e. year⁻¹ (Dehecq et al. 2019). The Karakoram glaciers of the Himalayas have shown significant anomaly, shifting from positive to negative mass balance (Kääb et al. 2012, 2015; Brun et al. 2017). Notably, the glaciers of HMA lost −16.3 ± 3.5 Gt year⁻¹ (gigatons per year) mass in 2000–2016 (Brun et al. 2017). These glacial changes may affect the entire glacier dynamics in terms of velocity, stored ice volume, and mass balance (Dehecq et al. 2019). A study by Zhou et al. (2018) showed that 8,814 Himalayan glaciers are experiencing slowdown in velocity owing to significant mass loss. A study by Dehecq et al. (2019) on surface glacier velocity of HMA shows a negative trend of −1.4 ± 0.4 m year⁻¹ decade⁻¹ or −9.8%± 2.9% decade⁻¹. These changes are associated with high ice-thinning rates in the last four decades (Zhang et al. 2022). Subsequently, glacier mass balance, thickness, and volume are critical for the prediction of water availability in future scenarios. A significant parameter in such investigations is the ice thickness distribution throughout the glacier, particularly for understanding the spatial patterns of ice thickness and mass changes for meltwater by which the storage capacity of the Himalayan glaciers could be estimated (Bandyopadhyay et al. 2019).

Multiple approaches have been employed to simulate the glacier thickness by incorporating mass-balance modelling and ice dynamics (Van Wyk de Vries et al. 2022). In practice, the modelling approaches provide an ice thickness for glaciers, except for small glaciers having a gentle topography, although with higher uncertainties (Linsbauer et al. 2012; Rabatel et al. 2018). Further, in situ observation based on intrusive/extrusive instrumentation such as hot water drilling, seismic or radar measurements, and gravimetry methods on glaciers in rugged terrain is challenging and often not feasible for entire glacier surfaces. Thus, extrapolated techniques are adopted to estimate the ice thickness by taking the traverse profile of glacier surfaces (Fischer 2009). The World Glacier Monitoring Service provides point ice thickness measurements of 2,000 glaciers, which renders calibration of the model parameters (Welty et al. 2020). Another impressive initiative is by the International Association of Cryospheric Sciences (IACS) in the Ice-thickness Models Intercomparison eXperiment, referred to as ITMIX by Farinotti et al. (2017), where they compared 17 different methods to estimate ice thickness based on basal shear stress, mass conservation, mass balance, ice flow velocity, artificial neural network, etc.

Studies by Bandyopadhyay et al. (2019), Singh et al. (2018), and Ramsankaran et al. (2018) utilize optical and microwave remote sensing to estimate ice thickness using techniques such as interferometric SAR (InSAR) for mass changes and thinning rate and some based on a model approach. Several ice thickness models have also been developed. A model based on the mass conservation of ice has been introduced by Farinotti et al. (2009), which uses digital elevation model (DEM) data, glacier boundary, and boundary of ice-flow catchments for estimating ice thickness. Some other models for ice thickness estimation include the slope-dependent ice thickness estimation (Haeberli & Hoelzle 1995), the glacier-flow-based ice thickness estimation (Gantayat et al. 2017), and the shear-stress-based model and its variants (Linsbauer et al. 2012).

In situ measurements of only 4,700 glaciers around the world are presently available (Millan et al. 2022). In the Himalayan ranges human accessibility is challenging owing to the harsh terrain and weather conditions, consequently in situ measurements are difficult to gather to produce a significant database of observations. The Himalayan glaciers are complex, and heterogeneity in their response pattern makes their assessment even more challenging (Fujita & Nuimura 2011; Kargel et al. 2011; Kääb et al. 2012; Gardelle et al. 2013; Harrison et al. 2021). Consequently, uniformity in the assessment of the glaciological properties of these glaciers also varies (Ji et al. 2022). Thus, it is essential to develop and maintain in situ datasets to provide a reference for calibration and continuous assessment of glacier dynamics models present today (Bolch et al. 2012). In this investigation, we combined two complementary methods, namely, the ice-flow-based methods and the basal-shear-stress-based methods, to estimate the glacier ice thickness and volume. We assess the dynamics of the Chhota Shigri Glacier in the northwest Indian Himalayas corresponding to the variations in velocity, ice thickness, and ice volume from 2017 to 2022. The approach follows the estimation of ice velocity based on the feature-tracking technique, and further modelling of ice thickness using an ensemble of different methods. The ensemble model of ice thickness consists of the five different methods that are derivatives of laminar-flow-based method and basal-shear-stress-based methods. Some of these methods in ice thickness modelling are known to overestimate and some to underestimate (Van Wyk de Vries et al. 2022). Subsequently, the spatial average from multiple methods ensures precise estimation of ice thickness. Combining the different methods also enables us to incorporate multiple parameters at different conditions yielding better estimates of ice thickness. A previous attempt at incorporating various methods to yield an overall ice thickness has been conducted previously on the Andes mountain range (Van Wyk de Vries et al. 2022). However, this study presents the first such implementation on a benchmark Himalayan glacier such as Chhota Shigri. The ice volume is determined as the product of ice thickness and glacier area derived from remotely sensed data and the volume–area scaling method, where ice thickness is derived from surface velocity and DEM.

In the present study, the boundaries of the Chhota Shigri Glacier were initially derived from the Randolph Glacier Inventory (RGI version 6.0) for clipping the remotely sensed data (Sentinel-2) and to determine the present extent of the glacier. For estimating the glacier velocities based on the feature-tracking technique, Sentinel-2 and Landsat true colour composites were used from 2017 to 2022. For additional comparisons with previous results, Landsat 4/5-Thematic Mapper (TM) images were also used for 2008–2009. A total of 214 multitemporal images have been used for surface glacier velocity at 50 m spatial resolution using glacier image velocimetry (GIV) (Van Wyk de Vries & Wickert 2021). The Advanced Land Observing Satellite–Digital Surface Model (ALOS DSM) at 30 m spatial resolution was used to measure ice thickness from 2017 to 2022. For the comparative experiments, the ALOS DSM is unavailable; hence, the Shuttle Radar Topography–Mission Digital Elevation Model (SRTM DEM) was used for the ice thickness modelling for 2008–2009. Further, MODIS (MOD11B3) Land Surface Temperature (LST) products were used to ascertain the corresponding variation of temperature with respect to the changes in glacier ice parameters. Table 1 summarizes the different datasets and tools used in this study. We used the codes developed by Van Wyk de Vries et al. (2022) and updated them per our requirements.

In addition, estimating the temperature for the study period was carried out using the MODIS-11 monthly data at 6 km resolution. The detailed usage of MODIS-11 data is given in Section 2.5.1 on ‘Temperature’. Along with the estimated parameters, i.e., surface velocity, slope, and temperature, the other parameters such as shape factor, basal sliding factor, and ice density have been used for the Monte Carlo Simulation model for the ice thickness ensemble model (Van Wyk de Vries et al. 2022). The primary methods use a shallow ice approximation (SIA) approach, which is a variant of the ice-flow-based method (Cuffey and Paterson 2010) and the basal-shear-stress-based method (Frey et al. 2014) for the estimation of ice thickness. Table 2 summarizes the methods used. A recent study published by Van Wyk de Vries et al. (2022) on the ice caps of the Andes Mountain ranges showed that using multiple methods that have a low correlation among them for the ice thickness calculation increases the accuracy of estimates using an ensemble model. We adopted the same approach in our investigations of Chhota Shigri Glacier. The estimated ice thickness was then used to calculate the glacier ice volume using the three methods illustrated in Figure 2.

The glacier boundary plays a significant role in ice volume estimation since changes in the area and spatial extent of glaciers are commonly observed owing to the rapid retreat in the Himalayas (Pfeffer et al. 2014). The RGI glacier boundaries are unable to represent the annual variation in these glacier boundaries and have high uncertainty leading to over- and underestimation of the spatial extents, which may not yield precise estimates of glacier ice volume (Van Wyk de Vries et al. 2022). Subsequently, precise mapping of glacier boundaries is important for ice masking and volume estimation. In our investigations, the RGI boundaries were observed to underestimate the snout position of the Chhota Shigri Glacier owing to the thick debris cover at the tongue (Zemp et al. 2023). Hence, we manually digitized the boundaries for 2017–2022 using Sentinel-2 satellite imagery having negligible snow cover. In addition, high-resolution Google Earth images have been utilized on an annual basis to map the temporarily snow-covered rocks and snout position, followed by manual preparation of ice masks for each year (2017–2022).

The glacier-surface velocity is a critical input in the SIA-based ice thickness estimation. In the present study, we used the GIV tool, which is based on the feature-tracking approach involving frequency domain multi-pass satellite image correlation for estimating the glacier-surface velocity (Van Wyk de Vries & Wickert 2021). The Sentinel-2 and Landsat 8 true colour composite images were used in the multi-pass image correlation for estimating the glacier-surface velocity. The GIV tool encompasses a temporal baseline, which is between nine and 800 days for generating image pairs, which is further segregated into short (9–90 days) and long (300–430 and 650–800 days) temporal baseline groups. We estimated the mean annual surface glacier velocity for a six-year period from 2017 to 2022 based on the short temporal baseline groups. The near anisotropic orientation filter (NAOF) was used for the pre-processing of the images to normalize the intensities in cloud-free, cloudy, and shadowed images (Van Wyk de Vries & Wickert 2021).

Displacement for each image pair is calculated using a multi-pass frequency domain cross-correlation algorithm. This technique gives the benefit of better feature-tracking combined with high resolution over a normalized image correlation algorithm (Van Wyk de Vries & Wickert 2021). The three-iteration technique reduces the window sizes with every lap in the multi-pass filter. The window sizes are reduced to 480, 240, and 120 m with a 50% overlap, and the displacement is finally calculated with a 75% overlap to ultimately generate a velocity map at 30 m resolution (Sveen 2004). The limiting values for the filters used in GIV (see Figure S1) to remove the erroneous pixel values were velocities exceeding the maximum threshold of 100 m year⁻¹, pixel values exceeding 50% of the surrounding pixels (four surrounding cells), and pixel values exceeding 200% of the surrounding pixels (ten surrounding cells). Also, the values with a signal-to-noise ratio less than 5, which provides less than 100% certainty for a strong signal and peak ratio more significant than 1.3, are excluded as well. The removed pixel values are then interpolated using the remaining valid cell values. Each velocity map is included in the time-series data generated by GIV. The mean velocity maps generated for the entire period of one year are sharpened through the Laplacian filter kernel, an edge-detection technique that is a 2D isotropic measure of the second spatial derivative. It particularly helps in finding the fine details in the image and correcting any rapid intensity change. Lastly, a velocity map is generated by averaging all individually processed velocity maps. The velocities out of 3σ bound were interpolated using a linear kernel and corrected with the surrounding pixels. Finally, the velocities for all six years were clipped to the modified RGI glacier boundary extent before estimating thickness and volume and are shown as Figure S1 in the Supplementary Material.

There are several methods to calculate glacier thickness, and each has its distinct advantages and limitations (Van Wyk de Vries et al. 2022). We adopted ensemble modelling to calculate ice thickness by taking a mean of widely used methods for glacial ice thickness approximation (Farinotti et al. 2017; Farinotti et al. 2021). In this analysis, three methods are based on ice flow velocity, and two are based on basal shear stress (Farinotti et al. 2017; Gantayat et al. 2017; Van Wyk de Vries et al. 2022).

In Equation (10a–10c), H_j,k is the ensemble mean value of each pixel in the ice thickness map, while r is the ground resolution of each pixel (assuming the square pixel is applicable to optical remote-sensing data), A_G represents the manually derived glacier area, represents the mean glacier ice thickness of the map, c is the coefficient of proportionality with a value of 0.2055, and γ is the exponential scaling parameter with a value of 1.36. The volume subscripts, namely, ‘p', ‘m', and ‘v–a’ represent the three approaches corresponding to the piece-wise, mean, and volume–area scaling, respectively. We used the values for these two scaling parameters, c and γ (Chen & Ohmura 1990). These values have been derived from the measurement of 63 European glaciers and give comparable results for the Himalayan glaciers (Frey et al. 2014).

A previous study by Azam et al. (2016) demonstrated that surface temperature is one of the sensitive parameters that influences the ice thickness and volume of a glacier, as also observed from Equations (4) and (8). Hence, a precise measure of temperature as input is imperative in the ice thickness retrieval model. For estimating the surface temperature, an inter-comparison between MOD-11 monthly LST data and MOD-11 eight-day LST data was carried out, to determine the better product. The comparison was based on remotely sensed temperature measurements with respective field measurements of air temperature for the years 2009–2010 to 2012–2013 given by Azam et al. (2016) for the meteorological conditions of the Chhota Shigri glacier air temperature from 2009–2010 to 2013–2014. The comparative study demonstrated that MOD-11 monthly LST data better approximate the measured air temperatures during that period. Consequently, we used the MODIS monthly LST products to estimate the temperature for our study period, 2017–2022, for the ice thickness estimation. The error in the temperature was calculated using the average deviation between the estimated LST and field-measured air temperature.

Errors in estimating glacier-surface velocity from the sub-pixel image correlation can emerge because of various factors such as erroneous image co-registration, false feature matching, errors in estimated sub-pixel displacement, and image distortions. Hence, for error analysis, velocity is estimated on the stable terrain (unchanged terrain) around the Chhota Shigri Glacier. Ideally, the surface velocity ds(x, y) for the stable terrain should be ∼0 m year⁻¹. Terrains that are deprived of snow for most years, landslides, and any stream flow are considered stable terrains, which comprise primarily bare rock surfaces in our study (Gantayat et al. 2014; Dehecq et al. 2019; Nela et al. 2022; Singh et al. 2023). Thus, this provides us with a means of analysing displacements s(x, y) at different points. We defined eight different regions of interest for stable ground, as shown in Figure S2. Any deviation in velocity from zero is considered an erroneous estimate. We calculated the mean and standard deviation of the displacements for these regions and used them to determine the uncertainty in ice velocity and in ice thickness estimation.

We evaluated the uncertainty of U_s by evaluating it in ice-free zones as discussed earlier. To define A, we vary temperature uniformly by ±0.45 K each year. Uncertainty in A varies from 19.2% in 2022 to 19.6% in 2019. We vary the lateral drag correction factor, f, uniformly between 0.7 and 0.9, maintaining the average at 0.8, i.e., an uncertainty of 12.5%, as discussed by Gopika et al. (2021). The basal sliding correction factor, β, is varied in this study between 0.1 and 0.4 (assume 25% overall basal sliding). Ice density, ρ, is assigned a constant value of 900 kg m⁻³, but the uncertainty in the matter is assumed to be 10% (Maanya et al. 2016). The error in slope comes from the DEM data and is assumed to be 10% of the mean gradient, which infers d sin α/sin α to be ±0.1. The uncertainty in the calculation of τ has been considered as 5%, i.e., ±7 kPa (Driedger & Kennard 1986; Haeberli & Hoelzle 1995; Ramsankaran et al. 2018).

The uncertainty in the glacier area for each year is calculated by estimating the percentage deviation between the area determined from the modified RGI glacier boundary used as the base polygon and the multiple-digitized glacier boundary from other sources, including Google Earth imagery and Sentinel-2 images at 10 m spatial resolution. Google Earth imagery is highly useful in identifying the snout position, rock surfaces, and other non-glaciated regions. Hence, the digitized glacier boundaries were derived from Sentinel-2 satellite imagery at 10 m spatial resolution (Ramsankaran et al. 2018) and high-resolution Google Earth imagery.

The results of glacier velocities indicate higher magnitudes in the accumulation and mid-ablation zones, as expected of a conventional mountain glacier (Benn & Evans 2010). The relatively lower velocities are observed near the snout and the edges, which follows the general concepts of glacier ice flow (Cuffey & Paterson 2010). Average velocities for each year are between 12.5 and 21.2 m year⁻¹. Subsequently, high values of ice thickness are observed in the accumulation area near the equilibrium line altitude (4,900–5,000 m a.s.l) (Chandrasekharan et al. 2018) and in the ablation zone (primary terminus) of the glacier (4,700–4,800 m a.s.l). The mean ice thickness estimates for all years are observed to be between 102 and 112 m. This section elaborates on the results primarily on the glacier ice velocity, thickness, and volume, from the proposed interventions from the methods discussed in the previous section. Further, insight on the uncertainty of the approach and the cross-sectional distribution of the ice thickness is also discussed here.

Another comparison between the ice thickness estimates from the present study was done with the ice depth data by Millan et al. (2022) for 2018 using statistical analysis. Millan et al. (2022) developed the ice thickness inventory for 98% of the 215,000 glaciers around the globe for the study period 2017–2018 using the laminar-flow-based method. The mean ice thickness for the Chhota Shigri Glacier in this global data inventory is approximately 98 m, and the mean depth in our study is approximately close for the period at about 104 m. It was observed that statistically within the 95% confidence interval the correlation between the estimates from this study and the ice depth from Millan et al. (2022) was 0.81. A study on the ice volume estimate of the Chhota Shigri Glacier based on the basal shear stress method revealed an ice volume of 1.74 ± 0.25 km³ as compared with the mean volume estimate of all the years in the present study of approximately 1.47 ± 0.26 km³ (Ramsankaran et al. 2018). This value is also comparable to the estimation of 1.4 km³ volume in 2015 given by Vashisht et al. (2017) using the volume–area scaling method, where thickness was estimated using the slope-dependent method. The glacier ice thickness and volume estimates given by Haq et al. (2021) for 2007 using artificial neural networks are 109 ± 17 m and 1.69 ± 0.26 km³, respectively, which are comparable to the estimations in the present study provided in Table S1. In comparison with these studies, the present approach is observed to slightly underestimate the glacier ice volume.

For the basal shear stress method, the same DEM was used for each year owing to the fact that these are unavailable on an annual time-scale. Consequently, the variation in the ice thickness and glacier volume are largely dependent upon the precise glacier boundaries. The changes in the glacier ice thickness can be compared with the mass-balance estimation, as given by Mandal et al. (2020) for the Chhota Shigri Glacier using remote-sensing techniques from 2017 to 2019. For 2017 and 2018, a negative trend in mass balance has been reported, where the mass loss is −0.28 and −0.4 m w.e. year⁻¹ for the years mentioned, respectively. For 2019, a positive trend in mass balance was reported, where the mass gain was +0.5 m w.e. year⁻¹. A consistency in mass-balance data and ice-flow-based estimates for ice thickness (Figures 4 and 5) was observed. The mean ice thickness based on the laminar flow model for the years 2017, 2018, and 2019 was 104.8, 104.4, and 112.1 m, respectively, showing a similar pattern to the mass-balance data of Mandal et al. (2020).

The velocities in the present study are marginally less than the estimates for the years before 2010 by Patel et al. (2019) and Tiwari et al. (2014), which also show a negative trend in the glacier velocities due to climatological impacts. A study by Sattar et al. (2019) shows that the velocity for the year 2016–2017 is 21 m year⁻¹, while in this study, the velocity for 2017 is comparable at 20.5 m year⁻¹.

Temperature is an important parameter that affects the ice thickness estimation in the laminar-flow-based method. We have noticed based on the MODIS LST products that although the temperature in 2018 was lower than in 2017, the mass balance was negative, and our results agree with the same. A possible explanation could be warming, which has been reported in some studies (Varade & Dikshit 2019; Varade & Dikshit 2020). Using MODIS data, we observed that the temperature difference between the accumulation and ablation zones for 2018 was higher than that observed in 2017. Daytime monthly temperature also indicated an early snow-melt in 2018. Also, during the peak melting season, the daytime temperature in 2018 was warmer than in 2017. These could be some of the possible reasons for the downtrend in the ice thickness in 2018. Also, resonating with the inferences by Van Wyk de Vries et al. (2022), we observed that the basin-divided approach gives an underestimated result using the basal-shear-stress-based method and the difference between ice thickness of individual basins can be seen around the edges of individual basins. This ice thickness pattern could be seen in the whole glacier and can be counted as a limitation of the basin-divided approach.

The ice thickness patterns, as shown in Figure 4, depict that the maximum accumulated ice is near the equilibrium line altitude at 5,050 m a.s.l. (Azam et al. 2016) and on the mid-ablation area, which is the only lobe in Chhota Shigri Glacier, as mentioned in the global inventory by the World Glacier Monitoring Services (WGMS) in 2021. From the lobe region, ice depth continuously decreases towards the snout. However, the disadvantage of not using a particular approach from the GlabTop2 method where the marginal cells are assigned an ice depth of zero is that an exaggerated value of ice depth can be encountered at the edges, especially at the tongue, owing to the lesser gradient. Further, averaging in the ice-flow-based method might not be effective for the edges, as the movement of surface particles in a non-glaciated region due to meltwater flux from the surrounding region (Nela et al. 2022) or technical errors, even for a minimal velocity, can derive significant ice depth, especially in an area with a lesser gradient. A terrain slope with even a decrement of 0.1° can have a staggering impact on ice thickness estimation. However, the value is not affected much by an increment of 0.1°, as shown in Figure S3 in the Supplementary Material. The ice thickness is observed to decrease by 1 m when the slope increases by 0.1° compared with the estimated thickness. Conversely, when the slope decreases by 0.1° for the whole glacier, the difference between the estimated thickness and the original estimation is observed to be 9 m.

Glaciers play a pivotal role in the Earth's climate system, and accurate ice thickness and volume estimates are crucial for understanding their behaviour and potential contributions to sea-level rise. This study combined two complementary methods, the ice-flow-based method and the basal-shear-stress-based method, to estimate glacier ice thickness and subsequently the glacier ice volume. The ice-flow-based method takes advantage of the precise ice velocity estimates obtained through state-of-the-art remote-sensing techniques to derive glacier ice thickness by utilizing the principles of ice flow dynamics. This method provides a spatially detailed view of ice thickness distribution across the glacier. In parallel, the basal-shear-stress-based methods focus on the interaction between a glacier's basal ice layer and its underlying substrate.

By analysing basal shear stress and incorporating information on glacier geometry, we refine our understanding of ice thickness. This approach offers insights into the glacier's mechanical behaviour at its base. Several studies in the past on the Chhota Shigri Glacier, one of the benchmark glaciers in India, have illustrated its significance from the perspective of the glaciohydrological processes in the region. In this study, we made an attempt to integrate ice velocity, basal shear stress, and glacier geometry to provide a comprehensive and accurate assessment of the Chhota Shigri Glacier ice thickness and volume.

The present study investigated by the ensemble approach, based on multiple methods for the estimation of Chhota Shigri Glacier ice thickness and subsequent derivation of the glacier ice volume by utilizing glacier area determined from manual digitization of the glacier boundaries from satellite images using the RGI reference. The surface velocity estimates based on the feature-tracking approach showed an alternate increase and decrease in the peak velocities between 2017 and 2022, where similar observations were made by Sahu & Gupta 2021. However, no significant changes have been observed in the ice thickness between 2017 and 2022 based on our study. The maximum variation in the peak ice thickness was observed to be around 14 m between 2017 and 2022. Similar observations can be made regarding the mean glacier volume, which varied by 0.14 km³ between 2017 and 2022. The VWDV and GT2B methods were, respectively, observed to significantly overestimate and underestimate the mean ice thickness and volume. The volume–area scaling approach for estimating glacier volume highly underestimated the results. Further, the results based on the GT2W method were observed to be the nearest to the mean among the various methods.

The objectives of this study were to investigate the spatial and volumetric dynamics of the Chhota Shigri Glacier in the state of Himachal Pradesh in India. The results reveal the potential of the ensemble approach for estimating the glacier ice parameters with remotely sensed data. While collecting in situ measurements on mountain glaciers is challenging owing to the rugged terrain and expensive logistics and field equipment, remote-sensing techniques provide a viable means to investigate glacier dynamics. The uncertainty analysis shows that substantial confidence can be placed in the estimates derived. It was observed that the quality of the retrieved glacier velocity, which is critical for the ice thickness estimation, depends upon the number of cloud-free satellite images available. In addition, high-resolution images with high spatiotemporal resolution are not easily available, which affects the accuracy of the glacier velocity estimates. One of the issues with ice thickness modelling is that there are methods that sometimes overestimate or underestimate ice thickness of the glaciers. For example, GlabTop2 underestimated the results in our study compared with other models that were based on laminar flow. This is one of the reasons that ensemble modelling is recommended for better results in investigations on ice thickness using remote sensing.

In the future, efforts can be directed at improving the glacier velocity estimates with fewer satellite images by incorporating velocity products or composite images from active microwave remote-sensing data. Further, incorporation of mass-balance-based approaches for the ensemble modelling of ice thickness estimation may also improve the estimates. We also recommend validating the thickness and volume estimations through different techniques with remote-sensing data like GRACE FO and geophysical data, which was not possible in this study due to the coarse resolution of the data and the corresponding smaller glacier area of the Chhota Shigri Glacier. The results of such studies may be used to study the implications of ice volume estimations on watershed management, climate change, policymaking, extreme events such as flash flooding, and so on.

The authors acknowledge the data support from the USGS. The authors also acknowledge the NFS seed grant SGT-100052 from IIT Jammu. The authors thank Dr Maximillian Van Wyk de Vries for providing the fundamental software codes for this work.

All relevant data are available from an online repository or repositories. The Sentinel-2 images are can be freely downloaded from the Copernicus dataspace https://dataspace.copernicus.eu/. The 30m resolution SRTM DEM is freely available through the Earth Explorer portal of the USGS at https://earthexplorer.usgs.gov/.

The authors declare there is no conflict.