3.1 Datasets used in this study

In order to estimate multitemporal glacier mass-balance results with the geodetic method, we used the X-band SRTM DEM and the generated glacier surface topographic data of Chinese ZiYuan-3 stereo optical images and TanDEM-X bistatic InSAR scenes. Two radar sensors of different wavelengths (X-band and C-band) were equipped on the SRTM mission. In this study, the SRTM-X DEM was employed, because both the horizontal and vertical accuracies of the SRTM-X DEM are better than those of the SRTM-C DEM (Rabus and others, Reference Rabus, Eineder, Roth and Bamler2003).

Two pairs of TanDEM-X bistatic InSAR images collected in February 2012 and February 2014 were employed to extract the glacier surface topography. We acquired the bistatic InSAR data from the German Aerospace Center (DLR) with the format of co-registered single look slant range complex (CoSSC). The spatial resolution of the used CoSSC data was ~2 m in both the azimuth and ground range. In this study, we processed the TanDEM-X InSAR images using the differential interferometric method developed by Liu and others (Reference Liu2019). Moreover, one pair of ZiYuan-3 stereo optical images acquired in December 2017 was used to generate the glacier surface DEM. The spatial resolution of the used ZiYuan-3 images was 2.1, 2.5 and 2.5 m for the nadir, forward and backward images, respectively (Zhou and others, Reference Zhou, Tang, Guo, Wang and Fan2018a). In addition, the ZiYuan-3 DEM was generated using Space Data Processor software from the Land Satellite Remote Sensing Application Center of the Ministry of Natural Resources of the People's Republic of China.

Glacier surface albedo was generated using the MODIS daily snow albedo products of MOD10A1/MYD10A1 (tile: h25v05). The daily snow albedo products provide level 3 global black-sky albedo for the cloud-free coverage at a 500 m pixel resolution, and are available for free from the NSIDC (Hall and Riggs, Reference Hall and Riggs2007). Moreover, the albedo products are in Hierarchical Data Format-Earth Observing Systems (HDF-EOS) format with a sinusoidal map projection, and are stored as integer type, ranging from 0 to 100 (%) (Schaaf and others, Reference Schaaf, Liu, Gao and Strahler2010). Considering the second MODIS instrument onboard the Aqua platform launched in 2002, the MOD10A1 product for 2000–20 and the MYD10A1 product for 2002–20 were employed in this study.

Given that glacial accumulation and ablation are mainly affected by snowfall and air temperature, respectively (Huintjes and others, Reference Huintjes, Neckel, Hochschild and Schneider2015; Li and others, Reference Li, Yao, Yang, Yu and Zhu2018), we used air temperature and snowfall data to evaluate the effect of regional climate variation on the annual glacier mass balance. The snowfall and 2 m air temperature obtained from the ERA5-Land monthly averaged data for 1979–2021 were used in this study.

3.2 Geodetic glacier mass-balance calculation

Geodetic glacier mass balance was calculated by converting the glacier-wide mean elevation change with the conversion factors suggested by Huss (Reference Huss2013). The glacier surface elevation changes were detected from a comparison of two sets of glacier surface topographic data from different times. Considering the penetration of radar waves into snow or ice, the X-band radar penetration depth needed to be corrected before the process of DEM differencing between TanDEM-X DEM and ZiYuan-3 DEM. In this study, the function between the X-band radar penetration depth and altitude in Liu and others (Reference Liu2019) was used to estimate the penetration depth for every pixel of the TanDEM-X DEM in 2014. In addition, in order to minimize the influence of geometric errors between each pair of DEMs, the universal co-registration method proposed by Nuth and Kaab (Reference Nuth and Kaab2011) was employed to accurately match the two used DEMs. Note that we did not fill the data gaps of the measured glacier surface elevation changes.

The glacier-wide mean value of surface elevation changes was calculated at a 50 m altitude band, according to the assumption that pixels at an altitude interval usually experience similar elevation changes (Berthier and others, Reference Berthier, Arnaud, Baratoux, Vincent and Rémy2004; Gardelle and others, Reference Gardelle, Berthier and Arnaud2012). First, the results of the DEM differencing were masked by the glacier boundaries to extract the surface elevation changes of the glacierized regions. Here, we employed the glacier outlines from the Randolph Glacier Inventory version 6.0 (RGI 6.0), which was released in July 2017 (RGI Consortium, 2017). Taking into consideration the pronounced glacier retreat during recent decades (Wei and others, Reference Wei2014), we modified the terminus locations in RGI 6.0 using cloud-free Landsat optical images acquired on a similar date to the older DEM for each result of the DEM differencing.

The annual rate of glacier-wide mean surface elevation change was computed by dividing the estimated mean value by the integer number of years for the study period. As the used DEM pairs were not always of the same or a similar date of the year, the systematic bias of the seasonal glacier surface elevation change needed to be corrected first. In this study, we employed the corrections for seasonal glacier elevation changes in Wang and others (Reference Wang, Yi, Chang and Sun2017), which were estimated by computing the elevation differences between the ICESat observations acquired from the approximately October campaigns and those obtained from the approximately March campaigns.

The conversion factor of 850 kg m⁻³ is generally used for calculating the geodetic glacier mass balance (Gardelle and others, Reference Gardelle, Berthier, Arnaud and Kääb2013; Ke and others, Reference Ke, Song, Yong, Lei and Ding2020; Liu and others, Reference Liu, Jiang, Wang, Ding and Xu2020). However, this conversion factor is recommended in the case of the time interval of the two used topographic datasets being larger than 5 years (Huss, Reference Huss2013). Here, we used this conversion factor for the time periods of 2000–12 and 2014–18, although the time interval for 2014–18 is 4 years. For the glacier-wide mean elevation change for 2012–14, the conversion factor of 784 kg m⁻³ suggested for the time interval of 2 years by Huss (Reference Huss2013) was employed to estimate the glacier mass balance.

3.3 Glacier-wide surface albedo extraction

Surface albedo, which is equal to the ratio of reflected solar radiation to incoming solar radiation, is a key part of glacier ice and snowmelt because the absorbed solar radiation is an important energy source for heating glaciers (Li and others, Reference Li, Yao, Yang, Yu and Zhu2018). Therefore, surface albedo is a significant parameter for estimating glacier energy balance and mass balance (Huintjes and others, Reference Huintjes, Neckel, Hochschild and Schneider2015). During recent years, many studies have proposed that the minimum value of glacier surface albedo in the melt season is an important indicator of annual glacier mass balance (Dumont and others, Reference Dumont2012; Brun and others, Reference Brun2015; Williamson and others, Reference Williamson, Copland, Thomson and Burgess2020; Banerjee and others, Reference Banerjee, Singh and Sheth2022). For glaciers in High Mountain Asia, the lowest surface albedo in a year is generally detected in the summer months (Ming and others, Reference Ming2015). Here, we selected the MODIS daily snow albedo products (MOD10A1 and MYD10A1) during the ablation season from 1 June to 30 September every year for the minimum glacier surface albedo calculation.

In order to calculate the averaged glacier-wide surface albedo for a certain day, the MOD10A1 and MYD10A1 albedo images were first subject to the processes of map reprojection, format conversion and spatial sub-setting using the MODIS Reprojection Tool. The surface albedo data of the Dongkemadi Glacier were extracted from the preprocessed images with the modified RGI 6.0 glacier outlines. Moreover, the MOD10A1 and MYD10A1 albedo products on the same day were merged to one scene by averaging these two surface albedo values for every pixel separately. If the MOD10A1 or MYD10A1 albedo data are absent for a pixel due to cloud coverage, the existed surface albedo value was given to the merged scene (Fig. 2). For example, the absence of albedo data was detected for some pixels in the MOD10A1 image (Fig. 2a) and the MYD10A1 image (Fig. 2b) on 1 June 2010. By merging the MOD10A1 and MYD10A1 albedo products, the absence of surface albedo data at the scale of pixel has been greatly alleviated (Fig. 2c).

Figure 2. Data void regions (white pixels) for the MOD10A1 image (a), the MYD10A1 image (b) and the merged image (c) on 1 June 2010.

It is noteworthy that, both the MOD10A1 or MYD10A1 albedo products were absent for some days in 2000 (14 d) and 2001 (18 d), because the original MODIS hyperspectral remote-sensing images were not available. Therefore, the averaged glacier-wide surface albedo cannot be directly calculated from surface albedo products for these days. In this study, these data gaps of the averaged glacier-wide surface albedo were filled by using the time-averaged values of a 10 d window. Moreover, the mean filter in time with a 10 d moving window was also employed for other days to decrease the influence of errors. Finally, the annual minimum-averaged glacier-wide albedo (AMGA) was obtained by detecting the minimum value of the filled and filtered daily glacier-wide mean albedo in a year.

3.4 Annual glacier mass-balance time-series estimation

Previous studies have found that a strong linear relationship exists between the AMGA and annual glacier mass balance (Dumont and others, Reference Dumont2012; Brun and others, Reference Brun2015; Davaze and others, Reference Davaze2018; Banerjee and others, Reference Banerjee, Singh and Sheth2022). Therefore, if the AMGA of an individual glacier has been extracted from the MODIS daily snow albedo products, annual glacier mass balance for the same year can be calculated with the albedo-based model, as follows:

(1)$$b = m \times \alpha + n$$

where m and n are constant parameters; and b and α are, respectively, the glacier mass balance and AMGA in a certain year. Furthermore, according to the characteristics of a linear equation with one variable, the average value of the annual glacier mass balance for the same time period was calculated as follows:

(2)$$\bar{b} = m \times \bar{\alpha } + n$$

Therefore, when the average value of the annual glacier mass balance and the mean AMGA were estimated for two different time periods, it was possible to calculate the two parameters of m and n. For example, we selected the geodetic glacier mass balance and mean AMGA for 2000–12 and 2012–14, and Eqn (2) could thus be rewritten as follows:

(3)$$\bar{b}_{t1} = m \times \bar{\alpha }_{t1} + n$$

(4)$$\bar{b}_{t2} = m \times \bar{\alpha }_{t2} + n$$

where $\bar{b}_{t1}$ and $\bar{b}_{t2}$ are the geodetic glacier mass balances for 2000–12 and 2012–14, respectively; and $\bar{\alpha }_{t1}$ and $\bar{\alpha }_{t2}$ are the mean AMGAs for the two time intervals. According to Eqns (3)–(4), parameters m and n can be calculated as follows:

(5)$$m = \displaystyle{{{\bar{b}}_{t1}-{\bar{b}}_{t2}} \over {{\bar{\alpha }}_{t1}-{\bar{\alpha }}_{t2}}}$$

(6)$$n = \displaystyle{{{\bar{b}}_{t2} \times {\bar{\alpha }}_{t1}-{\bar{b}}_{t1} \times {\bar{\alpha }}_{t2}} \over {{\bar{\alpha }}_{t1}-{\bar{\alpha }}_{t2}}}$$

In this study, by comparing the glacier surface topographic data for different times, we calculated the average glacier mass balance for the time periods of 2000–12, 2012–14 and 2014–18. Mean AMGAs of the three time periods were also obtained from the extracted AMGA time series. As listed in Table 1, three parameter scenarios were calculated with the retrieved geodetic glacier mass balance and time-averaged AMGA. Furthermore, with the extracted AMGA time series and the calculated values of constant parameters m and n, the annual glacier mass-balance time series were then computed using Eqn (1). Overall, the final glacier mass balance for a certain year was estimated by averaging the results obtained from all the three parameter scenarios (t12, t13 and t23).

Table 1. Parameter scenarios calculated with the geodetic glacier mass balance and average AMGA over the different time periods

3.5 Uncertainty assessment

The uncertainties of the estimated annual glacier mass-balance time series were assessed by applying the standard law of error propagation. According to Eqn (1), the uncertainty of the estimated glacier mass balance is caused by the error of the AMGA in the same year and the calculated constant parameters of m and n, and can be obtained with the following equation:

(7)$$\sigma _{b_i}^2 = \left({\displaystyle{{\partial b_i} \over {\partial \alpha_i}}} \right)^2\sigma _{\alpha _i}^2 + \left({\displaystyle{{\partial b_i} \over {\partial m}}} \right)^2\sigma _m^2 + \sigma _n^2 $$

The error of the AMGA for a certain year ($\sigma _{\alpha _i}$) was computed by dividing the error of the MOD10A1 or MYD10A1 data by the RMS of the number of glacier pixels (Dumont and others, Reference Dumont2012; Sirguey and others, Reference Sirguey2016). Here, the error of the MOD10A1 or MYD10A1 data was assumed to be ±10% of the maximum value in a year (Stroeve and others, Reference Stroeve, Box and Haran2006; Wu and others, Reference Wu2015). The errors of the extracted parameters (σ _m and σ _n) were calculated with the standard law of error propagation. Specifically, according to Eqns (5)–(6), the errors of parameter scenario t12 can be calculated by the following equations:

(8)$${\sigma _m^2 = \left({\displaystyle{{\partial m} \over {\partial {\bar{b}}_{t1}}}} \right)^2\sigma _{{\bar{b}}_{t1}}^2 + \left({\displaystyle{{\partial m} \over {\partial {\bar{b}}_{t2}}}} \right)^2\sigma _{{\bar{b}}_{t2}}^2 + \left({\displaystyle{{\partial m} \over {\partial {\bar{\alpha }}_{t1}}}} \right)^2\sigma _{{\bar{\alpha }}_{t1}}^2 + \left({\displaystyle{{\partial m} \over {\partial {\bar{\alpha }}_{t2}}}} \right)^2\sigma _{{\bar{\alpha }}_{t2}}^2 }$$

(9)$${\sigma _n^2 = \left({\displaystyle{{\partial n} \over {\partial {\bar{b}}_{t1}}}} \right)^2\sigma _{{\bar{b}}_{t1}}^2 + \left({\displaystyle{{\partial n} \over {\partial {\bar{b}}_{t2}}}} \right)^2\sigma _{{\bar{b}}_{t2}}^2 + \left({\displaystyle{{\partial n} \over {\partial {\bar{\alpha }}_{t1}}}} \right)^2\sigma _{{\bar{\alpha }}_{t1}}^2 + \left({\displaystyle{{\partial n} \over {\partial {\bar{\alpha }}_{t2}}}} \right)^2\sigma _{{\bar{\alpha }}_{t2}}^2 }$$

where $\sigma _{{\bar{b}}_{t1}}$ and $\sigma _{{\bar{b}}_{t2}}$ are the errors of the geodetic glacier mass balance during 2000–12 and 2012–14, respectively; and $\sigma _{{\bar{\alpha }}_{t1}}$ and $\sigma _{{\bar{\alpha }}_{t2}}$ are the errors of the mean AMGA in these time periods. According to the calculation of geodetic glacier mass balance, its uncertainty is mainly determined by the errors of the glacier area, mean surface elevation change and conversion factor (Liu and others, Reference Liu2019). Here, we assumed that the error of the glacier area was ±5% (Paul and others, Reference Paul2013) and the error of the conversion factor was ±7% for the geodetic estimate in 2000–12 (Huss, Reference Huss2013). Since period is shorter than 5 years, an error of ±14% was assumed for the conversion factors in 2012–2014 and 2014–2018. The errors of geodetic glacier mass balance and glacier-wide mean elevation change were estimated with the standard law of error propagation, and details can be found in Liu and others (Reference Liu2019). The error of the mean AMGA for a certain time period ($\sigma _{\bar{\alpha }}$) was calculated with the following equation:

(10)$$\sigma _{\bar{\alpha }}^2 = \displaystyle{{\sum {\sigma _{\alpha _i}^2 } } \over {k \times ( k-1) }}$$

where k is the number of years for the time period. It is noteworthy that the error of final glacier mass balance for a certain year was calculated as the average error of the three parameter scenarios, because the used geodetic results for calculating parameter scenario are partly the same.