remote sensing Geodetic Mass Balance of Haxilegen Glacier No. 51, Eastern Tien Shan, from 1964 to 2018

: The eastern Tien Shan hosts substantial mid-latitude glaciers, but in situ glacier mass balance records are extremely sparse. Haxilegen Glacier No. 51 (eastern Tien Shan, China) is one of the very few well-measured glaciers, and comprehensive glaciological measurements were implemented from 1999 to 2011 and re-established in 2017. Mass balance of Haxilegen Glacier No. 51 (1999–2015) has recently been reported, but the mass balance record has not extended to the period before 1999. Here, we used a 1:50,000-scale topographic map and long-range terrestrial laser scanning (TLS) data to calculate the area, volume, and mass changes for Haxilegen Glacier No. 51 from 1964 to 2018. Haxilegen Glacier No. 51 lost 0.34 km 2 (at a rate of 0.006 km 2 a − 1 or 0.42% a − 1 ) of its area during the period 1964–2018. The glacier experienced clearly negative surface elevation changes and geodetic mass balance. Thinning occurred almost across the entire glacier surface, with a mean value of − 0.43 ± 0.12 m a − 1 . The calculated average geodetic mass balance was − 0.36 ± 0.12 m w.e. a − 1 . Without considering the error bounds of mass balance estimates, glacier mass loss over the past 50 years was in line with the observed and modeled mass balance ( − 0.37 ± 0.22 m w.e. a − 1 ) that was published for short time intervals since 1999 but was slightly less negative than glacier mass loss in the entire eastern Tien Shan. Our results indicate that Riegl VZ ® -6000 TLS can be widely used for mass balance measurements of unmonitored individual glaciers.


Introduction
Glaciers store a large amount of fresh water and provide an indispensable supply of water to arid and semi-arid regions such as the Tien Shan [1]. The concept of "natural solid reservoir" is well represented in the Tien Shan, which is heavily glacierized [2,3]. Glaciers in the Tien Shan have been retreating since the end of the Little Ice Age in the context of climate warming [4]. The continuous retreat and subsequent disappearance of glaciers strongly affect the quantity and seasonal distribution of runoff in Central Asia's glacier-fed basins [4,5]. Glacier mass balance measurements are widely recognized as being key to studies of climate-glacier interactions, water resources and sea-level rise [5,6]. The annual and seasonal mass balance of individual glaciers can be measured using the direct glaciological method (i.e., measurements recorded from repeated snow pits observations and dense arrays of ablation stakes drilled into the glacier surface). However, there are only seven glaciers in the Tien Shan with a history of in situ mass balance measurements that extends more than five years, and two glaciers in the Tien Shan have glaciological mass balance records that cover more than 30 years (Tuyuksu glacier in Kazakhstan and Urumqi Glacier No. 1 in China) [3]. Thus, regional glacier change assessments are challenged by the small number and heterogeneous spatiotemporal distribution of in situ measurement series [7][8][9][10][11].
The geodetic mass balance (difference in surface elevation obtained from repeated geodetic survey) can be considered as complementary information to glaciological mass balance measurements [12]. Previous studies in the Tien Shan have focused on geodetic mass balance estimates for individual glaciers in the northern and central regions (e.g., [13][14][15][16][17][18]). Long-term geodetic mass balance estimates for individual glaciers in the eastern Tien Shan have not been systematically reported. To the best of our knowledge, only the mass balance of Urumqi Glacier No. 1 has been calculated using the geodetic method [19,20]. Although glacier mass change at regional and global scales can be estimated using remotely sensed geodetic measurements, such estimates need to be validated using reliable glacier mass balance measurements for individual glaciers [6,8]. The quality of estimates of geodetic mass changes depends strongly on the resolution and accuracy of the digital elevation models (DEMs) that are used [20,21]. Many available global DEMs, such as those from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM), Shuttle Radar Topography Mission (SRTM) DEM and TanDEM-X, have large vertical errors in mountainous areas with complex terrain [22][23][24]. Other global DEMs, such as the Advanced Land Observing Satellite (ALOS) World AW3D30 DEM, also contain a widespread elevation anomaly [25]. Such DEMs cannot resolve microtopographic variations of surface elevation for individual glaciers, especially in the accumulation area [11]. On the other hand, the penetration of the SRTM C-Band radar signal into snow and ice can introduce a bias to estimates of glacier surface elevation changes based on radar-derived DEM [11]. Data acquisition techniques, such as real-time kinematic global navigation satellite system (RTK-GNSS), unmanned aerial vehicle (UAV), and light detection and ranging (LiDAR) can generate high-resolution DEMs with sub-centimeter vertical accuracy [23]. Compared with RTK-GNSS, UAV and LiDAR are less time-consuming and can rapidly survey glacier surface terrain, and measurements from these instruments are increasingly being used to monitor glacier changes (e.g., [26][27][28][29]). The reliability of UAV surveys is dependent on ground control points (GCPs). If GCPs are unevenly distributed or have low precision, then the UAV-derived glacier mass balance may have higher uncertainty [28,29]. LiDAR is also called laser scanning, including airborne (ALS), mobile (MLS), and terrestrial laser scanning (TLS). MLS is not applicable in remote mountain areas because the instrument is mounted on a mobile platform, which is usually a ship, train, or car [30]. ALS can be used to rapidly map extensive areas but is limited to operating at flight altitude. Great topographic relief and high mountains often limit geodetic surveys made by aircraft because most ALS instruments have a limited operating flight altitude [20]. TLS is a flexible method with high measurement accuracy and has been increasingly used to monitor glacier mass changes at seasonal, monthly, and hourly scales (e.g., [27,31,32]).
A long-range Riegl VZ ® -6000 TLS was utilized to survey Haxilegen Glacier No. 51 in 2018. In situ measurements for Haxilegen Glacier No. 51 have been performed since 1999, making it one of the very few well-measured glaciers in the eastern Tien Shan. However, measurements stopped in 2011, and its mass balance has been reconstructed from 1999 to 2015 using a temperature index and an accumulation model [33]. Currently, the mass balance record for Haxilegen Glacier No. 51 has not been extended to the period before 1999, and this means that long-term mass changes in the glacier are poorly understood. In addition, to the best of our knowledge, previous studies of Haxilegen Glacier No. 51 have focused primarily on changes in its geometry (length, area, and thickness) (e.g., [34][35][36][37][38]). In this study, we extend the mass balance record of Haxilegen Glacier No. 51 using a 1:50,000-scale topographic map constructed from aerial photography and the Riegl VZ ® -6000 TLS-derived DEM. The main objectives of our study are (1) to give a detailed estimate of changes to surface elevation and mass for Haxilegen Glacier No. 51 over the past 54 years, (2) to compare these estimates with local and regional geodetic mass balance estimates, and (3) to infer potential drivers for the observed glacier mass changes.

Study Site
Tien Shan is the largest latitudinal mountain system in Central Asia, stretching across China, Kazakhstan, Kyrgyzstan, and Uzbekistan ( Figure 1a). It is far away from the ocean and covers a larger spatial area (the total length is about 2500 km and the average width from north to south is about 250-350 km). The total glacierized area in the Tien Shan is 11,864 km 2 , according to the latest Randolph Glacier Inventory V6.0 [39]. Haxilegen Glacier No. 51 (43.73 • N, 84.39 • E) is a small northeast-facing cirque glacier located in the headwaters of the Kuytun River in the eastern Tien Shan (Figure 1a). The glacier covers an altitudinal range of 3400 to 4000 m a.s.l., with a total area of around 1.48 km 2 and a length of 1.70 km in 1964 [36]. The glacier occupies a wide valley of~1.5 km at higher elevation parts and becomes narrower (<0.9 km) at lower altitudes ( Figure 1b). The surface terrain is relatively flat at the lower elevations, with an average slope of~10 • . In situ measurements for Haxilegen Glacier No. 51 include glaciological mass balance, surface velocity, thickness, and variations in the surface area and location of the front [33,34]. For glaciological mass balance measurements, approximately 20 well-distributed ablation stakes were drilled into the glacier in August 2017 using a stream drill (Figure 1b). Surveys made by groundpenetrating radar in 2010 showed that the maximum and mean thicknesses of the glacier were 73 m and 39 m, respectively [36]. The glacier surface velocity ranged from 1.5 to 3.1 m a -1 during 2000-2006, based on RTK-GNSS surveys [36].

Study Site
Tien Shan is the largest latitudinal mountain system in Central Asia, stretching across China, Kazakhstan, Kyrgyzstan, and Uzbekistan ( Figure 1a). It is far away from the ocean and covers a larger spatial area (the total length is about 2500 km and the average width from north to south is about 250-350 km). The total glacierized area in the Tien Shan is 11,864 km 2 , according to the latest Randolph Glacier Inventory V6.0 [39]. Haxilegen Glacier No. 51 (43.73°N, 84.39°E) is a small northeast-facing cirque glacier located in the headwaters of the Kuytun River in the eastern Tien Shan (Figure 1a). The glacier covers an altitudinal range of 3400 to 4000 m a.s.l., with a total area of around 1.48 km 2 and a length of 1.70 km in 1964 [36]. The glacier occupies a wide valley of ~1.5 km at higher elevation parts and becomes narrower (<0.9 km) at lower altitudes ( Figure 1b). The surface terrain is relatively flat at the lower elevations, with an average slope of ~10°. In situ measurements for Haxilegen Glacier No. 51 include glaciological mass balance, surface velocity, thickness, and variations in the surface area and location of the front [33,34]. For glaciological mass balance measurements, approximately 20 well-distributed ablation stakes were drilled into the glacier in August 2017 using a stream drill (Figure 1b). Surveys made by ground-penetrating radar in 2010 showed that the maximum and mean thicknesses of the glacier were 73 m and 39 m, respectively [36]. The glacier surface velocity ranged from 1.5 to 3.1 m a -1 during 2000-2006, based on RTK-GNSS surveys [36]. Haxilegen Glacier No. 51 is a summer-accumulation-type glacier (strong ablation and accumulation occur simultaneously during summer) in a continental climate setting [33,36]. The climate setting was confirmed by the annual climate records (1965-2017) of Haxilegen Glacier No. 51 is a summer-accumulation-type glacier (strong ablation and accumulation occur simultaneously during summer) in a continental climate setting [33,36]. The climate setting was confirmed by the annual climate records (1965-2017) of Usu Meteorological Station, which is located approximately 76 km northeast of Haxilegen Glacier No. 51 at an altitude of 1044 m a.s.l. and is approximately 2400 m a.s.l. below the glacier terminus in 2018. Based on climatic data from the station, the climate was characterized by very low air temperature in January (cold month) with a mean value of 15.9 • C and by high air temperatures during the summer months (May-August), with the warmest month being July (25.0 • C). More precipitation was observed in summer (May-September) than in winter, and winter precipitation accounted for only 13.7% of the total annual precipitation.

DEM Generated from Riegl VZ ® -6000 TLS
The novel long-range Riegl VZ ® -6000 TLS emits near-infrared laser signals (3B laser beams with a wavelength of 1064 nm), which enable high rates of reflection (>80%) for snowy and icy terrain [40]. The instrument can provide a maximum scanning range of 6 km with a spatial resolution better than 0.0005 • (1.8 arcsec). The scanned distances of Haxilegen Glacier No. 51 range from 80 to 1600 m. In order to prevent ground motion and enhance data acquisition quality, we performed TLS surveys at two fixed scan positions located in the stable terrain area with a wide view relative to the whole glacier surface ( Figure 1b). The instrument was mounted on a tripod for each scanning position. The viewing geometry has a significant influence on point cloud data quality because it can influence the working distance of the Riegl VZ ® -6000 and cause data voids. The instrument cannot receive the echo of the laser pulse when the angle of incidence is less than 10 • [40]. The glacier surface terrain is generally steep, relative to the two scan positions, which guaranteed the requirements of the angle of incidence, and therefore meant that most of the glacier surface could be surveyed. However, the left upper corner of the glacier was invisible from the two scan positions (Figure 1c), which resulted in data voids (the unscanned areas account for 7.3% of the whole glacier). In this study, these unscanned areas were not considered for the geodetic mass balance estimates. The accuracy of the point clouds relies on the accuracy of three-dimensional (3D) coordinates of the two scan positions. Prior to the laser scanning, we used the RTK-GNSS of Trimble R10 to survey 3D coordinates of two scan positions. The accuracy of the RTK-GNSS has been reported to be within ±1 cm horizontally and ±2 cm vertically [41]. The scan campaigns were performed at the end of the ablation season under an overcast sky to reduce the influence of large areas of fresh snow and atmospheric refraction [40,41]. We implemented the TLS surveys on 2 September 2018 with a scanning angle resolution of 0.02 • , resulting in a total of 58,684,617 points, with an average density of 18.67 points m −2 .
We generated a DEM based on point clouds of the Riegl VZ ® -6000 TLS survey. RiSCAN PRO ® v 1.81 software was used to process the cloud data. The procedures include direct georeferencing, data registration, data compression, data filtering and interpolation. We first converted the point cloud data from the scanner's own coordinate system into the global coordinate system using the RTK-GNSS surveyed 3D coordinates for the two scan positions. We then used multistation adjustment to complete data registration based on the iterative closest point algorithm, which iteratively adjusts the orientation between the two scans until the point cloud of the two scans is matched in space, at which point we merged the two scans into one layer. Then, we used an octree algorithm to compress the layer to generate points with equal spacing and to complete the data compression [20,41]. We used the terrain filter to remove noise and nonground data, and we checked the processed data and removed clear visual outliers through visual inspection. Finally, the processed point clouds were interpolated to generate a DEM using the coordinate system of the World Geodetic System (WGS) 84 datum and Universal Transverse Mercator (UTM) 45N projection, with a spatial resolution of 0.5 m. This was downscaled to 30 m resolution to match the spatial resolution of the DEM generated from a topographic map (see Section 3.2). The specific procedures of point cloud processing are described in detail in [41].

DEM Generated from the Topographic Map and Co-Registration
The 1:50,000-scale historical topographic map produced by the Chinese Military Geodetic Service in September 1964 was scanned at a resolution of 600 dpi to create a DEM for the past time. We used the one-kilometer gridlines to georeference the scanned map. The original coordinate system for the map was Beijing 1954, and the elevation system was the Yellow Sea system 1956. We reprojected the georeferenced map onto the WGS84 datum and UTM 45N projection to match the coordinate system for the TLS surveys. The topographic map and the Riegl VZ ® -6000 TLS use different surveying techniques, which may result in horizontal and vertical deviations between the two datasets [42]. Thus, we selected several clearly distinguishable stable features on the high-resolution TLS point clouds and used these to register the topographic map. Subsequently, we manually digitized 20 m elevation contours from the map. We used the digital contour lines to create a triangulated irregular network (TIN) via the Create TIN module from ArcGIS 10.2. Using the TIN To Raster module from ArcGIS 10.2, we generated a raster DEM with 30 m resolution by interpolating cell values from the elevations of the created TIN [20]. The nominal vertical accuracy (root mean square error of the elevations) of the topographic map was reported to be better than 5 m for flat and hilly areas with slopes of less than 6 • , and 8 m for mountainous areas with slopes of 6-25 • , according to the photogrammetric Chinese National Standard [43]. The vertical accuracy of the DEM generated from the topographic map should therefore be better than 8 m, as the mean slope was approximately 17.6 • in 1964.

Glacier Boundary Delineation
We manually delineated the glacier boundaries in 1964 and 2018 on the basis of the processed topographic map and TLS data. We used shaded reliefs based on the TLS-derived high-resolution DEM to delineate the glacier boundary in 2018 ( Figure 2). The processes are introduced in detail in [41]. The glacier boundary derived from a Landsat Thematic Mapper image (acquisition date: 8 August 2009) was taken from the second Chinese Glacier Inventory as an additional cross-reference and used to improve the digitizing accuracy for the unscanned areas [44].
topographic map (see Section 3.2). The specific procedures of point cloud processing are described in detail in [41].

DEM Generated from the Topographic Map and Co-Registration
The 1:50,000-scale historical topographic map produced by the Chinese Military Geodetic Service in September 1964 was scanned at a resolution of 600 dpi to create a DEM for the past time. We used the one-kilometer gridlines to georeference the scanned map. The original coordinate system for the map was Beijing 1954, and the elevation system was the Yellow Sea system 1956. We reprojected the georeferenced map onto the WGS84 datum and UTM 45N projection to match the coordinate system for the TLS surveys. The topographic map and the Riegl VZ ® -6000 TLS use different surveying techniques, which may result in horizontal and vertical deviations between the two datasets [42]. Thus, we selected several clearly distinguishable stable features on the high-resolution TLS point clouds and used these to register the topographic map. Subsequently, we manually digitized 20 m elevation contours from the map. We used the digital contour lines to create a triangulated irregular network (TIN) via the Create TIN module from ArcGIS 10.2. Using the TIN To Raster module from ArcGIS 10.2, we generated a raster DEM with 30 m resolution by interpolating cell values from the elevations of the created TIN [20]. The nominal vertical accuracy (root mean square error of the elevations) of the topographic map was reported to be better than 5 m for flat and hilly areas with slopes of less than 6°, and 8 m for mountainous areas with slopes of 6-25°, according to the photogrammetric Chinese National Standard [43]. The vertical accuracy of the DEM generated from the topographic map should therefore be better than 8 m, as the mean slope was approximately 17.6° in 1964.

Glacier Boundary Delineation
We manually delineated the glacier boundaries in 1964 and 2018 on the basis of the processed topographic map and TLS data. We used shaded reliefs based on the TLS-derived high-resolution DEM to delineate the glacier boundary in 2018 ( Figure 2). The processes are introduced in detail in [41]. The glacier boundary derived from a Landsat Thematic Mapper image (acquisition date: 8 August 2009) was taken from the second Chinese Glacier Inventory as an additional cross-reference and used to improve the digitizing accuracy for the unscanned areas [44].

Geodetic Mass Balance Calculations
The geodetic mass balance (B geod ), in units of m w.e. (meters water equivalent) was calculated for Haxilegen Glacier No. 51 by multiplying the glacier volume change (∆V) and a volume-to-mass conversion factor (ρ) over the glacier area: where S is the mean glacier area over the two acquisition years 1964 and 2018, assuming a linear change over the period [6]; ρ water is the density of water; and ρ is the average density of ∆V, for which we assume the widely recommended value of 850 ± 60 kg m −3 [45]. This assumption is valid because the geodetic surveys span more than 5 years [45]. We used the glacier boundary in 2018 to determine elevation changes due to the existence of a small glacial lake at the left fore of the glacier, which resulted in some data voids of the point cloud.

Uncertainties in Glacier Area
Uncertainties in the glacier area include systematic and random errors [46]. The uncertainties caused by snow cover, mountain shadow, side debris, and human misinterpretation can be regarded as systematic errors. Selecting high-quality remote sensing images can reduce these uncertainties, but the total systematic error cannot be accurately estimated [46]. Uncertainties related to remote sensing image resolution, spectral characteristics and debris cover can be considered as random errors [44]. We primarily considered the uncertainty caused by image resolution as the glacier surface was clean. The uncertainty (σ area ) can be expressed as where r is the spatial resolution of the topographic map and the TLS-derived DEM (30 m); N pixel is the total number of pixels passed by the glacier boundary [44]. Resulting values for σ area are given in Section 4.1.

Uncertainties in Geodetic Mass Balance
The uncertainties in glacier surface elevation and mass change depend mainly on the DEM differencing and volume-to-mass conversion processes [12,47]. Thus, the total uncertainty in the geodetic mass balance (σ geod ) over the period of record can be estimated as [48] where ∆h denotes the average glacier-wide change in surface elevation, which is equal to ∆V S ; σ ρ denotes the uncertainty in volume-to-mass conversion, which is set to 60 kg m −3 according to [44]. σ ∆h is the uncertainty in glacier surface elevation change, which can be evaluated using precise and well-distributed GCPs over the stable terrain [49]. Here, the standard deviation of the elevation changes over the stable terrain is used to estimate σ ∆h due to a lack of GCPs, and the corresponding error can be expressed as [50] where σ ∆hi is the standard deviation of elevation change over the stable terrain, and N eff is the standard deviation of the mean elevation change ∆h (Figure 3). N eff is the effective number of independent measurements and can be calculated as where N total is the total number of measurements; R is the spatial resolution of the applied DEM (30 m), and d is the distance of spatial autocorrelation (686 m), which is calculated using Moran's autocorrelation index [20,50].
DEM (30 m), and d is the distance of spatial autocorrelation (686 m), which is calculated using Moran's autocorrelation index [20,50]. Based on the law of error propagation, the average annual uncertainty in the geodetic mass balance ( . ) is calculated as [12] . = (6) where n is the number of years spanned by the geodetic surveys. The results are summarized in Table 1.

Changes in Glacier Area
The glacier area was 1.48 ± 0.09 km 2 in 1964 and 1.14 ± 0.08 km 2 in 2018; thus, the glacier lost 0.34 km 2 of its area over this period, at an average rate of 0.006 km 2 , or 0.4% per year. Most of the glacier area was lost from the lower elevations of the glacier ( Figure  1b; Table 1). A previous study reported that the glacier area was reduced by 0.12 km 2 between 1964 and 2006 at a rate of 0.003 km 2 or 0.2% per year [34], which suggests that glacier shrinkage has accelerated over the last ten years (Figure 4). Based on the law of error propagation, the average annual uncertainty in the geodetic mass balance (σ geod.a ) is calculated as [12] σ geod.a = σ geod n where n is the number of years spanned by the geodetic surveys. The results are summarized in Table 1. Table 1. Average annual glacier surface elevation changes (∆h), standard deviation (σ ∆h ) in surface elevation changes over stable terrain, average annual uncertainty in the geodetic mass balance ( σ geod.a ), and the corresponding annual geodetic mass balance (B geod.a ).

Changes in Glacier Area
The glacier area was 1.48 ± 0.09 km 2 in 1964 and 1.14 ± 0.08 km 2 in 2018; thus, the glacier lost 0.34 km 2 of its area over this period, at an average rate of 0.006 km 2 , or 0.4% per year. Most of the glacier area was lost from the lower elevations of the glacier (Figure 1b; Table 1). A previous study reported that the glacier area was reduced by 0.12 km 2 between 1964 and 2006 at a rate of 0.003 km 2 or 0.2% per year [34], which suggests that glacier shrinkage has accelerated over the last ten years (Figure 4).  The picture in 2006 is from [35]. The red frame indicates corresponding areas of the two years, where the exposed rock area increased significantly. A small glacial lake was observed at the fore of the glacier in 2018.

Changes in Glacier Surface Elevation and Mass Balance
Our results suggest that Haxilegen Glacier No. 51 experienced substantial thinning and mass loss during the period 1964-2018 ( Figure 5). The mean thinning of Haxilegen Glacier No. 51 was −23.09 ± 6.49 m or −0.43 ± 0.12 m a −1 , which was equivalent to 31.6% of the mean thickness of the glacier surveyed in 2010 [36]. The thickness reduced over almost the whole glacier area, and the most pronounced thinning occurred over the lower-elevation parts of the glacier, with the greatest reductions ranging from −60 to −50 m during the period 1964-2018 (Figures 1b and 5). Some clear surface lowering (approximately −30 to −20 m) was also observed in the middle-upper parts of the glacier, where the glacier was very thin, and where the exposed rock area increased significantly with glacier shrinkage (Figure 4). The picture in 2006 is from [35]. The red frame indicates corresponding areas of the two years, where the exposed rock area increased significantly. A small glacial lake was observed at the fore of the glacier in 2018.

Changes in Glacier Surface Elevation and Mass Balance
Our results suggest that Haxilegen Glacier No. 51 experienced substantial thinning and mass loss during the period 1964-2018 ( Figure 5). The mean thinning of Haxilegen Glacier No. 51 was −23.09 ± 6.49 m or −0.43 ± 0.12 m a −1 , which was equivalent to 31.6% of the mean thickness of the glacier surveyed in 2010 [36]. The thickness reduced over almost the whole glacier area, and the most pronounced thinning occurred over the lower-elevation parts of the glacier, with the greatest reductions ranging from −60 to −50 m during the period 1964-2018 (Figures 1b and 5). Some clear surface lowering (approximately −30 to −20 m) was also observed in the middle-upper parts of the glacier, where the glacier was very thin, and where the exposed rock area increased significantly with glacier shrinkage (Figure 4). The picture in 2006 is from [35]. The red frame indicates corresponding areas of the two years, where the exposed rock area increased significantly. A small glacial lake was observed at the fore of the glacier in 2018.

Changes in Glacier Surface Elevation and Mass Balance
Our results suggest that Haxilegen Glacier No. 51 experienced substantial thinning and mass loss during the period 1964-2018 ( Figure 5). The mean thinning of Haxilegen Glacier No. 51 was −23.09 ± 6.49 m or −0.43 ± 0.12 m a −1 , which was equivalent to 31.6% of the mean thickness of the glacier surveyed in 2010 [36]. The thickness reduced over almost the whole glacier area, and the most pronounced thinning occurred over the lower-elevation parts of the glacier, with the greatest reductions ranging from −60 to −50 m during the period 1964-2018 (Figures 1b and 5). Some clear surface lowering (approximately −30 to −20 m) was also observed in the middle-upper parts of the glacier, where the glacier was very thin, and where the exposed rock area increased significantly with glacier shrinkage (Figure 4). From 1964 to 2018, the total volume loss was approximately 0.034 ± 0.01 km 3 , and the corresponding water equivalent was approximately 28.90 ± 8.5 × 10 6 m 3 . The glacier wastage resulted in a negative cumulative geodetic mass balance of −19.63 ± 5.69 m w.e., with a mean annual change of −0.36 ± 0.12 m w.e. a −1 . A glacial lake was observed in 2018, which may be related to Haxilegen Glacier No. 51 having experienced mass loss over the period (Figure 4). The altitudinal distribution of glacier mass balance indicates that specific mass balance increased with elevation, and reduced glacier mass loss was observed at the glacier terminus ( Figure 6). This is a commonly observed pattern and agrees well with glaciological mass balance measurements for some "reference" glaciers (with >30 ongoing observations) [51]. From 1964 to 2018, the total volume loss was approximately 0.034 ± 0.01 km 3 , and the corresponding water equivalent was approximately 28.90 ± 8.5 × 10 6 m 3 . The glacier wastage resulted in a negative cumulative geodetic mass balance of −19.63 ± 5.69 m w.e., with a mean annual change of −0.36 ± 0.12 m w.e. a −1 . A glacial lake was observed in 2018, which may be related to Haxilegen Glacier No. 51 having experienced mass loss over the period (Figure 4). The altitudinal distribution of glacier mass balance indicates that specific mass balance increased with elevation, and reduced glacier mass loss was observed at the glacier terminus (Figure 6). This is a commonly observed pattern and agrees well with glaciological mass balance measurements for some "reference" glaciers (with >30 ongoing observations) [51].

Comparison with Previous Studies
Haxilegen Glacier No. 51 retreated at a rate higher compared to the global average rate for glacier retreat (0.35% a −1 ), greater than the average retreat rate for glaciers in China (~0.36% a −1 ), and higher than the average rate for glaciers in the Chinese Tien Shan (0.37% a −1 ) over the last 50 years [46,52,53]. Several studies have reported the glacier shrinkage of individual glaciers in the Tien Shan. For example, an area reduction rate of 0.41% a −1 for the Tuyuksu glacier was observed between 1975 and 2006 [14], and a retreat rate of 0.36% a −1 was observed during the period 1962-2018 for Urumqi Glacier No. 1 [54], which was slightly smaller than the retreat rate of Haxilegen Glacier No. 51. The area of Haxilegen Glacier No. 51 is relatively small, which may explain the large percentage value for the relative glacier retreat rate.
As mentioned in Section 1, available glaciological mass balance records are scarce in the Tien Shan. Previous studies reconstructed the mass balance of Haxilegen Glacier No. 51 from 1999 to 2015 based on glaciological mass balance measurements of seven discontinuous hydrological years between 1999 and 2011, and an average mass balance of −0.37 ± 0.22 m w.e. a −1 was found, which is in the same range as the mass balance estimates in our study [33]. Haxilegen Glacier No. 51 thinned by ~10 m (~0.28 m w.e. a −1 ) from the 1980s to 2010 [36]. Therefore, we can conservatively assume that glacier mass changes for the two periods 1964-1999 and 1999-2015 were similar, indicating that the mass loss rate for Haxilegen Glacier No. 51 had not increased over the past five decades. The nearest well-

Comparison with Previous Studies
Haxilegen Glacier No. 51 retreated at a rate higher compared to the global average rate for glacier retreat (0.35% a −1 ), greater than the average retreat rate for glaciers in China (~0.36% a −1 ), and higher than the average rate for glaciers in the Chinese Tien Shan (0.37% a −1 ) over the last 50 years [46,52,53]. Several studies have reported the glacier shrinkage of individual glaciers in the Tien Shan. For example, an area reduction rate of 0.41% a −1 for the Tuyuksu glacier was observed between 1975 and 2006 [14], and a retreat rate of 0.36% a −1 was observed during the period 1962-2018 for Urumqi Glacier No. 1 [54], which was slightly smaller than the retreat rate of Haxilegen Glacier No. 51. The area of Haxilegen Glacier No. 51 is relatively small, which may explain the large percentage value for the relative glacier retreat rate.
As mentioned in Section 1, available glaciological mass balance records are scarce in the Tien Shan. Previous studies reconstructed the mass balance of Haxilegen Glacier No. 51 from 1999 to 2015 based on glaciological mass balance measurements of seven discontinuous hydrological years between 1999 and 2011, and an average mass balance of −0.37 ± 0.22 m w.e. a −1 was found, which is in the same range as the mass balance estimates in our study [33]. Haxilegen Glacier No. 51 thinned by~10 m (~0.28 m w.e. a −1 ) from the 1980s to 2010 [36]. Therefore, we can conservatively assume that glacier mass changes for the two periods 1964-1999 and 1999-2015 were similar, indicating that the mass loss rate for Haxilegen Glacier No. 51 had not increased over the past five decades. The nearest well-monitored glacier is Urumqi Glacier No. 1, located approximately 200 km east of Haxilegen Glacier No. 51 (Figure 1a). Urumqi Glacier No. 1 experienced an annual mass loss of −0.37 ± 0.20 m w.e. a −1 from 1964 to 2018, based on glaciological mass balance measurements, which is similar to the mass balance estimates in our study [51].
The available mass balance data series for the central Tien Shan showed slightly more negative mass balances than for Haxilegen Glacier No. 51, especially Sary-Tor Glacier. The long-term glaciological measurements of the Tuyuksu glacier suggested that the glacier experienced mass loss at a rate of −0.43 ± 0.20 m w.e. a −1 , which is also slightly more negative than that experienced at Haxilegen Glacier No. 51; however, these mass balances agree well when error bounds are considered (within the 95% confidence interval, p = 0.017) ( Table 2). The differences in mass loss may be related to local climate and the aspect, size, geometry, and hypsometry of the investigated glaciers [55]. The wider accumulation zone of Haxilegen Glacier No. 51 may favor glacier accumulation [33]. Additionally, snowdrifts or avalanches can strengthen the accumulation of the small cirque glacier [33]. Table 2. Comparison of available mass balance series for glaciers in Tien Shan (TS). The locations of those glaciers are marked in Figure 1, and the unit of mass balance and its uncertainty is m w.e. a −1 . 1 The uncertainty of glaciological mass balance was estimated according to [18].

Glacier
A regional average mass balance rate of −0.50 ± 0.37 m w.e. a −1 was found for glaciers in the eastern Tien Shan for the period 1999-2000 to 2017-2018 [8]. Average mass loss rates of −0.32 ± 0.30 m w.e. a −1 and −0.41 ± 0.25 m w.e. a −1 were reported in the Borohoro range (a subregion in the eastern Tien Shan) for the years 1961 to 2012 and 2003 to 2009, respectively, according to glaciological modeling [3]. Mass budgets of −0.28 ± 0.20 w.e. a −1 and −0.40 ± 0.20 w.e. a −1 were estimated for glaciers in the whole Tien Shan and eastern Tien Shan, respectively, using satellite stereoimagery for the period 2000-2016 [57]. The aforementioned studies all estimated similar regional glacier mass balance for the last two decades in the eastern Tien Shan, which were more negative than the glacier mass changes for the whole Tien Shan. Mean mass losses of −0.40 ± 0.09 w.e. a −1 and −0.40 ± 0.07 w.e. a −1 were found for glaciers in the northern Tien Shan and in Ak Shirak (a subregion in the central Tien Shan) from 1964 to 2018, which were slightly more negative than the mass balance estimates for Haxilegen Glacier No. 51 [58]. The compared results only considered the absolute mass balance, and most of these estimates would be consistent if we considered the error bounds (0.07-0.37 w.e. a −1 ) ( Table 2). However, errors in mass balance estimates cannot be eliminated and advanced data acquisition technologies, such as UAV and LiDAR, should be used to improve the accuracy of mass balance measurements. Using these technologies may make it easier to compare glacier mass balance estimates for different regions, as the error bounds may become negligible.

Possible Drivers for Glacier Mass Changes
Glacier mass balance is deemed to be an essential climate variable and is acknowledged to be one of the most direct indicators of climate change [12,59]. Summer (May-August) air temperature and annual precipitation were the main contributors to mass balance fluctuation for Haxilegen Glacier No. 51 over the period 1999-2015 [33]. We analyzed data recorded at the Usu Meteorological Station to investigate possible drivers for the glacier mass changes from 1964 to 2018 (Figure 7). The annual average air temperature increased at a rate of 0.3 • C (10 a) −1 , which is approximately twice the average global warming rate over the same period [60]. The annual precipitation had been increasing gradually with an average rate of 12 mm (10 a) −1 . A particularly abrupt change in air temperature occurred in the 1990s, and the warming trend has been more significant since then (p < 0.01), with a decadal average increased value of 0.91 • C, according to data from the Usu Meteorological Station. Previous studies have suggested that mass balance for Haxilegen Glacier No. 51 is more sensitive to air temperature than precipitation [33,36]. We conservatively assumed that abnormal air temperature rise may be an important driver for glacier mass loss of Haxilegen Glacier No. 51. edged to be one of the most direct indicators of climate change [12,59]. Summer (May-August) air temperature and annual precipitation were the main contributors to mass balance fluctuation for Haxilegen Glacier No. 51 over the period 1999-2015 [33]. We analyzed data recorded at the Usu Meteorological Station to investigate possible drivers for the glacier mass changes from 1964 to 2018 (Figure 7). The annual average air temperature increased at a rate of 0.3 °C (10 a) −1 , which is approximately twice the average global warming rate over the same period [60]. The annual precipitation had been increasing gradually with an average rate of 12 mm (10 a) −1 . A particularly abrupt change in air temperature occurred in the 1990s, and the warming trend has been more significant since then (p < 0.01), with a decadal average increased value of 0.91 °C, according to data from the Usu Meteorological Station. Previous studies have suggested that mass balance for Haxilegen Glacier No. 51 is more sensitive to air temperature than precipitation [33,36]. We conservatively assumed that abnormal air temperature rise may be an important driver for glacier mass loss of Haxilegen Glacier No. 51. Some studies have suggested that glacier surface albedo reduction is an important driver for continuous and accelerated glacier mass loss in Tibetan Plateau and its surrounding areas [60][61][62]. The average albedo of glaciers in these areas decreases gradually as the ablation area increases, which enhances shortwave radiation heating and results in a further expansion of the ablation zone [60]. Surface albedo reduction can enhance glacier ablation by approximately 30-60%, and both annual and summer glacier surface albedo in the eastern Tien Shan presented an obvious decreasing trend from 2001 to 2018 [62]. This may be a possible driver for mass changes of Haxilegen Glacier No. 51, but the lack of observed surface albedo prevents a more detailed quantitative analysis.
For glaciers in the Tibetan Plateau and its surrounding areas, previous studies have also suggested that supraglacial rivers may influence energy exchange between the atmosphere and glacier, and thereby increase the shortwave radiation energy received by Some studies have suggested that glacier surface albedo reduction is an important driver for continuous and accelerated glacier mass loss in Tibetan Plateau and its surrounding areas [60][61][62]. The average albedo of glaciers in these areas decreases gradually as the ablation area increases, which enhances shortwave radiation heating and results in a further expansion of the ablation zone [60]. Surface albedo reduction can enhance glacier ablation by approximately 30-60%, and both annual and summer glacier surface albedo in the eastern Tien Shan presented an obvious decreasing trend from 2001 to 2018 [62]. This may be a possible driver for mass changes of Haxilegen Glacier No. 51, but the lack of observed surface albedo prevents a more detailed quantitative analysis.
For glaciers in the Tibetan Plateau and its surrounding areas, previous studies have also suggested that supraglacial rivers may influence energy exchange between the atmosphere and glacier, and thereby increase the shortwave radiation energy received by the glacier and expand the effective ablation area [60]. Based on field investigations in September 2018, visual interpretation found that supraglacial drainage developed in the ablation area of the glacier (Figure 4b), and this may increase ablation area. Glacier surface velocity in the Tien Shan decreased by −6.4 ± 1.0% (10 a) −1 over the period 2000-2017, a reduction in glacier velocity may reduce the supply of mass to the ablation area, and this may have contributed to the thinning at the glacier terminus [63]. Haxilegen Glacier No. 51 exhibited a slowing trend of glacier movement from 2000 to 2006 based on in situ measurements [36]. Slowdown for Haxilegen Glacier No. 51 may influence glacier mass loss. However, the possible drivers for glacier mass changes that are suggested in this paragraph are quite speculative. Further studies, including detailed observations and physics-based models, are needed to better understand these potential drivers for mass change.

Implications and Recommendations
Our results suggest that most of the glacier surface can be surveyed ( Figure 3) and demonstrate that the Riegl VZ ® -6000 TLS allows more detailed insight into the glacier's evolution when repeat surveys are used [46]. Therefore, we recommend the increased use of Riegl VZ ® -6000 TLS to monitor glacier mass balance. Nevertheless, there are both advantages and disadvantages to the TLS survey technique. It is not always easy to find suitable scan positions, where good visual angles can be obtained relative to the glacier surface. Thus, there are often data voids (unscanned areas), even for very small glaciers (e.g., [48]). Water vapor in the air (rain, cloud, and fog) can absorb the laser pulse and reduce the possible survey distance, which may result in further data voids. In addition, wind can influence the accuracy of the point cloud. These weather conditions are common in glaciated areas and the dry, fogless, and windless conditions that are ideal for TLS surveys are often not available [41,48]. In addition, the absence of substantial firn or fresh snow at the times of the TLS survey may increase the uncertainty of the resulting annual or seasonal geodetic mass balance estimates, because the recommended volume-to-mass conversion factor (density) is only appropriate where geodetic observations span more than 5 years [44]. We can select the end of the ablation season, when a small amount of fresh snow covers the glacier, to reduce the uncertainty in the assumed density. Ideally, the density of snow-firn in the accumulation area should be measured to establish an appropriate volume-to-mass conversion [32]. Our efforts at Haxilegen Glacier No. 51 included continuous measurements of glacier velocity because glacier emergence/submergence velocities have a significant impact on glacier surface evolution [41]. Foreseeable developments and future planned field studies include pasting several retroreflective targets to the surface of each ablation stake, increasing their reflectivity for incident laser radiation so that their 3D coordinates can be more easily surveyed. We hope that this will allow us to determine surface velocity more easily. Some recent studies propose the use of repeated snow line observations obtained from terrestrial automatic cameras to reconstruct glacier mass balance [16]. We plan to acquire snow line observations by installing terrestrial automatic cameras in a 30 m high glacier monitoring tower at the front of the glacier with the support of the "Key Scientific and Technological Infrastructure Construction Project for the Field Station Network of Chinese Academy of Sciences".

Conclusions
In this study, we used a 1:50,000-scale topographic map based on aerial photography acquired in September 1964, and long-range TLS data acquired in September 2018 to calculate area, volume, and mass changes for Haxilegen Glacier No. 51. We found that the glacier area reduced by 0.34 km 2 at a rate of 0.006 km 2 a −1 , or 0.42% a −1 , from 1964 to 2018. The relative shrinkage rate is slightly higher than global and regional average rates, which may be attributed to the relatively small area of the glacier. Our resulting estimates of surface elevation changes and geodetic mass balance are clearly negative, with surface lowering observed across almost the entire glacier. The calculated mean geodetic mass balance was −0.36 ± 0.12 m w.e. a −1 , which was in the same range as observed and modeled values for the mass balance published in recent studies for the period 1999-2015. However, our calculated value was slightly more positive than the glaciological mass balance of the Tuyuksu glacier over the past 50 years, and more positive than the glacier mass loss in the eastern Tien Shan over the last two decades.
Glacier mass balance measurements in the eastern Tien Shan are crucially needed. The Riegl VZ ® -6000 TLS is a well-established tool for monitoring the mass balance of individual glaciers. We recommend the increased use of this device to survey the glacier. Snow-firn density should be measured in the accumulation area to improve the accuracy of geodetic mass balance estimates. Further measurements and technologies such as terrestrial automatic cameras and glacier monitoring towers can be integrated.