Mass balance of lake terminating Gepang Gath glacier (western Himalaya, India) and the role of glacier–lake interactions

Abstract

The mass balance of lake-terminating glaciers responds to annual atmospheric variations, while calving-induced ice loss at the front is driven by local ice–water interactions. The current glaciological studies underestimate glacier response by neglecting the significant annual ice loss at the terminus through calving processes. This study integrates field measurements with remote sensing data to investigate the glaciological characteristics and proglacial lake evolution of the Gepang Gath glacier in the Chandra basin, Western Himalaya, India. Long-term observations reveal a continuous expansion of the proglacial lake from 0.21 ± 0.06 km² (1962) to 1.21 ± 0.05 km² (2023), along with terminus retreat of ∼2.76 km, attributed to calving at the ice–water interface. The glacier’s surface exhibits complex debris cover, with thicknesses up to 35 cm, creating significant spatial variations in surface mass balance. In-situ, glaciological measurements reveal a highly negative glacier-wide mass balance of −0.90 ± 0.30 m w.e. a⁻¹ between the years 2014 and 2023. The geodetic estimates also reveal a negative mass balance of −0.61 ± 0.1 m w.e. a⁻¹ over the past decade (2013–2023). The frontal area change (0.42 km²) and geodetic mass balance show a total volumetric ice loss of −21.77 × 10⁶ m³ w.e. during the same period. Overall, the yearly frontal ice loss exacerbates the mass loss by 17–22%. These findings suggest that the presence of proglacial lakes plays a significant role in intensifying ice mass loss from Himalayan glaciers, strongly regulating their overall evolution.

1. Introduction

The Himalaya is a critical global water tower, storing approximately 0.77% of the world’s freshwater in the form of glaciers and snow cover (Azam, 2021). These high-altitude mountain glaciers and ice caps, separated from the down valley settlement, provide essential freshwater resources for over a billion downstream inhabitants (Jansson and others, 2003). Several studies, primarily based on remote sensing techniques, have documented the continuous depletion of the Himalayan glaciers due to the ongoing warming, particularly over the past four decades (Bolch, 2012; Brun and others, 2017; Maurer and others, 2019; Shean and others, 2020; Hugonnet, 2021; Garg and others, 2022). However, field-based validations of remote sensing estimates are still limited in the Himalayan glaciers, with in-situ glacier observations (mass balance) restricted to only 35 glaciers (Azam and others, 2018; Azam, 2024). Glaciers across the Himalaya exhibit a consistent negative mass balance trend, although with spatiotemporal heterogeneity induced by regional climate and debris cover. For instance, Dokriani and Chorabari glaciers in the Garhwal Himalaya show an average loss of −0.72 and −0.32 m w.e. a⁻¹, while the Chhota Shigri glacier in the western Himalaya records −0.46 m w.e. a⁻¹. Mera and Pokalde Glaciers in Nepal experience negative mass balances of −0.23 and −0.79 m w.e. a⁻¹, respectively (Mandal, 2020; Vishwakarma, 2022). At regional scale, the average geodetic mass balance across the Himalaya reveals significant ice loss over the past 40 years (1975–2016), with ice loss doubling post-2000 (Maurer and others, 2019). Negative mass balance and continuous glacier retreat generate various surface morphological features such as supraglacial ponds, ice cliffs, and debris-covered areas over the ablation zone of glaciers (Zhang, 2024). Particularly, glacier retreat often leads to the formation of proglacial lakes as meltwater fills depressions left behind by receding ice (Basnett and others, 2013; Patel and others, 2017; Sattar, 2023). These proglacial lakes typically begin as small ponds near the glacier’s terminus, dammed by end moraine ridges (Carrivick and Tweed, 2013). Over time, with the glacier’s backwasting and favorable topographic condition, these ponds often expand to form larger water bodies in direct contact with the glacier stream and frontal ice cliffs (Richardson and Reynolds, 2000; Sato and others, 2022). Their formation is generally associated with glaciers lacking efficient drainage pathways for meltwater. The rate of proglacial lake growth is directly linked to the rate of glacier retreat and mass loss, with increasing lake area typically observed alongside a negative mass balance(Ives and others, 2010; Pandit and Ramsankaran, 2020; Gantayat and Ramsankaran, 2023). In some cases, the proglacial volume can even reach or surpass the volume of the remaining glacier ice (Pandit and Ramsankaran, 2020). Several studies have highlighted the potential for ongoing continuous mass loss of Himalayan glaciers and proglacial lake expansion to create new glacial lake outburst flood (GLOF) hazard zones (King and others, 2019; Wang, 2020; Sattar, 2023).

Glaciers that feed the proglacial lakes undergo a rapid change in surface morphology, such as the formation of a calving front, enhanced retreat rate and higher rate of ice mass loss (Shukla and others, 2018; Mohanty and Maiti, 2022; Peacey and others, 2023). Notably, calving processes can contribute significantly to annual ice ablation and have often been overlooked as a component of mass balance assessments in the Himalaya (Zhang, 2023). This phenomenon is highly site-specific, and few studies have specifically focused on its impact (Watson and others, 2020). A few regional-scale studies have successfully captured the general trend of mass loss associated with lake-terminating glaciers (Zhang and others, 2019; Zhang, 2023). This knowledge can be used to refine future studies and provide more accurate estimates of the current and future state of glaciers, as well as the response of proglacial lakes to climate change (Kaushik and others, 2020; Watson and others, 2020). Additionally, research has explored the role of lake width (lake–ice interface) in influencing mass balance dynamics of lake-terminating glaciers, investigating whether wider lakes promote more negative mass balance (Kirkbride and Warren, 1997).

Motivated by these observations, the present study focuses on the lake-terminating Gepang Gath glacier in the western Himalaya. Our primary aim is to assess the glacier mass balance trends, quantify retreat rates, and understand how proglacial lake expansion drives ice loss at the glacier front. The specific objectives of this study are: (a) to measure the in-situ (glaciological) mass balance of the Gepang Gath glacier during the last decade, (b) to compare the glaciological mass balance estimates with geodetic data for the corresponding period, (c) to quantify long-term changes in proglacial lake between 1962 and 2023 and (d) to assess the impact of glacier–lake interactions on the overall mass balance and state of the glacier. To achieve these, we employed an integrated approach of in-situ measurements with the multisensor and multitemporal remote sensing-based observations acquired between 2014 and 2023.

2. Glacier and Climate Setting

The Gepang Gath glacier is situated within the Chandra Basin of the western Himalaya (Himachal Pradesh, India) (Fig. 1). The total catchment area of the basin is approximately 2446 km². Notably, the Chandra Basin hosts a significant glacial coverage, with a total of 211 glaciers occupying an area of ∼631 km² (Oulkar and others, 2024). The meltwater generated from these glaciers contributes to the Chandra River. At an elevation of 2850 m a.s.l. (near Tandi), the Chandra River merges with the Bhaga River to form the Chandra-Bhaga River. The region falls in the higher Himalaya under the monsoon-arid zone (Bookhagen and Burbank, 2006).

Figure 1. Location map of the Gepang Gath catchment in the Chandra Basin, Himachal Pradesh, western Himalaya, India. The location of the main map (i.e. Gepang Gath glacier) is shown in the inset. A dense Differential Global Positioning System (DGPS) survey (yellow line) was done in 2021 and 2022, mainly over the lower ablation area and the peripheral of exiting proglacial lakes. Debris thickness was measured at 50 locations (blue dots), and the point surface mass balances at each site (red bar-ablation; yellow star-accumulation). For the presentation, only the locations of the 2018–2019 stakes are plotted. The background image is Sentinel 2 imagery of October 2018 at 10 m resolution. The glacier outline (green polygon) corresponds to the same image.

2.1. Gepang Gath glacier and proglacial lake

The Gepang Gath glacier is a compound-type valley glacier consisting of two tributary glaciers that merge with a main trunk glacier at an elevation of 4200–4300 m (Fig. 1). The main trunk glacier exhibits debris-free ice down to an elevation of 4200 m, while the debris-covered tributaries extend to a higher elevation of 4400 m. Details of the glacier characteristics are provided in Table (S1). Years of glacier recessions have created a topographic depression, now occupied by a proglacial lake having a total length of ∼2.76 km (as of 2023, Fig. 1). The lake’s spillway, located at an elevation of 4080 m, has maintained its position since at least 1962 (Rathore and others, 2015). Observations suggest that the right side of the terminal moraine is more prone to erosion by overtopping flow events (Worni and others, 2013; Sattar, 2023). The meltwater from Gepang Gath glacier directly feeds the lake (Fig. 1). Additionally, six small streams, mainly on the northern side of the glacier, contribute to the lake’s water volume (Fig. 1). Currently, the glacier terminus has transformed into a concave-shaped calving front with a total width of ∼900 m (Fig. 1). Because of calving processes, the lower ablation area of ∼1 km length is characterized by numerous transvers crevasses. The warm lake water deep into the ice, which may also afloat the glacier ice is yet to be confirmed.

3. Data and Methodology

Building on the previous studies (Patel and others, 2017, 2021), we digitized the area of Gepang Gath proglacial lake from 1962 (based on Survey of India toposheets) to 2017 (using Sentinel-2 Multispectral Instrument (MSI) satellite imagery). Subsequently, we obtained proglacial lake area data for 2020 from (Sattar, 2023) and have conducted a ground-based DGPS survey of the lake periphery in 2022. We have been monitoring the glacier-wide mass balance measurement since 2014, and this study presents glacier-wide mass balance data up to 2023 (Table S2). Our analysis focused on changes in the calving front, including cumulative and yearly retreat rates, as well as the change in the frontal glacierized area. The geodetic mass balance and the total volume loss at the calving front were derived by comparing the September 2013 and October 2023 digital elevation models (DEM) at 30 m resolution.

3.1. Kinematic DGPS survey

We made differential kinematic surveys using Trimble R9 carrier phase dual-frequency GPS receivers at a 1 s logging interval over the terminus area, the periphery of proglacial lakes and the lower ablation area of the Gepang Gath glacier during the summer field seasons of 2021 and 2022 (Fig. 1). Three rover receivers were installed on each of the wooden boxes that were carried as a backpack. Data were collected at 1 s intervals, giving the nominal data spacing of about 1 m, the data density was populated while the rover was stationary. The coordinate position was later processed in kinematic mode using a Canadian precise point-processing service (http://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php). The uncertainty in the kinematic survey over the glacier catchment, after post-processing, is ±10 cm in horizontal positioning and ±20 cm in vertical measurement.

3.2. Debris thickness measurement

Debris thickness across the ablation zone of the Gepang Gath glacier was measured manually at 50 points through debris pit measurements. Debris-free ice was assigned to zero debris thickness. Measurements were only conducted at locations where ablation stakes were installed every year. A handheld Garmin GPS unit is used to determine the location with a horizontal accuracy of ±3 m. The error in debris thickness is associated mainly with the height measurement from the ice surface to the debris surface. To minimize the error, we estimated depth on all sides and took the mean value as the final debris thickness. The standard deviation (std) of these values was assumed to be the error in debris thickness.

3.3. Frontal length and area

The annual retreat rate at the front was estimated using Sentinel-2 MSI imageries from 2014 to 2023. The images were obtained from the Copernicus Open Access Hub (https://scihub.copernicus.eu). The images are selected with a minimum cloud cover of <10%, mainly during September and October. The calving front of the glacier was identified and established as the interface between the glacier and the proglacial lake. The frontal outline was initially delineated for the year 2014 and then adjusted according to the temporal satellite images for the years 2015 to 2023.

3.4. Glaciological mass balance

For the mountain glacier, the glaciological mass balance method is well described (Østrem and Brugman, 1966; Cogley, 2010). We determined two components for Gepang Gath glacier-wide mass balance (i.e. surface ablation and surface accumulation). It represents the net difference between accumulation and ablation processes across the entire glacier and is expressed in meters of water equivalent (m w.e.). For surface ablation estimation, 6–14 bamboo stakes were drilled into the ice at several points over the glacier surface along with different altitudinal zones for every balance year (i.e. 1 October to 30 September) (Fig. 1 and Table S2). We followed the fixed date system, where 30 September was taken as the end of the mass balance year. Therefore, stake networking was done mostly in September. Surface ablation at a point location was calculated by the sum of the exposed length of stake measured over a balance year. To estimate ice ablation in meter water equivalent from the stake height change, a mean density of ice 870 ± 25 kg m⁻³ was considered (Pratap, 2019).

Surface accumulation was measured every year at the end of the ablation period (i.e. 15–30 September) via measuring snow density and thickness using a snow corer (Table S2). The thickness of snow accumulation was adjusted with density for total site-specific surface accumulation. The snow/firn density varied from 470 kg m⁻³ to 570 kg m⁻³. Due to sparse measurements and steep gradients in the accumulation area, we applied a linear increase between measured points and extrapolated up to 5000 m a.s.l. Beyond 5000 m a.s.l., particularly in steep slopes (covering 19% of the area), a constant extrapolation was applied to estimate total accumulation.

The mean annual altitudinal mass balances were estimated for each 50 m elevation band by averaging all point mass balance data within that band. Finally, the glacier-wide mass balance (Ba) is calculated by integrating these mean annual altitudinal mass balance values over the entire glacier surface area using the following equation:

(1)\begin{equation}B{\text{a}} = \frac{1}{S}\mathop \sum {b_{\text{z}}}{s_z}\end{equation}

where Ba is the glacier-wide mass balance (m w.e. a⁻¹) and bz is the mean altitudinal mass balance of the altitudinal area of Sz (m²). $S$ is the total glacier area (m²).

The uncertainty in glaciological mass balance depends on uncertainties of in-situ measurements of point surface ablation and accumulation, and the extrapolation of point mass balance to the entire glacier area. Following the refs (Kenzhebaev and others, 2017; Stumm and others, 2021; Pratap and others, 2023; Azam, 2025), we have estimated random and systematic errors for the glacier-wide annual mass balance.

The error in annual surface ablation (Δɑ _a) of each point measurement was derived from error in stake height change measurement (Δ_h), ice density used to convert height change into water equivalent (Δ_d), and error in debris thickness from the ice surface to debris surface (Δ_dt) owing to varying debris thickness (Dt). The final error associated with the surface ablation measurements of individual stakes during the mass balance year was calculated as follows (Pratap and others, 2023):

(2)\begin{equation}\frac{\Delta a_a}{a_a}=\text{ }\sqrt{\left(\frac{\Delta_h}H\right)^{2\text{ }}+\left(\frac{\Delta_d}D\right)^{2\text{ }}+\left(\frac{\Delta_{dt}}{Dt}\right)^{2\text{ }}}\end{equation}

This results in a mean uncertainty in net surface ablation of ±0.19 m w.e. a⁻¹ for the period from 2015 to 2023.

The error in surface accumulation measurements arises from snow core density variations, measurement errors and surface roughness effects (Thibert and others, 2008; Azam, 2025). While point accumulation measurements provide localized estimates, surface roughness and topographic variations introduce additional uncertainties. Furthermore, because total snow accumulation relies on sparse data limited to a specific elevation, this introduces uncertainties in the derived results, particularly the spatial distribution of snow accumulation. The overall uncertainly in net surface accumulation is estimated by

(3)\begin{equation}\Delta ac = \sqrt {{\Delta _{rgh}}{d^2} + \Delta {d^2}{\text{ }}{l^2}} \end{equation}

where ${\Delta _{rgh}}$ is surface roughness 0.35 m w.e. taken from (Thibert and others, 2008), $d$ is the mean density of snow/firn accumulation, $\Delta d$ is uncertainty in snow/firn densities along the length (l) of the core. The uncertainty in net surface accumulation was on average of ±0.22 m w.e. during 2015–2023.

Considering the uncertainties in net surface accumulation and net surface ablation, we estimated the glacier-wide mass balance uncertainty to range from ±0.20 to ±0.39 m w.e. a⁻¹ with a mean annual uncertainty of ±0.30 m w.e. a⁻¹ during 2014/15–2022/23. Uncertainty in mass balance measurements for individual Himalayan glaciers using the glaciological method varies between ±0.27 and ±0.53 m w.e. (Soheb and others, 2020; Angchuk, 2021; Wagnon, 2021). This range is consistent with the uncertainty estimated for the Gepang Gath glacier.

3.5. Geodetic mass balance

The geodetic elevation change of the Gepang Gath glacier was computed by comparing the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEMs of 17 September 2013 and 30 October 2023. The ASTER DEMs were acquired from NASA’s Earth Observing System Data and Information System (https://www.earthdata.nasa.gov). Before comparison, both the DEMs underwent pre-processing steps to ensure accuracy. First, any erroneous peaks and sinks identified in the DEMs were removed using the thin plate spline method. Then, the 2023 ASTER DEM was taken as the reference DEM and the 2013 ASTER DEM was co-registered to it. The horizontal congruence between both the DEMs was achieved by minimizing the standard deviation of the elevation differences over the stable terrain (Nuth and Kääb, 2011; Berthier and others, 2016). The stable terrain was achieved by masking the glaciers (using RGI7.0 glacier outlines), steep slopes (<4° and >45°) and elevation difference outliers (>±100 m) (Berthier and others, 2016; Garg and others, 2019). Subsequently, the DEMs were checked for vertical biases such as along and cross-track biases originating from satellite attributes. To correct them, the azimuth of ASTER’s ground track was used to rotate the coordinate system and then the elevation differences over stable ground were used to correct the biases (Gardelle and others, 2013). Then, considering the fact that the elevation changes are not homogeneous across the glacier surface and varies as a function of altitude (Pieczonka and Bolch, 2015), a method from (Zhou and others, 2018) was applied to identify outliers at each particular altitude:

(4)\begin{equation}{E_n} = {\text{ }}\frac{{{E_{mx}} - {E_{gl}}}}{{{E_{mx}} - {E_{mn}}}}\end{equation}

(5)\begin{equation}{{\Delta }}{h_{mx}} = A{\text{*}}{E_n}\end{equation}

where E_n is normalized elevation, E_mx and E_mn are the maximum and minimum glacier elevation, E_gl is the elevation of individual glacier pixels, Δh_mx is the maximum allowable thickness change at a particular elevation, and A is the empirical coefficient. The empirical coefficient A represents the maximum elevation change at the glacier front which was set to 65 m in order to preserve the true elevation change values while identifying the outliers (Zhou and others, 2018; Garg and others, 2022). The eliminated values using Δh_mx were replaced by the mean of corresponding 25 m altitude bins.

The volume change (i.e. area-weighted mean) for each 25 m altitude bin was computed by multiplying the area and elevation difference of respective bins. Finally, the geodetic mass balance (B _geo) was calculated using the ice density value of 870 ± 25 kg m⁻³ for the ablation zone and snow density of 570 kg m⁻³ for the accumulation zone as suggested for the study region (Pratap, 2019).

To quantify geodetic mass balance uncertainties, first, we used the standard deviation of elevation differences over non-glaciated stable terrain ( ${\delta _{ng}}$) to approximate the error in elevation change computations (δh_ng) by dividing it by the square root of a total number of effective pixels (N_eff) (Gardelle and others, 2013; Zhou and others, 2018). The N_eff was obtained using a decorrelation length of 500 m (Brun and others, 2017; Shean and others, 2020). Finally, the overall uncertainty in mass balance computation $\left( {{\delta _{Bgeo{\text{ }}}}} \right)$ was quantified by incorporating the error in density assumption $({\rho _{si}})$ as per the following:

(6)\begin{equation}{\delta _{Bgeo}} = {\text{ }}\sqrt {{{({\delta _{hng}})}^{2{\text{ }}}} + {{({\rho _{si}})}^{2{\text{ }}}}} \end{equation}

4. Results

4.1. Evolution of the proglacial lake and the surrounding area

The proglacial lake area of Gepang Gath glacier has changed significantly since the 1960s (Fig. 2). As the glacier has retreated over time, the proglacial lake has also expanded dramatically. We combined all published data with our kinematic DGPS survey around the lake perimeter in 2022 (Figs. 1 and 2). Over the past four decades (1979–2017), the lake area has expanded from 0.17 ± 0.06 km² to 0.84 ± 0.04 km². Incorporating the 2017 lake area of 0.84 ± 0.04 km², the DGPS survey conducted in October 2022, and Sentinel-2 imagery from October 2023, the total area is recorded to have expanded to 1.21 ± 0.05 km². A significant expansion rate in recent years was recorded between 2017-2023 when the area of the lake expanded by ∼0.06 km² a⁻¹. The trend in lake expansion shows a steady increase in lake area from 1962 to 2000, a moderate rise from 2000 to 2010, and a steep increase from 2010 to 2023, reflecting the accelerated expansion of the proglacial lake in recent times (Fig. 2). The higher expansion rate and steep slope post-2010 is linked to the larger width of the glacier terminus loss and the interaction of lake water with the glacier ice at the bottom.

Figure 2. Time series progression of the proglacial lake area from 1962 to 2023 based on the published and the year 2023 lake area data (area in m² provided in the literature was converted to km²). Dots show the lake area observed in a specific year. The solid color line shows the positive correlation of lake area expansion between the time intervals of (i) 1962–2000, (ii) 2000–2010, and (iii) 2010–2023. Standard errors (uncertainties) are included for each rate of change in the lake area.

4.2. Glacier surface morphology, frontal length and area change

The combined analysis of the satellite images and debris thickness estimation revealed a complex surface topography within the ablation zone of the Gepang Gath glacier. While the overall surface slope in this zone is notably gentle (less than 20°), there are small-scale undulations with steeper slopes (40°–60°) reaching up to 100–200 m in length. These undulations have an elevation difference of 10–15 m, primarily attributed to the formation of ice cliffs and supraglacial channels (Figs. 3a,b). Several surface features, such as the transverse crevasses, ice cliffs, and debris tills, are also evident over the ablation area (Figs. 3b,d). Two tributaries and a main centreline further characterize the glacier surface topography alongside the debris-covered region, where debris thickness reaches up to 35 cm, a portion of debris-free ice extending along the central line up to an elevation of 4250 m (Fig. 3a). The right and left lateral moraines are elevated by ∼100 m to the present glacier surface. The moraines are highly eroded and fragile. The supraglacial lakes are not prominent across this glacier’s ablation area. Furthermore, there has been no notable increase in supraglacial lakes over time, even in the context of an increase in debris coverage.

Figure 3. Characteristics of Gepang Gath glacier and its proglacial lake. (a) In-situ debris thickness measurement sites with the size of the orange color dots representing the debris thickness over the ablation area. Blue dots represent the debris-free sites based on the stakes installed during 2014–2023. Debris thickness measurements were mostly made in the year 2019, hence plotted over the glacier outline derived in the year 2019. (b) The frontal recession of the Gepang Gath glacier is shown over the Sentinel-2 MSI image of 2023. The frontal length change was observed on a yearly basis from 2014 to 2023. (c–d) field photos of the calving front and proglacial lake and the ablation and accumulation area of the glacier.

The past terminal retreat of this glacier was earlier reported by (Kumar and others, 2021) for four distinct periods vis-a-vis 1989–2000 (14 m a⁻¹), 2000-2009 (38 m a⁻¹), 2009–2013 (46 m a⁻¹) and 2013–2015 (50 m a⁻¹). We further examined the frontal length and area change from 2014 to 2023. Our study integrated with (Kumar and others, 2021) reveals a frontal retreat of ∼53 m a⁻¹ (Fig. 3b), indicating an accelerating retreat rate through time. Factors like increasing glacier widths (lake–ice-moraine interface) promote the loss of larger, thick-iced glacierized areas. The transverse profile (width) of the glacier’s front in 2005 was ∼760 m, while in 2014, the width increased to ∼840 m, and in 2023, it reached ∼940 m. This reflects the increased ice mass loss and the expansion of the lake with time (Fig. 3). The yearly length changes of the glacier correlate with the size of the calving ice and the width of the glacier’s frontal portion. Over the period from 2014 to 2023, the glacier has retreated about ∼482 m, and the total glacier area decreased by 0.40 km², which is accompanied by the expansion of the lake. As the glacier terminates into the lake water, it has a significant impact on the dynamics of the glacier ice, as several transverse crevasses are visible in that terminal area (Figs. 3c,d). This glacier is expected to shrink in the future, and conversely, the lake is expected to expand further, either by calving at the front or by the higher negative mass balance of the glacier.

4.3. Glacier-wide mass balance, gradient, ELA and AAR

The annual glacier-wide surface mass balance (Ba) of the Gepang Gath glacier has predominantly exhibited a negative trend, averaging at −0.90 ± 0.30 m w.e. a⁻¹, except for positive Ba of 0.10 ± 0.28 m w.e. recorded during the mass balance year 2018/19. These values are detailed in Table (1) and shown in Fig. 4. The highest negative Ba was observed during 2017/18, reaching −1.34 ± 0.22 m w.e., closely followed by the Ba − 1.30 ± 0.39 m w.e in 2021/22. The cumulative mass change from 2014/15 to 2022/23 was −7.66 ± 0.84 m w.e. In the lowest ablation area, which is mostly debris-covered and heavily crevassed, the point mass balance during 2014–2023 varies significantly from −5.00 to −1.00 to m w.e. a⁻¹ (Fig. 4). These differential point mass balance rates within similar elevation bands relate to the presence of spatially varying debris thickness. The high interannual variability in ablation rates results from stake replacement under varying debris thickness and differences in seasonal meteorological conditions such as seasonal snow cover and temperature. Also, the spatial distribution of debris cover/thickness over the glacier likely varies over time under the influence of widening crevasses, enlarging ice cliffs, and growing supraglacial ponds. The stakes along the central line (debris-free ice) area showed higher ablation (>−5.0 m w.e. a⁻¹). In comparison, stakes near the lateral moraines with an area under debris thickness of >20 cm had a significant reduction in ablation rate (<−2.0 m w.e. a⁻¹). To show the mass balance gradient, we have presented the annual point ablation and accumulation rates as a function of elevation for every studied year (Fig. 4). The mass balance gradient (i.e. db/dz) shows a linear trend with the mean value of 0.83 m w.e. (100 m⁻¹) (Table 1). The equilibrium line altitude (ELA) of the Gepang Gath glacier ranged from 4435 to 4640 m a.s.l. (average 4583 m a.s.l.). The accumulation area ratio (AAR) ranged from 42% to 62% (average 49%). Furthermore, our analysis revealed a balance budget ELA0 of 4445 m a.s.l. and an associated AAR0 of 62%. These values are consistent with the typical range of AAR0 (55–65%) for glaciers in the Himalaya.

Figure 4. (a) Hypsography, (b–i) point mass balance (ablation: orange dots, accumulation: yellow dots) and the mean annual altitudinal mass balances (black dots) at every 50 m elevation interval of Gepang Gath Glacier for the balance years between 2014/15 and 2022/23, with a gap year of 2019/20. The straight red represents the mass balance gradient (i.e. db/dz) as a linear fit with 95% confidence bond.

Table 1. Annual net surface ablation (Aa), net surface accumulation (Ca) and glacier-wide mass balances (Ba), ELA and AAR were calculated using the glaciological mass balance method. The ablation and accumulation area, and mass-balance gradients are also presented for the Gepang Gath glacier from 2014/15 to 2022/23. There is no observation for the year 2019/20

4.4. Geodetic mas balance

Geodetic mass balance at 30 m resolution was calculated for the Gepang Gath glacier (Fig. 5a). The surface elevation change analysis from 2013 to 2023 indicates significant thinning of the Gepang Gath Glacier, particularly in the lower ablation zone (Fig. 5a). Notably, the elevation difference on the glacier was highly heterogeneous on the spatial scale. Surface thinning (negative values of surface elevation changes are referred to as surface lowering) is observed on ∼50% of the glacier surface, with 100% thickness loss occurring primarily at elevations below 4150 m. This excessive surface lowering at the terminus reflects the influence of the proglacial lakes which causes the terminus portion to melt or break off through calving (Fig. 5a). The maximum volume loss (area-weighted surface lowering) is observed at 4100–4200 m altitude band while the maximum volume gain is observed at 4700–4775 m a.s.l. Fig. 5b illustrates the variability in geodetic mass balance across different elevations, including the average geodetic mass balance computed for every 25-m elevation. The geodetic mass balance gradient shows significant negative values at lower elevations, indicating substantial mass loss due to melting (Fig. 5b). As elevation increases, the mass loss decreases, reaching values near zero or slightly positive at higher elevations. The total geodetic mass balance for the 10 year monitoring period was calculated as −6.08 ± 1.2 m w.e., which translates to an annual mass balance of −0.61 ± 0.1 m w.e. a⁻¹.

Figure 5. (a) Surface elevation change profile on the Gepang Gath glacier between 2013 and 2023 and (b) average geodetic mass balance of every 25 m elevation band of Gepang Gath glacier.

5. Discussion

5.1 Comparison of in-situ and geodetic mass balances

The geodetic mass balance estimates for the Chandra basin are available in different time frames reflecting the mass loss of −0.68 ± 0.15 m w.e. a⁻¹ from 1999 to 2011 (Gardelle and others, 2013) and −0.65 ± 0.04 m w.e. a⁻¹ from 2000 to 2012 (Vijay and Braun, 2016). While there is no specific geodetic mass balance was observed for the Gepang Gath glacier. Using ASTER data, we have derived surface ice thickness change for the period from September 2013 to October 2023. Further, a spatial surface mass balance was generated by interpolating all yearly point-wise surface mass balance data from 2014 to 2023 to measure spatially variable surface mass balance with respect to surface elevation. The comparison between surface mass balance and ice thickness loss (geodetic mass balance) reveals key differences in ice mass loss (Fig. 6a and b). The surface mass balance ranges from −7.5 to 2 m w.e., showing a gradual transition across the glacier, while the geodetic mass balance varies from −7 to 2 m w.e, exhibiting localized thinning patterns with more spatial variability. The lower glacier tongue shows significant negative values in both cases, but geodetic mass balance captures more heterogeneous thinning, especially around the lower ablation zone with significant surface lowering of up to 90 m (Fig. 6b). The upper glacier regions show lower mass loss in both datasets, yet surface mass balance appears more evenly distributed, whereas geodetic mass balance reflects topographic influences.

Figure 6. Surface mass balance of Gepang Gath glacier: (a) Annual mass balance (points) extrapolated over a digital elevation map to illustrate the spatial surface mass balance pattern from 2014 to 2023. (b) Geodetic mass balance for the period from Sept 2013 to Oct 2023, a zoomed view of the rectangle is shown in (c–e) with three years of cumulative ice volume loss estimated based on the annual retreat rate. (f) Total geodetic ice volume loss from the glacier front during 2014–2023.

As noted, the average annual glacier-wide mass balance (−0.90 ± 0.30 m w.e. a⁻¹) of the Gepang Gath glacier is substantially high compared to the other valley glaciers in the Himalaya (Pratap and others, 2016; Mandal, 2020; Stumm and others, 2021; Romshoo and others, 2023). We explore this further by focussing on the impact of calving due to the proglacial lake on the annual glacier-wide mass balance. The frontal ice loss by calving, which refers to the breaking-off or detachment of ice chunks from the glacier’s terminus, plays a significant role in determining the overall mass balance of a glacier (Watson and others, 2020). When ice calves from the glacier to the lake, it contributes directly to the total ice loss and provides more surface area for the lake’s expansion. Recently, Zhang (2023) showed that the previous estimate of the total mass loss of lake-terminating glaciers in the Himalayas is underestimated by 6.5 ± 2.1%, emphasizing the need for proper accounting of proglacial lake–ice interaction. Figs. 6c–f shows the trend of ice loss and total volume loss at the frontal part of the glacier between 2014 and 2023. The total volumetric frontal ice loss at the glacier is −21.77 × 10⁶ m³ w.e., which includes a surface area loss of 0.42 km² during the period from 2014 to 2023. Initially, the frontal ice loss from 2014 to 2017 (3 years) was −6.95 × 10⁶ m w.e. a⁻¹, followed by −6.60 × 10⁶ m w.e. a⁻¹ during 2017-2020. Subsequently, there was an enhanced ice loss rate of −7.78 × 10⁶ m w.e. a⁻¹ during 2020–2023 (Fig. 6e). Frontal ice calving contributed to increased surface ice loss, affecting the estimation of glacier-wide surface mass balance. Overall, the response of the calving front to the annual glacier-wide mass balance of the Gepang Gath Glacier indicates an increase in ice loss by approximately 17–22% during the period from 2014 to 2023.

5.2 Significance of frontal calving on glacier ice loss and interdependency of glacier retreat and lake stability

An increase in proglacial lake volume confirms the retention of glacial meltwater in the glaciated basin. The rise in glacial lake volumes, in turn, significantly contributes to the acceleration of glacier melting through frontal calving and lake-glacier interface processes. Most of these lake-terminating glaciers in the Himalaya have shown an enhanced retreat rate in the last two decades (Patel and others, 2017; King and others, 2018; Sahu and Gupta, 2020). When glaciers terminate into the lake water, they lose glacierized area to the lake, leading to a more negative mass balance. As these glaciers experience accelerated melting and retreat due to climate warming, the ice mass loss from the total glacier volume causes the rapid formation of large proglacial lakes (Zhang and others, 2011; Furian and others, 2022) a phenomenon that is vividly reflected in the Gepang Gath glacier catchment. (Prakash and Nagarajan, 2018) inventoried 46 water bodies in the Chandra basin categorized as ice-dammed, bedrock-dammed, and moraine-dammed lakes. They identified seven lakes to be potentially dangerous sites in the Chandra Basin, with Gepang Gath and Samudra Tapu lakes classified as having high-outburst probability (Prakash and Nagarajan, 2017, 2018).

Our analysis of surface thickness change between 2013 and 2023 for the Gepang Gath glacier reveals significant changes, particularly in the lower part of the glacier. The surface changes of up to 90 m in 10 years indicate a notable total ice mass loss in the lower part of the glacier, coinciding with the expansion of a proglacial area. This clearly demonstrates a dynamic shift in the glacier’s morphology and volume change. Currently, the Gepang Gath lake encompasses about 1.21 km² of area from the origin of the proglacial lake phenomenon, traced back to 1962, when the lake area was approximately 0.2 km². This study clearly suggests that proglacial lakes significantly impact the overall ice mass loss. The presence of proglacial lakes at the glacier terminus also influences the glacier movement rates and its morphology, thus playing a significant role in overall glacier evolution (Benn, 2012; Pronk and others, 2021). Considering these dynamics and the glaciological aspects of losing ice, it is evident that the glacier ice mass loss over the Gepang Gath glacier will intensify in coming years, and the proglacial lake width may amplify in the coming future until the proglacial lake reaches the full longitudinal extent along the glacier terminus. This, along with the fact that it is categorized as a potentially hazardous lake, necessitates the establishment of an early warning system for Gepang Gath Lake.

6. Conclusions

Our comprehensive study on the lake terminating Gepang Gath glacier, incorporating field-based measurements of mass balance and debris thickness along with DGPS survey and geodetic estimates, provides a better understanding of the glacier response towards climate warming and its interactions with its proglacial lake. The analysis revealed that the proglacial lake associated with the glacier has expanded from 0.20 to 1.21 km² between 1962 and 2023, showing a nearly sixfold growth during this interval. This expansion is linearly correlated with time, showing accelerated expansion after 2010. The expanded lake area relates to glacier retreat, larger calving ice, and the frontal width of the glacier. The glacier exhibits predominantly negative glacier-wide mass balance values averaging at −0.90 ± 0.30 m w.e. a⁻¹ between 2014/15 and 2022/23. The ice loss at the frontal part due to calving adds approximately 17–22% increase in more negative mass balance. Moreover, the glacial lake expansion coincided with a total glacier retreat of 480 m, area loss of 0.42 km² and volume loss of −21.77 × 10 ⁶ m³ w.e in the terminus region of the glacier.

Our findings emphasize the significant contribution of the proglacial lake in intensifying the process of ice mass loss from the Gepang Gath glacier. This influence involves both surface mass balance alterations and the frontal ice loss attributed to calving. Notably, a visible disequilibrium between the annual balance and frontal ice loss becomes evident. At that same time, the retreat of the glacier also enabled the systematic expansion of its proglacial lake. All our findings consistently point towards a mutually dependent and compounding effect of the proglacial lake in exacerbating the overall ice mass loss faced by the Gepang Gath glacier. Considering these, it is expected that the Gepang Gath glacier will increasingly lose its mass in the coming years, making its proglacial lake highly vulnerable. Therefore, both the glacier and its proglacial lake need continuous and systematic monitoring at a high spatiotemporal scale and necessitates the establishment of an early warning system for potential cryospheric hazards like GLOF.