Dokriani Bamak Glacier (30°51′ N, 78°49′ E) is in the Garhwal range of the central Himalaya (Figure 1). It is a valley glacier measuring approximately 6 km in length and an area of 7.03 km², and elevation ranging from 4,050 to 6,632 m a.s.l. (Table 1) (Azam and Srivastava, 2020). The DBG has a north-west orientation and is guarded by three peaks: Jaonli (6,632 m a.s.l.) in the east, Draupadi Ka Danda I (5,716 m a.s.l.) in the south, and Draupadi Ka Danda II (5,670 m a.s.l.) in the west (Figure 1). The DBG tongue (4,050–4,900 m a.s.l.) is partially debris-covered (0.90 km², ~13% of DBG area) (Figure 1). The proglacial stream from DBG is called Din Gad which contributes to the Bhagirathi River of the Ganga River system. The DBG has extensively been investigated for its meteorological and mass balance conditions (Verma et al., 2018; Yadav et al., 2019, 2021; Azam and Srivastava, 2020; Dobhal et al., 2021; Garg et al., 2021).

Chhota Shigri Glacier (32°16′ N, 77°34′ E) is in the Lahaul-Spiti valley of the western Himalaya (Figure 1). This is a valley glacier with a length of ~9 km and an area of 15.5 km², and elevation ranging from 4,070 to 5,850 m a.s.l. (Table 1) (Azam et al., 2016). The CSG is a north-facing glacier, and its upper accumulation area is bounded by valley ridges, with Devachan Peak being the highest (6250 m a.s.l.) point. The terminus of CSG (<4,500 m a.s.l.) is covered with debris (~4% of CSG area) (Vincent et al., 2013). The CSG drains through a proglacial stream into the Chandra River, a tributary of the Indus River system (Figure 1). Since 2002, CSG is under continuous observations focusing on mass balances, SEB, ice thickness-volume-dynamics, and hydrology (Berthier et al., 2007; Wagnon et al., 2007; Soheb et al., 2017; Vashisht et al., 2017; Ramsankaran et al., 2018; Azam et al., 2019; Kumar et al., 2019; Mandal et al., 2020; Haq et al., 2021).

DBG and CSG have extensively been studied; hence different datasets are available from previous studies. On DBG Base Camp (BC, 3,774 m a.s.l.), an automatic weather station (AWS) logged the data over 2011–2016 while on CSG an AWS, mounted on a side moraine close to high camp (HC, 4,863 m a.s.l.), provided the meteorological data between 2009 and 2017. An automated precipitation gauge (Geonor T-200B) at CSG BC (3,850 m a.s.l.) provided the data since 2012. Supplementary Table S2 provides the logging details of meteorological data. The locations of AWSs and all-weather precipitation gauge are given in Figure 1. The measurement of B_a on DBG were conducted intermittently during 1992–2014 (−0.32 m w.e. a⁻¹; 1992–1995, 1996–2000, and 2007–2014) (Dobhal et al., 2021; Garg et al., 2021) while CSG represents the longest continuous B_a series since 2002 (−0.46 ± 0.40 m w.e. a⁻¹) in the Himalaya (Mandal et al., 2020). Altitudinal mass balances (b_a) are also available for 50-m bands over 2009–2013 for DBG (Pratap et al., 2015) and over 2002–2013 for CSG (Azam et al., 2016).

Daily reanalysis data from ERA5 was used to compute the surface energy fluxes and glacier-wide mass balances on DBG and CSG. ERA5 data is available since 1979 at 0.25° × 0.25° resolution [Copernicus Climate Change Service (C3S), 2017]. The ERA5 data was found to be readily accessible, consistent, and available over a long period, and it has already been used for mass- and energy-balance models in a few studies (Kumar et al., 2021; Patel et al., 2021). Daily incoming shortwave and net radiation (SWI and SWN), incoming longwave radiation (LWI), wind speed (WS), relative humidity (RH), air temperature (T_a), and precipitation (P) were downloaded for the nearest grids at DBG and CSG (Figure 1). The ERA5 raw data series for both the glaciers were bias-corrected using available in-situ meteorological data (Supplementary Table S2). The bias correction of T_a was done using a linear regression between mean monthly and daily raw ERA5 and in-situ T_a data for DBG and CSG, respectively. The bias correction of daily P, WS, RH, SWI, SWN, and LWI were performed using monthly factors derived from monthly in-situ and ERA5 raw data on DBG and CSG. LWI in-situ data were available from DBG; hence no bias correction was given. All the bias-corrected parameters showed a good coefficient of determination after the bias correction (R² > 0.90) (Supplementary Table S2). The details about the errors before and after bias correction and bias-correction factors are given in the Supplementary Material (Supplementary Figures S1–S13; Supplementary Table S3).

The mass- and energy-balance model (Figure 2) computes the SEB fluxes and b_a for each 50-m altitudinal range, and simulates snow accumulation, refreezing of rain/meltwater, surface melt, and sublimation/re-sublimation at daily time step using the long-term, bias-corrected daily ERA5 data between 1979 and 2020 (Section Study Area: Dokriani Bamak and Chhota Shigri Glaciers).

In the mass-, and energy-balance model, T_a is one of the most important parameters, as it decides the precipitation phase (snowfall or rain) (Hock, 2003; Shea et al., 2015a). In this study we calculated the extrapolated values of T_a using the temperature lapse rates (T_LR) developed using field observations (Azam et al., 2014a; Azam and Srivastava, 2020). Further, the T_p values were adopted from Jennings et al. (2018), where we used a T_p value of 0.7° and 1.1°C corresponding to 70–80% and 60–70% RH ranges for DBG and CSG, respectively, at which 90–100% precipitation was considered as snow.

The net solar radiation is also crucial in the surface mass, and energy-balance modeling, and the amount of insolation available for melt production largely depends on surface albedo (Azam et al., 2014b; Litt et al., 2019). Surface albedo values (α_s, α_i, and α_d) have high spatiotemporal variability over the glaciers. Deposition of dust and black carbon aerosols together with progressive snow metamorphism and compaction makes albedo values very uncertain (Oerlemans and Knap, 1998; Brock and Arnold, 2000). Moreover, the mass, and energy-balance models are highly sensitive to surface albedo (Kayastha et al., 1999; Acharya and Kayastha, 2019; Johnson and Rupper, 2020; Stigter et al., 2021); therefore, surface albedo values (α_s, α_i, α_d) are calibrated in the present study using the plausible ranges available from the study of Cuffey and Paterson (2010). Due to lack of information related to surface albedo evolution in the study area as well as to keep the model computationally simple, we have used static but separate calibrated albedo values for snow, ice, and debris surfaces, as also adopted in some other previous studies (Ragettli et al., 2013, 2015; Acharya and Kayastha, 2019). Energy, and mass-balance models are also sensitive to T_M, often unknown in the HK (Engelhardt et al., 2017; Azam et al., 2019; Azam and Srivastava, 2020). Further, the distribution of precipitation over glaciers is one of the biggest challenges in glaciological modeling, and it is spatially non-uniform in the HK region due to valley-specific precipitation gradients (P_G) (Maussion et al., 2014; Immerzeel et al., 2015; Sakai et al., 2015). Given that T_M, P_G, α_s, α_i, and α_d are highly sensitive parameters and least explored in the HK region; therefore, they are used for the calibration of the mass- and energy-balance model in this study.

For the calibration of the mass- and energy-balance model, Monte Carlo simulations are performed with 10,000 parameter sets where the parameters are varied over their plausible limits (Konz and Seibert, 2010; Rounce et al., 2020). The value, P_G is changed from 0 to 100% km⁻¹, T_M from −3°C to +3°C, α_s from 0.45 to 0.85, α_i from 0.35 to 0.55, and α_dfrom 0.1 to 0.2 following the study by Cuffey and Paterson (2010). The runs with minimum RMSE between modeled and in-situ B_a are selected for both the glaciers. The selected runs show an RMSE of 0.21 m w.e. a⁻¹ (1992–2014) and 0.32 m w.e. a⁻¹ (2002–2019) between modeled and in-situ B_a for DBG and CSG, respectively. The difference between the modeled and in-situ mean B_a are 0.01 and 0.06 m w.e. a⁻¹ on DBG and CSG, respectively (Table 2; Figure 3).

The model parameters are the key source of uncertainty in the modeled mass balances (Ragettli et al., 2013; Shea et al., 2015b). Parametric uncertainties are calculated by re-running the model while adjusting the parameters one-by-one within a reasonable range of their calibrated values and keeping the other model parameters unchanged (Table 2). The uncertainties in T_LR are taken as standard deviations of mean monthly values for both DBG and CSG. The uncertainties in other parameters (α_s, α_i, α_d, P_G, T_M and T_P) are unknown; hence these parameters are varied with the range of ±10% from their calibrated values (Anslow et al., 2008; Ragettli et al., 2013, 2015).

The total uncertainty in B_a is calculated by summing up all parametric uncertainties by applying the error propagation rule. The estimated mean uncertainties for B_a are 0.32 m w.e a⁻¹ and 0.38 m w.e a⁻¹ for DBG and CSG, respectively over 1979–2020. Among all the parameters, the highest uncertainty in B_a on both the glaciers is contributed by α_s (Table 2). The parametric uncertainty in the summer and winter mass balances are 0.38 and 0.01 m w.e a⁻¹ for DBG and 0.36 and 0.02 m w.e. a⁻¹ for CSG, respectively over 1979–2020.

In this study, a fixed glacier hypsometry was used on both the glaciers. For DBG, the hypsometry was manually delineated from high-resolution CNES-Airbus data from 10 July 2017 using a high-resolution Google Earth platform (Azam and Srivastava, 2020) while hypsometry for CSG was estimated using a digital elevation model (DEM) developed using Pléiades stereo pair from 18 August 2014 (Azam et al., 2016). This fixed area assumption introduces uncertainties to the modeled B_a, but these were found to be insignificant comparing the total estimated uncertainty in B_a in previous studies using T-index models (Azam et al., 2019; Azam and Srivastava, 2020), and hence they are ignored in this present study.

The HK Mountain range is located in the sub-tropical climate zone with a large annual temperature amplitude resulting in clear summer and winter seasons. The climate of HK is controlled by the Indian winter monsoon (IWM), embedding western disturbances, mainly during winters and the Indian summer monsoon (ISM), mainly during the summer (Gadgil et al., 2003; Dimri et al., 2016). The influence of IWM decreases eastwards; conversely, the intensity of the ISM decreases westwards along with the HK mountain range (Maussion et al., 2014). In-situ meteorological data from DBG and CSG are available for short periods (Section Available Field Data); therefore, long-term, bias-corrected ERA5 data over 1979–2020 are exploited to understand the mean seasonal characteristics on both the glaciers.

Mean monthly cycles of T_a and RH on both the glaciers followed roughly similar trends; however, the WS showed strong seasonality on DBG and moderate winds on CSG (Figure 5). The amplitudes of mean monthly T_a and RH (T_aDBG = −7.2°C and RH_DBG = 46% on DBG and T_aCSG = −6.1°C and RH_CSG = 47% on CSG) were sufficiently large to characterize the different seasons. A humid, warm, and less windy summer-monsoon from June through September and a less humid, cold, and windy winter season from December through March were demarcated on both the glaciers (Figure 5). A pre-monsoon over April–May and a post-monsoon over October–November were also defined (Supplementary Table S4). The same season demarcation was suggested on CSG by Azam et al. (2014a) using 3-years of AWS data (Figure 1).

The summer-monsoon was the warmest (T_aDBG = −1.5°C, T_aCSG = 1.8°C), least windy (WS_DG = 3.6 m s⁻¹, WS_CSG = 4.8 m s⁻¹), and the most humid (RH_DBG = 63%, RH_CSG = 61%) season of both the glaciers (Supplementary Table S4). Conversely, the winter season was the coldest, much below the freezing point (T_aDBG = −12.7°C, T_aCSG = −13.5°C), and windiest (WS_DBG = 8.8 m s⁻¹, WS_CSG = 5.7 m s⁻¹) on both the glaciers. Pre-monsoon and post-monsoon showed the moderate conditions for T_a, RH, and WS (Supplementary Table S4). Modeled surface temperature (T_s) was always negative on both the glaciers except for the summer-monsoon when it was close to 0°C, showing the melting at the surface due to higher T_a of the summer monsoon (Supplementary Table S4).

The mean monthly P cycles were remarkably different on both the glaciers (Figure 5). The ISM brought the major amount of annual P (74%) over DBG during the summer monsoon, while IWM brought the major amount of annual P (53%) over CSG during the winter. Therefore, these glaciers can be considered as summer accumulation-type and winter accumulation-type glaciers, respectively. The mean annual P on DBG was almost double that of CSG (Supplementary Table S4). Systematically, P amounts on DBG in all seasons were 1.5–2 times less than that of CSG except the summer monsoon when P on DBG was 9 times compared to CSG (Supplementary Table S4).

Despite the maximum solar angle in the summer monsoon, SWI was maximum during pre-monsoon on both the glaciers (SWI_DBG = 298 W m⁻², SWI_CSG = 418 W m⁻²) because monsoonal cloud cover impedes the SWI in the summer monsoon (Azam et al., 2014b; Litt et al., 2019). However, this effect was much stronger on DBG due to strong monsoonal influence (Figure 5). The reduced SWI (SWI_DBG = 294 W m⁻², SWI_CSG = 410 W m⁻²) during the summer monsoon was compensated by the highest LWI (LWI_DBG = 310 W m⁻², LWI_CSG = 272 W m⁻²) on both the glaciers—mostly emitted from warm, dense summer-monsoonal clouds, and surrounding valley walls. Post-monsoon and winter exhibited quite similar conditions, receiving lower SWI and LWI due to decreasing solar angle, T_a, and RH (Supplementary Table S4).

The modeled mass wastage was moderate and similar on both the glaciers with a mean wastage of −0.27 ± 0.32 m w.e. a⁻¹ (equivalent cumulative mass wastage of −11.16 ± 2.06 m w.e.) on DBG and −0.31 ± 0.38 m w.e. a⁻¹ (equivalent cumulative mass wastage of −12.61 ± 2.67 m w.e.) on CSG, over 1979–2020 (Figure 6). The years 1982/83 and 1988/89 showed the maximum B_a of 0.20 ± 0.33 m w.e. and 0.42 ± 0.23 m w.e., while the year 2000/01 showed the minimum B_a of −0.75 ± 0.32 m w.e. and −1.49 ± 0.76 m w.e. for DBG and CSG, respectively. The value of B_a was negative for 35 and 27 years and positive for 6 and 14 years on DBG and CSG, respectively. Though the mean mass wastage on both the glaciers was almost the same, the mass turnover on CSG (1.27 m w.e. a⁻¹) was higher than that of DBG (0.92 m w.e. a⁻¹) (Figure 6). This is due to the fact that CSG is a winter accumulation-type glacier and receives a lot of accumulation during winter and is melted out during the summer-monsoon, while for DBG, the accumulation and ablation seasons coincide during the summer monsoon (Section Meteorological Conditions and Seasonal Characteristics).

Modeled seasonal mass balances ranged from 0.07 to 0.69 m w.e. a⁻¹ and 0.17 to 0.72 m w.e. a⁻¹ for winter, and −1.00 to −0.16 and −1.69 to −0.06 w.e. a⁻¹ for summer on DBG and CSG, respectively. The mean summer and winter mass balances were −0.60 ± 0.30 and 0.32 ± 0.02 m w.e. a⁻¹ on DBG and −0.79 ± 0.36 and 0.48 ± 0.02 m w.e. a⁻¹ on CSG for the period 1979–2020, respectively.

Surface energy balance mainly depends on the seasons (Litt et al., 2019). In the summer-monsoon, SWN was the highest with mean values of 100 and 125 W m⁻² on DBG and CSG, respectively, with high daily variability from 71 to 151 W m⁻² and 76 to 171 W m⁻², respectively (Figure 7). The LWN was also maximum, with mean values of −4 and −41 W m⁻² on DBG and CSG, respectively, in the summer monsoon due to humid, warm, and dense cloud cover conditions that result in high values of LWI (Supplementary Tables S4, S5). Highest SWN and LWN resulted in the highest net radiation (R_n) at the surface during the summer monsoon with the mean seasonal values of 95 and 84 W m⁻² on DBG and CSG, respectively (Figure 7; Supplementary Table S5). Both the glaciers gained a small amount of energy through H. Conversely, a small amount of energy was released through LE—indicating some mass loss through sublimation (Supplementary Table S5). The resulting energy, Q, was the highest and positive during the summer monsoon on both the glaciers mainly due to the highest values of both SWN and LWN (Supplementary Table S5).

In winter, SWN was the least with 42 and 37 W m⁻² values while LWN remained moderate with −45 and −57 W m⁻² values on DBG and CSG, respectively (Supplementary Table S5). Winter SWN and LWN showed comparatively less variability on both the glaciers (Figure 7). In winter, DBG and CSG released some energy (−3 and −20 W m⁻², respectively) through R_n (Figure 7). Like CSG, a recent study also found negative R_n during winter on 8 glaciers in the Chandra valley, including CSG (Patel et al., 2021). Due to the higher temperature gradient and strongest WS in winter (Supplementary Table S4), H was the maximum and provided 13 and 23 W m⁻² energy at the surface of DBG and CSG, respectively (Supplementary Table S5). The LE was moderately negative during winter showing moderate glacier-wide sublimation (Supplementary Table S5). The net energy, Q, was also moderate but negative with −23 and −16 W m⁻² values, mainly due to the least values of winter SWN on both the glaciers (Supplementary Table S5).

In pre-monsoon and post-monsoon, the SWN values were moderate as 67 and 56 W m⁻² on DBG, and 87 and 61 W m⁻² on CSG, respectively (Supplementary Table S5). Energy loss through LWN was the highest in post-monsoon when compared to other seasons on both the glaciers (Figure 7). The value R_n was positive in pre-monsoon while negative in post-monsoon (Supplementary Table S5) due to the most negative LWN values in post-monsoon on both the glaciers. Both DBG and CSG gained more energy in the form of H in post-monsoon than in pre-monsoon (Figure 7). The most negative values of LE in pre-monsoon and post-monsoon (Supplementary Table S5) indicate the maximum mass loss through sublimation during these seasons. The Q value was the least in both pre-monsoon and post-monsoon due to higher negative values of LWN and LE, and moderate values of SWN (Supplementary Table S5).

Similar to previous studies on the HK region (Mölg et al., 2012; Azam et al., 2014b; Huintjes et al., 2015a,b; Johnson and Rupper, 2020; Patel et al., 2021), annual glacier-wide SWN contributed the maximum amount of energy to the total SEB on DBG and CSG. Further, both the glaciers lost energy through LWN at −35 and −57 W m⁻², respectively. The value of H brought some energy throughout the year with 9 and 15 W m⁻² values on DBG and CSG, respectively (Supplementary Table S5; Figure 7). Patel et al. (2021) also observed a similar mean annual value of H on 8 glaciers in the Chandra valley (Figure 7). LE remained negative throughout the year on both the glaciers indicating mass loss through sublimation, in line with other glacier-wide SEB studies in the HK region (Huintjes et al., 2015a,b; Patel et al., 2021). Annually, the highest SWN results in maximum R_n followed by H and LE on both the glaciers. Annual glacier-wide net energy, Q, was positive with a value of 10 W m⁻² on both DBG and CSG (Supplementary Table S5), indicating a net glacier-wide wastage (Section Glacier-Wide Annual and Seasonal Mass Balances).

The mean 50-m b_a varied from −4.62 to 2.20 m w.e. on DBG and −1.95 to 0.57 m w.e on CSG (Figure 8). On DBG, the terminus area less than 4,950 m a.s.l. showed less glacier wastage toward the valley walls compared to the middle of the glacier (Figure 8). This is due to the distribution of debris cover on DBG which is thicker toward valley walls (Pratap et al., 2015). Similarly, the CSG terminus area less than 4,400 m a.s.l. also showed less mass wastage (Figure 8) due to thick debris cover (Vincent et al., 2013). Despite the lowest albedo of debris cover which results in the highest SWN over the debris-covered glacier with mean values of 188 and 141 W m⁻² on DBG and CSG, respectively, the melting was the least due to a thick debris cover that protects the glacier from higher melt (Vincent et al., 2016; Banerjee, 2017). Going up on the glacier, the b_a increased and becomes positive in the accumulation areas on both the glaciers (Figure 8). The increase in b_a with altitude closely followed SWN which continuously reduced with altitude and achieved near-constant values of 56 and 59 W m⁻² at higher altitudes on DBG (>5,250 m a.s.l.) and CSG (>5,050 m a.s.l.) glaciers, respectively. This was probably due to the permanent snow cover in the accumulation area, having higher surface albedo values near stable SWN. Glaciers lost some energy through LWN which was highly negative over lower reaches (<4,500 m a.s.l.) compared to higher altitudes on both the glaciers (Figure 8). The value of R_n was higher over lower reaches (<4,850 m a.s.l.), and showed a reduction between 4,800 to 5,500 m a.s.l. and again increased slightly toward higher reaches (>5,500 m a.s.l.) due to the highest LWN at higher reaches on both the glaciers (Figure 8). The DBG showed high and positive H values at lower altitudes (<5,250 m a.s.l.) and slightly negative values at higher altitudes due to the negative air-surface temperature gradient (T_a – T_s) while it remained positive over the whole CSG (Figure 8). The LE, generally more negative at higher altitudes, showed altitudinal mass loss through sublimation equivalent to −6 to −42 W m⁻² on both DBG and CSG (Figure 8). The resulting energy, Q, was positive at lower altitudes (<5,000 m a.s.l.) and became negative at higher altitudes on both the glaciers (Figure 8).

On the DBG terminus area (<4,950 m a.s.l.), all the SEB components showed different behaviors toward valley walls as compared to the middle of the glacier due to the thick debris cover over the valley walls (Figure 8). Due to the low albedo of debris cover, SWN was the highest that resulted in maximum R_n and Q over those areas (Figure 8). H was slightly negative over debris-covered area compared to positive values at the middle of the glacier while LE was slightly less negative over the debris-covered area compared to more negative values at the middle of the glacier. This is due to the higher T_s than T_a over the debris-covered area due to the unavailability of snow or ice cover mainly during the summer-monsoon (Figure 8).

The mean glacier-wide sublimation was computed as −1.28 and −0.66 mm w.e. d⁻¹ over 1979–2020, with a strong spatial and temporal variability, on DBG and CSG, respectively (Figures 8, 9). The mean monthly sublimation was higher in May on both DBG (−1.93 mm w.e. d⁻¹) and CSG (−1.06 mm w.e. d⁻¹) glaciers, and sharply decreased over July–August as soon as monsoon arrived over these glaciers (Figure 9). Despite the lowest WS, the highest RH and T_a in the summer-monsoon months (Figure 9) reversed the specific humidity gradient. This reversal led to slightly positive values of LE, at least over the ablation area, in July–August indicating re-sublimation on both the glaciers (Section SEB in Ablation and Accumulation Zones). This provided the least glacier-wide values of LE and the least amounts of sublimation in the summer-monsoon on both the glaciers (Figure 9). The sign reversal of LE from the negative to positive values during humid and warmer conditions has also been observed from SEB studies on different mountain ranges, including the HK region (Wagnon et al., 1999, 2003; Oerlemans, 2000; Sicart et al., 2005; Azam et al., 2014b; Stigter et al., 2018; Litt et al., 2019). Even though the WS was the highest in winter, the mean monthly glacier-wide sublimations were moderate due to the lowest RH and T_a on both the glaciers, whereas sublimation was the highest during pre-monsoon and post-monsoon months corresponding to moderate WS, RH, and T_a (Figure 9). The computed glacier-wide sublimation losses account for significant amounts of 22 and 20% of total ablation on DBG and CSG, respectively; therefore, we stress that the inclusion of a simplified sublimation scheme in mass balance modeling using T-index models, have not yet included. Simplified scheme may parametrize the sublimation as a function of temperature, humidity, and wind speed (Azam et al., 2021).

Snow, clean ice, and debris-covered ice ablation contributed 56, 27, and 17% to the total ablation on DBG, while on CSG, these contributions were 65, 29, and 6%, respectively. In agreement to the highest glacier-wide snow ablation over both the glaciers, previous glacio-hydrological T-index modeling studies also suggested that the snowmelt contribution was the maximum on DBG and CSG (Engelhardt et al., 2017; Azam et al., 2019; Azam and Srivastava, 2020). A slightly higher percentage of snow ablation on CSG compared to DBG is probably due to the reason that it gets around 53% of its annual precipitation in winter months that melt out in the summer-monsoon months while DBG receives 74% of its annual precipitation in the summer-monsoon when T_a is the highest which might result in rainfall on DBG even up to 5,100 m a.s.l. (Pratap et al., 2015) (Section Meteorological Conditions and Seasonal Characteristics). Under similar mass wastage conditions (Section Glacier-Wide Annual and Seasonal Mass Balances), the percent contribution of melt from debris-covered ice on DBG was three times of CSG due to three-folds of debris-covered ice on DBG (~13%) compared to CSG (~4%) (Vincent et al., 2013; Pratap et al., 2015).

Most of the SEB studies in the HK have been performed at point-scale in the ablation zones (Azam et al., 2018; Litt et al., 2019). However, a few glacier-wide studies suggested that SEB is quite different in the ablation and accumulation zones of glaciers (Sun et al., 2014; Patel et al., 2021). To investigate the SEB in the ablation and accumulation zones of DBG and CSG, we estimated the mean annual ablation-wide and accumulation-wide SEBs separately, dividing the ablation and accumulation zones using the mean ELA from the literature (Table 1).

Mean monthly SWN budgets were similar over the ablation zones of both the glaciers: strong mean monthly cycles had the highest SWN in August during the summer monsoon and the lowest SWN in the winter, while relatively moderate values in the accumulation zone throughout the year except the winter when SWN were the lowest on both the glaciers (Supplementary Table S6; Figure 10). In winter, both the glaciers were completely covered by snow which resulted in higher surface albedo and similar ablation- and accumulation-wide SWN budgets (Supplementary Table S6). Negative LWN budgets showed higher loss of energy in the ablation zone compared to the accumulation zone throughout the year with most negative values in winter on both the glaciers (Supplementary Table S6; Figure 10), except slightly positive values in July–August on DBG most probably due to the highest LWI of heavy monsoonal cloud cover (Supplementary Table S4). Mean monthly H were positive as T_a was higher than T_s in the ablation zones of both the glaciers, while negative in the accumulation zones during the summer monsoon and pre-monsoon (Supplementary Table S6; Figure 10) as T_s becomes higher than T_a on both the glaciers (Supplementary Table S4). The LE was consistently negative in the accumulation zones of both the glaciers suggesting continuous mass loss through sublimation from higher altitudes; however it was slightly positive in the ablation zones over July–August on both the glaciers indicating re-sublimation during the core summer-monsoon (Figure 10). The resublimation was 3.50 and 2.04% on DBG and CSG respectively, compared to sublimation in the ablation zone. The value, Q remained negative throughout the year except for the summer monsoon in the ablation zones and July–August in the accumulation zones due to higher SWN on both the glaciers. Similar results have been discussed in other SEB studies on the Himalayan glaciers (Azam et al., 2014b; Litt et al., 2019; Patel et al., 2021).

Due to the scarcity of mass balance and meteorological data in the HK, the climatic drivers controlling the mass balances have been poorly discussed (Azam et al., 2014a, 2018; Shea et al., 2015a). To comprehend the major drivers controlling the glacier-wide seasonal and annual mass balances, the correlation coefficients (r) were developed amid annual and seasonal mass balances, bias-corrected mean annual ERA5 data, and surface energy fluxes over 1979–2020 on both the glaciers (Figure 11).

The value of B_a on DBG showed strong positive correlations (r = ~0.40–0.70) with P and surface albedo while moderately negative correlations (r = ~0.40–0.60) with SWN, R_n, and Q (Figure 11). Similarly, CSG also showed good correlations with P and surface albedo; however, the negative correlations with SWN, R_n and Q were stronger (r = ~0.80) (Figure 11). Due to their undersized role in total SEB (Section Seasonal and Annual Glacier-Wide Surface Energy Balance), H, LE, and LWN showed insignificant correlations with annual as well as seasonal mass balances on both the glaciers (Figure 11).

The value of B_a and summer mass balances on CSG showed moderate correlations with SWI while these correlations on DBG were insignificant (Figure 11). This is probably due to heavy monsoonal clouds that reduce the amount of SWI on DBG, resulting in low mean annual values (Supplementary Table S4). Winter mass balances on DBG showed weak positive correlation (r = 0.27) with P while a stronger positive correlation (r = 0.63) was observed on CSG (Figure 11). This is expected as DBG and CSG are summer and winter accumulation-type glaciers, respectively (Section Meteorological Conditions and Seasonal Characteristics). Further, summer mass balances on both the glaciers showed moderately positive but almost similar correlations (r = ~0.60) with P (Figure 11). Despite the fact that CSG receives its major annual precipitation during winter, almost similar correlation between summer mass balances and P occurs most probably due to sporadic summer-monsoonal snowfall events on CSG (Azam et al., 2019). A previous study on CSG investigated the critical role of summer-monsoon snowfalls in detail and concluded that these snowfalls often cover the whole or part of the ablation zone during peak melting months and abruptly reduce the SWI absorption and control the summer mass balances which further control the B_a (Azam et al., 2014b). The value of T_a showed weak and moderate negative correlations with B_a (r = −0.30) and summer (r = −0.52) mass balances, respectively on DBG, while showed moderate negative correlations with B_a (r = −0.44) and summer (r =-0.40) mass balances (Figure 11). The value of T_a was poorly correlated with winter mass balances on both the glaciers (Figure 11). However, as expected, T_a showed very strong correlations with T_s and LWO on both the glaciers (Figure 11).

The value of B_a on DBG and CSG showed moderately strong correlations with winter mass balances (r = 0.58 and r = 0.67, respectively) while showing very strong correlations with summer mass balances (r = 0.80, r = 0.97, respectively) (Figure 11). Higher dependency of B_a on summer mass balances suggests that both the glaciers have high vulnerability to regional warming; hence, expected to lose more mass in the continuation of warming (Banerjee and Azam, 2016; Kraaijenbrink et al., 2017; Krishnan et al., 2019; Mahto and Mishra, 2019).

Mass balance sensitivities were computed to understand the response of the glaciers to the changes in different model input parameters. These sensitivities were computed, one-by-one, by re-running the model with a unique set of each model parameter where h_H was the highest and h_L was the lowest value of parameter h, holding all the other parameters constant. Following Ragettli et al. (2013), the h_H and h_L were estimated by varying each parameter h by ±10% from its calibrated value except for T_m, T_p, T_a which were varied by 0.1, 0.1, and 1.0°C, respectively (Table 2). The mass balance sensitivities were estimated for the period 1979–2020 following (Oerlemans et al., 1998):

Where, B_a is the glacier-wide mass balance averaged over the period 1979–2020.

The estimated B_a sensitivities on DBG and CSG are given in Table 2. The B_a was the most sensitive to α_s, with the sensitivities of 0.29 and 0.37 m w.e. a⁻¹ on DBG and CSG, respectively (Table 2). Previous studies on other glaciers in the Alps and Himalaya also showed the maximum sensitivity of B_ato α_s (Klok and Oerlemans, 2004; Johnson and Rupper, 2020; Stigter et al., 2021). The modeled B_a showed moderate sensitivities to T_LR (DBG = 0.10 m w.e. a⁻¹; CSG = 0.12 m w.e. a⁻¹). Sensitivities were quite low to T_p, T_M, α_d, α_i and P_G for both the glaciers (Table 2).

The sensitivity of the modeled mean B_a to 1°C change in T_a was higher on DBG (−0.50 m w.e. a⁻¹) than CSG (−0.30 m w.e. a⁻¹) whereas the sensitivities to 10% change in P were roughly the same (DBG = 0.23 m w.e. a⁻¹, CSG = 0.13 m w.e. a⁻¹) (Table 2). Higher sensitivity to T_a on DBG is probably due to different precipitation regimes on both the glaciers. DBG receives its maximum of annual precipitation in the summer monsoon when T_a is the highest; hence, DBG had more sensitivity to T_a compared to CSG that receives its major precipitation in winters (Fujita, 2008; Azam et al., 2014a). Using a T-index model, a previous study on CSG computed higher sensitivity (−0.52 m w.e. a⁻¹) of mass balance to 1°C change in T_a and roughly similar sensitivity (0.16 m w.e. a⁻¹) to 10% change in P (Azam et al., 2014a). Another study on Zhadang Glacier in Tibet showed similar results using an energy balance model with the sensitivity of −0.47 m w.e. a⁻¹ to 1°C change in T_a and sensitivity of 0.14 m w.e. a⁻¹ to 10% change in P (Mölg et al., 2012). Our sensitivity results are quite comparable with these studies in the mountain glaciers.

In this section, we discuss the sublimation rates from different studies on HK glaciers. However, irrespective of our glacier-wide and round-the-year study, often the studies were (i) available at point-scale, (ii) from different months of the year, (iii) having different locations of automatic weather stations (on/off glacier), and (iv) installed on different surfaces (snow and ice) that hinders a direct comparison. In the present study, the mean glacier-wide sublimation was computed as −1.28 and −0.66 mm w.e. d⁻¹ over 1979–2020 on DBG and CSG, respectively. Previously, using in-situ AWS data from the middle of the ablation zone (4,670 m a.s.l.) on CSG, a point-scale SEB study computed a mean sublimation of −0.63 mm w.e. d⁻¹ over 2012–2013 (Azam et al., 2014a). Another recent point-scale SEB study on seasonal snow surface on a lateral moraine of CSG only (4,863 m a.s.l.) suggested that sublimation was −1.1 mm d⁻¹ during winters over 2009–2020 (Mandal et al., 2022). The mean glacier-wide sublimation was reported as −1.08 and −0.70 mm w.e. d⁻¹ on Zhadang Glacier (south-central Tibetan Plateau) and on Puruogangri ice (north-central Tibetan Plateau) between 2001 and 2011 (Huintjes et al., 2015a,b). In the central Himalaya, the sublimation rate was reported for short-term snow-cover at the Pindari Glacier AWS site (off-glacier; 3750 m a.s.l.) to be −0.3 mm d⁻¹ during winters (Singh et al., 2020). In the Nepal Himalaya, on Yala Glacier (5350 m a.s.l.) around −1.00 mm d⁻¹ mass was lost through sublimation during winters (Stigter et al., 2018). In line, another study showed significant sublimation rate of −7.1 and −1.9 mm d⁻¹ on Mera Glacier and −2.4 and −1.8 mm d⁻¹ on Yala Glacier during the post- and pre-monsoon season, respectively (Litt et al., 2019).

The in-situ glaciological method showed a mass wastage of −0.32 m w.e. a⁻¹ on DBG over 1992–1995, 1997–2000, and 2007–2013 (Dobhal et al., 2021), while our model showed slightly lesser mass wastage of −0.27 ± 0.31 m w.e. a⁻¹ over the same years (Figure 12A). On CSG, the glaciological method showed a mass wastage of −0.43 ± 0.40 m w.e. a⁻¹ over 2002–2019 (Mandal et al., 2020), while our model showed similar mass wastage of −0.42 ± 0.42 m w.e. a⁻¹ (Figure 12B).

Geodetic mass balances are also available on both the glaciers and are used here to validate the modeled mass-balance series. A recent study using high-resolution Cartosat-1 DEM and SRTM DEM estimated a mass wastage of −0.23 ± 0.10 m w.e. a⁻¹ on DBG over 1999–2014 (Garg et al., 2021). Over the same period, our model showed slightly higher mass wastage with a mean mass balance of −0.35 ± 0.41 m w.e. a⁻¹ on DBG. On CSG, our model showed higher mass wastage of −0.39 ± 0.43 m w.e. a⁻¹ against a mass wastage of −0.27 ± 0.13 m w.e. a⁻¹ over 2005–2014 derived using ASTER DEMs (Brun et al., 2017). Another geodetic study provided a mean mass wastage of −0.46 ± 0.34 m w.e. a⁻¹ over 2000–2012 (Vijay and Braun, 2016), while our model computed a mean mass wastage of −0.67 ± 0.54 m w.e. a⁻¹ over the same period. Another geodetic study combining SRTM and SPOT5 DEMs provided a mean mass wastage of −1.12 m w.e. a⁻¹ using density assumption 1 and −1.02 m w.e. a⁻¹ using density assumption 2 which is in close agreement with our model value (−1.01 m w.e. a⁻¹) over 1999–2004 (Berthier et al., 2007) (Figures 12A,B).

Our modeled mean mass balance of −0.27 ± 0.39 m w.e. a⁻¹ on DBG and −0.31 ± 0.38 m w.e. a⁻¹ on CSG were in better agreement with the modeled mean mass balance of −0.25 ± 0.37 m w.e. a⁻¹ and −0.26 ± 0.29 m w.e. a⁻¹ from a simplified T-index model over 1979–2020 (Srivastava et al., 2021). The present study, showing a moderate mean mass balance on both the glaciers over 1979–2020, is in close agreement with most of the previous studies and nicely captures the inter-annual variability with other modeled MBs (Figures 12A,B).