1. Introduction
The Tibetan Plateau accommodates a large number of glaciers, and the meltwater from these glaciers provides valuable freshwater for the people living downstream (Pritchard, Reference Pritchard2019; Immerzeel and others, Reference Immerzeel2020; Yao and others, Reference Yao2022). Most of the glaciers across the Tibetan Plateau have experienced shrinkage and mass loss in recent decades (Yao and others, Reference Yao2012; Brun and others, Reference Brun, Berthier, Wagnon, Kääb and Treichler2017; Zhou and others, Reference Zhou, Li, Li, Zhao and Ding2018; Shean and others, Reference Shean, Bhushan, Montesano, Rounce, Arendt and Osmanoglu2020; Zemp and others, Reference Zemp2020), thereby highlighting the pressing need to accurately predict the timing of ‘peak water’ as global warming continues to have an ever-increasing impact on these crucial freshwater reservoirs (Gao and others, Reference Gao, Li, Duan, Ren, Meng and Pan2018; Huss and Hock, Reference Huss and Hock2018). Glacier mass balance and ice volume are two key components for predicting the timing of peak water (Huss and Hock, Reference Huss and Hock2018; Welty and others, Reference Welty2020).
Glacier mass balance can be estimated from glaciological or geodetic methods or on the basis of numerical modeling. The glaciological method is based primarily on stake and pit measurements, which can quantify mass balance at a high temporal resolution (Zemp and others, Reference Zemp2020). The geodetic method is an effective means for estimating the glacier mass balance at multiple scales (spatial and temporal) using remote-sensing-based elevation data (Brun and others, Reference Brun, Berthier, Wagnon, Kääb and Treichler2017; Hugonnet and others, Reference Hugonnet2021; Wu and others, Reference Wu2021); however, field observations are important for verifying the glacier-scale mass-balance estimates, given the uncertainties in remote-sensing data (Cogley, Reference Cogley2009; Xu and others, Reference Xu, Li, Li, Wang and Zhou2019). To reconstruct past and project future changes to glacier, one certainly needs a model. A variety of mass-balance modelling approaches are present, each with limitations, but that most require careful calibration and evaluation (Hock, Reference Hock2005; Huss and Hock, Reference Huss and Hock2015; Oerlemans and Reichert, Reference Oerlemans and Reichert2000; Yang and others, Reference Yang, Yao, Guo, Zhu, Li and Kattel2013). Statistical models presume a direct relationship between glacier-wide mass balance and climate variables and can be readily applied to glaciers with long observational records of mass balance, such as Xiao Dongkemadi Glacier (XDG). Only long-term conventional observations (temperature and precipitation) are available at the high altitude meteorological stations across the Tibetan Plateau, the statistical models used to correlate the glaciological mass balance with air temperature and precipitation is effective in reconstructing glacier mass balance in this region (Wang and others, Reference Wang, He, Pu, Jiang and Jing2010, Reference Wang, Yao, Tian and Pu2017; Zhang and others, Reference Zhang, Li, Wang and Wang2014). However, the large financial and logistical requirements of obtaining in situ field observations on the remote glaciers across the central Tibetan Plateau has resulted in continuous mass-balance observations of <20 glaciers to date (Yao and others, Reference Yao2012; WGMS, 2021).
Glacier thickness measurements are essential for estimating glacier volume. The ice thickness can be measured directly at full-depth boreholes to the glacier bed or via geophysical methods, including radar, seismic, geoelectric and electromagnetic methods (Welty and Contributors, Reference Welty2020). Ground-penetrating radar (GPR), which is based on the transmission, reflection and subsequent detection of radio waves, is the most common method (Schroeder and others, Reference Schroeder2020). The Glacier Thickness Database v3 (GlaThiDa v3) indicates that in situ ice thickness observations are available for only ∼3000 glaciers worldwide (GlaThiDa Consortium, 2020) and <40 of these being located across the Tibetan Plateau (Liu and others, Reference Liu, Sun, Shen and Li2003; Cao and others, Reference Cao, Pan, Guan, Wang and Wen2017; Che and others, Reference Che, Wang, Liang and Chen2022). Gurenhekou and Zhadang glaciers are the only glaciers in the central Tibetan Plateau with GPR-measured ice thickness constraints to date (Ma and others, Reference Ma, Tian, Wei and Wei2008; Zhu and others, Reference Zhu, Yao, Yang and Tian2014). Due to the lack of in situ ice thickness measurements, the modelled ice volumes of these glaciers vary widely across this region (Frey and others, Reference Frey2014; Bahr and others, Reference Bahr, Pfeffer and Kaser2015; Farinotti and others, Reference Farinotti2019; Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022).
The Tanggula Mountains in the central Tibetan Plateau form the northern boundary of the South Asian monsoon, as defined by isotopic evidence (Tian and others, Reference Tian, Masson-Delmotte, Stievenard, Yao and Jouzel2001; Yao and others, Reference Yao2013). Systematic field investigations in the Tanggula Mountains began in 1989, with in situ glaciological measurements acquired at regular intervals across XDG (Fujita and Ageta, Reference Fujita and Ageta2000; Pu and others, Reference Pu2008). The long-term mass-balance record and field meteorological observations of XDG make it a key glacier for studying glacier–climate interactions across the central Tibetan Plateau (Fujita and Ageta, Reference Fujita and Ageta2000; Zhang and others, Reference Zhang, He, Ye and Wu2013; Ke and others, Reference Ke, Ding and Song2015; Shi and others, Reference Shi, Duan, Liu, Yang, Zhang and Sun2016; Liang and others, Reference Liang, Cuo and Liu2018). Furthermore, shallow ice core (Li and others, Reference Li2015), albedo (Wu and others, Reference Wu, Wang, Lu, Pu, Guo and Zhang2015; Zhang and others, Reference Zhang, Jiang, Liu, Sun and Wang2018), mercury deposition (Paudyal and others, Reference Paudyal2017), ice flow velocity (Wu and others, Reference Wu, Liu and He2016), trace elements (Li and others, Reference Li, Huang, Li and Zheng2020; Li and others, Reference Li2011), carbonaceous matter (Li and others, Reference Li2023), bacterial diversity (Xie and others, Reference Xie, Wang, Pu and Chen2009) and glacier runoff (Gao and others, Reference Gao, He, Ye and Gao2011) studies have been conducted on XDG. We have previously estimated the geodetic mass balance of the entire Tanggula Mountains (including XDG) between ∼1969 and ∼2015 using remote sensing data (Chen and others, Reference Chen, Wang, Li, Wu, Zhang and Guo2017).
In this study, we sought to improve on the understanding of XDG volume and historic volume change. We first extended the geodetic volume change record to 2020, making use of multispectral stereoscopic satellite data, to evaluate long-term area and volume changes. Second, we conducted in situ measurements of ice thickness using a GPR, enabling us to constrain the glacier’s contemporary and historic volume. Finally, we reconstructed long-term mass-balance changes using a statistical model based on three decades of measurements. The field observations and remote sensing results presented in this study will enhance our understanding of the glacier variations across the central Tibetan Plateau.
2. Study area
XDG (33°10′ N, 95°08′ E) is located within the Dongkemadi Ice Field, Tanggula Mountains (Fig. 1). It is a small valley glacier that separated from Dongkemadi Glacier (DG) between 2000 and 2007 due to continued glacier shrinkage (Shi and others, Reference Shi, Duan, Liu, Yang, Zhang and Sun2016). The glacier area was 1.77 ± 0.08 km2 in 1969, based on a topographic map, with the glacier extending from 5380 m above sea level (a.s.l.) to 5926 m a.s.l (Pu and others, Reference Pu2008). Glaciological mass-balance measurements have been obtained continuously since 1988/89, with the annual mass balance varying from 525 mm water equivalent (w.e.) in 1988/89 to −1066 mm w.e. in 2009/10 (Pu and others, Reference Pu2008; Zhang and others, Reference Zhang, He, Ye and Wu2013). Brief weather observations from XDG during 2008–12 indicated that the annual air temperature and precipitation near the equilibrium line were −8.6°C and 680 mm, respectively (Zhang and others, Reference Zhang, He, Ye and Wu2013).
Location of XDG. (a) Geographical location of Dongkemadi Ice Field in the Tibetan Plateau, with the red star is the location of Anduo Meteorological station; (b) location of XDG in the Dongkemadi Ice field; (c) stake network on XDG in 2015.
3. Data and methods
3.1. Mass balance during 1969–2020
3.1.1. Glaciological method
The glaciological mass balance of XDG has been monitored continuously since 1989 (Pu and others, Reference Pu2008). By 2015, a network of 24 stakes were established to cover the entire glacier surface ranging from 5400 to 5750 m a.s.l. (Fig. 1c). We determined the net mass balance according to the glaciological method, by measuring changes in stake heights and snowpit characteristics in both the ablation and accumulation areas. At each stake, we recorded the height of the stake above the glacier surface, snow density and thickness, presence of superimposed ice layers and the structure of the snowpit profile (Pu and others, Reference Pu2008). The mass balance at each stake is assigned to a corresponding altitude range, and the specific mass balance is depicted as a function of altitude. These measurements are then extrapolated to estimate the mass balance for the entire glacier (Chen and other, Reference Chen, Wang, Li, Wu, Zhang and Guo2017). The mass-balance records from 1989 to 2010 were published by Yao and others (Reference Yao2012), while subsequent observations have extended this record to 2016.
The mass balance derived using the glaciological method involves several sources of uncertainty, including uncertainties in field measurements at specific locations, spatial averaging of these measurements across the entire glacier and changes in glacier area and elevation (Zemp and others, Reference Zemp2013). Thibert and others (Reference Thibert, Blanc, Vincent and Eckert2008) suggested that contributing to mass-balance uncertainty—such as measurement, extrapolation and glacier change uncertainties—can be calculated separately. Following Xu and others (Reference Xu, Li, Li, Wang and Zhou2019), we estimate the uncertainty of the mass balance by 0.21 and 1.53 m w.e. for the glacier-wide and cumulative mass balance over the period from 2010 to 2016.
3.1.2. Geodetic method
The prior study of Chen and others (Reference Chen, Wang, Li, Wu, Zhang and Guo2017) assessed the glacier surface-elevation changes for the entirety of the Tanggula Mountains between ∼1969 and ∼2015. Surface elevation measurements of the glacier surface at different periods were obtained from a topographic map, Shuttle Radar Topography Mission (SRTM) 1 data, and Landsat Thematic Mapper/Operational Land Imager (TM/OLI), Terra Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and ZiYuan-3 (ZY-3) satellite images (Table 1). A Topo Digital Elevation Model (DEM) (generated from topographic map) and ASTER DEM of XDG have previously been generated (Chen and others, Reference Chen, Wang, Li, Wu, Zhang and Guo2017). We included the ZY-3 stereopsis (spatial resolution of 3.5 m) acquired on 9 October 2020 in this study to generate the ZY-3 DEM in 2020.
Details of multi-source remote sensing data
Notes: *Date acquired given in dd/mm/yyyy format, noting that the specific dates for the Topographic map and SRTM were unavailable.
** The ID of topographic map is confidential. ***The spatial resolutions of the orthophotos and stereopsis are 2.1 and 3.5 m, respectively.
The Topo DEM was generated from the 1969 topographic map using the ANUDEM 5.23 software through vectorized contour lines and elevation points. The ASTER and ZY-3 DEMs were generated using the ‘DEM Extraction’ module in ENVI 5.2. The coordinate system was defined using the datum WGS1984 UTM zone 46, and the spatial resolution was resampled to 30 m via bilinear interpolation. The SRTM approximates the surface elevation in 1999 by considering the seasonality correction (Chen and others, Reference Chen, Wang, Li, Wu, Zhang and Guo2017).
The automated, iterative implementation of the co-registration algorithm for different DEMs was carried out following the method described by Nuth and Kääb (Reference Nuth and Kääb2011) and subsequently generalized by Shean and others (Reference Shean2016) using Python and shell scripts. The relative uncertainties of DEMs were evaluated using the non-glaciated terrain, which remained stable throughout the study period. The uncertainty in glacier elevation changes (
${\varepsilon _{{\text{DEM}}}}$) was calculated using Eqn. (1), incorporating the radar wave penetration accuracy (
${\varepsilon _{\text{p}}}$) and normalized median absolute deviation (NMAD, ∆σ) for the glacier-free terrain after DEM co-registration (Chen and others, Reference Chen, Wang, Li, Wu, Zhang and Guo2017):
\begin{equation}{\varepsilon _{DEM}}{\textrm{ = }}\sqrt {{{\left( {\Delta \sigma } \right)}^2}{{\left( {{\varepsilon _p}} \right)}^2}} \end{equation} The relative uncertainty between Topo DEM and ZY-3 DEM can be directly described by NMAD, as both DEMs are derived from optical images. The radar wave penetration accuracy across Dongkemadi Ice Field was determined to be 3.5 m, based on comparisons of the ICESat GLA14 footprints and SRTM elevation data (Chen and others, Reference Chen, Wang, Li, Wu, Zhang and Guo2017). The overall mass-balance uncertainty (
${\varepsilon _{\text{M}}}$) was estimated using Eqn. (2), following the methods outlined by Braun and others (Reference Braun2019):
\begin{equation}{\varepsilon _{\text{M}}} = \sqrt {{{\left( {\frac{{\Delta M}}{{\Delta t}}} \right)}^2}\left( \!\! {{{\left( {\frac{{{\varepsilon _{{\text{DEM}}}}}}{{\Delta h}}} \right)}^2} {+} {{\left( {\frac{{{\varepsilon _{\text{A}}}}}{A}} \right)}^2} {+} {{\left( {\frac{{{\varepsilon _\rho }}}{\rho }} \right)}^2}} \right) {+} {{\left( {\frac{{{\varepsilon _{\text{p}}}}}{{\Delta t}} {\times} \rho } \right)}^2}} \end{equation} where ΔM/Δt represents the mass-balance estimate; 
$\Delta h$, A and ρ (850 kg m−3, Huss, Reference Huss2013) denote the mean elevation difference (MED) of the glacier, the glacier area and the density for volume-to-mass conversion, respectively; 
${\varepsilon _{\text{A}}}$ and 
$\,{\varepsilon _\rho }$ (60 kg m−3, Huss, Reference Huss2013) represent the uncertainties of the glacier area and conversion density, respectively.
3.2. Glacier area
XDG glacier boundaries were manually delineated using a November 1969 topographic map, a Landsat-5 satellite image from 30 August 2000 and a ZY-3 satellite image from 9 October 2020. Aerial photographs from 1969 were used to produce a 1:100 000 topographic map in 1973. The map was scanned, re-projected to the Universal Transverse Mercator coordinate system and referenced to the World Geodetic System 1984 datum. The Landsat-5 and ZY-3 satellite orthophotos have spatial resolutions of 30 and 2.1 m, respectively. The uncertainty in the glacier boundary delineated from the topographic map was set to 5%, and the uncertainties (
${\varepsilon _{\text{A}}}$) in the boundaries that were delineated from the satellite images were estimated by applying a half-pixel buffer to each analyzed image, as follows (Paul and others, Reference Paul2013):
\begin{equation}{\varepsilon _{\text{A}}} = N \times \frac{{{S^2}}}{2}\end{equation} where N is the number of pixels of the glacier boundary and S is the spatial resolution of the image. The uncertainty in the glacier area change (
${\varepsilon _{{\text{AC}}}}$) for each pair of images was calculated based on error propagation (Guo and others, Reference Guo, Liu, Wei and Bao2013), as follows:
\begin{equation}{\varepsilon _{AC}} = \,\sqrt {\varepsilon _{A1}^2 + \varepsilon _{A2}^2} \end{equation} where 
${\varepsilon _{{\text{A}}1}}$ and 
$\,{\varepsilon _{A2}}$ represent the uncertainties of the glacier area of each image.
3.3. GPR measurements and glacier volume
A GPR survey was conducted on XDG in August 2020 using a pulseEKKO PRO system to constrain the ice thickness of the glacier. A total of 122 GPR measurements were obtained along one longitudinal profile and six transverse profiles (Fig. 2a), with measurement spacings of ∼30 and 40 m, respectively. The GPR operated at a frequency of 100 MHz, with the transmitting and receiving antennas separated by 1.5 m at each measurement point. The radargrams clearly showed the glacier bedrock interface for each profile (Fig. 2b), and the ice thickness was calculated based on the two-way travel time of the GPR reflection from glacier bedrock interface, using a radio-wave velocity of 0.169 m ns−1 for glacier ice. The GPR measurement locations were surveyed with a Unistrong MG868 global positioning system (GPS), which provides single point positioning accuracy of 1.2 m. The same GPR parameters were employed at Bayi Glacier, where a 1% uncertainty was reported between the measured ice thickness and ice core length through the glacier (Wang and Pu, Reference Wang and Pu2009). The ice thickness data were interpolated using the Ordinary Kriging method and constrained by glacier outline to create an ice thickness raster. This raster was subsequently used to calculate the mean ice thickness and ice volume of XDG over the five-decade analysis period.
GPR measurements on XDG: (a) distribution of GPR measurement points; (b) GPR profile along transect AA′ in (a).
We created a glacier bed DEM for XDG by subtracting the 2020 GPR derived ice thickness raster from the 2020 ZY-3 DEM. Using this glacier bed DEM, we calculated the elevation differences between the glacier bed DEM and the 1969 topographic map, the 2000 SRTM and the 2015 ASTER DEM to estimate ice thicknesses for 1969, 1999 and 2015, respectively. By combining these ice thickness estimates with the glacier area at different times, we reconstructed the corresponding ice volumes for XDG. The DEM projections were transformed to the WGS1984/UTM zone 46 coordinate system, and their spatial resolutions were resampled to 30 m.
The uncertainties in both the glacier area (
${\varepsilon _{\text{A}}}$) and mean ice thickness (
${\varepsilon _{\overline{H}}}$) are needed to assess the glacier volume uncertainty (
${\varepsilon _{\text{V}}}$) as follows:
\begin{equation}{\varepsilon _{\text{V}}} = \sqrt {{{\left( {\bar H \times {\varepsilon _{\text{A}}}} \right)}^2} + {{\left( {A \times {\varepsilon _{\overline{H}}}} \right)}^2}} \end{equation} The uncertainty in the glacier area is detailed in Section 3.2, A is glacier area. The uncertainty in ice volume change (
${\varepsilon _{{\text{VC}}}}$) for each period was calculated using error propagation methods:
\begin{equation}{\varepsilon _{{\text{VC}}}} = \sqrt {{{\left( {{\text{MED}} \times {\varepsilon _{\text{A}}}} \right)}^2} + {{\left( {{A_{\text{F}}} \times {\varepsilon _{{\text{DEM}}}}} \right)}^2}} \end{equation} where 
${A_{\text{F}}}$ and 
${\varepsilon _{\text{A}}}$ represent the glacier area and the uncertainty of glacier area at former time in each period. The MED and 
${\varepsilon _{{\text{DEM}}}}$ are detailed in Section 3.1. The uncertainties of reconstructed ice thicknesses (
${\varepsilon _{\overline{H}}}$) are the uncertainty in glacier elevation changes (
${\varepsilon _{{\text{DEM}}}}$) between glacier bed DEM and multitemporal DEMs. The uncertainty in the ice thickness (
${\varepsilon _{\overline{H}}}$) derived from GPR measurements is divided into two components: the GPR measurement uncertainty (
${\varepsilon _{{\text{Hdata}}}}$) and the uncertainty associated with the Ordinary Kriging method results (
${\varepsilon _{{\text{HRMSE}}}}$) (Che and others, Reference Che, Wang, Liang and Chen2022):
\begin{equation}{\varepsilon _{\overline{H}}} = \sqrt {\varepsilon _{{\text{Hdata}}}^2 + \varepsilon _{{\text{HRMSE}}}^2} \end{equation} The GPR measurement uncertainty (
${\varepsilon _{{\text{Hdat}}{{\text{a}}_i}}}$) is influenced by the pulsed-radar measurement uncertainty (
${\varepsilon _{{\text{HGP}}{{\text{R}}_i}}}$) and the horizontal positioning uncertainty (
${\varepsilon _{{\text{Hx}}{{\text{y}}_i}}}$). Consequently, the ice thickness uncertainty at each measurement point can be estimated as follows (Lapazaran and others, Reference Lapazaran, Otero, Martin-Espanol and Navarro2016):
\begin{equation}{\varepsilon _{{\text{Hdat}}{{\text{a}}_i}}} = \sqrt {\varepsilon _{{\text{HGP}}{{\text{R}}_i}}^2 + \varepsilon _{{\text{Hx}}{{\text{y}}_i}}^2} \end{equation} We estimate the uncertainty in the mean ice thickness (
${\varepsilon _{Hdata}}$) of entire glacier as the average of the uncertainties from all the GPR measurement points:
\begin{eqnarray}\varepsilon_{Hdata}=\frac{1}{n}{\sqrt{\sum_{i=1}^{n}\varepsilon^{2}_{Hdata_{i}}}} \end{eqnarray} The horizontal positioning uncertainty is determined by the GPS accuracy (1.2 m) because the GPS was not directly linked to the GPR in the field. The pulsed-radar measurement uncertainty (
${\varepsilon _{{\text{HGPR}}}}$) is derived from the uncertainties in the radio wave velocity in ice (
${\varepsilon _{{{\text{H}}_{\text{c}}}}}$) and the observed travel time (
${\varepsilon _{{{\text{H}}_{\text{t}}}}}$), as follows:
\begin{equation}{\varepsilon _{{\text{HGPR}}}} = \sqrt {\varepsilon _{{{\text{H}}_{\text{c}}}}^2 + \varepsilon _{{{\text{H}}_{\text{t}}}}^2} \end{equation} where 
${\varepsilon _{{{\text{H}}_{\text{c}}}}} = ({\varepsilon _{\text{c}}} \times t)/2$ and 
${\varepsilon _{{{\text{H}}_{\text{t}}}}} = ({\varepsilon _{\text{t}}} \times c)/2$ are the components of the thickness uncertainty due to uncertainties in the radio wave velocity in ice and the observed travel time 
$({\varepsilon _{\text{t}}})$ respectively. Here, t and c represent the time and velocity of the radar signal in ice, respectively. We assume 
$c = 0.169\,m\cdot{\text{n}}{s^{ - 1}}$, 
${\varepsilon _{\text{c}}} = 0.02c$ and 
${\varepsilon _{\text{t}}} = 1/{f_{\text{r}}}$, where 
${f_{\text{r}}}$ is the frequency of the GPR (100 MHz in this study).
We used 85% of the measured points for interpolation and reserved the remaining 15% for cross validation. These cross-validation points were randomly distributed across the entire glacier to quantitatively assess uncertainties in the interpolation method (Fig. 3a). The uncertainty in the Ordinary Kriging method results (
${\varepsilon _{{\text{HRMSE}}}}$) is represented by the root-mean-square error (RMSE), calculated from the interpolated (
$\widehat {{y_i}}$) and measured (
${y_i}$) ice thickness values:
\begin{equation}{\varepsilon _{{\text{HRMSE}}}} = \sqrt {\frac{{{{\left( {{y_1} - \widehat {{y_1}}} \right)}^2} + {{\left( {{y_2} - \widehat {{y_2}}} \right)}^2} + \ldots + {{\left( {{y_i} - \widehat {{y_i}}} \right)}^2}}}{i}} \end{equation}Cross-validation of the ordinary kriging method: (a) distribution of the verification points; (b) comparison between interpolated and measured ice thickness values.
The cross-validation test revealed an RMSE of 3.31 m between interpolated and measured ice thicknesses (Fig. 3b). The high consistency with an R 2 value of 0.99 indicates that the ice thickness estimates for XDG using the Ordinary Kriging method are reliable.
3.4. Meteorological data
Anduo meteorological station (32°21′ N, 91°06′ E, 4800 m) is located on the southern slope of the Tanggula Mountains (Fig. 1a) and is the closest station (121 km) to XDG. We selected meteorological data from Anduo station to investigate the potential drivers behind the observed variations in XDG’s mass balance. Temperature and precipitation data were provided by the Chinese National Meteorological Center (http://data.cma.cn/site/index.html). The average annual temperature and precipitation over the long-term period from 1967 to 2020 are −2.4°C and 447.6 mm, respectively.
To better understand the meteorological conditions at XDG, we assessed the climatic elements at the median elevation of 5600 m (Fig. 4). This analysis is based on observational data from Anduo station and the altitudinal gradients for the warm season: −0.63°C/100 m for temperature (Shi and others, Reference Shi, Duan, Nicholson, Han, Klaus and Yang2020) and 23 mm/100 m for precipitation (Zhang and others, Reference Zhang, He, Ye and Wu2013).
Variations in warm season temperature (red) and precipitation (blue) at the median height (5600 m a.s.l.) of XDG from 1967 to 2020.
Considering the differences in the extent to which glaciological mass balance is influenced by temperature and precipitation in different periods of the year, it is necessary to select a reasonable period of the year when there is a significant correlation between the temperature and/or precipitation and glacier mass balance. The climate data were therefore divided into several seasonal periods: cold season (October–April), warm season (May–September), winter (December–February) and summer (June–August) (Pu and others, Reference Pu2008). We will identify the climatic factors most strongly correlated with glacier mass balance to serve as the key drivers in constructing the statistical model.
3.5. Statistical model of glacier mass balance
Given that the statistical model is effective for reconstructing mass balance using only conventional meteorological data (Wang and others, Reference Wang, He, Pu, Jiang and Jing2010), we select this approach to reconstruct the mass-balance history of XDG over the past 50 years. A correlation model between glacier mass balance and climatic factors for XDG is established by using the 28 year time series of XDG mass-balance observations and the long-term climate data at the median elevation of XDG. This method enabled us to reconstruct past mass-balance variations and assess the sensitivity of the glacier mass balance to climate change. We constructed a multivariate regression based on the correlation analysis results in Section 4.4, with the temperature and precipitation as independent variables and glacier mass balance as the dependent variable. We chose to assess the uncertainty of reconstructed mass balances (
${\varepsilon _{{\text{RM}}}}$) by constructing confidence intervals of them (He and Liu, Reference He and Liu2015), as follows:
\begin{equation}{\varepsilon _{{\text{RM}}}} = {z_{\alpha /2}} \times {\text{SE}}\end{equation} where 
$\alpha $ is the significance level, and the value of the statistical model is used in this paper, 
${z_{\alpha /2}}$ is the two-sided quantile of the standard normal distribution, which can be obtained by looking up the standard normal distribution table. 
${\text{SE}}$ is the standard error of the reconstructed mass balance.
4. Results
4.1. Geodetic and glaciological mass balances of XDG
Figure 5 shows the annual (blue) and cumulative (red) mass balance for 1989–2016, as calculated using the glaciological method. The cumulative glaciological mass balance of XDG was −8.47 ± 1.53 m w.e., and the annual mass balance varied between +0.52 ± 0.21 and −1.07 ± 0.21 m w.e. a−1, with a mean of −0.30 ± 0.21 m w.e. a−1 over the 28 year period. We extended the XDG mass-balance time series to 2020 by combining the results of this study with those from our previous research (Chen and others, Reference Chen, Wang, Li, Wu, Zhang and Guo2017, Fig. 5, Table 2). XDG experienced continued mass loss from 1969 to 2020, with a cumulative geodetic mass balance of −9.69 ± 3.06 m w.e., equivalent to an annual rate of −0.19 ± 0.02 m w.e. a−1. The mass loss varied over different periods: limited loss of −0.06 ± 0.10 m w.e. a−1 from 1969 to 1999, an increased rate of −0.41 ± 0.14 m w.e. a−1 from 1999 to 2015 (which is close to the glaciological value of −0.42 ± 0.21 m w.e. a−1) and −0.32 ± 0.14 m w.e. a−1 from 2015 to 2020.
Annual glaciological mass balance (blue), cumulative mass balance (red), mean glaciological value for 1999–2015 (yellow line) and annual rate of geodetic mass balances (black lines, with the orange regions indicating the uncertainty of geodetic results) of XDG over the investigated periods, 5 yearly geodetic mass balances (purple lines, with the light blue regions indicating the uncertainty of geodetic results) provided by Hugonnet and others (Reference Hugonnet2021) for DG.
Changes in area, surface elevation and geodetic mass balance of XDG during the 1969–2020
Note: The geodetic results for the period before 2015 are sourced from Chen and others (Reference Chen, Wang, Li, Wu, Zhang and Guo2017), but the mass-balance uncertainty is estimated using a new equation presented in this study; MED denotes the mean elevation difference and SMB stands for specific mass balance.
4.2. Glacier area changes during 1969–2020
The changes in the ice-covered area of XDG from 1969 to 2020, based on comparisons between the 1969 topographic map and the 1999, 2015 and 2020 satellite images (Fig. 6, Table 2), show a decrease from 1.77 ± 0.08 to 1.48 ± 0.02 km2, amounting to a reduction of 0.29 ± 0.08 km2 (16.38 ± 4.66%). XDG has experienced continuous shrinkage throughout the different investigated periods, with the highest annual shrinkage rate of −0.64 ± 0.52% · a−1 occurring between 1999 and 2015. The loss is primarily due to the retreat of the glacier terminus and surface lowering at the western headwall (indicated by the orange circle in Fig. 6).
XDG boundaries in 1969 (red, derived from a topographic map); 1999 (yellow, derived from Landsat 5 satellite imagery); 2015 (purple, derived from Landsat 8 satellite imagery); and 2020 (blue, derived from ZY-3 satellite imagery). The background image is a ZY-3 satellite imagery acquired on 9 October, 2020. The dashed circle highlights an area of extensive surface lowering along the western headwall of XDG.
4.3. Glacier thickness in 2020 and ice volume changes during 1969–2020
We interpolated measurements from 104 GPR points to calculate the thickness of XDG. In 2020, the glacier had a mean ice thickness of 54.78 ± 3.69 m, with the thickest ice (146.88 ± 1.69 m) found in the upper firn basin. The thickness distribution indicates a gradual increase from the glacier edges toward the central part (Fig. 7a). Along this central part, the ice thickness exceeds 120 m for approximately half of the glacier’s length, with the maximum thickness occurring in the upper firn basin, where it surpasses 135 m (Fig. 6a). The glacier bed’s topography resembles that of a typical valley (Fig. 7b), showing an increasing breadth–depth ratio with higher elevations.
Glacier thickness in 2020 (a) and morphology of ice bed (b) for XDG.
In 2020, XDG had an ice volume of 0.0811 ± 0.0056 km3, calculated from a mean ice thickness of 54.78 ± 3.69 m and a glacier area of 1.48 ± 0.02 km2. The reconstructed ice volumes for XDG were 0.1175 ± 0.0066 in 1969, 0.1073 ± 0.0104 in 1999 and 0.0848 ± 0.0047 km3 in 2015, as determined from topographic map, SRTM, and ASTER DEM, respectively. These volumes corresponded to mean ice thicknesses of 66.40 ± 3.44, 64.28 ± 3.73 and 56.56 ± 2.13 m, respectively (Table 3). Surface elevation change maps (Fig. 8) indicate that XDG has experienced ice loss throughout all investigated periods, with a total loss of 0.0364 ± 0.0051 km3 (31.01 ± 4.59%) over from 1969 to 2020. The ice loss rate was 0.29 ± 0.18% · a−1 from 1969 to 1999, which is much lower than the rates of 1.00 ± 0.37 and 0.89 ± 0.75% · a−1 during 1999–2015 and 2015–20, respectively.
XDG surface-elevation differences for the (a) 1969–99, (b) 1999–2015 and (c) 2015–20 periods.
Changes in mean ice thickness and ice volume of XDG during the 1969–2020
Note: The uncertainty of the mean ice thickness in 2020 was estimated from GPR measurements, and the other uncertainties were based on via geodetic method.
4.4. Statistical mass-balance reconstruction from 1967 to 2020
To identify the optimal climate factors for driving the statistical model of glacier mass balance, we analyze the correlation coefficients between the XDG mass balance and temperature and precipitation data from Anduo station from the period 1989–2016 (Table 4). The correlation coefficients between mass balance and temperature were approximately −0.79 (p < 0.01) for both the warm season and summer. Given that glacier melting at XDG typically begins in May, we selected the temperature during the warm season as a key climatic factor influencing the glacier’s mass balance. Significant correlations were found between the mass balance and precipitation during both the annual and warm seasons (Table 4). Over 90% of the precipitation in the Dongkemadi basin falls during the warm season (He and others, Reference He, Ye and Ding2009), and recent field observations show minimal mass change during the cold season (October–April). Consequently, we consider warm season precipitation as an additional key climatic factor affecting glacier mass balance.
To reconstruct and predict the mass balance of XDG, we used multivariate regression to determine the relationship between the mass balance (M b, mm w.e.), warm season temperature (T w, ℃) and precipitation (P w, mm) at the median glacier height as follows:
\begin{equation}{{\text{M}}_{\text{b}}}{\text{ = }} - {\text{498}}{\text{.1 \times }}{{\text{T}}_{\text{w}}}{\text{ + 1}}{\text{.1 \times }}{{\text{P}}_{\text{w}}} - {\text{399}}{\text{.6}}\;\left( {{R^2} = 0.65,{\text{ }}n = 28} \right)\end{equation}Correlation coefficients between the XDG mass balance, and temperature and precipitation at Anduo meteorological station for different seasonal periods
Note: *significant correlation with less than 0.05 probability of random chance for the relationship,
** Significant correlation with less than 0.01 probability of random chance for the relationship.
The variance test of the regression formula showed that the significance level was within 0.01, indicating a strong relationship between the parameters. The uncertainty of the reconstructed mass balances, calculated using Eqn. (12), was 0.12 m w.e. The reconstructed mass balance of XDG ranged from −0.82 to +0.44 m w.e., with a long‐term average of −0.20 ± 0.12 m w.e., over the period from 1967 to 2020. The cumulative mass balance during this time was −10.60 m w.e., indicating mass loss during this period (Fig. 9a). The mass-balance reconstruction shows that XDG experienced two distinct patterns of change during 1967–2020, the mass balance exhibited weak fluctuations from 1967 to 1993 (with the cumulative mass balance of −0.15 m w.e.) and overall mass losses have occurred since 1994 (with the cumulative mass balance of −10.45 m w.e.). This pattern aligns with similar observations reported by Shangguan and others (Reference Shangguan, Liu, Ding, Zhang, Du and Wu2008), and Zemp and others (Reference Zemp2020) documented a comparable trend of negative glacier mass balance on a global scale.
Variations in the mass balance of XDG from 1967 to 2020 are illustrated as follows: (a) triangles denote observed values (black), closed circles represent reconstructed values (blue) and the reconstructed cumulative mass balance is shown in red; (b) comparison between observed and reconstructed values.
5. Discussion
5.1. Validation of the mass-balance estimates from different methods
Mass balance is a basic parameter for monitoring glacier changes. We estimated the mass balance of XDG by using glaciological methods, geodetic methods and statistical models. The glaciological method is a widely accepted approach for measuring mass balance, and the glaciological mass balance is used to verify the mass balance acquired from other methods. We had verified the good agreement between glaciological and geodetic mass balances of XDG, and the absolute differences in mass balances between the two methods were <5% for the different periods between 1999 and 2015 (Chen and others, Reference Chen, Wang, Li, Wu, Zhang and Guo2017). For the period of 1989–2016, the reconstructed cumulative mass balance of −8.17 m w.e. closely matches the glaciological value of −8.47 ± 1.53 m w.e., and the correlation coefficient between reconstructed and glaciological values was 0.80 (RMSE = 0.25, p < 0.01, Fig. 9b). The reconstructed mean mass balances were −0.21 ± 0.12, −0.08 ± 0.12, −0.34 ± 0.12 and −0.58 ± 0.12 m w.e. a−1 for the periods of 1969–2020, 1969–99, 1999–2015 and 2016–20, respectively. These values are approximate to the geodetic mass balances of −0.19 ± 0.02, −0.06 ± 0.10, −0.41 ± 0.14 and −0.32 ± 0.14 m w.e. a−1 for the corresponding periods. The verification results show that the mass balances of XDG estimated by different methods have good agreement, which gives us the confidence to compare with others.
5.2. Comparison of glacier variations with other studies
Remote sensing data have documented significant glacier shrinkage across the Tibetan Plateau, particularly in its central part. Previous studies have focused on regional-scale glacier changes, noting a 28.94% reduction in glaciers across the Tanggula Mountains between 1980 and 2020 (Zhang, Reference Zhang2023) and a surface elevation change of −0.36 ± 0.06 m a−1 between ∼1969 and 2015 (Chen and others, Reference Chen, Wang, Li, Wu, Zhang and Guo2017). A GPS survey conducted on XDG in October 2007 indicated a surface lowering rate of 0.21 ± 0.14 m a−1 between 1969 and 2007 (Shangguan and others, Reference Shangguan, Liu, Ding, Zhang, Du and Wu2008). This rate falls between the surface-lowering rates observed for the 1969–99 and 1999–2015 periods in this study. Hugonnet and others (Reference Hugonnet2021) provided mass-balance time series of DG, which includes Dadongkemadi Glacier and XDG (Fig. 1 and Table 5). Their study calculated a mass balance of −0.87 ± 0.46 m w.e. a−1 for 2015–20, which is higher than the geodetic result of −0.32 ± 0.14 m w.e. a−1 for the same period in this study. There were three possible reasons for the discrepancy between these two studies. First, the calculation methods differ. We computed glacier elevation changes by differencing DEMs from two distinct epochs, while Hugonnet and others (Reference Hugonnet2021) derived elevation changes from time series data and interpolated elevation observations into continuous time series using Gaussian process regression. Secondly, the specific time period of the geodetic mass balance from Hugonnet and others (Reference Hugonnet2021) is from 1 January 2015 to 1 January 2020, corresponding mass-balance year is from 2014/15 to 2018/19 after considering the seasonality correction. For our geodetic mass balance, the corresponding mass-balance year is from 2015/16 to 2019/20. The reconstructed records indicated that the mass loss in 2014/15 (−0.80 m w.e. a−1) was approximately twice that in 2019/20 (−0.42 m w.e. a−1). The third possible reason is different pixel resolutions between Hugonnet and others (Reference Hugonnet2021, 100 m) and our study (30 m). The comparison of annual mass balance revealed an RMSE of 0.30 m w.e. between the glaciological method of XDG and that of Hugonnet and others (Reference Hugonnet2021) for the period of 2000–16 (Fig. 10). The weak consistency (R 2 = 0.34) of two datasets between the datasets suggests that the regional ASTER data have limitations for estimating annual mass balance, but better consistency observed over 5 year and decade scales (Table 5).
Comparison of annual mass balance between glaciological method for XDG and those provided by Hugonnet and others (Reference Hugonnet2021) for DG.
Mass-balance values for XDG and DG during different periods, as obtained from the glaciological method and Hugonnet and others (Reference Hugonnet2021)
5.3. Comparison between gridded and observed ice thickness
Accurate assessments of a glacier’s ice volume and thickness distribution are critical for evaluating its water resource potential and predicting glacier variations. Based on multi-source DEMs, Farinotti and others (Reference Farinotti2019) and Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) estimated the global ice thickness using a weighted ensemble of five models and shallow-ice approximation, respectively. They shared the gridded products, which significantly enhanced the understanding of the distribution of global ice thickness. However, the gridded products of central Asia released by Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) were most likely calculated from the GDEM, but the acquisition time of GDEM is unclear. Farinotti and others (Reference Farinotti2019) estimate the ice volume of XDG as 0.1061 km3, based on SRTM data for the Tibetan Plateau. This value is close to the reconstructed ice volume of 0.1073 ± 0.0104 km3 for 1999 in this study. However, the glacier thickness distribution released by Farinotti and others (Reference Farinotti2019) shows thinner ice in the middle and thicker ice near the edges compared to the reconstructed and observed ice thicknesses in 1999 (Fig. 11a) and 2020 (Fig. 7a), respectively. Additionally, the thickest ice of XDG in the gridded products of Farinotti and others (Reference Farinotti2019) was recorded as 104.37 m in 1999, which is ∼42 m less than the observed maximum thickness in 2020. These discrepancies in ice thickness distribution and the maximum ice thickness suggest that additional glacier measurements in the Tibetan Plateau are necessary to refine and improve the accuracy of existing glacier thickness products.
The distribution of glacier thickness for XDG in 1999: (a) reconstructed values of this study and (b) estimated values released by Farinotti and others (Reference Farinotti2019).
Ice thickness measurements for glaciers across the Tibetan Plateau are limited due to their remote locations, and the available measurements carry a high degree of uncertainty. The uncertainty in ice thickness is divided into two components: GPR measurement uncertainty and interpolation method uncertainty. GPR measurement uncertainty is a systematic error and can be quantitatively assessed using specific equations (Lapazaran and others, Reference Lapazaran, Otero, Martin-Espanol and Navarro2016). In contrast, uncertainties in estimating the thickness of the entire glacier are primarily influenced by the interpolation methods used, which depend on the number and distribution of measurement points. Langhammer and others (Reference Langhammer, Grab, Bauder and Maurer2019a) employed elements of sequential optimized experimental design for determining cost-optimized GPR survey layouts, but their results were difficult to offer general recommendations. When high-precision ice thickness maps are required, it is therefore advisable to acquire as much data as can be afforded (Langhammer and others, Reference Langhammer, Grab, Bauder and Maurer2019a). However, ground-based GPR measurements on glacier often struggle to capture ice thicknesses at the upper regions and edges due to accessibility issues, yet these areas are crucial for estimating the overall glacier thickness. The advancement of helicopter-borne GPR presents a significant opportunity for more comprehensive ice thickness surveys, greatly enhancing the quantity and coverage of measurements. Since the first study reported by Kennett and others (Reference Kennett, Laumann and Lund1993), helicopter-borne ice thickness measurements of several alpine glaciers have been completed (Blindow and others, Reference Blindow, Salat and Casassa2012; Langhammer and others, Reference Langhammer2019b). Pritchard and others (Reference Pritchard, King, Goodger, McCarthy, Mayer and Kayastha2020) developed the Bedmap Himalayas airborne radar system mainly based on an airborne ice-sounding radar, with which successful experiments were conducted on several glaciers in Nepal. It is believed that the airborne ice-sounding radar will be applied to the ice thickness measurement in the central Tibet Plateau in the near future, which will provide an effective means for obtaining accurate ice thickness information.
5.4. Mass-balance models of XDG
XDG has the longest record in the central Tibetan Plateau, making it an excellent test location for various glacier mass-balance models. A distributed temperature-index mass-balance model coupled with volume-area scaling method has previously been applied for predicting variations of XDG and Dongkemadi Ice Field in the future (Shi and others, Reference Shi, Duan, Liu, Yang, Zhang and Sun2016, Reference Shi, Duan, Nicholson, Han, Klaus and Yang2020). Liang and others (Reference Liang, Cuo and Liu2018) developed a modified distributed surface energy-balance model (DSEBM) with high spatial-temporal resolution, which has been employed to reconstruct the mass balance of the DG from 1960 to 2009 and to analyze its influencing factors (Liang and others, Reference Liang, Cuo and Liu2019, Reference Liang, Cuo and Liu2020, Reference Liang, Cuo and Liu2022). Recently, a hybrid modelling approach has been introduced, coupling an empirical glacier evolution model with a temperature-index model to project response of DG to climate change through 2100 (Han and others, Reference Han, Long, Zhao and Slater2023). Compared with the sophisticated models mentioned above, the statistical mass-balance model offers advantages such as lower input requirements (limited to temperature and precipitation) and ease of calculation, making it an effective tool for reconstructing glacier mass balance on the Tibetan Plateau (Wang and others, Reference Wang, He, Pu, Jiang and Jing2010, Reference Wang, Yao, Tian and Pu2017; Zhang and others, Reference Zhang, Li, Wang and Wang2014). The reconstructed cumulative mass balance of XDG, calculated using the statistical mass-balance model, was −4.65 m w.e. for the period of 1967–2009, which was similar to the value for DG (−5.45 m w.e., calculated using the DSEBM) over the same period (Liang and others, Reference Liang, Cuo and Liu2019). However, the model does not account for the dynamic evolution of glacier morphology, limiting its ability to project future glacier variations.
The statistical model indicated a mass-balance sensitivity of XDG to warm season temperature and precipitation of −498 mm w.e./°C and +1.1 mm w.e./1 mm, respectively. A similar analysis result for Urumqi Glacier No. 1 suggested that a 1°C increase in summer temperature led to a 441 mm w.e. decrease in glacier mass balance, while a 1 mm increase in annual precipitation led to a 1.2 mm w.e. increase in glacier mass balance (Zhang and others, Reference Zhang, Li, Wang and Wang2014). Wang and others (Reference Wang, Yao, Tian and Pu2017) reported that the mass-balance sensitivity to air temperature for Qiyi Glacier in the Qilian Mountains was −239 mm w.e./°C, while the sensitivity to the precipitation was +1.1 mm w.e. per 1 mm increase. These studies reveal that while mass-balance sensitivity to air temperature varies across regions on the Tibetan Plateau, the sensitivity to precipitation remains consistent. Previous studies have reported that the effects of a 1°C warming on the mass balance were equivalent to a >25% increase in precipitation (Oerlemans, Reference Oerlemans2005; Zhang and others, Reference Zhang, Li, Wang and Wang2014; Wang and others, Reference Wang, Yao, Tian and Pu2017). Linear-fitting results reveal a significant increase in the warm season temperature (0.27°C 10 a−1, p < 0.01) and a slight increase in the warm season precipitation (11.15 mm 10 a−1; equivalent to 1.87% 10 a−1) at the median height of XDG during 1967–2020 (Fig. 4), which explained the ongoing shrinkage and ice loss of the glacier since the 1960s.
Han and others (Reference Han, Long, Zhao and Slater2023) predicted that by 2100, the ice volume of DG will be reduced to just 3% of its 2018 level under the SSP585 scenario, which aligns with the projected ∼95% ice loss across the entire Tanggula Mountains under the RCP8.5 scenario using the Python Glacier Evolution Model (Rounce and others, Reference Rounce, Hock and Shean2020). Both studies utilized ice thickness data from Farinotti and others (Reference Farinotti2019). However, the overestimated ice thickness near the edges of XDG in the gridded products of Farinotti and others (Reference Farinotti2019) may lead to a slower rate of recent glacier shrinkage, while the underestimated thickness in the middle part of glacier suggests that it will retain less ice by 2100. Nevertheless, with anthropogenic warming already at 1.25°C above the 1850–1900 baseline and expected to exceed 1.5°C within the next decade (Matthews and Wynes, Reference Matthews and Wynes2022), the Tibetan Plateau climate is likely to follow the SSP585 scenario for the rest of this century. This indicates that XDG is at risk of disappearing before 2100. Future work needs to undertake a more sophisticated glacier evolution projection using the refined understanding of XDG ice thickness and mass balance.
6. Conclusions
As one of the alpine glaciers with the longest continuous mass-balance observational record in the central Tibetan Plateau, XDG serves as a key indicator for understanding how glaciers in this region may evolve in a warming global climate. This study documented the overall shrinkage and ice loss of XDG from 1969 to 2020 using a combination of field observations and remote sensing techniques. During this period, the ice-covered area of XDG decreased from 1.77 ± 0.08 km2 in 1969 to 1.48 ± 0.02 km2 in 2020, with a corresponding ice volume loss of 0.0364 ± 0.0051 km3, equivalent to 31.01 ± 4.59% of the 1969 ice volume. By 2020, the mean ice thickness of XDG was measured at 54.78 ± 3.69 m, translating to an ice volume of 0.0811 ± 0.0056 km3. Although the ice volume of XDG in GlaThiDa v3 closely matched the reconstructed value, there were obvious discrepancies in the distribution of ice thickness and location of thickest ice. Enhancing the accuracy of the glacier thickness estimates requires increasing the number and coverage of measurements, which underscores the need for broader adoption of airborne ice-sounding radar technology.
Continuous field observations at XDG indicated that the mean annual mass balance was −0.30 ± 0.21 m w.e. for the period from 1989 to 2016, with records extended to the 2015–20 period through geodetic methods, showing a value of −0.32 ± 0.14 m w.e. Long-term mass-balance reconstructions using statistical models indicate a cumulative mass balance of −10.60 m w.e., averaging −0.20 ± 0.12 m w.e. from 1967 to 2020. The year 1993 marked a turning point in mass-balance trends, with relatively stable conditions up to that point, followed by significant ice loss thereafter. Warming temperatures have been the dominant driver of shrinkage and ice loss of XDG since the 1960s. XDG faces the risk of disappearing before 2100 if the climate scenario of Tibetan Plateau toward to SSP585 scenario.
Data availability statement
The mass-balance observations and ice thickness measurements were provided by Qiyi Glacier Station and access can be granted upon request by contacting the corresponding author.
Acknowledgements
We thank the Scientific Editor Evan Miles, the Associate Chief Editor Nicolas Eckert and the anonymous reviewers for their constructive comments. This work is supported by the Second Tibetan Plateau Scientific Expedition and Research Program (No. 2019QZKK020102), the National Natural Science Foundation of China (Nos. 42130516, 42171139, 41871053), Natural Science Basic Research Program of Shaanxi Province (No. 2023-JC-QN-0300) and the Qin Chuangyuan Innovation and Entrepreneurship Talent Project (No. QCYRCXM-2023-096).
Author contributions
Ninglian Wang and An’an Chen designed the study, An’an Chen led data analysis and drafted the manuscript. Zhen Li, Yuwei Wu and An’an Chen carried out field investigations and data collection. Xi Jiang, Zhongming Guo, Xuejiao Wu and Xuewen Yang were involved in the discussions. All authors discussed results, implications, provided feedback and approved of the final manuscript.
Open access
