Changes in the Surface Elevation of the Laohugou Glacier No. 12 in Western Qilian Mountains

As the largest valley glacier in the Qilian Mountains, the Laohugou glacier No. 12 (LHG12) has shrunk significantly since 1957. In this study, two topographic maps and a WorldView-2 satellite stereopair image data were used to assess the volume and cumulative mass balance of LHG12 located at the western Qilian Mountains during 1957–2015. During the study period, the LHG12 exhibited changes in two processes: slightly ablation and stability in a brief period during 1957–1989 and strong melting and accelerated ablation during 1989–2015. During 1957–2015, the volume of LHG12 decreased by 0.38 km3, the average thickness decreased by 17.23 m, the cumulative mass balance (MB) was −14.69 ± 3.00 m w. e., and ablation was found glacier-wide. By comparing the previous MB simulation and digital elevation model (DEM) differencing results, it was found that the MB simulation results underestimated the strong melting trend of LHG12 since the 1990s. Temperature rose, especially in autumn and winter, and could cause the ice temperature of LHG12 to increase, and LHG12 may become more sensitive to climate change.


INTRODUCTION
Due to global warming, glaciers have shrunk significantly worldwide, and much of this loss was not reversed (Moon, 2017;Liu et al., 2020;Shean et al., 2020). Glacier mass balance (MB), which is a key component in glaciology, is an important factor for studying the changes in the climate, water resource, and sea level (Zemp et al., 2015;Sold et al., 2016). In recent centuries, glacier mass balances have been considered to be sensitive indicators of climate change (Oerlemans and Fortuin, 1992). Combining traditional observations with satellite altimetry and gravimetry, glacier mass budgets were reconciled in order to obtain an estimate of the glacier contribution to sea level change in the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, but the analysis was only possible over a short time period (Gardner, 2013). Long series glacier mass balance data are important and useful for investigating climate trends and for numerical simulations of glacier dynamics, but only 33 glaciers worldwide have an annual mass-balance series longer than 40 years (Dyurgerov and Meier, 1999;Vincent, 2002;Le Meur and Vincent, 2003). Glacier mass balance can be monitored using traditional glaciological or geodetic methods. Traditional glaciological methods provide in situ observations of the annual and sometimes seasonal mass balance (Zemp et al., 2013), stake measurements and snow pits provide ablation and accumulation data, and point observations can be extrapolated to a glacier-wide mass balance using the contour-line or profile methods (Østrem and Brugman, 1991). Geodetic methods use multi-temporal digital elevation models (DEMs) generated by repeated mapping or stereo image pairs to calculate a glacier-wide mass balance Zemp et al., 2013). The geodetic method has been widely used in glacier mass balance research (Ruiz et al., 2017;Xu et al., 2017a;He et al., 2020). This technique can be used to calibrate long-term glaciological mass balance series (Berthier et al., 2007;Huss et al., 2009).
The Qilian Mountains, located on the northeastern part of the Qinghai-Tibetan Plateau, northwestern China, are an important freshwater resource in the Hexi-Corridor. According to the Second Chinese Glacier Inventory, the Qilian Mountains contains 2684 glaciers, covering an area of 1597.81 ± 70.30 km 2 with an ice volume of 84.48 km 3 . Over the past half-century, the area and volume of the glaciers have decreased by 420.81 km 2 (−20.88%) and 21.63 km 3 (−20.26%), respectively (Sun et al., 2018). As the largest valley glacier in the Qilian Mountains, the Laohugou glacier No. 12 (LHG12) has shrunk significantly, with reduction in the terminus, area, and volume of 402.96 m (3.99%), 1.54 km 2 (7.03%), and 0.1816 km 3 , respectively. Between 1957 and 2015, the reduction rate accelerated . Observations of LHG12 began in 1958. However, the data series is not continuous, as the mass balance observations were interrupted by a summer flood in 1962 and were restarted in 2005. Because of the lack of mass balance measurements, two mass balance simulation data series were obtained using the degree-day factor method (DDF) Chen et al., 2020), but some differences were found between the two sets of simulation results. Thus, the aims of this study are (1) to generate authentic data for mass balance simulations; (2) to generate the glacier surface elevation changes using two topographic maps and a WorldView2 image and to calculate the glacier volume change; and (3) to obtain the glacier wide net mass loss based on the mass balance conversion using the DEM differencing algorithm.

STUDY AREA
LHG12 (glacier number: 5Y448D0012) is located in the Daxueshan region, western Qilian Mountains, northern Tibetan Plateau (Figure 1). It is the largest valley glacier in the Qilian Mountains. LHG12 is 9.7 km long and covers an area of 20.37 km 2 . During 1960-2015, the glacier terminus retreated by about 400 m, and the glacier area decreased by 1.54 km 2 (7.03% in total) . The average equilibrium line altitude (ELA) was 4830 m a.s.l. during 1958-1977(Kang and Ding, 1981, and the ELA0 was 5015 m a.s.l. during 2010-2012 (Chen et al., 2017). The main observations were conducted in the east branch, the confluence, and the terminus regions. LHG12 is a typical continental valley glacier. The annual mean air temperature of LGH12 recorded by an automatic weather station (AWS) at 5040 m a.s.l. was −11.8°C. The annual precipitation was about 443 mm water equivalent (w.e.) and was highly concentrated from May to September (85%) at 4990 m a.s.l. (Sun et al., 2012;Du et al., 2016).

DATA AND METHODS
In this study, two topographical maps derived from aerial photographs acquired in 1957 and 1989 and a WorldView-2 stereo image acquired in 2015 were used to calculate the changes in the glacier's surface elevation. The coordinate system of all images and vector layers processed using QGIS (https://www. qgis.org/en/site/) are WGS84/UTM47N. The climate background was analyzed using the meteorological data recorded at Tuole station (Tuole AWS) from 1957 to 2015 and LHG12 automatic weather station (LHG12 AWS).

Topographic Map
The 20 m interval contours of two geographical maps with scales of 1:50,000 were digitized in order to generate DEMs based on the Beijing54 coordinate. The seven-parameter datum transformation model was used to re-project the DEMs onto the WSG84 coordinate. The transformation error was less than 0.002 m (Wang et al., 2003).

Remote Sensing Image
A Landsat5 TM image acquired on January 5, 1989, was used to measure the boundary and glacier area, and the path and row of the image was p136r33. The glacier boundary was extracted using artificial vectorization, which afforded a satisfaction precision. The error of the artificial vectorization was less than 2% (Liu et al., 2013).

Meteorological Data
The climate background was analyzed using the temperature and precipitation recorded at the Tuole meteorological station during 1957-2015. The Tuole meteorological station, which is a national weather station, is located about 200 km from LHG12. The LHG12 AWS located at confluence region (4550 m a. s. l) was set in 2008, and temperature data were obtained during 2010-2015 ( Figure 1).

Mass Balance Calculations
Glacier DEM differencing generated by repeated mapping can be used to determine the volume and surface changes. The results of the glacier volume change (ΔV) can be converted into a specific mass balance over a period of record (PoR) in units of metre water equivalent (m w.e.) (Zemp et al., 2013): where ρ is the average density of ΔV, assuming no change in the bulk glacier density over the balance period.
S is the average glacier area during the two surveys at times t0 and t1, assuming a linear change with time, i.e., S S t0 + S t1. 2 (2) ρ is a key parameter in mass convention. In this study, 850 ± 60kg · m −3 was used, which is recommended for periods of longer than 5 years, with stable mass balance gradients, the presence of a firn area, and volume changes that significantly differ from zero (Huss, 2013).

Uncertainty Analysis
Due to technical limitations, the mapping era, and other factors, the error of the topographic map is unknown. In order to evaluate the error of the glacier volume change, it was necessary to assume that the non-glacier region was stable terrain. However, in the non-glacier region, the topography may be modified by freeze-thaw action, runoff erosion, fluvial-glacial erosion, and so on. Therefore, the spatial distribution of the check points should be on different slopes and rivers should be avoided.
In this study, the uncertainty was calculated using 39 points in the non-glacial region. The uncertainty of the check point at different times was described by the Root Mean Square Error (RMSE): where H A and H B are the elevation check points on the stable terrain on two topographic maps; n is the number of check points. The DEMs produced from the topographic map and WorldView-2 image were re-sampled to a 30 m ground resolution (GSD). The RMSE calculation results show that the accuracy of the DEMs in 1957 and 1989 were lower ( Table 1). The precisions of the glacier volume change are ±2.0403 × 10 −5 km 3 , ±1.0422 × 10 −5 km 3 , and ±1.8738 × 10 −5 km 3 for the periods of 1957-1989, 1989-2015, and 1957-2015, respectively. The mass-balance convention uncertainty σ in the homogenized results from the errors in the geodetic mass change can be calculated as follows (Xu et al., 2017b): where σ ρ is the uncertainty of the density. Δh is the mean of the geodetic elevation changes, and the related uncertainty depends on the accuracy of the two DEMs. σ Δh is the elevation uncertainty in the non-ice region, and it can be calculated as follows: where σ Δh is the standard deviation of the elevation changes for the two DEMs for the stable terrain; N is the number of check points. The MB convention uncertainty is shown in Table 2.

Changes in Glacier Volume and MB
During 1957-1989, the glacier volume increased above 4800 m a.s.l. and decreased below 4800 m a.s.l. At the same time, glacier area above 4800 m a.s.l. except at 5100-5200 m a.s.l. decreased (Figure 2), which means that the ice thickness at those elevation zones were increased obviously. The maximum regions of increase and decrease were 4900-5000 m a.s.l. and 4500-4600 m a.s.l., respectively. During 1989-2015, the entire glacier shrank, and the mean surface elevation in every zone decreased. Overall, the trend of the glacier's surface elevation change in 1957-2015 was the same as in 1989-2015 ( Figure 3A). At 4200-4300 m a.s.l. and 5400-5500 m a.s.l., the terminus and top region of the LHG12 had a small area and thin thickness (Wu et al., 2009;Wu et al., 2011;Wang et al., 2016), so the surface elevation changes had a low magnitude in those areas. During 1957-1989, the maximum surface elevation decreased at the confluence region of LHG12. During 1989-2015, the maximum surface elevation decreased at 4800-5200 m a.s.l., and at the firn basin region, it had significant ablation. Glacier surface elevation changes can be converted to MB using Eq. 1. During 1989-2015, MB was negative in every elevation zone. The most intense ablation region was the terminus area (4300-4400 m a.s.l.). From 1957 to 2015, the change in the MB with altitude was the same as during 1989-2015 ( Figure 3B).
In 1957-2015, the DEM differencing result showed that the trend of the LHG12 volume reduction was accelerating. In 1957-1989, the glacier volume decreased by 0.01 ± 2.0403×10 −5 km 3 , and the average glacier thickness decreased by 0.39 m. In 1989-2015, the glacier volume FIGURE 2 | Distribution of area in elevation zone (1957, 1989, and 2015) with 100 m interval.

Comparison of Simulation and DEM Difference
Two MB simulation results were provided by Chen et al. (2020) (Data C) and Zhang et al. (2018) (Data Z), and there are some differences between them. The correlation between Data C and Data Z is 0.7, the mean averages of Data C and Data Z are −126.79 mm w. e. and −307.86 mm w. e., and their standard deviations are 221.66 and 229.75, respectively. Data C and Data Z were both obtained using the DDF method. For Data C, air temperature data from six stations were used to interpolate the air temperature of the LGH12 using the ordinary kriging interpolation method, and the precipitation data for LHG12 was simulated from the data collected at Tuole station. For Data Z, four national weather stations were used to reconstruct the daily air temperature and precipitation, and the precipitation was reconstructed using precipitation gradients and the inverse distance weight (IDW) method. The DEM differencing result contained the total mass balance of the glacier over a period of time. The results of Data Z show that the MB of LHG12 was −5.52 m w. e. in 1959-1989, −12.25 m w. e. in 1989-2015, and −17.55 in total in 1959-2015. The results of Data C show that the MB of LHG12 was 0.44 m w. e. in 1960-1989, −7.57 m w. e. in 1989-2015, and −7.10 m in total in 1960-2015. For 1957-1989, the results of Data C are similar to the DEM differencing, and the simulation and DEM differencing results show a weak positive and negative MB, respectively. However, the results of Data Z are 1564.25% lower than the DEM differencing result. For 1989-2015, the results of both Data Z and Data C are higher than the DEM differencing result by 14.71 and 47.27%, respectively. For 1957-2015, the results of Data Z are 19.47% lower than the DEM differencing result, but the results of Data C are 51.66% higher than the DEM differencing ( Table 3).
Before the 1980s, the LHG12 experienced a slight acceleration in the rate of retreat, but it stabilized in the 1980s. Since the 1990s, the LHG12 has shrunk rapidly, and the trend has accelerated . The MB simulation results did not show the intense melting process after the 1980s, but Data C produces a more accurate result for 1957-1989. The cumulative mass balance produced using Data Z is close to the DEM difference result, but it is not as accurate as the process simulation, and it overestimates the amount of ablation before 1989 and underestimates the amount of ablation after 1989. The cumulative mass balance obtained using Data C underestimates the amount of ablation, especially the strong melting after the 1980s, but the simulation results are accurate before 1989.

Climate Background
The automatic weather station was located at the LHG12 confluence region (4550 m a.s.l.) with continuous observation data in 2010-2015. The daily temperature between LHG12 and Tuole in the same period is closely associated (R 2 = 0.921) ( Figure 5). Therefore, the temperature and precipitation data of Tuole meteorological station were used to analyze climate background. Based on the 1957-2015 Tuole meteorological station data analysis, in western Qilian Mountains, the temperature and precipitation fluctuations increased; the rates of increase were 0.34°C 10 a −1 and 14.00 mm 10 a −1 , respectively ( Figure 6).
The number of positive and negative temperature days and positive and negative accumulated temperatures in each year were counted respectively. In 1957-2015, the number of positive temperature days increased at a rate of 0.27 day year −1 , and correspondingly, the number of negative temperature days   Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 832701 6 decreased at the same rate; the positive accumulated temperature increased at a rate of 4.89°C year −1 , while the negative accumulated temperature decreased at a rate of 7.70°C year −1 , which was higher than the increase of the positive accumulated temperature (Figure 7).
Temperature variation in season was analyzed, and in spring and summer temperature increased at a rate of 0.02°C a −1 , 0.03°C a −1 , respectively, in autumn the temperature increased at a rate of 0.04°C a −1 , and in winter temperature increased at a rate of 0.05°C a −1 . The temperature rise in winter was significantly higher than the other seasons ( Figure 8). Temperature rise, especially in autumn and winter, could cause the ice temperature of LHG12 glacier to increase and become more sensitive to climate change. Analysis of ice temperature in-situ monitoring data from LHG12 showed a significant increase in ice temperature (Zhu et al., 2019).

CONCLUSION
In this study, two topographic maps and a WorldView-2 satellite image data were used to assess the volume change and cumulative mass balance of the LHG12 located at the western Qilian Mountains in 1957-2015. The conclusions are as follows: 1) During the study period, the changes of the LHG12 included two processes: slight ablation and stability for a brief period in 1957-1989, and strong melting and accelerated ablation in 1989-2015; 2) Due to temperature increase, the changes in region of LHG12 at 4500-4600 m a. s. l and 4800-5200 m a. s. l were the most obvious. In 1957-2015, the volume of LHG12 decreased by 0.38 km 3 and the average thickness thinned to 17.23 m. The cumulative mass balance was −14.69 ± 3.00 m w. e., and ablation was found glacier-wide; 3) Comparing DEM differencing result with previous MB simulation results, the MB simulations underestimated the strong melting trend of LGH12 since the 1990s. It is necessary to analyse the variations in the input parameters of the MB simulation in further research; 4) The increase of temperature, especially in autumn and winter, leads to positive accumulated temperature increase and negative accumulated temperature decrease in a glacier region, and the sensitivity of glaciers to climate change is increased.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
YL, DQ, and XQ: Designed study and final edit. YL: Analyzed meteorological data. ZL and JW: Processed topographic map and generated DEMs. ZJ and LW: Analyzed mass balance data and modified figure. Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 832701