Glacier thickness and volume estimation in the Upper Indus Basin using modeling and ground penetrating radar measurements

Abstract

In the Himalaya, ice thickness data are limited, and field measurements are even scarcer. In this study, we employed the GlabTop model to estimate ice reserves in the Jhelum (1.9 ± 0.6 km³) and Drass (2.9 ± 0.9 km³) sub-basins of the Upper Indus Basin. Glacier ice thickness in the Jhelum ranged up to 187 ± 56 m with a mean of ~24 ± 7 m, while the Drass showed ice thickness up to 202 ± 60 m, with a mean of ~17 ± 5 m. Model results were validated using Ground Penetrating Radar measurements across four profiles in the ablation zone of the Kolahoi glacier in the Jhelum and nine profiles across the Machoi glacier in the Drass sub-basin. Despite underestimating ice-thickness by ~10%, the GlabTop model effectively captured glacier ice-thickness and spatial patterns in most of the profile locations where GPR measurements were taken. The validation showed high correlation coefficient of 0.98 and 0.87, low relative bias of ~ −13% and ~ −3% and a high Nash–Sutcliffe coefficient of 0.94 and 0.93 for the Kolahoi and Machoi glaciers, respectively, demonstrating the model's effectiveness. These ice-thickness estimates improve our understanding of glacio-hydrological, and glacial hazard processes over the Upper Indus Basin.

1. Introduction

The Himalaya has the largest glacier area outside the poles (Bolch and others, 2012; Sakai, 2019) sustaining lives and livelihood of millions of people downstream (Immerzeel and others, 2010; Tuladhar and others, 2021). Glaciers in the Himalaya and elsewhere in the world exert a complex influence on land surface and climate processes (Milner and others, 2017; Johnson and Rupper, 2020) and are anticipated to affect the regional hydrological regimes under the projected climate change (Romshoo and others, 2020a; Chen and Yao, 2021). The assessment of land system changes, changes in the local, and regional climate and hydrological regimes and glacial hazards requires an accurate estimate of glacier volume and ice thickness distribution (Huss and Hock, 2018). However, despite far-reaching implications, glacier volume and ice thickness distribution estimates over the Himalayan region are limited largely due to technological limitations, remoteness, and challenging topography (Bolch and others, 2012) and the consequent limited field observations (Wagnon and others, 2013; Zhang and others, 2022). As a result, knowledge about the amount of water stored in these glaciers and their response to changing climate is limited. It is important to note that accelerated glacier melting and the consequent impacts on various dependent sectors have attracted the attention of researchers from all over the world to understand the response and behavior of the Himalayan glaciers (Cogley and others, 2010; Bhambri and others, 2011; Gardelle and others, 2013; Gardner and others, 2013; Kääb and others, 2015; Brun and others, 2017; Salerno and others, 2017; Maurer and others, 2019; Abdullah and others, 2020). However, most of these studies have investigated glacier retreat (Kamp and others, 2011; Pandey and others, 2011), mass balance (Ghosh and Pandey, 2013), glacier elevation changes (Abdullah and others, 2020; Romshoo and others, 2022a), climate change impacts (Rashid and others, 2017) and only a few have studied glacier ice thickness or volume (Linsbauer and others, 2009; McNabb and others, 2012; Frey and others, 2014; Gantayat and others, 2014; Linsbauer and others, 2016; Farinotti and others, 2017; Farinotti and others, 2019; Sattar and others, 2019; Pandit and Ramsankaran, 2020; Millan and others, 2022; Nela and others, 2023). It is pertinent to mention that direct ice thickness measurements over the Himalayan region are available for only about 15 glaciers (Mishra and others, 2022).

The glacier thickness and volume are the basic and most important parameters for projecting the future glacier evolution (Le Meur and others, 2007; Kaser and others, 2010; Gabbi and others, 2012; Immerzeel and Bierkens, 2012; Farinotti and others, 2019; Liang and Tian, 2022), future water availability (Huss and others, 2008), and estimation of future sea-level rise (Gabbi and others, 2012). Information about glacier thickness, besides being required for glacier, volume estimation, is also important for various glacio-hydrological studies (Huss and others, 2008), regional and local climate modeling (Kotlarski and others, 2010) and assessment of glacier hazards (Frey and others, 2014). Due to the paucity of sufficient information, particularly about glacier thickness and volume (Haq and others, 2021), it is difficult to assess glacier dynamics and future glacier projections under climate change in the Himalaya (Jacob and others, 2012).

The existing glacier volume estimates are largely based on a simple volume–area relation (Chen and Ohmura, 1990; Bahr and others, 1997) and slope-dependent glacier thickness estimation (Hoelzle and Haeberli, 1995). It is noteworthy that volume estimates based on volume–area scaling are inherently unstable, nonunique and sensitive to small changes in glacier boundaries (Bahr and others, 2012) and, the estimates quite often vary considerably between these methods. For example, Frey and others (2014) estimated the glacier volume over the Himalaya–Karakoram region to range from 2955 to 4737 km³ based on various approaches and the glacier mean thickness varies from 94–158 m in the Karakoram and 54–83 m in the Himalayan region (Frey and others, 2014). Similarly, Bolch and others (2012) reported glacier ice volume estimates over the Himalaya ranging from 2300 to 6500 km³ depending on the approach employed. On the other hand, Ohmura (2009) estimated glacier volume between 3800 and 4850 km³ over the Himalaya encompassing the parts of Pakistan, Nepal, Bhutan and India. Depending on the glacier inventory used, Cogley (2011) estimated glacier volume in the range of 3600–7200 km³ for the Karakoram–Himalaya region. Recently, various spatially distributed glacier ice-thickness estimation models have been used to estimate glacier volume in the Himalaya and elsewhere (Linsbauer and others, 2009; McNabb and others, 2012; Gantayat and others, 2014; Farinotti and others, 2017, 2019; Millan and others, 2022) but there are significant differences in these estimates (Farinotti and others, 2017; Farinotti and others, 2019). Based on an ensemble of five distributed ice-thickness models, Farinotti and others (2019), for example, estimated glacier volume of the Jhelum (2.2 km³) and Drass (3.4 km³) sub-basins of the Upper Indus Basin. These simulated estimates usually lack validation from ground observations, particularly in the Upper Indus Basin (Haq and others, 2021).

Against this background, we simulated ice-thickness and volume of all the glaciers in the two sub-basins of the Upper Indus Basin; Jhelum and Drass using GlabTop distributed ice thickness model (Linsbauer and others, 2012) and validated the model estimates with field-based thickness estimates observed from several transects across two glaciers, one each from the Jhelum and Drass sub-basins of the Upper Indus Basin, employing 8 MHz GPR. We also used a different approach for ice-thickness estimates but the results did not look promising in comparison to field-based observations. In order to reduce the errors in the estimation of water stored in the Himalayan cryosphere, the purpose of this work was to establish a reliable approach based on simple parameterization for predicting ice thickness and volume estimates over a large region such as Himalaya. The knowledge of glacier thickness and volume is crucial for a wide range of applications, including predicting future water availability, assessing glacier hazards, glacio-hydrological investigations, estimating future sea-level rise and predicting future status and evolution of glaciers, that are rapidly melting due to climate change.

2. Study area

The Upper Indus Basin has an area of ~ 17 000 km² and extends from the HinduKush through Karakoram to the western margins of the Tibetan Plateau (Immerzeel and others, 2010). With elevation ranging from 200 to 8600 m a.s.l., the Upper Indus Basin has a mean annual temperature of −1.3°C and mean annual precipitation of ~644 mm (Zou and others, 2021). Around 10% of the Upper Indus Basin is covered by 11 711 glaciers and 80% of these glaciers are distributed in sub-basins of Hunza, Shigar and Shyok (Bajracharya and Basanta, 2011). The study area in this research comprises of the Jhelum and Drass sub-basins of the Upper Indus Basin. The Jhelum sub-basin has ~ 150 glaciers covering an area of ~ 85 km² (Romshoo and others, 2020b). The basin has ~ 0.7% of its area covered by glaciers mainly confined to the Lidder and Sind watersheds (Romshoo and others, 2021). The Jhelum sub-basin has temperate type of climate dominated by the western disturbances (Dimri and Mohanty, 2009). The Drass sub-basin has a cold semiarid type of climate which is also dominated by the western disturbances and hosts 190 glaciers with an area of ~167 km² (Romshoo and others, 2022b). Both sub-basins receive precipitation in the form of snow largely during winter (Zaz and others, 2019; Romshoo and others, 2022b). Most of the glaciers in both sub-basins are located in the elevation range 4500–5000 m a.s.l.

We selected the Kolahoi and Machoi glaciers in the Jhelum and Drass sub-basins, respectively, for the validation of the modeled ice thickness. The Kolahoi glacier (KG) with an area of ~11 km² forms the largest glacier by area in Jhelum basin (Rashid and others, 2017). It has northerly aspect and lies between 34°07′ to 34°12′N latitude and 75°16′ to 75°23′E longitude (Fig. 1). The melt-waters emanating from the Kolahoi glacier in the form of west Lidder River, join the east Lidder River at Pahalgam, ~35 km from the present snout of the glacier and thereafter flows as Lidder River before discharging into the Jhelum River, one of the major tributaries of the Indus River. The glacier has mean altitude of ~ 4450 m a.s.l. and the slope varies from 7° to 60° with mean value of ~20°. The Machoi glacier (MG) is located in the Drass sub-basin, about 26 km from Sonamarg, the major tourist attraction in the Kashmir valley (Fig. 1). The glacier has an area of ~6 km² and is the largest glacier in the Drass sub-basin (Romshoo and others, 2022b). The Drass River originating from the melt-waters of the Machoi glacier joins the Suru River at the Kargil town (Pall and others, 2019). The glacier, facing north, has a mean altitude of ~4600 m a.s.l. and the slope ranges from 1°–60° with a mean value of 21°.

Fig. 1. Location map of the study region. (a) Jhelum and Drass sub-basins of the Upper Indus Basin. Location of various GPR profiles used for validation of the modeled ice thickness on (b) the Kolahoi glacier and (c) Machoi glacier.

3. Methods and data

3.1. Data sets

The glacier outlines and the branch lines were delineated manually from the snow- and cloud-free Landsat OLI Satellite imagery of October, 2016, complimented by the ASTER Global Digital Elevation Model (GDEM) V2 (Tachikawa and others, 2011). The Landsat OLI image has the spatial resolution of 30 m. The minimal seasonal snow-cover in the study area during the month of October facilitated easy discrimination of various glacier features on the satellite images (Rashid and Abdullah, 2016). The glacier topographic characteristics like elevation, slope and aspect were derived from the ASTER GDEM V2, with horizontal resolution of 30 m and vertical accuracy of 17 m (Tachikawa and others, 2011). The ASTER GDEM has been widely used in glaciological studies at different spatial scales worldwide to derive various topographic parameters such as elevation, slope and aspect (Vignon and others, 2003; Bhambri and others, 2011; Kamp and others, 2011; Wu and others, 2014; Wang and Kääb, 2015; Haireti and others, 2016; Lu and others, 2021). Therefore, we did not see the need to test different DEMs for determining the ice thickness and instead relied on the existing literature where ASTER GDEM has been used. GPR measurements (point measurement) were carried out over the Kolahoi and Machoi glaciers during 15–30 August 2018, along several transects.

3.2. Methodology

3.2.1. Ice thickness modeling

The distributed glacier ice thickness was estimated using the GlabTop model (Linsbauer and others, 2012; Paul and Linsbauer, 2012). The model uses glacier boundaries, glacier branch lines and a DEM as inputs to simulate glacier ice-thickness. The glacier boundaries and branch lines were manually digitized from the satellite images and the Digital Elevation Model, whereas the mean glacier slope was estimated from the ASTER GDEM. The GlabTop model is based on the parameterization scheme presented by Hoelzle and Haeberli, (1995) using a constant value for the basal shear stress (τ_b). The GIS implementation of the GlabTop model (Paul and Linsbauer, 2012) used in the present study estimates the value for the basal shear stress τ_b as a function of elevation relief to derive ice-thickness along branch lines as follows:

(1)$$h = \displaystyle{{\tau _b} \over {\,f\rho gsin\alpha }}$$

where, h is the ice-thickness along the central branch lines, the shape factor f accounts for the friction of glacier body with the valley walls and is computed as the ratio of glacier cross sectional area and the product of its perimeter and the centre-line thickness, ρ is density of glacier ice (900 kg m⁻³), g is acceleration due to gravity (9.8 m s⁻²) and α is the mean glacier surface slope along the branch lines (α ≠ 0). In the present study we used a constant value of f = 0.8 and τ _b = 120 kPa (Zou and others, 2021). The technical details of the GlabTop model and its implementation in GIS are described in detail by Paul and Linsbauer (2012) and Linsbauer and others (2012), and are therefore only briefly summarized here. The GlabTop model calculates the ice-thickness at the points along the branch lines. The ice-thickness values between the base points were interpolated to produce a continuous thickness along the branch lines. The ice-thickness was calculated along 26 and 10 branch lines for the Kolahoi and Machoi glaciers, respectively. Spatial interpolation of the ice-thickness values along the branch lines was performed in ArcGIS (ESRI, 2008) using the ANUDEM interpolation scheme (Hutchinson, 1989) to generate spatially distributed ice thickness maps of both glaciers.

3.2.2. Field measurements

The field-based glacier ice-thickness measurements were carried out using 8 MHz Ground Penetrating Radar, a turnkey system developed by Blue System Integration Ltd (Mingo and Flowers, 2010; https://www.bluesystem.about/ice-penetrating-radar.html). The system comprises of radar receiving unit, a set of antenna (Transmitter and Receiver), rigging components, and collapsible antenna protection tubing together with a ski-based sled system (Mingo and Flowers, 2010). The IceRadar Embedded Processing Unit (EPU) is equipped with data acquisition electronics, GPS unit and a touch-screen embedded computer. The radar and GPS data are acquired and logged continuously at user defined intervals, besides, the radar data acquired is displayed real time in 1-D and 2-D radargrams. All the components are fitted within a water- and dust-proof rugged enclosure (Fig. 2).

Fig. 2. Ground Penetrating Radar (GPR) operation during glacier field survey; (a) Deployment of the ground penetrating ice-radar for surveys on the Kolahoi glacier; (b) The receiver, digitizer and embedded computer system are housed in a rugged water- and dust-proof enclosure. The receiving antenna is connected to the digitizer through a port drilled into the back of the case and (c) Transmitter and battery are housed in another case and mounted on skis during field survey. The transmitting antenna is threaded through ports drilled into case.

The mono pulse transmitter (Narod and Clarke, 1994) delivers 1100 V (±550 V) into a 50Ω resistively loaded dipole antenna at a rate of 512 Hz (Mingo and Flowers, 2010). In the present study, 4 m half-length antenna i.e., 8 m each for receiving and transmitting, with a central frequency of 8 MHz, was used for glacier ice thickness measurements. With a wave speed of 168 m μs⁻¹ in glacier ice and air speed of 300 m μs⁻¹, the ice depth D is calculated as follows:

(2)$$D = \displaystyle{{\;1} \over 2}\left[{{168}^2\left({t + \displaystyle{d \over {300}}} \right)-d^2} \right]^{( 1/2) }\;, \;$$

where, d is the antenna separation and t is the two-wave travel time.

The IceRadar uses the National Instruments NI-5133 digitizer with a sampling rate up to 250 MSs⁻¹, 100 MHz bandwidth, and 512 Hz PRF with 12-bit resolution (Mingo and Flowers, 2010). A small, waterproof, rugged and USB-powered GARMIN NMEA 18 ×  GPS on-board is used for recording the coordinates at each measurement point. The GPR measurements were recorded once every five seconds. In this study, several glacier transects on the two glaciers were surveyed for the ice thickness (Fig. 1). IceRadar Analyzer version 4.2.6 software, as part of the IceRadar system, was used for data viewing, picking, analysis and export for generation of the ice-thickness data from the radargrams. The software allows applications of Dewow, Detrend, gain control and filtering of the radargram data.

3.2.3. Validation of the simulated ice thickness

The GPR measurements are stored as point features with the associated geolocation coordinates. Using the ‘Extract Multi Values to Points’ function in ArcGIS, ice-thickness values from the spatially distributed rasters were extracted for every GPR point measurement.

In order to ascertain the accuracy of the simulated ice-thickness estimates, statistical tests including Correlation Coefficient (CC, Benesty and others, 2009), Relative Bias (RB) (Gumindoga and others, 2016) and the Nash–Sutcliffe Coefficient (NSC, Nash and Sutcliffe, 1970) were used. A CC value of +1 indicates a perfect positive fit, −1 value indicates a perfect negative fit whereas 0 value indicates no correlation at all. The CC is calculated as follows:

(3)$$CC = \displaystyle{{\mathop \sum \nolimits^ \left({Obs-\mathop {Obs}\limits^{{\_}{\_}{\_}{\_}{\_}{\_}} } \right)\left({Sim-\mathop {Sim}\limits^{{\_}{\_}{\_}{\_}{\_}} } \right)} \over {\sqrt {\mathop \sum \nolimits_i^n \left({Obs-\mathop {Obs}\limits^{{\_}{\_}{\_}{\_}{\_}{\_}} } \right)} \sqrt {\mathop \sum \nolimits_{i = 1}^n \left({Sim-\mathop {Sim}\limits^{{\_}{\_}{\_}{\_}{\_}{\_}} } \right)} }}\;, \;$$

The performance of a model, on the basis of RB, is generally categorized into three basic classes: underestimation (bias ≤10%), overestimation (bias >10%), and approximately equal (−10 to 10%) (Brown, 2006). RB is calculated by using the following equation:

(4)$$RB = \left[{\displaystyle{{\mathop \sum \nolimits_{i = 1}^n ( {Sim_i-Obs_i} ) } \over {\mathop \sum \nolimits_{i = 1}^n Obs_i}}} \right]\times 100\;, \;$$

The Nash–Sutcliffe Coefficient (NSC) (Nash and Sutcliffe, 1970), also known as the coefficient of efficiency, is a good indicator for predicting efficiency of a model and varies from minus ∞ to 1, with 1 being the perfect skill and values close to 1 indicate high model efficiency. The NSC is calculated as follows:

(5)$$NSC = 1-\displaystyle{{\mathop \sum \nolimits_{i = 1}^n {( {Obs_i-Sim_i} ) }^2} \over {\mathop \sum \nolimits_{i = 1}^n {( {Obs_i-\overline {Obs} } ) }^2}}\;, \;$$

where, n is the total number of pairs of simulated (Sim) and observed (Obs) ice thickness values, i is the i ^th value of the simulated and observed ice thickness values, $\overline {Sim}$ and $\overline {Obs\;}$ are mean values of modeled and observed glacier ice thickness, respectively.

3.2.4. Error estimation

We estimated uncertainty in the glacier-area estimates, field observed ice-thickness measurements, simulated ice-thickness estimates and glacier-volume estimates using various approaches. The uncertainty in glacier area was found to be equal to 0.2 and 0.1 km² for Kolahoi and Machoi glaciers, respectively, obtained using the following equation (Braun and others, 2019):

(6)$$\delta A = \displaystyle{{Rp/A} \over {Rp/A_P}}\;\;\times \;0.03\;, \;$$

where, Rp/A is the glacier perimeter–area ratio and Rp/A _P is a constant equal to 5.03 km^–1 (Paul and others, 2013). Several other studies have used the same approach to determine uncertainty in glacier area delineation from remote sensing data (Malz and others, 2018; Farías-Barahona and others, 2020). The uncertainty observed in glacier area in this study is in line with the error estimates reported in previous studies over the region (Abdullah and others, 2020; Romshoo and others, 2020b, 2022a).

For determining uncertainty of the GlabTop simulated ice thickness, we relied on the error estimation analysis done by Linsbauer and others (2012) who calculated the model uncertainty by systematic variation in the values of the input parameters including shear stress (τ: ± 30%), glacier slope (α: ± 10%) and shape factor (f: ± 12.5%). The study reported an uncertainty of ±30%, and the same has been considered in other studies like Paul and Linsbauer (2012) and Ramsankaran and others (2018).

Similarly, we considered an overall uncertainty of ±5% for the GPR measurements (Fischer and Kuhn, 2013), which is due to the uncertainty of the signal velocity in a glacier, uncertainty in the antenna separation, uncertainty in the interpretation of multiple reflections, the accuracy of the oscilloscope readings (Fischer and Kuhn, 2013), and errors caused by off-nadir reflections as we did not perform off-nadir reflection corrections (Plewes and Hubbard, 2011; Björnsson and Pálsson, 2020).

The uncertainty in the total glacier volume u(V)/V depends on the ice-thickness (h) and glacier area (A) and was calculated using the error propagation method (Ramsankaran and others, 2018), which is described as follows:

(7)$$\displaystyle{{u( V ) } \over V} = \sqrt {{\left[{\displaystyle{{u( A ) } \over A}} \right]}^2 + {\left[{\displaystyle{{u( h ) } \over h}} \right]}^2} , \;$$

where, u(A) and u(h) are uncertainties in glacier area and ice thickness, respectively.

4. Results

4.1. Ice thickness simulations of the Kolahoi and Machoi glaciers

The estimated, spatially distributed ice thickness of the Kolahoi and Machoi glaciers are presented in Figure 3. The modeled ice-thickness of the Kolahoi Glacier ranges up to 158 ± 47 m with a mean value of 79 ± 24 m. Ice-thickness up to 135 ± 41 m was modeled in the flatter parts of the glacier ablation zone. With an area 10.9 ± 0.2 km², the total ice volume was determined to be 0.9 ± 0.3 km³ based on the GlabTop model.

Fig. 3. GlabTop simulated spatially distributed ice thickness of (a) the Kolahoi glacier and (b) the Machoi glacier. Dashed lines show the central branch lines used in the ice thickness modeling.

Mean ice thickness of the Machoi glacier was 81 ± 24 m ranging up to 162 ± 48 m and the maximum ice-thickness was found in the ablation zone near the Equilibrium Line. The ablation zone of the Machoi glacier has in general thicker ice compared to the accumulation zone. The ice thickness drops from ~60 ± 18 m to ~45 ± 13 m above the ice falls, which coincides with the ELA of the glacier. Based on the GlabTop model, the Machoi glacier stores ~0.44 ± 0.2 km³ of ice within an area of 5.9 ± 0.1 km².

4.2. Validation of ice thickness simulations

The simulated glacier thickness over the two glaciers was validated with the observed GPR ice thickness measurements carried out on the two glaciers.

4.2.1. Kolahoi glacier

The modeled ice thickness across all the four profiles matches quite well with the observed GPR measurements. A representative radargram of the Kolahoi glacier GPR measurements is provided in the supplementary figure (Fig. S1). The GlabTop model in general underestimated the ice-thickness when compared with the GPR measurements (Fig. 4) with an average RB of −10%, CC of 0.96 and NSC of 0.94 across the four profiles. The lowest RB was observed in the Profile 1 (KP1) in the lower ablation zone (Fig. 4) with a mean observed ice thickness of 92 ± 5 m against the mean simulated ice thickness of 86 ± 26 m. Comparison of the observed and simulated ice thickness along KP1 revealed a good agreement, indicated by the high correlation coefficient of 0.97 and a low relative bias of −6%. The NSC value of 0.93 also indicates effectiveness of the GlabTop for simulating glacier ice thickness.

Fig. 4. GlabTop Simulated and GPR observed ice-thickness estimates of four profiles of the Kolahoi glacier. The line graphs have been smoothed using the exponential smoothing filter with smoothing factor of 0.3 in the Excel software. The profiles are numbered from terminus towards accumulation zone of the glacier and were taken across the width of the glacier.

Despite RB of −13% with KP4, highest among the four profiles, the model performed well in capturing the spatial patterns of the ice thickness transect, which is evident by high values of CC and NSC for KP4 (Table 1).

Table 1. Simulated and GPR-observed ice-thickness data along transects on the Kolahoi and Machoi glaciers and statistical evaluation of the relationship between simulated and GPR-observed ice thickness estimates.

4.2.2. Machoi glacier

For Machoi Glacier, modeled ice thickness across all the profiles matches well with the GPR measurements. The simulated ice-thickness of the glacier was validated against nine GPR profiles taken in the ablation zone of the glacier (Fig. 1). A representative radargram of the Glacier is provided as supplementary figure (Fig. S2). The RB between the field-observed and simulated ice thickness was found to be between +6% and ~ −11% with a mean value of ~ −3% (Table 1), indicating good agreement between the two. A high CC value 0.92 and NSC of 0.90 for the nine profiles further corroborates the effectiveness of the GlabTop model (Fig. 5, Table 1).

Fig. 5. GlabTop simulated and GPR-observed ice-thickness estimates of nine profiles on the Machoi glacier. The lines graphs have been smoothed using the exponential smoothing function with a smoothing factor of 0.3 in Excel software. The profiles are numbered from terminus towards accumulation zone of the glacier and were taken across the width of the glacier.

The thickest ice was observed in Profile MP4 with a mean simulated ice-thickness of 91 ± 16 m, which is in good agreement with the mean observed ice-thickness of 96 ± 5 m with a low RB of −6% only. With a slight overestimation of ~ 2%, the lowest bias was observed for the Profile MP5. For Profile MP7, the RB of −11%, was the highest among the nine profiles.

4.3. Ice-reserves of the Jhelum and Drass sub-basins

Given the good agreement between the simulated and observed ice thickness, we estimated the glacier ice-reserves of the Jhelum and Drass sub-basins in the Upper Indus Basin using the GlabTop model. The model simulations revealed that the ice thickness of glaciers in the Jhelum basin ranges up to 187 ± 56 m with a mean ~24 ± 7 m (Fig. 6). The glaciers in the Jhelum basin cover an area of ~ 85 ± 11 km², and hold a glacier volume of ~2 km³ corresponding to an ice-reserve of 1.8 ± 0.6 Gt. The investigation further revealed that 85% of the ice-reserves in the Jhelum are concentrated in the elevation zone 4000–4800 m a.s.l., whereas only 10% of the total ice volume is located between 3300 m and 4000 m a.s.l. (Fig. 6). The simulated ice thickness as well as the altitude distribution of glacier area and ice mass in the Jhelum basin are shown in Figure 6.

Fig. 6. Spatial distribution of the simulated ice thickness of glaciers in the Jhelum basin. The inset graph depicts glacier area and ice storage in different elevation zones of the sub-basin. (a) and (b) provide zoomed-in views of the glacier ice-thickness distribution.

Similarly, with a glaciated area of 167 ± 17 km², the Drass sub-basin has an ice reserve of 2.7 ± 0.9 Gt (Fig. 7). The ice thickness in the Drass sub-basin ranges up to 202 ± 60 m with a mean of ~17 ± 5 m. Around 96% of the total ice mass in the Drass sub-basin is stored in the elevation range 4200–5100 m a.s.l. The simulated ice thickness as well as the altitude distribution of glacier area and ice mass in the Drass sub-basin are shown in Figure 7.

Fig. 7. Spatial distribution of the simulated ice thickness of the glaciers in the Drass sub-basin. The inset graph depicts glacier area and ice storage in different elevation zones of the sub-basin. (a) and (b) provide zoomed-in views of the glacier ice-thickness distribution.

5. Discussion

A comparison of the simulated and field-observed ice-thickness revealed an underestimation of 16% (averaged over four profiles) and ~3% (averaged over nine profiles) for the Kolahoi and Machoi Glaciers, respectively (Table 1). These differences are well within the uncertainty range (±30%) of the GlabTop model (Paul and Linsbauer, 2012). Similar differences have been reported in previous studies using GlabTop-based ice-thickness simulations. For instance, Petrakov and others (2016) reported a 16% underestimation in GlabTop ice-thickness estimates over the Tien Shan region. Similarly, Ramsankaran and others (2018) found an uncertainty of ±14% in GlabTop model ice-thickness estimates when compared to field measurements on the Chotta Shigri glacier in the western Himalaya. The average difference between the field measurements and GlabTop-based ice-thickness estimates found in this study are consistent with the ±15 m range reported by Azam and others (2012).

The robustness and accuracy of the GlabTop model in simulating glacier ice thickness are substantiated by various statistical tests employed in this study. The small deviation of the simulated ice-thickness with respect to the field-based measurements from the two glaciers is supported by high correlation coefficients, low relative biases and a robust Nash–Sutcliffe Efficiency for both the Kolahoi and Machoi glaciers. These findings indicate that the model is efficient for simulating glacier ice thickness (Table 1). The model consistently captured the general ice-thickness pattern, and the magnitudes across all surveyed profiles over the two glaciers. The widely-available input parameters required by GlabTop compared with other available models, which required data such as mass balance and velocity (Gardner and others, 2013; Vincent and others, 2013; Ramsankaran and others, 2018), is an advantage of the GlabTop model for computing distributed ice-thickness over the data-scarce Himalaya.

Furthermore, the reliability of the GlabTop model for estimating glacier volume in the Himalayan region has been previously reported by Linsbauer and others (2014) and Frey and others (2014). A recent study by Zou and others (2021) further demonstrated the reliability of GlabTop for ice-thickness estimation over a part of the Upper Indus Basin. It is noteworthy that we also explored a velocity-based approach to estimate ice thickness, where we used the by laminar flow equation to relate glacier surface velocity and ice thickness (Cuffey and Paterson, 2010; Gantayat and others, 2014). However, in comparison to the GPR measurements, this velocity-based model performed poorly in simulating glacier ice thickness, showing an average relative bias of −28%. Consequently, we have not included the results from this approach in this study.

The average ice-thickness in the Jhelum (24 m) and Drass (17 m) sub-basins is lower than the reported average ice thickness for the entire Upper Indus Basin (~75 m) by Zou and others (2021). This difference can be attributed to the distribution of ice thickness at high altitudes, probably influenced by steep headwalls that hinder snow accumulation and glacier formation at these elevations (Nagai and others, 2016; Zou and others, 2021). However, both this study and Zou and others (2021) found a good agreement between modeled and GPR-based ice thickness. As expected, the larger glaciers have thicker ice in both sub-basins, while ice tends to thinner on steeper slopes at higher altitudes compared to ice on more gentle slopes at lower altitudes, consistent with the findings by Salerno and others (2017) and Hoelzle and Haeberli (1995). It is noteworthy that the glacier volume estimates for the Jhelum and Drass sub-basins presented in this study are respectively 13 and 17% lower than those provided by Farinotti and others (2019). The volume estimates of the Jhelum sub-basin by Millan and others (2022) are 2.3 km³, or ~17% higher than results presented in this study. Millan and others (2022), furthermore, estimated the volume of the Drass sub-basin glaciers to be 5.5 km³, nearly double the volume estimate of this study. Additionally, Zou and others (2021) estimated the ice volume in the Drass sub-basin as 7.6 km³, or 150% higher than the result of this study. The differences between our study and previous studies can be attributed to the different approaches employed. For instance, Farinotti and others (2017, 2019) averaged their ice-thickness estimates from multiple models, while Millan and others (2022) relied on a velocity-based approach. Additionally, disparities may arise from differences in model input parameters, including the glacier boundaries, delineation errors, and the extent and number of glacier branch lines. These differences probably account for the discrepancies observed in the ice-thickness estimates between our study and that of Zou and others (2021).

Additionally, we computed the glacier volume across the study area using the volume–area scaling approach (Bahr and others, 1997) and found 3.8 km³ for the Jhelum and 7.5 km³ for Drass sub-basins. These values are considerably higher (97% for Jhelum and 150% for Drass) than our results (Supplementary Fig S3).

6. Conclusions

The lack of ice-thickness measurements has led to significant uncertainty regarding the total glacier ice volume in the Upper Indus Basin and consequently, projections of changes in glacier mass and volume under a changing climate. This study is the first reporting of ice- thickness measurements for two glaciers in the Jhelum and Drass sub-basins of the Upper Indus Basin. The study found that the GlabTop model generally performs well in modeling the spatial distribution of glacier ice thickness, albeit typically underestimating the thickness by on the order of 10% on average. The model proved robust in modeling the ice thickness along almost all selected field transects for GPR-based measurements on the two glaciers, as demonstrated by the statistical evaluation of the simulated ice-thickness profiles. The relative bias between observed and simulated ice thickness for the Kolahoi and Machoi glaciers was −11% and ~ 2%, respectively, falling within acceptable error limits.

These improved estimates of glacier ice thickness and volume are essential for accurate projections of glacier-melt contributions to sea-level rise, the impact on streamflow, and informed water resource management under a changing climate. This is particularly important for policymaking for food, energy and water security in the Indus basin, where waters are shared among neighboring countries across political boundaries.