Annual mass-balance time series of Dongkemadi Glacier, 2000–20, from a linear albedo-based model using geodetic data and validated with the glaciological method

Abstract

Long-term time series of annual glacier mass balance is important for revealing a glacier's response to regional climate variations. However, for the Tibetan Plateau, time series of annual glacier mass balance with more than 10 consecutive years remains scarce due to the inaccessibility and harsh climate conditions. In this study, we established an albedo-based model to reconstruct annual glacier mass balance for 2000–20 over the Dongkemadi Glacier, based on multitemporal geodetic estimates and annual minimum glacier-wide mean surface albedos. Geodetic glacier mass-balance estimates for 2000–12, 2012–14 and 2014–18 were obtained by comparing glacier surface topographic data. Minimum surface albedos were obtained from the moderate-resolution imaging spectroradiometer daily snow albedo products for 2000–20. The estimated results are supported by the field measurements of annual glacier mass balance. During the early 2000s, we detected a relatively balanced state of glacier mass change, whereas a pronounced mass loss of more than −0.5 m w.e. was found for most years in the 2010s. By analyzing the regional climate variations with the ERA5-Land monthly averaged data, we found that the accelerated glacier mass loss can be attributed to both decreased annual snowfall and increased summer air temperature.

1. Introduction

Glacier mass balance is one of the key indicators of glacier status (WGMS, 2008; IPCC, 2013). During recent decades, an apparent glacier mass loss has been generally reported for glaciers on the Tibetan Plateau and its surroundings (Gardelle and others, 2013; Brun and others, 2017; Zhou and others, 2018b; Ke and others, 2020; Liu and others, 2020). However, most of the existing glacier mass-balance estimates are simply averaged values of specific mass change over a time frame of several years to several decades. In addition, because direct field measurements are usually time consuming and labor intensive, to date, time series of annual glacier mass-balance observations are scarce for the Tibetan Plateau and its surroundings.

Field measurement is the traditional method for obtaining annual glacier mass balance (WGMS, 2022). However, in situ observation over high-altitude glacierized regions is very difficult because of the sheer inaccessibility and severe weather conditions. For the tens of thousands of glaciers on the Tibetan Plateau and its surroundings, only tens of glaciers have field annual mass-balance measurements since the 1960s (Yao and others, 2012). In particular, annual glacier mass balance has been consecutively measured over only a few glaciers during a time period of more than 10 years (WGMS, 2022). The gravity recovery and climate experiment (GRACE) mission has provided monthly gravity data for the whole of the Earth with a spatial resolution of ~300 km since 2002 (Chen and others, 2008). However, the observations of the GRACE mission combine the mass changes of hydrology, the solid Earth and the cryosphere (Wang and others, 2013). On the Tibetan Plateau, glacier mass balance is almost impossible to distinguish from lake water changes when using GRACE observations.

By comparing glacier surface topographies in different times, glacier mass balance can be calculated by the geodetic method (Berthier and others, 2004; Lei and others, 2012; Gardelle and others, 2013; Liu and others, 2019). However, almost all of the estimated geodetic results are time-averaged, with a time interval of several years to several decades on the Tibetan Plateau. Zhou and others (2018b) calculated the mean glacier mass balance between the mid-1970s and 2000 for 11 glacierized regions on the Tibetan Plateau and its surroundings by subtracting the DEMs generated by KH-9 optical scenes from the shuttle radar topography mission (SRTM) DEM. Moreover, time-averaged glacier mass balances for 2000–2010s have been estimated for a large number of glaciers on the Tibetan Plateau through a comparison between the SRTM DEM and more recent DEMs (Gardelle and others, 2013; Liu and others, 2016; Brun and others, 2017; Berthier and Brun, 2019; Ke and others, 2020). Recently, by applying the DEMs generated from TanDEM-X bistatic interferometric synthetic aperture radar (InSAR) data, Liu and others (2019) calculated the annual glacier mass balance over the Puruogangri ice field. However, long-term time series of annual glacier mass balance cannot be obtained from the TanDEM-X DEMs, due to the limited data availability (Krieger and others, 2007).

Glacier surface albedo, which is equal to the ratio of reflected solar radiation to incoming solar radiation, is an important parameter for estimating glacier energy balance and mass balance (Henneman and Stefan, 1999; Sicart and others, 2008; Li and others, 2018). Over the Saint-Sorlin Glacier in the alps, Dumont and others (2012) detected a linear relationship between the field-measured annual glacier mass balance and the minimum glacier-wide averaged surface albedo in the summer season. Furthermore, by employing this linear correlation, Zhang and others (2018) successfully predicted annual mass fluctuations for glaciers with similar climate regimes. However, the albedo-based regression model is not suitable for most of the glacierized regions on the Tibetan Plateau, due to the paucity of field measurements of glacier mass balance.

In this study, we aimed to establish an albedo-based linear model using multitemporal geodetic glacier mass-balance estimates and corresponding time-averaged minimum glacier-wide mean surface albedos. The SRTM DEM and three recent DEMs generated from TanDEM-X bistatic InSAR data and Chinese ZiYuan-3 stereo images were employed to estimate the geodetic glacier mass balance. Glacier surface albedo was extracted from the moderate-resolution imaging spectroradiometer (MODIS) daily albedo products released by the National Snow and Ice Data Center (NSIDC). Annual glacier mass balance for 2000–20 over the Dongkemadi Glacier was estimated using the established albedo-based linear model and was validated using field measurements. The influence of regional climate variations on the glacier mass balance was also investigated using annual snowfall and summer air temperature data retrieved from the ERA5-Land monthly averaged data.

2. Study area

The Dongkemadi Glacier is located on the interior region of the Tibetan Plateau, with the central coordinates of 33°05′ N and 92°04′ E (Fig. 1). The area covered by the Dongkemadi Glacier has been ~16 km² during recent years (Guo and others, 2015). This glacier is composed of the Da Dongkemadi and Xiao Dongkemadi (labeled as XD in Fig. 1) glaciers. The glacier mass balance of the Xiao Dongkemadi Glacier has been measured using the glaciological method since 1989, and continuous measurements are still being recorded (Yao and others, 2012; Zhang and others, 2013). Like other glaciers on the inner Tibetan Plateau, the climate of this glacier is mainly determined by continental climatic conditions (Yao and others, 2012). Therefore, the Dongkemadi Glacier is of continental type (Zhang and others, 2013; Shi and others, 2016).

Figure 1. Geographic map of the Dongkemadi Glacier (dark red line) and its surroundings. Contours (m a.s.l.) are drawn at intervals of 50 m. Field measurements of the glacier mass balance have been conducted over the Xiao Dongkemadi Glacier (labeled as XD). Thirty-five ablation stakes are positioned in regions A and B (Shi and others, 2016).

3. Data and methods

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, 2003).

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 (2019). 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, 2018a). 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, 2007). 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, 2010). 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, 2015; Li and others, 2018), 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 (2013). 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 (2019) 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 (2011) 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, 2004; Gardelle and others, 2012). 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, 2014), 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 (2017), 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, 2013; Ke and others, 2020; Liu and others, 2020). 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, 2013). 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 (2013) 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, 2018). Therefore, surface albedo is a significant parameter for estimating glacier energy balance and mass balance (Huintjes and others, 2015). 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, 2012; Brun and others, 2015; Williamson and others, 2020; Banerjee and others, 2022). For glaciers in High Mountain Asia, the lowest surface albedo in a year is generally detected in the summer months (Ming and others, 2015). 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, 2012; Brun and others, 2015; Davaze and others, 2018; Banerjee and others, 2022). 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, 2012; Sirguey and others, 2016). Here, the error of the MOD10A1 or MYD10A1 data was assumed to be ±10% of the maximum value in a year (Stroeve and others, 2006; Wu and others, 2015). 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, 2019). Here, we assumed that the error of the glacier area was ±5% (Paul and others, 2013) and the error of the conversion factor was ±7% for the geodetic estimate in 2000–12 (Huss, 2013). 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 (2019). 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.

4. Results

4.1 Glacier elevation change and geodetic glacier mass balance

By comparing the DEMs for different times, the results of surface elevation changes for 2000–18 were extracted for the Dongkemadi Glacier and its surroundings (Fig. S1). The histograms of the elevation differences over the off-glacier regions demonstrate a nearly Gaussian distribution (Fig. S2), indicating that errors generally are random. A pronounced surface thinning was detected over the glacier tongue region during the studied period (Fig. 3). Moreover, a slight surface thickening was observed for most of the upper zone for 2000–12 (Fig. 3a) and 2014–18 (Fig. 3c); however, this high-altitude zone generally experienced slight surface thinning during 2012–14 (Fig. 3b). In general, the spatial pattern of the glacier elevation changes over the Da Dongkemadi Glacier is similar to that over the Xiao Dongkemadi Glacier. Overall, an upward decrease in surface thinning and an upward increase in surface thickening was found for the tongue region and upper zone of the Dongkemadi Glacier, respectively.

Figure 3. Measured glacier elevation change (m) over the Dongkemadi Glacier during the time periods of 2000–12 (a), 2012–14 (b) and 2014–18 (c). Black arrows indicate the direction of ice flow.

Glacier-wide mean surface elevation change values of −0.36 ± 0.11, −1.48 ± 0.43 and −0.77 ± 0.29 m a⁻¹ were retrieved over the Dongkemadi Glacier for the time periods of 2000–12, 2012–14 and 2014–18, respectively (Table 2). By employing the conversion factors suggested by Huss (2013), we estimated a serious glacier mass loss of −1.16 ± 0.38 m w.e. a⁻¹ for 2012–14 on the Dongkemadi Glacier (Table 2). Glacier mass loss was relatively modest during the period of 2014–18 (−0.65 ± 0.27 m w.e. a⁻¹). The least glacier mass loss during the studied time period was extracted for 2000–12 (−0.31 ± 0.09 m w.e. a⁻¹). Through a comparison of the estimated geodetic glacier mass balance for 2000–12 to that for 2012–14 and 2014–18, we found that the Dongkemadi Glacier likely experienced much more serious mass loss in the 2010s than in the first decade of the 21st century.

Table 2. Observed glacier elevation change and geodetic mass balance in the time periods of 2000–12, 2012–14 and 2014–18

Numbers listed in this table are rounded to two decimal places.

4.2 Temporal variation of the AMGA

The AMGA of the Dongkemadi Glacier for 2000–20 is illustrated in Figure 4. The maximum value of the AMGA (0.555) was retrieved for 2008, and an AMGA of larger than 0.5 was also calculated for 2001 (0.513) and 2003 (0.521). Moreover, we extracted the minimum value of the AMGA (0.290) for 2013, and the AMGA was ~0.3 for 2006 (0.307), 2010 (0.311) and 2016 (0.300). The mean value of the AMGA was 0.446 for 2000–10, whereas the average AMGA dropped to 0.374 in the time period of 2011–20. In general, we detected a decreasing trend for the AMGA over this glacier for 2000–20 (Fig. 4). This likely indicates an accelerating rate of glacier melt during the studied period, because more and more solar radiation should be absorbed for heating a glacier with the decreased surface albedo. Furthermore, the accelerating rate of glacier melt is supported by the temporal variation of the estimated geodetic glacier mass balance for 2000–18.

Figure 4. Extracted AMGA over the Dongkemadi Glacier for 2000–20.

4.3 Annual glacier mass-balance time series

The annual glacier mass-balance time series for the Dongkemadi Glacier for 2000–20 were extracted by applying the AMGA time series and the parameter scenario which was calculated with the estimated geodetic glacier mass balance and mean AMGA during the same time period. As listed in Table 3, three parameter scenarios were calculated in this study, and thus the annual glacier mass-balance time series were calculated using parameter scenarios t12, t13 and t23. In general, the glacier mass-balance estimates for these three parameter scenarios are similar, especially for the years with a moderate glacier mass loss (Table 4). For example, we obtained glacier mass-balance values of −0.28 ± 0.28, −0.29 ± 0.21 and −0.27 ± 0.49 m w.e. in 2017 for parameter scenarios t12, t13 and t23, respectively.

Table 3. Calculated parameters of the linear albedo-based model

Numbers listed in this table are rounded to two decimal places.

Table 4. Estimated AMGA and glacier mass-balance time series during 2000–20

Numbers listed in this table are rounded to two decimal places.

The final glacier mass balance for a certain year was estimated by averaging the three parameter scenarios (t12/t13/t23). The results of the annual glacier mass-balance time series indicate that a relatively balanced state of glacier mass change was mainly extracted for the early 2000s, while a pronounced glacier mass loss of more than −0.5 m w.e. was estimated for most years in the 2010s (Table 4). Furthermore, the Dongkemadi Glacier experienced negative mass changes for all the years in the 2010s, and a positive glacier mass change was only detected for three years in the 2000s (2001, 2003 and 2008). Consequently, the modeled annual glacier mass-balance time series clearly reveal an accelerating trend of glacier mass loss in 2000–20.

5. Discussion

5.1 Comparison with field measurements of glacier mass balance

In order to validate the geodetic estimates and modeled annual glacier mass-balance time series, we compared the estimated and modeled results to field measurements of glacier mass balance. The traditional glaciological method has been used to measure the annual glacier mass balance of the Xiao Dongkemadi Glacier since 1989 (Yao and others, 2012; Zhang and others, 2013). Currently, the most recent Fluctuations of Glaciers Database (version 2022) released by the World Glacier Monitoring Service (WGMS) only provides the measured annual glacier mass balance of the Xiao Dongkemadi Glacier for 1989–2010 (WGMS, 2022). Taking into consideration the studied period (2000–20), we employed the field-measured glacier mass balance for 2000–10 to validate the estimated and modeled results in this study.

As illustrated in Figure 5, both the geodetic estimate and the modeled annual glacier mass-balance time-series estimates are supported by the field measurements. The mean value of field-measured glacier mass balance in 2000–10 is −0.36 m w.e. a⁻¹, which is generally in agree with our geodetic estimate in 2000–12 (−0.31 ± 0.09 m w.e. a⁻¹). A subtle difference between the estimated results and field measurements is apparent for the years with moderate glacier mass loss (~−0.50 m w.e.) and some years with a relatively balanced glacier mass change. For example, a glacier mass balance of −0.59 m w.e. was obtained from field measurements in 2007, and the modeled glacier mass balance is −0.54 ± 0.32 m w.e. for this year. Furthermore, a quantitative comparison of the modeled glacier mass-balance time series with the field measurements was conducted using the correlation coefficient (R) and the significance level (P). A strong positive correlation (R = 0.92) is apparent between modeled results and field measurements at the 99% confidence interval (P < 0.01). In addition, mean absolute error of 0.16 m w.e. was estimated by averaging the differences between the modeled results and in situ observations.

Figure 5. Modeled annual glacier mass balance (m w.e.) for 2000–20. The red dotted line represents field measurements of glacier mass balance (m w.e.) in 2000–10. Dark red lines represent the annual mean results of geodetic estimates (m w.e. a⁻¹).

5.2 Application of previously published geodetic estimates

During recent years, some regional or global geodetic glacier mass-balance estimates have been published (Brun and others, 2017; Shean and others, 2020; Hugonnet and others, 2021). Over the Dongkemadi Glacier, glacier mass balances of −0.24 m w.e. a⁻¹ (2000–09) and −0.74 m w.e. a⁻¹ (2010–19) were estimated from the elevation change map in Hugonnet and others (2021). By using these geodetic estimates, we calculated parameter scenario (scenario Hugonnet) of 4.94 (m) and −2.51 (n) for the linear albedo-based model, and estimated annual glacier mass-balance time series (Table S1). Moreover, we re-calculated the final glacier mass-balance estimates by averaging the results of parameter scenarios t12, t13, t23 and Hugonnet (Table 4). The correlation coefficient between re-calculated estimates and field measurements is still 0.92, because the correlation coefficient is determined by the relationship between the AMGAs and field measurements. The mean absolute error of the re-calculated estimates (scenario t12/t13/t23/Hugonnet) was estimated to be 0.14 m w.e., which is smaller than that of scenario t12/t13/t23 (0.16 m w.e.). This indicates that a more accurate annual glacier mass-balance time series is probably estimated by adding the results modeled with the geodetic estimates in Hugonnet and others (2021).

It is noteworthy that, only the geodetic estimates of the average value of annual glacier mass balance can be used for the construction of the linear albedo-based model. In this study, we didn't employ the results published by Brun and others (2017) and Shean and others (2020), because these geodetic estimates were obtained by calculating the linear trend of elevation time series.

5.3 Limitations of the MODIS daily snow albedo data

As shown in Eqn (1), in addition to the calculated parameter scenario, the accuracy of the modeled annual glacier mass-balance time series was also determined by the used MODIS daily snow albedo data. The glacier mass balance in a certain year was calculated using the albedo-based model and the AMGA, which was obtained by detecting the minimum daily glacier-wide albedo between June and September. In order to accurately extract daily glacier-wide surface albedo, the MOD10A1 and MYD10A1 snow albedo products on the same day were merged to one image using the principles proposed in Zhang and others (2018). As illustrated in Figure 2, the coverage of the merged albedo image is nearly 100%, whereas void regions were clearly detected for both the MOD10A1 and MYD10A1 albedo data. However, the MYD10A1 daily snow albedo product is absent for 2000 and 2001, due to the subsequently launched Aqua satellite (Hall and Riggs, 2007). It means that the MOD10A1 snow albedo data were utilized as a substitute for the merged imagery for all the days in 2000 and 2001.

The absence of the MYD10A1 snow albedo data likely lead to the apparent error for the extracted daily glacier-wide average albedo in 2000 and 2001. As shown in Figure 2, if only the MOD10A1 data are available, the calculated average value of glacier-wide albedo should be underestimated for this day, because of a pronounced data gap for the MOD10A1 image over the upper accumulation zone (Fig. 2a). Moreover, when the void regions of the MOD10A1 image are mainly detected over the lower ablation zone, the glacier-wide average snow albedo will be overestimated. The underestimation or overestimation of daily glacier-wide average albedo likely causes the pronounced error for the extracted AMGAs. Consequently, the apparent difference between the modeled glacier mass balances and field measurements in 2000 and 2001 (Fig. 5) is mainly attributed to the absence of the MYD10A1 albedo products for these two years. However, the linear relationship between the AMGA and glacier mass balance is nearly not impacted by the pronounced errors of the AMGAs in 2000 and 2001. A strong positive linear relationship was detected between the AMGA and the field measured glacier mass balance for the three time periods of 2000–10, 2001–10 and 2002–10 (Fig. 6).

Figure 6. Relationship between the AMGA and the field measured glacier mass balances for the time periods of 2000–10 (a), 2001–10 (b) and 2002–10 (c).

For the years of both the MYD10A1 and MOD10A1 albedo products, data gap is possibly detected for the merged image (e.g. the white pixels in Fig. 2c), which also possibly lead to an error for the extracted AMGA and thus for the modeled annual glacier mass balance. As proposed by Banerjee and others (2022), snowline altitude and normalized difference snow index are also effective remote-sensing data for estimating annual glacier mass balances. Therefore, the synergistic use of these satellite remote-sensing products will be a promising method for improving the accuracy of the albedo-based modeled results.

5.4 Influence of regional climate variations

Glacier mass balance is equal to the total sum of all the accumulation and ablation across the entire glacier. For a certain year, glacier accumulation is made up of snowfall and refreezing water, and glacier ablation mainly consists of snow/ice melt and sublimation. In general, an apparently temporal change of annual snowfall and snow/ice melt is apparent for glaciers on the interior Tibetan Plateau, whereas the levels of annual refreezing of water and sublimation are relatively stable (Huintjes and others, 2015; Li and others, 2018). For the Dongkemadi Glacier, ablation occurs mainly in the summer months. Therefore, in order to analyze the effect of regional climate variations on the glacier mass balance during the studied period, we extracted summer air temperature (June, July and August) and annual snowfall for 1979–2021 from the ERA5-Land monthly averaged data.

As illustrated in Figure 7a, a pronounced decreasing trend is apparent for the annual snowfall during 1979–2021. During 1979–2000, the number of years with snowfall of more than 350 mm w.e. is 14; however, this drops sharply to 4 for 2001–21. The decrease in annual snowfall is partly related to the observed accelerating trend of glacier mass loss, because glacier accumulation is mainly determined by snowfall. Furthermore, an apparent increasing trend is apparent for the monthly mean air temperature in June, July and August (Fig. 7b). During the early 1980s, the mean air temperature during the summer season was ~3°C, whereas the average air temperature quickly increased to ~5°C in the late 2010s. Importantly, the rapid increase in summer air temperature likely caused the pronounced rise in snow/ice melt for 1979–2021.

Figure 7. Annual snowfall (a) and monthly mean air temperature in summer (b) for 1979–2021, obtained from the ERA5-Land monthly averaged data.

Simple linear regression was further used to evaluate the influence of regional climate variations on glacier mass balance by computing the correlation coefficient (R) and the significance level (P). A strong positive correlation (R = 0.75) is apparent between glacier mass balance and annual snowfall at the 99% confidence interval (P < 0.001) (Fig. 8a). In addition, it is also apparent that glacier mass balance is negatively correlated with summer mean air temperature (R = 0.76) at the 99% confidence interval (P < 0.001) (Fig. 8b). Overall, both annual snowfall and summer mean air temperature were significant controls on glacier mass balance for 2000–20. Consequently, the observed accelerating trend of glacier mass loss was likely due to both the decreased annual snowfall and the increased air temperature in the summer months during the studied period.

Figure 8. Correlation between annual glacier mass balance and regional climate variations (annual snowfall (a) and summer mean air temperature (b)) for the Dongkemadi Glacier for 2000–20. R is the correlation coefficient, and P is the significance level.

6. Conclusions

In this study, we extracted annual glacier mass-balance time series for 2000–20 with an albedo-based linear model over the Dongkemadi Glacier on the Tibetan Plateau. The linear model was established by applying bitemporal geodetic estimates of glacier mass balance and the corresponding time-averaged AMGA. The geodetic glacier mass balance was calculated through the comparison of glacier topographic data in different times for the periods of 2000–12, 2012–14 and 2014–18. The AMGA time series for 2000–20 was estimated from the MODIS daily surface albedo products. Three parameter scenarios for the linear model were calculated with the retrieved geodetic glacier mass balance and the time-averaged AMGA. The final glacier mass balance for a certain year was estimated by averaging the results obtained from all the three parameter scenarios. Overall, the modeled annual glacier mass-balance time series is generally supported by field measurements.

In general, we detected an apparent accelerating rate of glacier mass loss over the Dongkemadi Glacier during 2000–20. This glacier experienced a balanced mass change (slight mass loss or gain) in the early 2000s, whereas pronounced glacier mass loss of more than −0.5 m w.e. was found for most years in the 2010s. By employing the ERA5-Land monthly averaged data, we found decreased annual snowfall and increased summer air temperature for 1979–2021. Moreover, both the annual snowfall and summer mean air temperature were found to be significant controls for glacier mass balance during 2000–20. Consequently, the accelerating rate of glacier mass loss during the studied period was determined by this temporal feature of the regional climate variations.