Reconstruction of Annual Glacier Mass Balance from Remote Sensing-Derived Average Glacier-Wide Albedo

Abstract

Annual mass balance is an important reflection of glacier status that is also very sensitive to climate fluctuations. However, there is no effective and universal albedo-based method for the reconstruction of annual mass balance due to the scarcity of field observations. Here, we present an improved albedo–mass balance (IAMB) method to estimate annual glacier surface mass balance series using remote sensing techniques. The averaged glacier-wide albedo derived with the MODImLab algorithm during the summer season provides an effective proxy of the annual mass change. Defined as the variation in the albedo as a function of elevation change, the altitude–albedo gradient ($\partial z/\partial \alpha$) can be obtained from a glacier digital elevation model (DEM) and optical images. The Chhota Shigri glacier situated in the western Himalayas was selected to test and assess the accuracy of this method over the period from 2003 to 2014. Reconstructed annual mass budgets correlated well with those from the observed records, with an average difference and root mean square error (RMSE) of −0.75 mm w.e. a⁻¹ and 274.91 mm w.e. a⁻¹, respectively, indicating that the IAMB method holds promise for glacier mass change monitoring. This study provides a new technique for annual mass balance estimation that can be applied to glaciers with no or few mass balance observations.

1. Introduction

Glacier mass balance, which is primarily controlled by mass ablation and accumulation, is a direct indicator of glacier status [1] and a natural record of climate change [2]. As the atmospheric air temperature increases [3,4], obvious recessions can be observed in global mountain glaciers [5,6]. The unabated shrinkage and further demise of glaciers worldwide will lead to multiple impacts for social-ecological systems, resulting not only in an increase in global sea level and changes in water availability but also affecting agriculture, tourism, ecosystems and human livelihoods [7]. For glaciers with rugged mountainous terrain, long-term changes in cryospheric (e.g., glaciers, permafrost) conditions are likely to lead to serious natural hazards. Recently, a large rock and ice avalanche descended into the Ronti Gad valley, causing widespread devastation and damage among the downstream inhabitants [8]. Therefore, studying glacier change is very important, especially to understand glacier characteristics and future trends at regional and global scales. The well-established glaciological method is a traditional and widely used method for annual mass balance monitoring. However, field observations of annual glacier mass balance remain scarce due to the inaccessibility and harsh environmental conditions of the glaciated regions [9]. Thus, it is essential to provide an efficient method that can be used to quantitatively evaluate the annual glacier mass balance, as well as its response to changes in climatic controls. Glacier albedo is defined as the ratio between the reflected and incident flux density at the glacier surface [10,11]. This variable determines the net shortwave radiation absorbed by the glacier and, in turn, the energy available for melting [12,13,14]. Similarly to the accumulation–area ratio (AAR) [15] and the equilibrium line altitude (ELA) [16,17], the minimum glacier albedo has been proved to be significantly correlated with annual glacier surface mass balance and used for mass balance reconstruction [18,19,20,21]. Additionally, Zhang et al. (2018) utilized the MODIS-derived annual minimum-averaged glacier-wide albedo to predict specific mass budgets for glaciers with similar climate regimes, implying that the minimum albedo is potentially promising for deriving large-coverage and long-term annual glacier mass fluctuations [22]. Nevertheless, the existing albedo-based method requires adequate glaciological measurements for calibration [20,22], which prevents its wider application with glaciers with no field mass balance observations.

In addition to the minimum glacier albedo, the averaged albedo for an entire glacier over the summer period—i.e., the average summer albedo—is also used as a proxy of net annual mass budgets. Abundant prior research has confirmed that the average albedo is highly correlated with the net surface mass balance records, such as work on the Svalbard glaciers [23,24] and Greenland glaciers [25,26]. Recently, Williamson et al. (2020) undertook an exhaustive comparison of MODIS Terra albedo and net annual mass balance for five Canadian Arctic glaciers over the period 2002–2016. The results indicated that the interpolated average albedo was a more accurate proxy than the minimum or raw average values for the glacier net mass change, as use of the interpolated average albedo could reduce the impacts of inter-annual climate variability, cloud cover and statistical accidents [27]. However, there is no robust and universal albedo-based method that can be used to reconstruct annual mass balance at the regional scale, and there is no research to date that evaluates the usage of the average summer albedo to calculate the annual mass balance for Tibetan Plateau (TP) glaciers.

In this study, we aimed to reconstruct the annual mass balance series without in situ mass balance measurements using the averaged summer surface albedo. We first present an improved albedo–mass balance (IAMB) method, including the details for the establishment of the algorithm, parameter estimation, sensitivity analysis and uncertainty analysis. Then, the estimation of the mass balance–albedo gradient ($\partial b/\partial \alpha$) and altitude–albedo gradient ($\partial z/\partial \alpha$) using remote sensing techniques is introduced, and the effect of elevation on their estimation is also analyzed. Finally, we tested the accuracy of our IAMB method with the Chhota Shigri glacier and reconstructed the annual mass balance fluctuations over the period 2003–2014 by employing the average summer glacier-wide albedo derived with the MODImLab algorithm [18,28].

2. Chhota Shigri Glacier

Situated in the Chandra-Bhaga river basin of Lahaul valley, India, the Chhota Shigri glacier (32°17′N, 77°35′E) is a valley-type glacier situated over the western Himalayas (Figure 1). It covers an area of about 15.48 km², with an elevation ranging from 4050 to 5830 m a.s.l. [29,30] and a mean ELA of 4960 m a.s.l. [31]. This glacier is oriented roughly south to north in its ablation area and the tributaries have a variety of orientations in the accumulation area (Figure 1b). In general, the Chhota Shigri glacier is a non-surging, temperate and clean glacier, except for a proportion of about 12% comprising debris-covered areas in the lower parts (<4500 m a.s.l.) [32]. Moreover, the Chhota Shigri glacier has continuous in situ annual mass balance records from 2002 [33,34], and these long and consecutive time series of field measurements were used to facilitate the testing and validation of the established method.

3. Data and Methodology

3.1. Datasets Used in This Study

3.1.1. MODImLab Products and Glacier-Wide Albedo

MODImLab, developed by Sirguey et al. (2009), is an algorithm that can derive subpixel snow fractions and snow and ice albedos from the Moderate Resolution Imaging Spectroradiometer (MODIS) at 250 m spatial resolution [28]. It retrieves the surface albedo based on atmospheric and topographic correction of MODIS Level 1B swath top-of-atmosphere (TOA) calibrated radiances data, SRTM DEMs, atmospheric profiles and cloud mask data [18,28,35,36]. MODImLab provides estimates of both white-sky and blue-sky albedo, and the overall accuracy of the retrieved broadband albedo under clear-sky conditions was found to be about 10% based on field measurements [18,20]. Therefore, the MODImLab albedo product can accurately depict the characteristics of ice and snow surfaces, and it is preferred to the coarser MODIS daily snow product (MOD10A1), especially for mountainous glaciers [19,20,36].

White-sky albedo was considered here due to its reduced sensitivity to illumination conditions [18]. Near-daily albedo maps from the summer season (1 June to 30 September) over the period 2003–2014 were processed for Chhota Shigri glacier. Additionally, a glacier-wide albedo map for the melt season was produced by averaging all available values (cloud-free pixels) for each grid cell. The annual average glacier-wide albedo was then averaged into a single value for each year using an equally weighted mean. We used 10% of the average glacier-wide albedo value as the value for the uncertainty of albedo retrieval, following previous studies [18,20].

3.1.2. Landsat Imagery

Landsat TM/ETM+/OLI data (level 1T-corrected) obtained from the USGS website (http://glovis.usgs.gov/, accessed on 16 January 2022) with a spatial resolution of 30 m were utilized in this study to estimate the albedo variation across the glacier surface. Since no clear-sky Landsat images of the Chhota Shigri glacier during the ablation period from 2002 to 2007 were available, we expanded the span of the study period to 2017. Moreover, the influence of low sun elevation angles (<40°) on terrain shielding is well-known and results in lower albedo values [37]. Therefore, the Landsat imagery with lower sun elevation angles (<45°) was excluded to obtain more accurate glacier surface albedo estimates in this study. Finally, 16 cloud-free Landsat images from over the Chhota Shigri glacier acquired between 2000 and 2017 (between 1 June and 30 September in each year) were selected and utilized for albedo estimation. Details of the Landsat scenes used are summarized in the following

Table 1. The Landsat data for the Chhota Shigri glacier used in this study.

Satellite	Sensor	Imagery Date (YYYY-MM-DD)	Sun Elevation/°
Landsat-7	ETM+	2000-09-29	49.94
Landsat-7	ETM+	2001-06-28	66.06
Landsat-5	TM	2001-09-24	50.07
Landsat-5	TM	2008-06-23	65.93
Landsat-5	TM	2008-09-27	49.75
Landsat-5	TM	2009-09-30	49.64
Landsat-5	TM	2010-07-15	64.95
Landsat-8	OLI	2013-09-25	52.59
Landsat-8	OLI	2014-08-11	63.27
Landsat-8	OLI	2014-09-28	51.54
Landsat-8	OLI	2015-07-29	65.20
Landsat-8	OLI	2015-08-30	59.50
Landsat-8	OLI	2015-09-15	55.42
Landsat-8	OLI	2016-08-16	62.31
Landsat-8	OLI	2017-09-04	58.22
Landsat-8	OLI	2017-09-20	53.87

.

3.1.3. Field Measurements

In situ mass balance was first measured at the Chhota Shigri glacier in 1987 using the direct glaciological method, and measurements were restarted to obtain continuous observations in 2002 [33,34]. About 22 bamboo stakes were set up in the ablation area at sites ranging from 4300 to 4900 m a.s.l. to estimate ablation, and four to six snow pits were dug in the accumulation zone to estimate net accumulation between 5150 and 5550 m a.s.l. [29,31,38]. The distribution of ablation and accumulation sites and the details of the direct mass balance measurements were provided by Wagnon et al. (2007) [31]. In this study, the field measurements of annual glacier-wide mass balance series from between 2003 and 2014 [29,38] were used for calculation of the mean mass balance and validation of the method.

In addition, the vertical mass–balance gradient ($\partial b/\partial z$), which represents the time- and space-averaged mass balance gradient with altitude, was also utilized for annual mass balance estimation. Many previous studies have demonstrated that the $\partial b/\partial z$ has been relatively stable in the Chhota Shigri glacier over the last twenty years [29,31,38]. Therefore, the vertical mass–balance gradient of 0.66 ± 0.09 m w.e. (100 m)⁻¹ from Azam et al. (2016) [29] was used in our method, considering that the observed time period was the most consistent with our study.

3.1.4. Glacier Outlines

The glacier boundaries derived from the Randolph Glacier Inventory (RGI 6.0, released in July 2017) were employed to identify the glacierized regions of the Chhota Shigri glacier. In order to avoid discrepancies in glacier outlines caused by different definitions and delineations [39], following Dumont et al. (2012) [18] and Brun et al. (2015) [19], we manually defined the glacier mask using a Landsat OLI image acquired on 28 September 2014. Moreover, we also removed the pixels corresponding to debris-covered areas (<4600 m a.s.l.) to reduce their influence on the satellite-derived snow and ice albedo. Adjusted outlines of the entire Chhota Shigri glacier (black solid line) and of the debris-free regions (red dashed line) are shown in Figure 1b. This process resulted in 204 pixels in the MODImLab 250 m resolution albedo maps of the Chhota Shigri glacier, which is about twice (89 pixels) that in Brun et al. (2015) [19]. The larger glacier mask can help us to achieve more accurate estimations of averaged glacier-wide albedo when cloud cover is unavoidable.

3.2. The Improved Albedo–Mass Balance (IAMB) Method

With similarities to the ELA and AAR, a simple linear model can be used to describe the relationship between annual surface mass balance and the average summer albedo $\left({\alpha }_i\right)$. In this study, we used an improved albedo–mass balance (IAMB) method inspired by the ELA method [40]. Specifically, over the n years of the study period, the annual mass balance $b_i$ for year i can be expressed as:

$$b_i=\frac{\partial b}{\partial \alpha }{\alpha }_i-p$$

(1)

where $p$ is a constant and $\partial b/\partial \alpha$, not identically zero, represents the mass balance–albedo gradient across the equilibrium line. It results that

$$\overline{B}=\frac{1}{n}{\displaystyle \sum _{i=1}^n{b}_i}=\frac{1}{n}\frac{\partial b}{\partial \alpha }{\displaystyle \sum _{i=1}^n{\alpha }_i}-p$$

(2)

where $\overline{B}$ is the mean mass balance over the study period. Defining ${\alpha }_{eq}=\frac{p}{\partial b/\partial \alpha }$ yields:

$$\frac{\overline{B}}{\partial b/\partial \alpha }=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\alpha }_i}-{\alpha }_{eq}$$

(3)

Substituting into Equation (1), the estimated annual mass balance $b_i$ can be obtained as:

$$b_i=\overline{B}+\frac{\partial b}{\partial \alpha }({\alpha }_i-\frac{1}{n}{\displaystyle \sum _{i=1}^n{\alpha }_i})$$

(4)

Provided that the average summer albedo for each year (${\alpha }_i$) is known, the estimated annual mass balance $b_i$ depends on the determination of $\overline{B}$ and $\partial b/\partial \alpha$. $\overline{B}$ can be computed from field measurements, using photogrammetry or with other geodetic methods [41,42,43,44]. The estimation of $\partial b/\partial \alpha$ over a glacier is described in the next section.

3.3. Estimation of the Mass Balance–Albedo Gradient $\partial b/\partial \alpha$

$\partial b/\partial \alpha$ is the mass balance–albedo gradient in the vicinity of the ELA, which represents the trend in the mass balance as a function of albedo for a specific glacier. Since glacier albedo changes over time and across space, $\partial b/\partial \alpha$ cannot be readily obtained in field observations. Instead, $\partial b/\partial \alpha$ is derived from the altitude–albedo gradient ($\partial z/\partial \alpha$) and the mean vertical mass–balance gradient ($\partial b/\partial z$) as follows:

$$\frac{\partial b}{\partial \alpha }=\frac{\partial b}{\partial z}\frac{\partial z}{\partial \alpha }$$

(5)

$\partial b/\partial z$ can be obtained from field measurements. $\partial z/\partial \alpha$ is defined as the variation in the albedo with altitude for a glacier and can be acquired using remote sensing through the following three steps:

• Estimation of the glacier albedo. Cloud-free Landsat imageries were utilized to calculate the broadband glacier albedo with the method proposed by Knap et al. (1999) [45]. Data pre-processing, including calibration, Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) atmospheric correction, topographic correction [46] and spatial clipping, were undertaken prior to the albedo estimation;
• Extraction of the DEM and albedo data. We used the C-band Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) with approximately 30 m spatial resolution to acquire the glacier elevation information. Co-registration of the glacier DEM and Landsat imagery was completed to ensure the spatial consistency and extract the corresponding glacier elevation and albedo for each pixel;
• Calculation of $\partial z/\partial \alpha$. Assuming that the glacier elevation ranges from $H_{\min }$ (the lowest elevation) to $H_{\max }$ (the highest elevation), the glacier can be divided into $s$($s=\frac{H_{\max }-H_{\min }}{10}$) 10 m altitude intervals. The average glacier albedo (${\overline{\alpha }}_{ak}$) for altitude band $k (1\le k\le s)$ was calculated, providing a curve representing the variation in albedo with altitude, from which $\partial z/\partial \alpha$ could be obtained.

The averaged glacier albedo in the summer season increases with glacier elevation, and $\partial z/\partial \alpha$ varies with different years, glacier altitudes, locations and other local parameters. In attempting to obtain the glacial fluctuation over the studied period, we only calculated the $\partial z/\partial \alpha$ near the glacier ELA, which has been shown to be an important proxy of glacier annual mass balance [16,17]. Moreover, in this study, $\partial z/\partial \alpha$ values acquired in different years were aggregated into a mean value for the studied period. The detailed processing steps for the IAMB method for annual mass balance reconstruction are shown in Figure 2, and estimation of $\partial b/\partial \alpha$ is depicted in the central part of this figure.

3.4. Uncertainty Analysis

Using the above IAMB model, the annual glacier mass balance $b_i$ in Equation (4) can be simplified as a linear equation:

$$b_i=\overline{B}+g{t}_i$$

(6)

where $g=\partial b/\partial \alpha$, and $t_i={\alpha }_i-\frac{1}{n}{\displaystyle \sum _{i=1}^n{\alpha }_i}$ is the summer albedo anomaly in year $i (1\le i\le n)$.

The overall uncertainty in the annual glacier mass balance can be computed from the uncertainty for the three separate variables in Equation (6); i.e., ${\sigma }_{\overline{B}}^{}$, ${\sigma }_g^{}$ and ${\sigma }_{t_{{}_i}}^{}$. According to the error propagation method, the error variance for the estimated annual mass balance can be written as:

$${\sigma }_{b_i}^2={\sigma }_{\overline{B}}^2+t_i^2{\sigma }_g^2+g^2{\sigma }_{t_i}^2$$

(7)

• The uncertainty in ${\sigma }_{\overline{B}}^2$ associated with $\overline{B}$ is related to the glacier mass balance data used in the study. If $\overline{B}$ is averaged from the observed mass balance, ${\sigma }_{\overline{B}}^2$ is determined with in situ data from each year. If $\overline{B}$ is obtained from the geodetic mass balance over several years, then ${\sigma }_{\overline{B}}^2$ depends on the uncertainties associated with each glacier DEM, the geo-referencing and co-registration errors of the DEMs and the glacier density used to convert the change in elevation into the mass balance [47]. In this study, $\overline{B}$ was calculated from field annual mass balance series, and the uncertainty range for the measured data was fixed at ±0.40 m w.e. [29]. Therefore, the uncertainty in ${\sigma }_{\overline{B}}^2$ was also a constant number over the studied period;
• The uncertainty in $g$ can be computed using the error propagation approach [48], considering the uncertainty in each term of $\partial b/\partial z$ and $\partial z/\partial \alpha$. Given that $\partial b/\partial z$ is different from year to year and glacier to glacier, we used a variability of ±0.09 m w.e. (100 m)⁻¹ as an estimate of uncertainty [29]. The uncertainty in $\partial z/\partial \alpha$ results from (1) the glacier albedo retrieval from optical imagery, (2) the image registration of Landsat and SRTM DEM data and (3) the accuracy of the DEM. Combining these different sources of error, we estimated that the uncertainty of $\partial z/\partial \alpha$ was equal to the standard deviation, which was 12.69 m/(0.01 albedo) for the Chhota Shigri glacier (see

  Table 2. The variation in the altitude–albedo gradient ($\partial z/\partial \alpha$) within elevations of 200 m of the equilibrium line altitude (ELA) for the Chhota Shigri glacier.

  Date	$\partial \mathit{z}/\partial \mathit{\alpha } (\mathbf{m}/0.01\mathbf{albedo})$	Correlation Coefficient (R)
  2001-06-28	12.60	0.995
  2001-09-24	40.20	0.987
  2008-06-23	28.76	0.963
  2010-07-15	11.90	0.975
  2014-08-11	14.83	0.991
  2014-09-28	17.29	0.990
  2015-07-29	17.53	0.977
  2015-08-30	23.41	0.985
  2015-09-15	23.57	0.985
  2016-08-16	53.75	0.967
  2017-09-04	12.26	0.997
  2017-09-20	18.52	0.990
  Mean	22.88
  Standard deviation	12.69

  in Section 4.2);
• The uncertainty in the glacier summer albedo anomaly $t_i$ is not only affected by the summer albedo in year i but also those of all years. Based on the glacier-wide average albedo error referred to in Section 3.2, the uncertainty in $t_i$ can be estimated with the following expression:

  $${\sigma }_{t_i}^2=\frac{{\sigma }_{{\alpha }_1}^2+{\sigma }_{{\alpha }_2}^2+\dots +{(n-1)}^2{\sigma }_{{\alpha }_i}^2+\dots +{\sigma }_{{\alpha }_n}^2}{n^2}$$

  (8)

4. Results

4.1. The Averaged Glacier-Wide Albedo

Summer albedo maps retrieved from MODImLab products for the Chhota Shigri glacier over the 2003–2014 period are illustrated in Figure 3. The pixel cells in the albedo maps represent the annual mean albedo values over the melt season, and the change in colors from red to blue indicates a positive trend.

Although the spatial distribution of the average summer albedo varied each year, the distribution of albedo followed a similar pattern in general. Higher albedo was found in the accumulation area in the southeast part of the glacier, whereas lower albedo occurred in the ablation area over the main glacier body and the western flank, respectively. This also indicated that the surface albedo for the Chhota Shigri glacier is correlated with the glacier elevation. Furthermore, there were several discrete speckles on the edge of the glacier boundary in the albedo maps, which were probably due to the complex glacial terrain. Overall, Figure 3 shows that MODImLab performs well in mapping and characterizing the surface albedo at 250 m resolution over a rugged glacier topography.

The annual average summer albedos and field annual mass balances between 2003 and 2014 are illustrated in Figure 4. It can be seen that the annual variation in the average summer albedo was relatively stable, with the maximum, minimum and average values being 0.594, 0.464 and 0.528 respectively. Compared with the in situ mass balance, it can be observed that the most positive mass balance was recorded in 2009 and 2010 and the most negative in 2003 and 2004, which is in good agreement (R = 0.90, p < 0.01) with the trend for the average summer albedo. This indicates that the annual average summer albedo is a potentially good proxy of the surface mass balance.

4.2. Mass Balance–Albedo Gradient $\partial b/\partial z$

Based on the 16 cloud-free Landsat imageries, the glacier surface albedo and its variation with elevation at the Chhota Shigri glacier is given in Figure 5. Considering that the satellite-derived albedo may have been affected by dirty ice, bedrock and other debris in the lower parts, only the albedo above 4600 m a.s.l. (red dashed line in Figure 1b) was considered. The ELA of the Chhota Shigri glacier is also indicated by a gray dotted line in the Figure 5.

Note that the albedo results are divided into two categories, the first delineated by solid lines and the other by dashed lines. For the glacier albedo expressed by the solid lines, although the extent of albedo variation fluctuated from year to year, the shape of the change in the albedo with altitude remained similar. Specifically, a significant albedo increase was detected with the glacier elevation, which was accompanied by a larger and smaller albedo gradient below and above the glacier ELA, respectively. In addition, the surface albedo reached its maximum value around 5300 m a.s.l., with a relative difference (the absolute difference between the maximum and minimum divided by the magnitude of the maximum value) greater than 46%. The ablation and accumulation regions of the glacier situated under and above the glacier ELA, respectively, corresponded well with lower and higher albedo values, respectively, which is consistent with Ming et al.’s (2015) description of the evolution of glacier albedo in the summer season [13].

Nevertheless, the four albedo datasets in Figure 5 delineated by dashed lines (2000-09-29, 2008-09-27, 2009-09-30 and 2013-09-25) exhibit shapes that contrast with those of the solid lines. In particular, all of the albedo values at different altitudes were large and did not change significantly (<20%), with the minimum value being above 0.5. In terms of the spatial distribution, the change was less pronounced compared to the former, which was observable in a gradual albedo increase below the glacier ELA that was subsequently interrupted by a slight decrease. This was caused by early snow fall, as shown in the Figure 6, leading to a relatively high and homogeneous surface albedo. Therefore, we could not determine the glacier accumulation and ablation areas or the albedo changes with altitude from these remote sensing images alone. Finally, the albedo data shown by the solid lines in Figure 5 were used for the estimation of the altitude–albedo gradient ($\partial z/\partial \alpha$).

Table 2 displays the calculated $\partial z/\partial \alpha$ and the correlation coefficient (R) between the elevation and albedo in the vicinity of the ELA (±200 m) at the Chhota Shigri glacier. The $\partial z/\partial \alpha$ ranged from 11.90 to 53.75 and was probably associated with the melting state of the glacier surface at the Landsat acquisition time. For the Chhota Shigri glacier, the arithmetic mean of 22.88 m/ (0.01 albedo) was used to quantify the $\partial z/\partial \alpha$ over the entire study period.

Using the field $\partial b/\partial z$ measurement (see Section 3.1.3), the mass balance–albedo gradient ($\partial b/\partial \alpha$) for the Chhota Shigri glacier during 2003–2014 was estimated as 151 mm w.e./ (0.01 albedo).

4.3. Reconstructed Annual Mass Balance Series

Using the average summer albedo (${\alpha }_i$) presented in Section 4.1, the mean mass balance ($\overline{B}$) from Azam et al. (2016) [29] and the mass balance–albedo gradient ($\partial b/\partial \alpha$) presented in Section 4.2, the annual glacier-wide mass balance for the Chhota Shigri glacier was estimated for the period 2003–2014 with Equation (4). Figure 7a shows a comparison between the reconstructed annual mass balance values obtained with the IAMB method and the values from the field measurements. Considering annual data for the entire study period, the overall uncertainties in the derived mass balance for each year were almost constant over the entire period, with a mean value of 321 mm w.e. a⁻¹. In addition, the majority of the mass balance uncertainty resulted from the uncertainty in the $\partial b/\partial \alpha$ (see Section 3.3). Figure 7a shows a good correlation between the observed yearly glacier mass balance series and those computed using the IAMB method with R² = 0.81 (p < 0.01). The average difference and the root-mean-square error (RMSE) for the series were −0.75 mm w.e. a⁻¹ and 274.91 mm w.e. a⁻¹, respectively. Figure 7b compares the cumulative in situ mass balance [29,38] and that estimated with the IAMB method for the Chhota Shigri glacier. It can be noted that the IAMB glacier mass balances are in good agreement with the in situ cumulative mass balance.

Nevertheless, the year of 2005 exhibited a relatively large discrepancy of about 600 mm w.e. This may have been caused by the incomplete glacier mass balance observations in 2005. There were no mass balance measurements for part of the accumulation areas in this year [31]; hence, the real mass balance was lower than that shown by the measured values. Moreover, the inaccurate estimate of the average summer albedo for this year is another possible reason for the large difference. Although the estimation of mean glacier-wide albedo should be accurate when all albedo data for the summer season are averaged, it can be complicated by (1) cloud-cover shadows over the glacier surface, (2) recent snowfall concealing ablation areas and (3) the steep and rugged glacier terrain influencing the shortwave radiation budget in the Chhota Shigri glacier. Fortunately, increases in the number of satellite sensors and improvements in sensor performance should make it easier to acquire more precise glacier albedos in the near future.

It is worth mentioning that, in our method, the reconstructed annual glacier-wide mass balance did not represent the net mass balance but rather the mass balance fluctuation over the studied period. However, the results for the Chhota Shigri glacier show that the reconstructed mass balance from the albedo reproduced the observed glacier-wide mass balance well. This was mainly because ${\alpha }_i$ and $\partial b/\partial \alpha$ exhibited relatively large spatiotemporal variability across the glacier; therefore, our IAMB method could be used to estimate the mass balance fluctuation over the studied period. To calculate the net mass balances, an independent method capable of determining the total mass volume over the whole period should be considered; e.g., photogrammetry and glaciological monitoring.

5. Discussion

5.1. Sensitivity of the IAMB Mass Balance to $\partial b/\partial \alpha$

In this study, the mass balance–albedo gradient ($\partial b/\partial \alpha$) was calculated by multiplying the altitude–albedo gradient ($\partial z/\partial \alpha$) by the mean vertical mass–balance gradient ($\partial b/\partial z$). This means that $\partial b/\partial \alpha$ varies on a yearly basis and is specific to each glacier. Its variation is assumed to be due to various local controlling factors, such as the air temperature, precipitation, glacier type, glacier size, ELA, aspect, slope, exposure, etc. Therefore, it is essential to understand its effect on the reconstructed annual glacier-wide mass balance estimates ($b_i$). A sensitivity analysis of $\partial b/\partial \alpha$ with regard to the annual glacier mass balance estimates was carried out as follows.

Provided that $\overline{B}$ and ${\alpha }_i$ are invariable but $\partial b/\partial \alpha$ increases by $\Delta \partial b/\partial \alpha$, the newly calculated glacier mass balance for the $ith$ year ($b_i^{\prime }$) can be expressed as:

$$b_i^{\prime }=(\frac{\partial b}{\partial \alpha }+\Delta \frac{\partial b}{\partial \alpha })\ast ({\alpha }_i-{\alpha }_{eq})$$

(9)

The change in the annual glacier mass balance ($\Delta b_i^{\prime }$) resulting from $\Delta \partial b/\partial \alpha$ can be easily obtained as:

$$\Delta b_i^{\prime }=b_i^{\prime }-b_i=\Delta \frac{\partial b}{\partial \alpha }\ast ({\alpha }_i-\frac{1}{n}{\displaystyle \sum _{i=1}^n{\alpha }_i})$$

(10)

For the Chhota Shigri glacier, a systematic change of ±30% was applied to $\Delta \partial b/\partial \alpha$ to identify its effect on the annual mass balance from 2003 to 2014 (Figure 8). It can be observed that the newly estimated mass balance implies a different variation in each year when the $\partial b/\partial \alpha$ changes. Specifically, a high absolute value for $\Delta b_i^{\prime }$ exists for the years when the glacier albedo anomaly is relatively large, while a low absolute value of $\Delta b_i^{\prime }$ can be found for the years when the glacier albedo anomaly is relatively small. Moreover, positive or negative departures of $\Delta b_i^{\prime }$ depend on the product of $\Delta \partial b/\partial \alpha$ and the glacier albedo anomaly. For the largest value for the glacier albedo anomaly in 2010, the change in the newly estimated mass balance was merely −300.19 mm w.e., which was smaller than that for the derived mass balance uncertainties. This reveals that the impact of $\partial b/\partial \alpha$ is limited in cases where its variation is as high as ±30%. However, variations in $\partial b/\partial \alpha$ across different glaciers also need intensive investigations.

5.2. Sensitivity of Estimated Mass Balance to $\overline{B}$

The mean mass balance ($\overline{B}$) is another important input variable in the IAMB method that can be obtained based on traditional glaciological measurements or, more desirably, using a remotely based geodetic approach, such as airborne or satellite photogrammetry, radargrammetry or the gravimetric method. Accordingly, we carried out a simple sensitivity analysis with respect to $\overline{B}$ to evaluate its effect on the glacier mass balance.

Assuming that $\Delta \partial b/\partial \alpha$ and ${\alpha }_i$ are invariable but $\overline{B}$ increases by $\Delta \overline{B}$, the newly calculated glacier mass balance for the $ith$ year ($b_i^{″}$) becomes:

$$b_i^{″}=(\overline{B}+\Delta \overline{B})+\frac{\partial b}{\partial \alpha }\ast ({\alpha }_i-\frac{1}{n}{\displaystyle \sum _{i=1}^n{\alpha }_i})$$

(11)

The change in the annual glacier mass balance ($\Delta b_i^{″}$) resulting from $\Delta \overline{B}$ can be easily obtained as:

$$\Delta b_i^{″}=b_i^{″}-b_i=\Delta \overline{B}$$

(12)

In the above equation, it can be observed that the variation in the mean mass balance leads to an equal change ($\Delta \overline{B}$) in the glacier mass balance for every year. To the best of our knowledge, many existing methods can provide accurate and reliable estimates for $\overline{B}$ approximation, such as geodetic and gravimetric methods [49,50,51,52]. Therefore, the reconstructed annual glacier mass balance can mitigate the impact of $\overline{B}$ on our IAMB estimation.

5.3. Effect of Elevation on The $\partial z/\partial \alpha$

$\partial z/\partial \alpha$ was used to calculate the mass balance–albedo gradient ($\partial b/\partial \alpha$) and represents the transfer parameter from an ELA variation to annual mass balance fluctuation. Its value, therefore, can influence the estimated mass balance results. Figure 9 gives the temporal variation in $\partial z/\partial \alpha$ in the vicinity of the ELA (ELA ± 100 m, ELA ± 150 m, ELA ± 200 m, ELA ± 250 m, ELA ± 300 m) between 2001 and 2017. In general, the $\partial z/\partial \alpha$ values of these five altitude intervals are similar to each other, with a relative difference of less than 35% except for two special cases on 30 August 2015 (38%) and 15 September 2015 (43%). The average $\partial z/\partial \alpha$ estimates (R > 0.95) were 20.78 ± 12.88, 22.59 ± 13.33, 22.88 ± 12.69, 22.88 ± 11.80 and 23.35 ± 10.85 m/ (0.01 albedo), respectively. These low discrepancies indicate that the choice of altitude interval does not significantly affect the $\partial z/\partial \alpha$ estimation. Nonetheless, considering the practical ELA fluctuations at the Chhota Shigri glacier [38], we used altitude ranges of 200 m around the ELA for the $\partial b/\partial \alpha$ calculation.

Furthermore, in order to quantify the impact of the chosen $\partial z/\partial \alpha$ on the annual glacier mass balance reconstruction, we also applied our method with different $\partial z/\partial \alpha$ values. Considering the small discrepancies in absolute terms, the reconstructed mass balance results were not significantly different: the RMSEs of the calculated annual mass balance series were 281.79, 275.32, 274.91, 274.91 and 274.62 mm w.e. a⁻¹ for the five respective $\partial z/\partial \alpha$ values above. Therefore, the average altitude–albedo gradient across the ELA for a glacier appears to be an appropriate estimate in light of our results.

5.4. Why Not Choose the Minimum Albedo?

Previous studies have proven that the minimum melt season albedo for a whole glacier can realistically approximate the annual or seasonal mass balance for glaciers in the European Alps [18,21], New Zealand [20], the Himalayas [19] and the Tibetan Plateau [22]. However, the albedo method based on the minimum averaged albedo requires glaciological measurements with sufficient length for robust model calibration [20] and merely generalizes them to the nearby glaciated regions under similar climate regimes [22]. In other words, it remains challenging to find a universal model for the quantification of the relationship between minimum averaged albedo and surface glacier mass budget. This constrains the application of mass balance estimation for temporal and spatial extrapolations, especially for large-scale glaciers in remote areas without adequate in situ observations. Furthermore, cloud cover, which can hide the glacier surface in any season, is another inevitable obstacle for the calculation of the true minimum albedo value [27]. As a consequence, the application to large-scale glacier studies of the minimum albedo is inferior to the use of the ELA or AAR.

Differently from the minimum albedo, the averaged glacier-wide albedo can be obtained easily via remote sensing techniques, and the annual fluctuation in the average albedo is more accurate for larger glacier regions. Many experiments have confirmed that the average albedo is highly correlated with the surface glacier mass balance [24,25,26,27], thereby expanding the potential applications of the average albedo in annual mass balance reconstruction. Additionally, except for the usually adopted advanced very high resolution radiometer (AVHRR) and MODIS data, other new sensors, such as the Visible Infrared Imaging Radiometer Suite (VIIRS) instrument [53] and the Ocean and Land Color Instrument (OLCI) [54], can all feasibly provide global and consecutive albedo observations comparable with those of MODIS. This indicates that the average albedo can be further and sustainably applied for glacier mass balance estimation when the adaptation of the IAMB method to other sensors has been tested.

5.5. The IAMB Method: Future Prospects and Limitations

Long-term and continuous annual mass balance data are indispensable to understand the relationship between climate and glacier changes, as well as the contributions of glaciers to water resources. Unlike the existing albedo-based method, the IAMB method uses a model based on mean mass balance ($\overline{B}$) rather than the in situ annual mass balance. $\overline{B}$ can be quantified in terms of the glacier volume from the difference between glacier elevation data at different timescales, which is the so-called geodetic method. A wide variety of datasets, such as direct topographic measurements; topographical maps; aerial photographs; and satellite data from optical, radar or laser sensors, have been utilized for glacier elevation estimation [41,42,43,44]. The ability to leverage $\overline{B}$, which is widely available from geodetic methods, ensures that the IAMB method is not limited by the in situ mass balance observations. Therefore, this method has great potentialities for annual mass balance reconstruction with other glaciers with few or no field measurements. Moreover, the improved method has the advantage of reconstructing the glacier annual mass balance with a lower logistical burden and without additional precipitation or temperature information being required.

Nonetheless, the IAMB method needs the field mean vertical mass–balance gradient ($\partial b/\partial z$) for mass balance–albedo gradient ($\partial b/\partial \alpha$) estimation. Furthermore, $\partial b/\partial z$ is not fixed and changes with regional climate circumstances in different years and different glaciers. Fortunately, many previous studies have found that $\partial b/\partial z$ differs slightly for glaciers located under the same climate conditions [55]; therefore, $\Delta \partial b/\partial \alpha$ may not be larger than ±30% for neighbor glaciers lying in the same climatic zone, making it possible to apply the method in regional glacier monitoring. Furthermore, the IAMB method is built on the foundation of using the average albedo to capture the annual mass balance fluctuation signals, so it is not applicable to glaciers that are not suitable for average albedo estimations, such as glaciers with no accumulation zone in the ablation season, glaciers that are too small in size (e.g., less than one MODIS pixel) and surge glaciers. For debris-covered glaciers, it is best to exclude the debris areas before estimation of the average albedo and mass balance.

6. Conclusions

We proposed an improved albedo–mass balance (IAMB) method inspired by the ELA method to quantify the time series of annual glacier surface mass balance. We used the method to reconstruct the annual glacier mass balance based on the mean mass balance ($\overline{B}$) over the study period, the average summer albedo derived with the MODImLab algorithm and the mass balance–albedo gradient ($\partial b/\partial \alpha$). Estimation of the $\partial b/\partial \alpha$ from remote sensing data was also presented in detail. By using $\partial b/\partial \alpha$ and $\overline{B}$, the IAMB method greatly reduces reliance on observed mass balance measurements, thus underlining the high potential for the application of this method for glaciers without in situ mass balance data. In other words, our IAMB method can efficiency quantify the annual glacier-wide mass balance using remote sensing data if the average mass balance can be determined using geodetic or other methods.

Given the glacier size and the continuity of the field annual mass balance observations, this method has to date only been applied to the Chhota Shigri glacier for annual mass balance estimation over the period 2003–2014. The reconstructed mass budget correlated well with the observed yearly mass balance, with an average difference and root-mean-square error (RMSE) of −0.75 mm w.e. a⁻¹ and 274.91 mm w.e. a⁻¹, respectively, suggesting that the IAMB method is viable for glaciers with no or inadequate field measurements and promising for its ability to considerably increase the number of glacier mass balance series. Testing and verification of the IAMB method with other glaciers are still needed.

Nonetheless, it is worth noting that further work is needed to validate, calibrate and develop this method as follows: (1) further work on this method should preferentially be conducted with other glaciers situated in middle latitudes, polar areas and equatorial or monsoon-regime regions that have many measured or geodetic mass balance estimates to examine its applicability; (2) research should be undertaken on the spatial distribution of and change in $\partial b/\partial \alpha$ at regional scales (e.g., the Himalayas or the Tibetan Plateau) to verify the application potential in large regional glaciers; and (3) the method should be improved through model calibration or by adding seasonal surface albedos to simultaneously quantify the seasonal (summer and winter) mass balance components.