Reanalysis of the longest mass balance series in Himalaya using a nonlinear model: Chhota Shigri Glacier (India)

The glacier-wide mass balance (MB) series on Chhota Shigri Glacier has been reanalysed by combining the traditional MB reanalysis framework and a nonlinear MB model. The nonlinear model is preferred over the traditional glaciological method to compute the glacier-wide MBs, as the former can capture the spatiotemporal variability in point MBs from a heterogeneous in situ point MB network. Further, the nonlinear model is also used to detect erroneous measurements from the point MB observations over 2002–2023. ASTER and Pléiades stereo imagery show limited areal changes but negative mass balances of −0.38 ± 0.05 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ during 2003–2014 and −0.51 ± 0.06 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ during 2014–2020. The nonlinear model outperforms the traditional glaciological method and agrees better with these geodetic estimates. The reanalysed mean glacier-wide MB over 2002–2023 is −0.47 ± 0.19 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$, equivalent to a cumulative loss of −9.81 ± 0.87 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ Our analysis suggests that the nonlinear model can also be used to complete the MB series if for some years the field observations are poor or unavailable. With this analysis, we revisit the glacier-wide MB series of Chhota Shigri Glacier and provide the most accurate and up-to-date version of this series, the longest continuous ever recorded in the Himalaya. We recommend applying the nonlinear model on all traditional glaciological mass balance series worldwide whenever data are sufficient, especially in the Himalaya, where in situ data are often missing due to access issues.

1 Introduction

Glaciers are excellent indicators of changing climate; therefore, long-term glacier mass changes are observed to understand the impacts of climate change (Oerlemans, 2001; Zemp et al., 2019). Glacier monitoring is also essential to understand possible glacial hazards (Harrison et al., 2018; Shukla et al., 2018; Shugar et al., 2021; Gantayat and Ramsankaran, 2023), regional hydrology (Azam et al., 2021; Yao et al., 2022; Nepal et al., 2023), and sea level rise (Gardner et al., 2013; Rounce et al., 2023). The glacier mass balance (MB) can be estimated from satellite data, through modelling approaches, or measured using the field-based traditional glaciological method (Cogley, 2009; Zemp et al., 2015; Kumar et al., 2018; Miles et al., 2021; Berthier et al., 2023).

Over the last decade, rapid development has been made through satellite geodetic MB estimates covering almost all glacierized areas in the Himalaya (Brun et al., 2017; Bolch et al., 2019; Shean et al., 2020; Hugonnet et al., 2021; Jackson et al., 2023). These geodetic estimates are primarily available at a multiannual scale and thus cannot be used to study the inter-annual variability in glacier MB. Conversely, field measurements using standard methods (Østrem and Stanley, 1969) yield data on the seasonal/annual response of glacier MB to local meteorological conditions (Zemp et al., 2015). Field MB observations remain scarce in the Himalaya compared to the other mountain ranges (Azam et al., 2018) and have been limited to only 35 glaciers (Vishwakarma et al., 2022). Most observations are available from easily accessible and small glaciers for short periods that are generally shorter than 10–15 years. The ongoing MB series include Chhota Shigri, Hōksar, Kolahoi, and Sutri Dhaka glaciers in the western Himalaya (Oulkar et al., 2022; Mandal et al., 2020; Romshoo et al., 2022; 2023); Mera, Pokalde, Rikha Samba, Trambau, West Changri Nup, and Yala glaciers in the central Himalaya (Sunako et al., 2019; Wagnon et al., 2021; Stumm et al., 2021); and Ganju La and Thana glaciers in the eastern Himalaya (Tshering and Fujita, 2016).

For annual glacier-wide MB estimation, the traditional glaciological method has been used in the Himalaya (Azam et al., 2018). This method involves interpolation/extrapolation of point MB measurements from fixed locations to the whole glacier area by applying different approaches, including contouring, profiling, and kriging (Østrem and Brugman, 1991; Zemp et al., 2013), or by the application of observed MB gradients to the glacier hypsometry (Funk et al., 1997; Wagnon et al., 2021). The selected point measurement sites may not be representative of the surrounding areas because (1) ablation stakes are often inserted away from the steep slopes towards the valley walls for safety reasons, and thus, snow avalanche inputs onto valley glaciers are not included; (2) crevassed areas are not sampled; (3) snow accumulation is site-specific and largely depends on local topography that controls snow drift and deposition; and (4) harsh weather sometimes restricts access to accumulation measurement sites. Almost all the MB series are victims of one issue or other such issues; therefore, the estimated glacier-wide MBs often carry systematic biases (Thibert et al., 2008). These biases can be corrected by calibrating the MB series using satellite-derived geodetic mass estimates generally over 5–10 years (Zemp et al., 2013; Wagnon et al., 2021).

Furthermore, it is practically difficult to keep the position fixed for point measurements due to accessibility issues, stake displacement due to glacier dynamics, use of different surveying equipment (GPS, dGPS, total station, theodolite, etc.), and different researchers' involvement for decades of monitoring. Hence, the measurement network differs in space and time. In this situation, heterogeneous in situ measurements do not always allow us to catch the large spatiotemporal variability in point MBs within the same elevation range (Funk et al., 1997; Vincent and Six, 2013); consequently, the point MB–elevation relationship is insufficient to investigate the changes in glacier-wide MBs (Kuhn, 1984; Huss and Bauder, 2009; Thibert et al., 2013).

To include the spatiotemporal variability in point MB measurements, Lliboutry (1974) proposed a linear statistical model and tested it over the small ablation area of Saint-Sorlin Glacier (France), assuming similar temporal changes in the MB over the whole area. Vincent et al. (2018) suggested that the linear model of Lliboutry (1974) was valid over a limited elevation range but ignored the decreasing spatiotemporal variability in point MBs with elevation (Oerlemans, 2001). To address this issue, they proposed a nonlinear model that considers the decreasing spatiotemporal changes in point MBs over the large elevation range and successfully tested their model on four different glaciers from different climate regimes, including Chhota Shigri Glacier (India).

The MBs on Chhota Shigri Glacier were estimated using the nonlinear model over 2002–2016 and then calibrated using geodetic MB over 2005–2014 (Vincent et al., 2018). In the present study, we extended the MB series on Chhota Shigri Glacier up to 2023 using the traditional method, estimated the areal changes and geodetic MBs over the 2003–2014 and 2014–2020 periods, estimated the debris cover as of September 2020, and reanalysed the annual MB series since 2002 using a novel reanalysis framework that combines the Vincent et al. (2018) nonlinear model and the reanalysis framework proposed by Zemp et al. (2013). Additionally, we assessed areal changes and geodetic MBs of neighbouring glaciers Hamtah and Sichum over the same periods based on available satellite stereo images.

Since 2002, the MB series of Chhota Shigri Glacier has been continuously monitored, making it the longest series in the Himalaya. Azam (2021) highlighted the importance of Chhota Shigri as a reference glacier for large-scale MB and hydrological studies; therefore, the main aim of the present study is to produce the most accurate glacier-wide MB series in this region. First, the nonlinear model of Vincent et al. (2018) was used to detect the erroneous point MB measurements in the data set. Second, the nonlinear model was applied using the observed point MBs to estimate the glacier-wide MB at an annual scale. Third, homogenization of the glacier-wide MB series accounting for glacier areal changes was performed, and fourth, the glacier-wide MB series was calibrated using geodetic MBs. Additionally, we compared the performance of the nonlinear model with the traditional method for estimating glacier-wide MB. We also assessed the nonlinear model's ability to estimate glacier-wide MB using end-of-season snow line data when field measurements were unavailable in a particular year.

2 Study area

Chhota Shigri Glacier (32.28° N, 77.58° E) is in the Chandra River basin, a tributary of Upper Indus Basin, Lahaul–Spiti valley of the western Himalaya (Fig. 1). Chhota Shigri flows from 5830 to 4100 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, with a length of ∼ 9 km and an area of 15.47 km² (in 2020). Based on the most updated map obtained in September 2020, 12 % of its total surface area is covered with debris between the snout and 4500 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, including medial and lateral moraines from 4100 to ∼ 4900 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ and a debris-covered eastern tributary glacier (Fig. 1). Debris thickness ranges from less than a few centimetres of thin debris to a few metres of boulders. Valley walls bound its accumulation area with the highest Devachan peak (6250 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). The accumulation area has two east- and west-oriented tributaries that feed to the main ablation area (< 5070 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), having a north aspect and divided into two parallel flows by a medial moraine.

Figure 1 Chhota Shigri Glacier, showing the location of ablation and accumulation point measurement sites. Orange strips show the debris-covered glacier area. The background image is a Pléiades satellite image taken on 12 September 2020 (© CNES 2020; Distribution Airbus Defence and Space). The glacier extent corresponds to 12 September 2020. The coordinates are in the Universal Transverse Mercator (UTM) north zone 43.

Chhota Shigri is a well-studied glacier for various aspects, including traditional MBs, energy balance, dynamics, ice thickness, and hydrology (Wagnon et al., 2007; Azam et al., 2012; Ramsankaran et al., 2018; Haq et al., 2021; Srivastava and Azam, 2022a; Mandal et al., 2020, 2022). Several studies have also observed its geodetic MBs (Berthier et al., 2007; Vincent et al., 2013; Brun et al., 2017; Mukherjee et al., 2018). Long-term annual MBs have been reconstructed over 1950–2020, applying a temperature index model (Srivastava et al., 2022), and over 1979–2020, using an energy balance model (Srivastava and Azam, 2022b). Due to recent glacier wastage on Chhota Shigri Glacier, the western tributary (WT) glacier got disconnected in the summer of 2012 (Srivastava et al., 2022). The fragmented tributary is now clearly visible in the high-resolution Pléiades image from 12 September 2020 (Fig. 1).

In this study, we focus on Chhota Shigri Glacier, but the available satellite stereo images also cover neighbouring Hamtah and Sichum glaciers; therefore, we also estimated the areal changes and geodetic MBs for these two glaciers (Sects. 3.4 and 3.5). Hamtah Glacier has been studied for its MBs and avalanche contribution (Vincent et al., 2013; Laha et al., 2017). Furthermore, for all three glaciers, we also delineated the debris cover corresponding to 2020 (Table 1).

3 Methods

3.1 Traditional mass balance method

Glacier-wide annual MBs (B_a) have been estimated using a network of 22–25 ablation bamboo stakes (inserted up to 10 m into the glacier) distributed over 4300–4900 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ along the main axis of the glacier (Fig. 1) and 4–6 accumulation pits/cores over 5160–5550 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ distributed over the eastern and western tributaries of the glacier (Wagnon et al., 2007). The traditional glaciological profile method was used to estimate the glacier-wide MB from the observed point MBs (Østrem and Stanley, 1969). First, using the observed point MBs, the mean altitudinal MBs were estimated for each 50 m elevation band from available point MBs within each elevation band (Fig. 1). In case no measurements were available (due to loss of stakes or missing accumulation measurements) the MBs were estimated using linear interpolation/extrapolation of neighbouring bands. Second, the B_a (in $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$) was estimated as follows:

$$\begin{array}{}\text{(1)}& B_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{S}}{\sum }_{z=\text{min}}^{z=\text{max}}{b}_z{s}_z,\end{array}$$

where b_z is the mean altitudinal MB (in $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$) of a given elevation band, z, for area s_z (m²); and S is the total glacier area (m²). In the ablation area, emergence changes at each ablation stake were converted to the point MB using a fixed density of 900 kg m⁻³ for ice and 350 kg m⁻³ for snow (Wagnon et al., 2007; Cogley et al., 2011), while in the accumulation area, the varying snow, firn, or ice densities (350–900 kg m⁻³) were measured in the field. The hydrological year for MB calculations is defined from 1 October to 30 September of the following year; however, the exact measurement dates on site varied from a couple of days to a week. Following Thibert et al. (2008), an overall uncertainty of ± 0.40 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ for glacier-wide MB was estimated by incorporating the errors in point measurements and their distribution over the glacier (Azam et al., 2012).

Due to access difficulties, such as snowstorms on 22–24 September 2018, or logistical or budget issues, a limited number of point MB measurements could be carried out in some years. This was the case for October 2015, when only two accumulation measurements could be performed, or 2018, when measurements were done early in the season before the storm. For those 2 years, point MB data in the accumulation zone, where no measurements had been taken, were estimated using previous years with a similar ablation pattern (Mandal et al., 2020). In 2020, only two in situ point MB data were available, preventing the traditional method from being applied. Furthermore, no measurements could be performed in 2021; hence, no MB could be estimated. Table S1 in the Supplement provides all information about the point MBs and field expeditions since 2002.

3.2 Nonlinear mass balance model

The nonlinear MB model suggests that the observed point MB, b_i,t, at any site i for year t can be decomposed into (1) a spatial effect term, α_i, and (2) a temporal term, β_t, combined with a spatial effect, γ_i, and can be written as (Vincent et al., 2018)

$$\begin{array}{}\text{(2)}& b_{i,t}={\mathit{\alpha }}_i+{\mathit{\beta }}_t{\mathit{\gamma }}_i+{\mathit{\epsilon }}_{i,t},\end{array}$$

where α_i, the spatial effect at location i, is the average point MB at the site over the whole study period; β_t is the annual deviation from the average point MB (thus Σβ_t = 0); and ${\mathit{\gamma }}_i={\mathit{\sigma }}_i/{\mathit{\sigma }}_{\text{max}}$ is a scaling factor defined as the ratio of the standard deviation (SD) of annual MB at site i and the maximum standard deviation (σ_max) observed from the point MB measurements over a long period. The ε_i,t term represents residuals resulting from measurement errors and inconsistencies between the model and observed data. The spatiotemporal decomposition proposed in Eq. (2) assumes that β_t is the same at all point locations for any given year (t) and thus has a glacier-wide significance, while γ_i term accounts for nonlinear effects with elevation (Vincent et al., 2018).

To compute the scaling factor, γ_i, on Chhota Shigri Glacier, standard deviations were computed from the point MBs available for each 50 m elevation band as the point MBs are not available each year from the same fixed locations (Fig. 2). The standard deviations were computed only for 50 m elevation bands, where mean annual MBs were available from in situ measurements over a minimum of 10 years, and it was assumed that the computed standard deviations are representative of the whole period of investigation (2002–2023). This resulted in 16 standard deviation values over the whole glacier, with a maximum standard deviation of 1.17 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ at 4525 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (4500–4550 band) and a minimum standard deviation of 0.40 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ at 5325 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (5300–5350 band). The decreasing magnitude of the standard deviation with elevation indicates the decreasing sensitivity of the annual MB to temperature and precipitation (Fig. 2), as already suggested by several studies on glaciers worldwide (Kuhn, 1984; Soruco et al., 2009; Basantes-Serrano et al., 2016; Vincent et al., 2018; Wagnon et al., 2021). The measurements are poor in the accumulation area, and no measurement was available above 5325 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$; therefore, after some trials, we adjusted the standard deviation at 6000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ to be zero (Fig. 2). A decreasing trend in standard deviation values below 4525 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 2) is due to the presence of debris cover over the tongue of Chhota Shigri Glacier (Fig. 1) that undermines the glacier's sensitivity to climate (Vincent et al., 2013; Banerjee and Shankar, 2013). The scaling factor, γ_i, at each point MB location, was computed from the second-degree polynomial function fitted over the standard deviation vs. elevation scatter plot (Fig. 2).

Figure 2 Standard deviations of the annual MBs versus elevation. The black line corresponds to a polynomial fit (degree of freedom = 2). The standard deviations were estimated for those 50 m elevation bands, where a minimum of 10 years of point measurements were available at each site, and it is assumed to be zero at 6000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (above the glacier top at 5830 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$).

The nonlinear model was run at 200 m × 200 m spatial resolution over 2002–2023 using all available point MBs (413 point measurements, excluding the erroneous measurements; Sect. 3.3) and polynomial equation (Fig. 2; details can be found in the supplement of Vincent et al., 2018). The MB is assumed to be spatially constant over each 200 m × 200 m grid for a given year. If there is more than one observation in a grid in a given year, then the mean MB of the available observations was used for MB computation. The grid size is a compromise between the spatial variability and the density of the available point measurements.

Figure 3 Distribution of all 413 point MB measurements (yellow dots) available over 2002–2023 on Chhota Shigri Glacier. The grids (in light blue) show a spatial resolution of 200 m × 200 m of the nonlinear model. For 2020/2021, no field measurement was conducted; hence, two grids (shown with green colour outline) corresponding to zero MB were selected on the delineated SLA to run the model. The background is a Sentinel image from 6 September 2021, which is used to delineate the SLA.

Figure 4 (a) The residuals between measured and modelled point MBs from the nonlinear model using all available 423 point MBs as a function of distance from the glacier snout for each hydrological year between 2002 and 2023. The shaded dark blue and light blue envelopes represent the 1 and 2 SD values, respectively. (b) The probability density function (normal distribution curve) of all point MB residuals between 2002 and 2023.

Field measurements were unavailable in the 2020/2021 year (Sect. 3.1); hence, the nonlinear model cannot be run for this hydrological year. To run the model, at least one point MB measurement is required each year (Vincent et al., 2018). We assumed the snow line altitude (SLA) at the end of the ablation season to be equivalent to the equilibrium line altitude (ELA) (Rabatel et al., 2005; Brun et al., 2015; Davaze et al., 2020; Barandun et al., 2021). The SLA was delineated on a 6 September 2021 Sentinel image, and zero MBs (MB at ELA = 0 $\mathrm{m}\mathrm{w}.\mathrm{e}.$) were assumed for two 200 m × 200 m grids for which MB observations were available from other years (Fig. 3). It should be noted that there was no other cloud-free image from September 2021. The MB estimation from SLA using the nonlinear model is discussed in detail in Sect. 5.3.

The model output provides the mean α_i and mean γ_i for each point location over 2002–2023 and β_t for each year (Eq. 2). A total of 54 values for α_i and γ_i and 21 values for β_t (corresponding to each hydrological year) were computed (Tables S2 and S3 in the Supplement). For the calculation of glacier-wide MB, a spatial distribution of α_i over the whole surface area of the glacier is needed. First, for each 50 m elevation range (e), mean α_e was estimated from all available α_i, by taking a simple arithmetic mean, and γ_e from all available γ_i from respective elevation bands (Eq. 2). The modelled point MBs were available over the 4355–5512 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ elevation range, and beyond this range, the mean α_e and γ_e from the lowest (4300–4350 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) and highest (5500–5550 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) ranges were used to cover the lowest (0.15 km²; 0.97 % of total area) and highest (0.68 km²; 4.40 % of total area) parts of the glacier. Second, applying α_e, γ_e, and β_t from all elevation bands in Eq. (1), along with corresponding elevation areas, the annual glacier-wide MBs over 2002–2023 were estimated.

3.3 Tracking the erroneous in situ point mass balances

The nonlinear model computes the residuals (difference between the measured and theoretical values) of each measured point MB and can detect errors in in situ point MB data (Vincent et al., 2018). The distribution of residuals over the glacier as a function of distance from the snout showed no spatiotemporal pattern (Fig. 4a), indicating that the nonlinear model does not provide any apparent bias for any specific year. As expected, the residuals followed a normal distribution with a standard deviation (SD) of 0.35 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ (Fig. 4b). To detect the measurement errors in the point MBs in the Chhota Shigri measurement network over 2002–2023, we assumed all the point MBs have residuals > 2 SD (0.70 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$) to be suspicious. Of 423 point MB measurements, 15 such point MBs were found and investigated further. Five point MBs had been wrongly reported from the notebooks and thus have been corrected. We could not find any reason for the rest of the suspicious points. Therefore, they have been considered erroneous and discarded in the final model run. The erroneous data were collected in different years (five ablation point measurements from 2009, 2012, 2018, and 2022 and five accumulation point measurements from 2009, 2011, 2014, and 2022) (Fig. 4). The standard deviation of the residuals from the nonlinear model was reduced from 0.35 to 0.30 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ after correction/removal of suspicious point MB measurements.

3.4 Areal changes and debris cover estimation

The areal changes and debris cover were estimated on Chhota Shigri, Sichum, and Hamtah glaciers by manual delineation following the Global Land Ice Measurements from Space (GLIMS) guidelines from the available ASTER (8 October 2003) and Pléiades images (26 September 2014 and 12 September 2020) (Raup et al., 2007). We have preferred manual delineation as it was considered the most accurate method for delineating glacier outlines (Stokes et al., 2007; Garg et al., 2017; Shukla and Qadir, 2016). The ice divides were interpreted using the Pléiades digital elevation model (DEM). The changes were estimated for the ablation area for 2014 and 2020 as the changes in the accumulation area were insignificant. The generated glacier outlines (2003, 2014, and 2020) were used to estimate the glacier area changes during 2003–2020. The uncertainties associated with the glacier area were calculated using the buffer method (Bolch et al., 2010; Chand and Sharma, 2015). The buffer size was half the pixel value (Bolch et al., 2010; Andreassen et al., 2022).

3.5 Geodetic mass balances

The geodetic MBs were estimated over two periods (2003–2014 and 2014–2020) for Chhota Shigri, Sichum, and Hamtah glaciers using satellite stereo images from ASTER (15 m resolution) acquired on 8 October 2003 and Pléiades (0.70 m resolution) acquired on 26 September 2014 and 12 September 2020, respectively. The ASTER October 2003 stereo pair was preferred to other ASTER or SPOT5 stereo pairs acquired in late summer 2002, 2004, and 2005 because it resulted in the smallest uncertainties. The stereo images were acquired close to the end of the hydrological year, reducing the impact of any seasonal offset. The DEM generation, co-registration, and MB calculation procedure is the same as in Falaschi et al. (2023). Uncertainties for the glacier-wide geodetic MB were estimated using the patch method (Wagnon et al., 2021). This method aims to empirically determine the uncertainty associated with the mean elevation change by sampling patches of stable terrain of various sizes to measure the decay of the error with the averaging area.

Geodetic MBs were estimated over 10.97 years (from 8 October 2003 to 26 September 2014) and 5.96 years (from 26 September 2014 to 12 September 2020) and linearly scaled to estimate the geodetic MBs over 11- and 6-year periods, respectively, to make a direct comparison with the in situ MBs (estimated from end of September to the end of September in the next year). Furthermore, the geodetic MBs included both the WT glacier, which fragmented around 2012 (Srivastava et al., 2022), and the main Chhota Shigri (area-weighted) (Table 1) for a direct comparison with the traditional and nonlinear MBs that include the WT glacier.

Figure 5 Conceptual diagram of the overall methodology: homogenization, uncertainty estimation, validation, and calibration steps.

3.6 Homogenization of glacier-wide mass balances

In initial studies (Wagnon et al., 2007; Azam et al., 2012), a fixed hypsometry (glacier area and elevation) from SPOT5 2005 DEM was used, while in follow-up studies (Azam et al., 2014; Mandal et al., 2020) a fixed hypsometry from the Pléiades August 2014 DEM was used to estimate the traditional MBs on Chhota Shigri Glacier. These fixed hypsometries insert bias in the MB series (Cogley et al., 2011; Zemp et al., 2013). Here, the Chhota Shigri Glacier annual MBs (from the traditional method and nonlinear model) are homogenized with the linearly changing annual hypsometries from ASTER and Pléiades DEMs over 2003–2014 and Pléiades DEMs over 2014–2020 (Sect. 4.1). We adopted the approach suggested by Zemp et al. (2013) that assumes a linear area change over a record period (N years) and estimates the area (s) of an elevation band (e) for each year (t) as follows:

$$\begin{array}{}\text{(3)}& s_{\mathrm{e},t}=s_{\mathrm{e},\mathrm{0}}+{\displaystyle \frac{t}{N}}\cdot (s_{\mathrm{e},N}-s_{\mathrm{e},\mathrm{0}}),\end{array}$$

where s_e,0 and s_e,N are the elevation bin areas from the first and the second geodetic survey, respectively, and the time t is zero in the year of the first survey. The homogenization process of both traditional and nonlinear MB series changed the annual glacier-wide MBs at most by 0.02 $\mathrm{m}\mathrm{w}.\mathrm{e}.$, reflecting the negligible impact of areal changes over the 2003–2020 period on Chhota Shigri Glacier (Sect. 4.1). Post-2020, the hypsometry of the 2020 year was used to estimate the MBs until 2023. Figure 5 summarizes the overall methodology step by step, including homogenization, validation/calibration, and error estimation (Sect. 3.7 and 3.9).

3.7 Validation and calibration of glacier-wide mass balances

Previously, we validated the traditional MBs with geodetic MB available over 2005–2014 (Azam et al., 2016). The systematic biases were within the uncertainty ranges of traditional and geodetic MBs; hence, no calibration was done. In this study, we repeated this validation over two periods for which the geodetic MBs were calculated (Sect. 4.2).

The traditional, as well as nonlinear, MBs over 2003–2014 were not statistically different from the geodetic MB, and the null hypothesis H₀ (the cumulative glaciological MB is not statistically different from the geodetic MB) was accepted at 95 % and 90 % levels (Zemp et al., 2013). However, over 2014–2020, both traditional and nonlinear MBs were statistically different from the geodetic MBs, and the null hypothesis H₀ was rejected at 95 % and 90 % levels. This showed that the systematic biases were significant over 2014–2020 (Table 2). Even though we did not observe a significant bias over 2003–2014, we decided to calibrate the traditional and nonlinear MBs over both periods, as suggested in previous studies (Thibert et al., 2008; Huss et al., 2009; Andreassen et al., 2016; Wagnon et al., 2021).

In the calibration procedure, the annual relative variability in glacier-wide MBs is taken from the MB series and the series was fitted to the multi-annual geodetic MB, B_g, as follows:

$$\begin{array}{}\text{(4)}& B_{\mathrm{a},\text{cal}}=B_{\mathrm{a}}+{\displaystyle \frac{\left(B_{\mathrm{g}}-{\sum }_N{B}_{\mathrm{a}}\right)}{N}},\end{array}$$

where B_a,cal is the annual calibrated glacier-wide MB, and N is the number of years over which the geodetic MB has been estimated. It should be mentioned that the MBs obtained from traditional method or nonlinear model refer only to the surface MB, whereas the geodetic MBs also integrate the internal and basal MBs, assumed to be small compared to the surface MB (Cuffey and Paterson, 2010).

3.8 Calibration of mean altitudinal mass balances

The mean altitudinal MBs (b_e,t) for each 50 m elevation band (e) and each year (t) were computed using Eq. (1), exploiting the values of α_i, β_t, and γ_i obtained from the nonlinear model. These altitudinal mean MBs were adjusted to fit the calibrated annual glacier-wide MBs, following Zemp et al. (2013). The calibrated altitudinal mean MB ($b_{\mathrm{e},t,\text{cal}}$) for each year is estimated as

$$\begin{array}{}\text{(5)}& b_{\mathrm{e},t,\text{cal}}=b_{\mathrm{e},t}-B_{\mathrm{a}}+B_{\mathrm{a},\text{cal}},\end{array}$$

where B_a is the uncalibrated annual nonlinear MBs, and B_a,cal is the calibrated annual nonlinear MBs. The equilibrium line altitude (ELA_cal) and MB gradient for each year (t) are also estimated by plotting the linear regression over the calibrated annual mean altitudinal MBs ($b_{\mathrm{e},t,\text{cal}}$) over an elevation range of 4375–5225 m. Finally, using the calibrated ELAs, the calibrated accumulation area ratios (AARs) were estimated each year (Table 3).

Table 1 The areal and geodetic mass changes on Chhota Shigri, Sichum, and Hamtah glaciers over the 2003–2014 and 2014–2020 periods.

3.9 Random error estimation in nonlinear mass balances

The random error (${\mathit{\sigma }}_{{\mathrm{B}}_{n,\text{cal}}}$) in calibrated nonlinear glacier-wide MB is estimated following:

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{{\mathrm{B}}_{n,\text{cal}}}=\pm \sqrt{{\displaystyle \frac{{\mathit{\sigma }}_{B_{\mathrm{g}}}^{\mathrm{2}}}{N}}+\mathrm{\Sigma }{s}_i^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}}.\end{array}$$

${\mathit{\sigma }}_{B_{\mathrm{g}}}$ is the error in the geodetic MBs (${\mathit{\sigma }}_{B_{\mathrm{g}}}$ = 0.57 and 0.36 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over 2003–2014 and 2014–2020, respectively), N is the number of years for geodetic MB estimation (Sect. 3.3), s_i terms represent the relative areas of each 50 m elevation band (except for the 5400–5850 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ range that has been treated as a single band) compared to the total glacier area (therefore, Σs_i = 1), and σ_ε = 0.30 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ is the standard deviation of the residual term of Eq. (2) obtained with the nonlinear model (Sect. 3.2). Equation (6) is valid for the hydrological years within calibration periods (2003–2014 and 2014–2020). The random errors in nonlinear glacier-wide MBs for 2002/2003 and 2020–2023 hydrological years were estimated following the procedure described in Wagnon et al. (2021). The mean annual random error, ${\mathit{\sigma }}_{{\mathrm{B}}_{n,\text{cal}}}$, of the calibrated nonlinear glacier-wide MB was estimated to be ± 0.19 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over 2002–2023, with slightly higher random errors for the years outside the calibration period (Table 3).

Figure 6 Glacier area change in Chhota Shigri, Sichum, and Hamtah glaciers between 2003 and 2020 (background image is Pléiades satellite imagery of 12 September 2020; CNES 2020, Distribution Airbus Defence and Space).

4 Results

4.1 Glacier area changes since 2003

Chhota Shigri, Sichum, and Hamtah glaciers showed limited areal changes since 2003, mostly restricted to the snout area (Table 1; Fig. 6). The estimated debris cover, corresponding to September 2020, was 12 %, 22 %, and 79 % of the total area on Chhota Shigri, Sichum, and Hamtah glaciers, respectively (Table 1). During 2003–2020, the total area change for each glacier was very small, with a deglaciation rate of −0.07 ± 0.22 % a⁻¹, −0.07 ± 0.22 % a⁻¹, and −0.03 ± 0.19 % a⁻¹ for Chhota, Sichum, and Hamtah, respectively (Table 1).

Figure 7 The thickness changes for Chhota Shigri, Sichum, and Hamtah glaciers differencing the ASTER 2003 (8 October 2003) and Pléiades (26 September 2014) DEMs over 2003–2014 and Pléiades DEMs (26 September 2014 and 12 September 2020) over 2014–2020.

Table 2 Cumulative MBs (in parenthesis; mean annual MBs) from the traditional method, nonlinear model, and geodetic estimates over available periods. The balance year 2002/2003 is not included here as it is not covered in the geodetic estimate available over 2003–2014. The cumulative traditional MB over the 2014–2020 period has been estimated by adding the modelled annual MB for 2019/2020 (Srivastava and Azam, 2022b). All units are in $\mathrm{m}\mathrm{w}.\mathrm{e}.$ ($\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$).

^∗ Estimated from traditional MBs (2014–2019) and modelled MB (2019/2020).

4.2 Geodetic mass balances

The maps of elevation changes for 2003–2014 and 2014–2020 periods indicate a general pattern of thinning for the glacier tongues and limited changes in the upper reaches of the glaciers (Fig. 7). The area-weighted geodetic MB of Chhota Shigri Glacier (including WT) was −0.43 ± 0.08 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over 2003–2020 (Table 1), with a higher annual wastage of −0.51 ± 0.06 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over 2014–2020 compared to −0.38 ± 0.10 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over 2003–2014 (Table 2). Sichum and Hamtah glaciers showed slightly stronger annual mass wastage of −0.57 ± 0.08 and −0.51 ± 0.08 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$, respectively, over 2003–2020, with similarly an increased mass wastage over the recent period (2014–2020) (Table 1). The slightly more negative glacier-wide MBs on all these glaciers during 2014–2020 agree with a recent study, suggesting an increased wastage over the recent decade in the Himalaya (Hugonnet et al., 2021).

The mean annual geodetic mass wastage of −0.43 ± 0.08 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ on Chhota Shigri Glacier over 2003–2020 is in good agreement with the region-wide mean glacier mass wastage of −0.37 ± 0.15 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over the whole Lahaul–Spiti region (glacierized area = 7960 km²) during a slightly different period (2000–2016) from multiple ASTER DEMs (Brun et al., 2017). Hence, Chhota Shigri is a reference glacier in the Himalaya (Azam, 2021) and a representative glacier for the whole Lahaul–Spiti region, as already suggested (Vincent et al., 2013).

Figure 8 Calibrated nonlinear annual glacier-wide MBs (with random errors) over 2002–2023, traditional cumulative MBs over 2002–2023, nonlinear cumulative MBs over 2002–2023, calibrated nonlinear cumulative MBs over 2002–2023, calibrated traditional cumulative MBs over 2002–2023, and geodetic MBs over 2003–2014 and 2014–2020 (with estimated uncertainties). The cumulative traditional MB series (2002–2019) is completed until 2023 by adding the modelled MB of 2019/2020 and 2020/2021 from Srivastava and Azam (2022b).

4.3 Annual and cumulative glacier-wide mass balances since 2002

Table 2 and Fig. 8 show the traditional and nonlinear MBs (before and after calibration) and geodetic MBs over available periods. The traditional MBs were not available for 2019/2020 and 2020/2021 (Sect. 3.1); therefore, to calibrate these MBs and to cover the geodetic observations, the modelled MBs (2019/2020 = 0.07 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ and 2020/2021 = −1.17 $\mathrm{m}\mathrm{w}.\mathrm{e}.$) from surface energy balance approach (Srivastava and Azam, 2022b) were added to the series.

Table 3 Calibrated nonlinear MBs ($B_{{\mathrm{a}}_{n,\text{cal}}}$), calibrated traditional MBs ($B_{{\mathrm{a}}_{t,\text{cal}}}$), MB gradients ($\mathrm{d}b/\mathrm{d}z$), ELA_cal, and AAR_cal on Chhota Shigri Glacier between 2002 and 2023.

^∗ The calibrated traditional MBs for 2019/2020 and 2020/2021 years are originally from the model (Srivastava and Azam, 2022b).

Compared to uncalibrated traditional MB series, uncalibrated nonlinear MB series showed much smaller biases, with a slightly negative bias of −0.03 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ (against a bias of −0.10 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ in traditional MBs) over 2003–2014 and of −0.17 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ (against a bias of 0.33 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ in traditional MBs) over 2014–2020 (Table 2; Fig. 8). Therefore, following Eq. (4), the nonlinear annual MBs were systematically increased by 0.03 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over 2003–2014 and by 0.17 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over 2014–2020, while traditional MBs were systematically increased by 0.10 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over 2003–2014 and decreased by 0.33 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over 2014–2020 to match the geodetic estimates (Fig. 8). The hydrological years 2002/2003 and 2020–2023 are outside the calibration periods, but these years were also calibrated by the mean values of biases observed over 2003–2014 and 2014–2020, respectively. To avoid confusion, we discussed only the calibrated nonlinear glacier-wide MBs in the paper, although the calibrated traditional MBs are given in Table 2 and 3 for reference.

The annual calibrated glacier-wide MB from the nonlinear model varied from 0.53 ± 0.16 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ in 2018/2019 to −1.71 ± 0.24 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ in 2021/2022, with a standard deviation of 0.65 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ during 2002–2023 (Table 3). In the 21-year-long MB series, 6 hydrological years (2004/2005, 2008/2009, 2009/2010, 2018/2019, 2020/2021, and 2022/2023) showed positive/near-steady-state MBs. The mean annual glacier-wide MB was estimated to be −0.47 ± 0.19 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$, equivalent to a cumulative loss of −9.81 ± 0.87 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ over 2002–2023 (Table 3). The uncertainty in cumulative mass loss comes from error propagation law.

4.4 Equilibrium line altitude and accumulation area ratio

Using the calibrated mean altitudinal MBs (Sect. 3.8), the equilibrium line altitude ELA_cal, accumulation area ratio AAR_cal and MB gradients ($\mathrm{d}b/\mathrm{d}z$) were also estimated. The maximum ELA_cal was 5248 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ corresponding to the most negative MB of −1.71 ± 0.24 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ and minimum AAR_cal of 19 % in 2021/2022, while the minimum ELA_cal was 4911 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, corresponding to a positive MB of 0.49 ± 0.19 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ and a maximum AAR_cal of 67 % in 2004/2005. The mean ELA_cal was 5070 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, corresponding to a mean mass wastage of −0.47 ± 0.19 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ and mean AAR_cal of 44 % over 2002–2023.

Figure 9 The ELA and AAR as a function of annual glacier-wide MB.

The annual ELA_cal and AAR_cal showed good correlations with annual glacier-wide MBs (r² = 0.98 and 0.97, respectively) over 2002–2023 (Fig. 5). The ELA_cal for a zero glacier-wide MB (ELA₀) was also computed from the regression between glacier-wide MBs and ELA_cal over 2002–2023 and calculated as ∼ 5001 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 9). Similarly, AAR₀was computed as ∼ 54 % for steady-state glacier-wide MB.

5 Discussion

5.1 Biases in glacier-wide mass balances and performance of nonlinear model

A total of 358 annual ablation and 65 annual accumulation point measurements were collected on Chhota Shigri Glacier over 2002–2023 to estimate the glacier-wide MBs (five ablation and five accumulation point MB measurements were removed before final model run; Sect. 3.3). Figure 10 shows the temporal evolution of the number of these point measurements, and Table S1 provides the details about these point MBs. In general, the point MB measurement network (especially the accumulation points) has been poor after 2014 (Sect. 3.1; Fig. 10). The eastern accumulation site at 5550 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 1) could only be accessed five times (2003, 2004, 2005, 2009, and 2011) over the 2002–2023 period, while no accumulation measurements were done in 2018, 2020, and 2021 (Sect. 3.1). Occasionally, the ablation measurements were also missing due to missing stakes (heavy ablation or destroyed stakes). In the traditional method, these missing measurements were filled with extrapolated values from nearby ablation/accumulation MB measurements or the previous years' point MB measurements to estimate the glacier-wide MBs (Azam et al., 2016; Mandal et al., 2020; Table S1).

Figure 10 Number of available ablation, accumulation, and total point MBs for each hydrological year between 2002 and 2023.

The systematic biases in glacier-wide annual MB series with the same monitoring network are expected to be of the same sign throughout the observation period, and the series is systematically adjusted to match the geodetic MBs available over one or more periods (Zemp et al., 2013; Wagnon et al., 2021). Nonlinear MB series on Chhota Shigri Glacier showed negative biases (−0.03 and −0.17 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over the 2003–2014 and 2014–2020 periods, respectively), suggesting that the nonlinear model can reasonably estimate the glacier-wide MBs with the existing monitoring network. Conversely, the traditional MB series showed a negative bias (−0.10 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$) over the 2003–2014 period and a large, positive bias (0.33 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$) over the 2014–2020 (Fig. 8; Table 2). The major disagreement between the cumulative nonlinear and traditional MB curves after 2017 (Fig. 8) is likely due to a degradation of the quality of field observations due to harsh weather, field surveys that are too short, or observers not being sufficiently experienced (Fig. 10; Table S1; Sect. 3.1). Wagnon et al. (2021) performed a thorough analysis on Mera Glacier (Dush Koshi basin, Nepal) and identified the precise source of systematic bias in the glacier-wide MB by comparing the surface-specific mass balance calculated using the traditional glaciological method of a specific zone on the glacier with that derived from the ice flux method (based on the mass conservation equation). Unfortunately, we could not conduct such an analysis in the current study due to insufficient data availability. However, future research will include this comparative analysis to uncover any systematic biases in the glacier-wide MB data series for the Chhota Shigri Glacier.

To further investigate the performance of the nonlinear model compared to the traditional MB method, we calibrated both the MB series with the geodetic MB estimated using ASTER (8 October 2003) and Pléiades (12 September 2020) DEMs (details in SI) and used the geodetic MB over 2003–2014 (Sect. 4.2) to validate both the calibrated series. The calibrated nonlinear MB series showed a good agreement with the available geodetic MB (−3.88 against −4.18 $\mathrm{m}\mathrm{w}.\mathrm{e}.$), while the traditional MB showed very strong deviation from the geodetic MB over 2003–2014 (−6.13 against −4.18 $\mathrm{m}\mathrm{w}.\mathrm{e}.$) (Fig. S1 in the Supplement). This good agreement between nonlinear and geodetic MBs over 2003–2014 shows the robustness of the nonlinear model for the glacier-wide mass balance estimation. Furthermore, this comparison also highlights the importance of using short-duration geodetic MB estimates for the calibration process, as with two calibration periods, the calibrated traditional MB is in better agreement with the geodetic MB (Fig. S1).

The nonlinear model shows a much better agreement with geodetic MBs than the traditional method (Fig. 8; Table 2) mainly due to the (i) capability of the nonlinear model to better capture the spatial variability in surface MB from a heterogeneous, discontinuous, and limited point MB data series than the traditional method (Vincent et al., 2018); (ii) correction/exclusion of erroneous measurements (Sect. 3.3); and (iii) exclusion of the extrapolated ablation/accumulation points in the nonlinear model that might have introduced biases in traditional MB (Fig. S2 in the Supplement). The extrapolated point MBs in the accumulation area showed a difference ranging from −1.98 to 1.74 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ between modelled and extrapolated values, especially post-2014 (Figs. S2 and S3). The better performance of the nonlinear model suggests that the extrapolation of point accumulations (in case of missing point measurements) to estimate the glacier-wide MB using the traditional method is risky.

Figure 11 Glacier-wide MBs (a), temporal (β_t) (b) and spatial terms (α_e and γ_e) (panels (c) and (d), respectively) obtained with the nonlinear model, following two different scenarios as a function of their original values obtained with the full data set. In the first scenario (2022/23_2020), we remove all the data from 2022/2023 (24 point MBs), except for 2 located at the observation points in 2019/2020 (see Sect. 5.2). In the second scenario (2022/23_SLA), we remove all the data from 2022/2023 and keep only two point MB data (= 0 $\mathrm{m}\mathrm{w}.\mathrm{e}.$) obtained along the SLA (see Sect. 5.3). The filled dots highlight the test year of 2022/2023.

5.2 2019/2020 glacier-wide mass balance from two point mass balances

The spatial and temporal terms in Eq. (2) are computed from a data sample available from the whole series; therefore, MB computation is expected to be affected by missing data from any single year (or, in general, from all years whenever data are missing). The glacier-wide MB for 2019/2020 was estimated using only two point MB observations (Sect. 3.2; Table S1); therefore, it might have biases (Lliboutry, 1974; Vincent et al., 2018).

To investigate the additional error, we selected the year 2022/2023 to test the performance of the nonlinear model. The 2022/2023 year was selected because it is among the years with the maximum of point MB observations, and they were performed at their original locations. The nonlinear model was re-run over the 2002–2023 period, keeping only 2 point MB data (out of 26) for 2022/2023 year, corresponding to the locations of the 2 point MB measurements in 2019/2020. With only two point MBs, the glacier-wide MB for 2022/2023 was recomputed to be 0.13 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ against the original MB of 0.04 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$, with a difference of 0.09 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$, while all other years' glacier-wide MBs were changed by a maximum of ± 0.01 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ (Fig. 11a). As expected, the changes in the temporal term, β_t, having a glacier-wide significance, showed significant deviation from 0.93 to 1.06 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ for 2022/2023 year, while for other years it changed by a maximum up to ± 0.04 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ (Fig. 11b). Conversely, the deviations in mean altitudinal spatial terms α_e and γ_e were very small (maximum up to ± 0.06 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ and ± 0.005, respectively) (Fig. 11c and d). Therefore, the temporal term (β_t) in Eq. (2) mainly controls the annual glacier-wide MB.

The deviation of 0.09 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ in glacier-wide MB estimated with only two point MBs is less than the estimated random error of 0.27 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ in the 2022/2023 glacier-wide MB in the original model run; therefore, it is assumed that the error in 2019/2020 glacier-wide MB due to restricted number of MB measurements is also less than the estimated random error of 0.16 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ (Table 3). Unlike the traditional MB method, the nonlinear model can fill the gaps in glacier-wide MB for which some point MB observations are missing and can provide a consistent series of temporal fluctuations.

5.3 2020/2021 glacier-wide mass balance from nonlinear model–SLA method

The glacier-wide MB for 2020/2021 year was estimated by inferring two point MB input from end-of-summer SLA, assuming it to be equivalent to ELA (i.e. MB = 0 $\mathrm{m}\mathrm{w}.\mathrm{e}.$) (Sect. 3.2; Fig. 3). Due to only two point MB input data, the modelled glacier-wide MB for 2020/2021 may also have additional errors.

To quantify this error, we repeated the same exercise as in Sect. 5.2 for the year 2022/2023, this time again keeping two point MB data of 2022/2023 but at the two sites where point MB data have been assessed to be zero in 2020/2021. The resulting 2022/2023 glacier-wide MB is 0.26 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$, 0.22 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ higher than the original value (Fig. 11a), mainly explained by the β_t term (Fig. 11b). This difference is still lower than the estimated random error of 0.27 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ in 2022/2023 (Table 3). However, there are still possible biases in glacier-wide MB of 2020/2021 year as the SLA was delineated from a Sentinel image from 6 September 2021 (Sect. 3.2; Fig. 3) that is not exactly from the end of ablation season (30 September) on Chhota Shigri Glacier. The surface energy balance model estimated a MB of −0.19 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ over the 6–30 September 2021 (Srivastava and Azam, 2022a). However, this seasonal offset correction in SLA-derived annual MB may be given, but it was avoided as the differences are within the estimated random error of 0.20 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ (Table 3). Our analysis shows that the glacier-wide MB can also be estimated from SLA using the nonlinear model if the field measurements cannot be carried out for some specific years.

However, the nonlinear model-SLA method has several limitations: (i) the delineated SLA must pass through a grid(s) having a previous point MB observation(s) (Fig. 3), as at least one previous measurement is required to run the model; (ii) the delineated SLA must be from the end of ablation season to consider it ELA; (iii) SLA delineation has its challenges, and often it is difficult to find the cloud-free image for delineation at the end of ablation season (Brun et al., 2015; Racoviteanu et al., 2019); and (iv) SLA is severely affected by recent snowfall and hence must be checked with in situ precipitation data before using SLA in nonlinear model. This latter point implies that the ELA can be inferred from the end-of-ablation-season SLA, which is not always possible over glaciers, especially in monsoon-dominated regions (Brun et al., 2015).

5.4 Recommendation: apply the nonlinear model on other glaciers

This study demonstrates that the nonlinear model outperforms the traditional method for estimating glacier-wide MB (Sect. 5.1). Apart from the present research, the nonlinear model has been applied only to the Mera Glacier in the Himalaya (Wagnon et al., 2021) and on Argentière, Saint-Sorlin, Mer de Glace, and Gébroulaz (France, Alps), Zongo (Bolivia, Andes), and Nigardsbreen (Norway, Scandinavia) glaciers (Vincent et al., 2018).

Equation (1) includes a spatial effect term (γ_i) that accounts for the standard deviations in point MBs across elevation. This term typically requires around 10 years of point MB observations to be accurately estimated (Vincent et al., 2018). Therefore, applying the nonlinear model wherever MB observations are available for around 10 years is advisable, especially in the Himalaya, where data accessibility issues often lead to gaps in observations (Azam et al., 2018). We recommend extending the application of the nonlinear model to other Himalayan glaciers that have consistent MB observations spanning approximately 10 years, such as Kolahoi, Hōksar, and Sutri Dhaka glaciers in the western Himalaya and Chorabari, Dokriani Bamak, Pokalde, Rikha Samba, Yala, and West Changri Nup glaciers in the central Himalaya. However, the estimated glacier-wide MBs may contain systematic biases due to the distribution of point measurements across the glacier. Therefore, they should be verified and, if necessary, reanalysed using geodetic estimates.

6 Conclusions

This work reanalyses glacier-wide MBs by combining the traditional reanalysis framework (Zemp et al., 2013) and the nonlinear MB model (Vincent et al., 2018). Previously, the annual glacier-wide MBs had been estimated on Chhota Shigri Glacier since 2002, when the traditional glaciological method was applied using heterogeneous in situ point MB measurements. The heterogeneous measurement network does not always catch the large spatiotemporal variability in point MBs; hence, the point MB–elevation relationship is insufficient to investigate the changes in glacier-wide MBs. Therefore, we applied the nonlinear model to compute the glacier-wide MBs of Chhota Shigri Glacier as it enables the computation of the glacier-wide MB from a heterogeneous in situ point MB network. The nonlinear model was used to detect the measurement errors. Of the 423 point measurements, 5 were corrected from field notebooks, and 10 were recognized as wrong observations and discarded before running the final model.

ASTER and Pléiades DEMs were used to estimate the geodetic MBs over 2003–2014 and 2014–2020, and they have been used to reanalyse the nonlinear MBs. Nonlinear MBs agreed well with the geodetic estimates available over 2003–2014 and 2014–2020, unlike traditional MBs that showed large differences, especially over the 2014–2020 period. The reanalysed nonlinear MBs showed a large annual variability ranging from 0.53 ± 0.16 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ in 2018/2019 to −1.71 ± 0.24 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ in 2021/2022. The Chhota Shigri Glacier is imbalanced, with a mean mass wastage of −0.47 ± 0.19 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$, equivalent to a cumulative loss of −9.81 ± 0.87 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ over 2002–2023.

With the 21-year-long MB observations, the Chhota Shigri Glacier MB series is the longest in the Himalaya. This work has enabled the data set to be extended, optimized, and corrected to provide the best possible mass balance series for this benchmark glacier. We plan to monitor this glacier over a long period, with repeated satellite image acquisitions by the Pléiades Glacier Observatory to regularly validate/calibrate the glacier-wide MB, typically every 5 years.

Our detailed analysis suggests that the nonlinear model performs better in calculating the glacier-wide MB than the traditional method as (i) the nonlinear MBs are in much better agreement with the geodetic MB estimates, (ii) it can detect erroneous measurements, (iii) it provides better glacier-wide MBs than those of the traditional method when the observational network is very limited, and (iv) glacier-wide MB can be computed using SLA if the ablation-end SLA passes through a grid cell that contains point MB observations from previous years. Therefore, the application of the nonlinear model is suggested on all monitored glaciers whenever data are sufficient. It becomes even more relevant in the Himalaya, where data are sometimes missing due to access issues. However, the estimated glacier-wide MBs may contain systematic bias (arises from the distribution of point measurements over the glacier) and, therefore, should be checked and, if necessary, reanalysed with geodetic estimates.