Spatio-Temporal Changes of Mass Balance in the Ablation Area of the Muz Taw Glacier, Sawir Mountains, from MultiTemporal Terrestrial Geodetic Surveys

: The glaciers in the Sawir Mountains are an important freshwater resource, and glaciers have been experiencing a continuing retreat over the past few decades. However, studies on detailed glacier mass changes are currently sparse. Here, we present the high-precision evolution of annual surface elevation and geodetic mass changes in the ablation area of the Muz Taw Glacier (Sawir Mountains, China) over the latest three consecutive mass-balance years (2017–2020) based on multitemporal terrestrial geodetic surveys. Our results revealed clearly surface lowering and negative geodetic mass changes, and the spatial changing patterns were generally similar for the three periods with the most negative surface lowering (approximately − 5.0 to − 4.0 m a − 1 ) around the glacier terminus. The gradient of altitudinal elevation changes was commonly steep at the low elevations and gentle in the upper-elevation parts, and reduced surface lowering was observed at the glacier terminus. Resulting emergence velocities ranged from 0.11 to 0.86 m a − 1 with pronounced spatial variability, which was mainly controlled by surface slope, ice thickness, and the movement of tributary glaciers. Meanwhile, emergence velocities slightly compensated the surface ablation at the ablation area with a proportion of 14.9%, and dynamic thickening had small contributions to glacier surface evolution. Limited annual precipitation and glacier accumulation may result in these weak contributions. Higher-resolution surveys at the seasonal and monthly scales are required to get insight into the mass balance processes and their mechanism. To our knowledge, glacier thickness makes the most important contribution to ice velocity, which may result in higher velocities at the regions [52]. The movement of tributary glaciers can slightly increase the ﬂow of mass; thus, the ablation stake of E3 was characterized with higher velocity. However, those were only quite qualitative descriptions, and we needed more observations to validate the interpretations.


Introduction
Glacier mass balance is known to be one the most direct and immediate indicators of climate variations [1,2], which is of crucial importance for the assessment of local water resources, sea-level rise, and disasters [3][4][5]. West China and its surroundings host the world's largest mid-latitude glaciers, but the extensive in situ measurements of glaciological method contribute to the small number and heterogeneous spatiotemporal distribution of in situ measurement series; only one continuous glaciological mass balance measurement (Urumqi Glacier No. 1, eastern Tien Shan) is available for more than 30 years, and there are several discontinuous observations in the Tibetan Plateau [6][7][8]. This can result in the spatiotemporal pattern of glacier mass change in West China being strongly heterogeneous [6,9]. Consequently, many studies have focused on glacier mass changes in the Altai Mountains, Tien Shan, and Tibetan Plateau by using the geodetic method, which calculated glacier mass changes based upon the digital elevation model (DEM) differencing (e.g., [9][10][11][12][13]). Some studies suggest that the spatial variability of glacier mass balances in the Tibetan Plateau is related to the westerly and monsoon atmospheric circulation processes [6,14]. The Sawir Mountains are controlled by the prevailing westerlies interacting with the Asian anticyclone and polar air mass in winter [15]. Glacier meltwater in the Sawir Mountains is an important freshwater resource for populations and hydroeconomies in Jeminay and its surrounding counties. Hence, studying the glacier changes in the Sawir Mountains has important scientific significance. However, to our knowledge, the present study mainly focuses on glacier area, length, and thickness in the Sawir Mountains, detailed glacier mass changes had not been reported so far [16,17].
The geodetic method can overcome the time-consuming and labor-intensive drawbacks of glaciological measurements, but the method is limited by the accuracy and spatiotemporal resolution of available DEM [18,19]. Recently, cutting-edge technologies such as unmanned aerial vehicles (UAVs) and three-dimensional (3D) laser scanning have been widely used on geodetic mass balance estimates [20][21][22][23]. Airborne laser scanning (ALS) provides measurements over large areas with easy and efficiency, but large topography undulation and high-altitude rock outcrops around glaciers reduce the ability of aerial survey, since most ALS instruments have limited operating flight altitude, and the economic cost of aircraft is also very high [21]. Terrestrial laser scanning (TLS) generally covers limited areas, but it can provide information on surface evolution and mass changes of individual glaciers with high spatiotemporal resolution [23]. Note that glacier emergence/submergence velocities play a significant role in glacier surface evolution [24]. The emergence velocities are upward and result in dynamic thickening in the ablation area; submergence velocities are downward and result in dynamic thinning in the accumulation area [25]. The emergence/submergence velocities correspond to the net mass balance for any given point of the glacier exactly when a glacier is in steady state (i.e., no geodetic glacier volume change) [24]. For an imbalanced glacier, surface elevation changes are the arithmetic sum of the net balance and the emergence/submergence velocities [24,25]. Thus, the geodetic method needs to consider the effects of glacier dynamics to quantify glacier surface mass balance [23]. Many studies have analyzed emergence/submergence velocities for glaciers in the Alps [23,26,27], Svalbard [28], and the Himalayas [20,29,30]. However, no studies have focused on the estimates of dynamic thinning and thickening in the Sawir Mountains by combining in situ glaciological measurements and high-resolution glacier surface elevation changes.
The novel Riegl VZ ® -6000 TLS offers a long measurement range of 6 km and is exceptionally well suited for glacier mapping, since the device uses class 3B laser beams, theoretically enabling high rates of reflection from snow-and ice-covered terrains (>80%) [31]. Thus, in the aforementioned context, the aim of this paper is to study the spatiotemporal evolution of mass balance in the ablation area of the Muz Taw Glacier (Sawir Mountains, China) over the latest three mass-balance years (2017-2020) by using the repeated Riegl VZ ® -6000 TLS surveys. The main contents include (1) analyzing and discussing the spatiotemporal pattern of glacier mass changes and (2) quantifying the influences of surface ablation and emergence velocity to glacier surface lowering.

Study Area
The Sawir Mountains stretch across China and Kazakhstan and are the transitional section between the Tien Shan and the central Altai Mountains (Figure 1a). The Mountains are not only a very obvious watershed between the inland Xinjiang and the Arctic Ocean but also the highest mountain range in the northernmost part of the western mountains of the Junggar Basin. The mountains run east to west with a steep south slope and a gentle north slope. The highest peak is Mustau Mountain with an altitude of 3875 m a.s.l. The Chinese Sawir Mountains host 21 glaciers with a total area of 16.84 km 2 according to the first glacier inventory of China (GIC) [32]. The glaciers in the Sawir Mountains are characterized by high latitude and low altitude, which are significantly different from glaciers in Tien Shan and Tibetan Plateau. Muz Taw Glacier (47°04′ N, 85°34′ E) is a northeast-faced valley glacier that is situated on the northern side of the Sawir Mountains ( Figure 1), with an area of 3.13 km 2 and length of 3.2 km in 2016 [33]. Glaciological investigations of Muz Taw Glacier were initi ated in 2013, and subsequently, an observation and research station (Altai Observatio and Research Station of Cryospheric Science and Sustainable Development, Chines Academy of Sciences) was set up for the long-term measurements. According to the lates observations, the glacier area decreased continuously with a value of 10.51 km 2 (45.72% over the period 1977-2017, and the mean annual length change (front variation) reached up to 11.5 m a −1 from 1989 to 2017 [17]. The glacier is controlled by the prevailing wester lies in a continental climate setting, Jimunai Meteorological Station (984 m a.s.l., located about 46 km northeast of the glacier) showed that the average annual air temperature and precipitation were 4.3 °C and 212 mm, respectively, during the period 1961-2016. Bot annual air temperature and precipitation presented an increasing trend with rates of 0. °C (10 a) −1 and 12 mm (10 a) −1 , respectively. High air temperature mainly occurred from May to August, and air temperature was very low during winter. The annual distributio of precipitation was approximately homogeneous and precipitation in July was relativel high, and the winter precipitation accounted for 10-30% of the annual total ( Figure 2). Muz Taw Glacier (47 • 04 N, 85 • 34 E) is a northeast-faced valley glacier that is situated on the northern side of the Sawir Mountains ( Figure 1), with an area of 3.13 km 2 and a length of 3.2 km in 2016 [33]. Glaciological investigations of Muz Taw Glacier were initiated in 2013, and subsequently, an observation and research station (Altai Observation and Research Station of Cryospheric Science and Sustainable Development, Chinese Academy of Sciences) was set up for the long-term measurements. According to the latest observations, the glacier area decreased continuously with a value of 10.51 km 2 (45.72%) over the period 1977-2017, and the mean annual length change (front variation) reached up to 11.5 m a −1 from 1989 to 2017 [17]. The glacier is controlled by the prevailing westerlies in a continental climate setting, Jimunai Meteorological Station (984 m a.s.l., located about 46 km northeast of the glacier) showed that the average annual air temperature and precipitation were 4.3 • C and 212 mm, respectively, during the period 1961-2016. Both annual air temperature and precipitation presented an increasing trend with rates of 0.4 • C (10 a) −1 and 12 mm (10 a) −1 , respectively. High air temperature mainly occurred from May to August, and air temperature was very low during winter. The annual distribution of precipitation was approximately homogeneous and precipitation in July was relatively high, and the winter precipitation accounted for 10-30% of the annual total ( Figure 2).

TLS Surveys and Data Processing
As mentioned above, we utilized the Riegl VZ ® -6000 TLS to acquire the glacier surface terrain. The TLS was produced by Riegl Laser Measurement Systems GmbH (Austria) in 2012 with a weight of 14.5 kg, and it was dust-and water-proof. The continuous operating voltage and air temperature range from 11 to 32 V and from −10 to 50 °C, respectively. The device can offer a maximum measurement range of 6 km for natural targets and highspeed data acquisition of 222,000 points per second. The angle measurement resolution is better than 0.0005° (1.8 arcsec), and laser divergence is 0.12 mrad (equivalent to an increase of 12 mm of beam diameter per 100 m distance). The scanning range ranges from 0 to 360° horizontally and from 60-120 °vertically (from zenith). The Riegl VZ ® -6000 TLS is a wellestablished tool for measuring snowy and icy terrain even at long distances, since the device operates at the wavelengths in the near-infrared band (class 3B laser beams, 1064 nm), which enables high rates of reflection (>80%) [31]. Note that the volume-to-mass conversion factor is a fundamental issue for geodetic mass balance calculations (see Section 3.2); the TLS surveys are ideally carried out at the end of the ablation season to reduce the influence of fresh snow. Here, point clouds of the ablation area were acquired by the TLS surveys implemented on 13 September 2017, 4 September 2018, 13 September 2019, and 15 September 2020, and the ablation area was mainly covered with bare ice on the days of the TLS surveys ( Figure 3). Specific measurements included three procedures.
(1) Placement of scan positions. According to [34] and our field experience, the placement of scan positions adhere to the following basic criteria: (a) scan positions must be placed in the stable terrain with broad perspectives and easily accessible; (b) scan positions should be distributed at different elevations around mountain glaciers, and the number of scan positions should be as small as possible to reduce error related to point cloud processing; (c) the distance between the scanner and the target should be within the scanning range, and the angle between the laser pulse and the target should be more than 10° to receive the echo of the laser pulse; (d) the overlap percentage of two neighboring scan positions should be more than 30% to guarantee the quality of data registration. Based on the aforementioned criteria, two scan positions were used for the TLS surveys, and their spatial distribution was presented in Figure 1c. We fixed each scan position on the stable bedrock by using a reinforced concrete structure; then, we combined and installed stainless steel measuring nails and protection boxes to build a standard GNSS-leveling (GNSS, Global Navigation Satellite System) point. Despite the visual angles of the two scan positions not being very ideal relative to the entire glacier surface, the TLS surveys completely covered the glacier tongue. In addition, we did not place scan positions at the right side

TLS Surveys and Data Processing
As mentioned above, we utilized the Riegl VZ ® -6000 TLS to acquire the glacier surface terrain. The TLS was produced by Riegl Laser Measurement Systems GmbH (Austria) in 2012 with a weight of 14.5 kg, and it was dust-and water-proof. The continuous operating voltage and air temperature range from 11 to 32 V and from −10 to 50 • C, respectively. The device can offer a maximum measurement range of 6 km for natural targets and high-speed data acquisition of 222,000 points per second. The angle measurement resolution is better than 0.0005 • (1.8 arcsec), and laser divergence is 0.12 mrad (equivalent to an increase of 12 mm of beam diameter per 100 m distance). The scanning range ranges from 0 to 360 • horizontally and from 60-120 • vertically (from zenith). The Riegl VZ ® -6000 TLS is a well-established tool for measuring snowy and icy terrain even at long distances, since the device operates at the wavelengths in the near-infrared band (class 3B laser beams, 1064 nm), which enables high rates of reflection (>80%) [31]. Note that the volume-tomass conversion factor is a fundamental issue for geodetic mass balance calculations (see Section 3.2); the TLS surveys are ideally carried out at the end of the ablation season to reduce the influence of fresh snow. Here, point clouds of the ablation area were acquired by the TLS surveys implemented on 13 September 2017, 4 September 2018, 13 September 2019, and 15 September 2020, and the ablation area was mainly covered with bare ice on the days of the TLS surveys ( Figure 3). Specific measurements included three procedures. of Figure 1c, since the region was covered with sand and gravel, and the terrain of the right side was unstable and steep, which will be detrimental to TLS surveys.
(2) Three-dimensional (3D) coordinate of the two scan positions. We used the realtime kinematic (RTK) GNSS of Trimble R10 to acquire coordinate, which was used for point cloud direct georeferencing (see data processing). Note that it is better to place some retro-reflective targets surveyed with the GNSS at much greater distances relative to scan positions, since the measuring accuracy will reduce with the increase of scanning range. However, this study mainly focused on the relative changes in glacier surface elevation between two consecutive scan campaigns, and the requirement of absolute accuracy was not very strong. Thus, we had not implemented this type of survey.
(3) Point cloud acquisition. We mounted the TLS on a sturdy tripod for each scan position after placement and 3D coordinate of scan positions. Considering that the environmental conditions around glaciers vary significantly, here, the input environmental parameters (air temperature, atmospheric pressure, relative humidity, and height above sea level) of the TLS surveys were provided by portable meteorological instruments and Trimble R10 GNSS. All of the four scan campaigns were carried out under an overcast sky to reduce the influence of atmospheric refraction ( Figure 3) [31]. Note that the altitude drop and scanning range were relatively small, since this study only focused on the ablation area of the Muz Taw Glacier, which may result in unremarkable variations of atmospheric refraction. In addition, dry, fogless, and windless weather conditions were also considered to enhance measurement accuracy (see Section 3.5). After these steps, we configured the laser pulse repetition rate as 50 kHz to complete coarse scanning with vertical and horizontal angles range of 60-120° from zenith and 0~360°, respectively; then, the glacial regions of the coarse scanning images were selected to finish fine scanning with the laser pulse repetition rate of 30 kHz. Detailed surveying parameters are listed in Table  1.   (1) Placement of scan positions. According to [34] and our field experience, the placement of scan positions adhere to the following basic criteria: (a) scan positions must be placed in the stable terrain with broad perspectives and easily accessible; (b) scan positions should be distributed at different elevations around mountain glaciers, and the number of scan positions should be as small as possible to reduce error related to point cloud processing; (c) the distance between the scanner and the target should be within the scanning range, and the angle between the laser pulse and the target should be more than 10 • to receive the echo of the laser pulse; (d) the overlap percentage of two neighboring scan positions should be more than 30% to guarantee the quality of data registration. Based on the aforementioned criteria, two scan positions were used for the TLS surveys, and their spatial distribution was presented in Figure 1c. We fixed each scan position on the stable bedrock by using a reinforced concrete structure; then, we combined and installed stainless steel measuring nails and protection boxes to build a standard GNSS-leveling (GNSS, Global Navigation Satellite System) point. Despite the visual angles of the two scan positions not being very ideal relative to the entire glacier surface, the TLS surveys completely covered the glacier tongue. In addition, we did not place scan positions at the right side of Figure 1c, since the region was covered with sand and gravel, and the terrain of the right side was unstable and steep, which will be detrimental to TLS surveys.
(2) Three-dimensional (3D) coordinate of the two scan positions. We used the real-time kinematic (RTK) GNSS of Trimble R10 to acquire coordinate, which was used for point cloud direct georeferencing (see data processing). Note that it is better to place some retro-reflective targets surveyed with the GNSS at much greater distances relative to scan positions, since the measuring accuracy will reduce with the increase of scanning range. However, this study mainly focused on the relative changes in glacier surface elevation between two consecutive scan campaigns, and the requirement of absolute accuracy was not very strong. Thus, we had not implemented this type of survey.
(3) Point cloud acquisition. We mounted the TLS on a sturdy tripod for each scan position after placement and 3D coordinate of scan positions. Considering that the environmental conditions around glaciers vary significantly, here, the input environmental parameters (air temperature, atmospheric pressure, relative humidity, and height above sea level) of the TLS surveys were provided by portable meteorological instruments and Trimble R10 GNSS. All of the four scan campaigns were carried out under an overcast sky to reduce the influence of atmospheric refraction ( Figure 3) [31]. Note that the altitude drop and scanning range were relatively small, since this study only focused on the ablation area of the Muz Taw Glacier, which may result in unremarkable variations of atmospheric refraction. In addition, dry, fogless, and windless weather conditions were also considered to enhance measurement accuracy (see Section 3.5). After these steps, we configured the laser pulse repetition rate as 50 kHz to complete coarse scanning with vertical and horizontal angles range of 60-120 • from zenith and 0~360 • , respectively; then, the glacial regions of the coarse scanning images were selected to finish fine scanning with the laser pulse repetition rate of 30 kHz. Detailed surveying parameters are listed in Table 1. We processed point cloud by using RiSCAN PRO ® v 1.81 software and the overall workflow used for data processing was given in Figure 4, which included direct georeferencing, data registration, compression, filtering, multi-temporal registration, and interpolation. The specific processes were as follows: (1) direct georeferencing, we used the 3D coordinates of the two scan positions to transform the point cloud data from the scanner's own coordinate system into a global coordinate system [35]; (2) data registration, we used multistation adjustment (MSA) to finish data registration of each scan position according to the iterative closest point (ICP) algorithms [36] and subsequently merged the two scans in one layer; (3) compression, an octree algorithm was used to compress the merged point cloud with a equal space of 0.5 m; (4) filtering, the terrain filter was used to filter out the noise and non-ground data; additionally, we used visual interpretation to check the data individually and deleted outliers [37]; and (5) multi-temporal registration, the attitude angles of each scan campaign are different as the orientation of each scan position is continually adjusted to complete data registration, which can lead to horizontal and vertical shifts between two consecutive scan campaigns and increase the error of TLS-derived DEM differencing. Thus, multi-temporal registration directly determined the accuracy of glacier surface elevation change [38,39]. Since the coarse scanning and fine scanning had the same attitude angles, for each scan campaign, we first removed the glacierized regions of the point clouds from coarse scanning, and then we set the processed point clouds of 2017 as the reference layer and used ICP algorithms to implement the alignment of scan campaigns in 2018, 2019, and 2020 based on the point clouds over stable non-glacier terrain. Hence, the point clouds of fine scanning between two consecutive scan campaigns can coincide well after these processes. Finally, interpolation of the processed point clouds from fine scanning created high-resolution DEMs (0.5 m) after multi-temporal registration. to the iterative closest point (ICP) algorithms [36] and subsequently merged the two scans in one layer; (3) compression, an octree algorithm was used to compress the merged point cloud with a equal space of 0.5 m; (4) filtering, the terrain filter was used to filter out the noise and non-ground data; additionally, we used visual interpretation to check the data individually and deleted outliers [37]; and (5) multi-temporal registration, the attitude angles of each scan campaign are different as the orientation of each scan position is continually adjusted to complete data registration, which can lead to horizontal and vertical shifts between two consecutive scan campaigns and increase the error of TLS-derived DEM differencing. Thus, multi-temporal registration directly determined the accuracy of glacier surface elevation change [38,39]. Since the coarse scanning and fine scanning had the same attitude angles, for each scan campaign, we first removed the glacierized regions of the point clouds from coarse scanning, and then we set the processed point clouds of 2017 as the reference layer and used ICP algorithms to implement the alignment of scan campaigns in 2018, 2019, and 2020 based on the point clouds over stable non-glacier terrain. Hence, the point clouds of fine scanning between two consecutive scan campaigns can coincide well after these processes. Finally, interpolation of the processed point clouds from fine scanning created high-resolution DEMs (0.5 m) after multi-temporal registration. . Schematic illustration of the workflow used for point cloud processing. "GNSS" is the abbreviation of "Global Navigation Satellite System" and "RTK" is the abbreviation of "real-time kinematic".

Geodetic Mass Balance
Geodetic mass balance is calculated as a multiplication between glacier volume change and volume-to-mass conversion factor [5]. Volume change ΔV over the TLS survey period between t0 and t1 can be determined as: where r is the pixel size of TLS-derived DEM (0.5 m), △ ℎ is the glacier surface elevation change of the two consecutive DEMs at pixel k, and K denotes the number of total pixels covering the extent of the TLS surveyed area. Glacier volume change ΔV is converted to geodetic mass balance (m water equivalent (w.e.)) using: where is the mean glacier area of the two survey times, (1000 kg m −3 ) is the density of water, and is the volume-to-mass conversion factor for converting glacier volume change to mass balance. Here, we conservatively assumed that equaled the density of ice (900 kg m −3 ), since the monitoring area was mainly covered with bare ice on the days of the TLS surveys ( Figure 3). . Schematic illustration of the workflow used for point cloud processing. "GNSS" is the abbreviation of "Global Navigation Satellite System" and "RTK" is the abbreviation of "real-time kinematic".

Geodetic Mass Balance
Geodetic mass balance is calculated as a multiplication between glacier volume change and volume-to-mass conversion factor [5]. Volume change ∆V over the TLS survey period between t 0 and t 1 can be determined as: where r is the pixel size of TLS-derived DEM (0.5 m), h k is the glacier surface elevation change of the two consecutive DEMs at pixel k, and K denotes the number of total pixels covering the extent of the TLS surveyed area. Glacier volume change ∆V is converted to geodetic mass balance (m water equivalent (w.e.)) using: where S is the mean glacier area of the two survey times, ρ water (1000 kg m −3 ) is the density of water, and ρ is the volume-to-mass conversion factor for converting glacier volume change to mass balance. Here, we conservatively assumed that ρ equaled the density of ice (900 kg m −3 ), since the monitoring area was mainly covered with bare ice on the days of the TLS surveys ( Figure 3).

Glaciological Mass Balance
The glaciological mass balance of the Muz Taw Glacier had been implemented by measuring ablation stakes and snow pits; more than 20 stakes were drilled into the glacier using a steam drill and distributed in eight rows at different altitudes with roughly three stakes in every row (Figure 1b). The measurements include the stake vertical height above the glacier surface, superimposed ice thickness, and the density and thickness of each snow/firn layer at snow pits. The mass balance of each measured point can be calculated as the arithmetic sum of snow, glacier ice, and superimposed ice mass balance. Point values are extrapolated to glacier-wide mass balance using the contour-line or profile method [40,41].

Emergence Velocity
The geodetic surveys of glacier elevation change mainly include surface mass balance and vertical velocity [25]. The glaciological method measures the surface mass balance [42]. Therefore, conservation of mass at a point (ablation stake) on the surface of a glacier can be expressed as: where .
h is glacier surface elevation changes, .
b is the glaciological mass balance, ρ is density, and ∇· → Q is ice flux. When Equation (3) is integrated over the entire glacier surface, ∇· → Q becomes zero, and the glaciological mass balance equals the elevation change multiplied by the volume-to-mass conversion factor (i.e., geodetic mass balance). However, ∇· → Q is not zero for a given point, and the vertical velocity (w s ) plays a confounding role [43].
As shown in Figure 5, vertical velocity refers to the upward or downward flow of ice relative to the glacier surface at a fixed coordinate, and it is called "emergence velocity" in the ablation area and "submergence velocity" in the accumulation area [25]. Emergence and submergence velocity result in dynamic thickening and thinning of a glacier, respectively [44]. Previous studies assume that emergence velocity w s depends on the kinematic boundary condition at the glacier surface [25,43] and can be calculated as: where u s and v s are the horizontal components of ice velocity at the glacier surface S [25,43]. We neglect the advection of the glacier surface topography induced by horizontal ice flux since short time intervals of the TLS surveys [45], which yields: and for emergence velocity:

Uncertainty Estimates
The uncertainty in TLS-derived glacier surface elevation and mass changes depen on point cloud data acquisition and processing and volume-to-mass conversion [21]. Po cloud data acquisition corresponds to the stability of the two scan position and atm pheric conditions of the TLS surveys, and therefore, we fixed scan positions on stable ro surfaces to minimize the uncertainty. Additionally, point clouds of the glacier surf were acquired in dry, fogless, and windless weather conditions to reduce wind-and mo ture-induced errors, since the wind will result in the vibration of the TLS instrument a moisture can absorb laser pulse. In addition, the systematic error of the device can a induce related uncertainty in data acquisition, but it varies greatly and cannot be quan fied well, for instance, with the scanning distance to the target. So, we only discussed i Section 5.1. For point cloud data acquisition processing error, the standard deviation errors (σMSA) obtained from registering the point cloud can be used to evaluate multi-te poral registration quality [46], the resulting errors were in the range of 0.07-0.10 m a generally less than the error of TLS-derived glacier surface elevation ( Table 2). We o discussed the results in Section 5.1. The octree algorithm built the topological relationsh of scattered points, which can preserve terrain information as much as possible. The er related to DEM (with a resolution of 0.5 m) creation can be neglected because of the h average point density (1.90-5.21 points m −2 , Table 1). For volume-to-mass conversion, uncertainty of glacier ice density was assumed to be ± 17 kg m −3 according to [47].
The standard deviation of the stable non-glacier terrain ( ∆ ) was used to calcul the uncertainty of the TLS-derived glacier surface elevation ( ∆ ) due to the lack of p cise and well-distributed stable ground control points, but the spatial auto-correlat must be removed when we calculate the uncertainty of the glacier-wide elevation chang Here, ∆ was quantified using the geostatistical analysis of [48] and expressed as: where is the correlation area. We assume that = because of the high po density (>1 point m −2 ) [49]. Uncertainty in geodetic mass balance ( ) can be estima using: Figure 5. Conceptual representation of the different velocity components, which was redrawn base on [25].
As described in Section 3.3, here, we used the point mass balance from glaciological measurements in 2018-2019 to estimate the emergence velocity for each individual ablation stake.

Uncertainty Estimates
The uncertainty in TLS-derived glacier surface elevation and mass changes depends on point cloud data acquisition and processing and volume-to-mass conversion [21]. Point cloud data acquisition corresponds to the stability of the two scan position and atmospheric conditions of the TLS surveys, and therefore, we fixed scan positions on stable rock surfaces to minimize the uncertainty. Additionally, point clouds of the glacier surface were acquired in dry, fogless, and windless weather conditions to reduce wind-and moisture-induced errors, since the wind will result in the vibration of the TLS instrument and moisture can absorb laser pulse. In addition, the systematic error of the device can also induce related uncertainty in data acquisition, but it varies greatly and cannot be quantified well, for instance, with the scanning distance to the target. So, we only discussed it in Section 5.1. For point cloud data acquisition processing error, the standard deviation of errors (σ MSA ) obtained from registering the point cloud can be used to evaluate multitemporal registration quality [46], the resulting errors were in the range of 0.07-0.10 m and generally less than the error of TLS-derived glacier surface elevation ( Table 2). We only discussed the results in Section 5.1. The octree algorithm built the topological relationship of scattered points, which can preserve terrain information as much as possible. The error related to DEM (with a resolution of 0.5 m) creation can be neglected because of the high average point density (1.90-5.21 points m −2 , Table 1). For volume-to-mass conversion, the uncertainty of glacier ice density was assumed to be ± 17 kg m −3 according to [47]. Table 2. Error (σ MSA ) of multi-station adjustment (MSA) and the number of points (n) used for multi-temporal registration of two consecutive TLS surveys, the uncertainty (σ ∆h TLS ) of TLS-derived glacier surface elevation, and the corresponding uncertainty of geodetic (σ geod ) and glaciological (σ glac.p ) mass balance. The standard deviation of the stable non-glacier terrain (σ ∆h TLS ) was used to calculate the uncertainty of the TLS-derived glacier surface elevation (σ ∆h TLS ) due to the lack of precise and well-distributed stable ground control points, but the spatial auto-correlation must be removed when we calculate the uncertainty of the glacier-wide elevation changes. Here, σ ∆h TLS was quantified using the geostatistical analysis of [48] and expressed as:

Mass Balance Year
where S cor is the correlation area. We assume that S = S cor because of the high point density (>1 point m −2 ) [49]. Uncertainty in geodetic mass balance (σ geod ) can be estimated using: where ∆h TLS is the mean value of the TLS-derived glacier elevation change. The calculated uncertainties are shown in Figure 6 and listed in Table 2.
Remote Sens. 2021, 13, x 9 of 19 Figure 6. Spatial distribution of surface elevation changes over stable non-glacier terrain (a-c) and corresponding frequency distributions of these changes (d-f). The median (μ), the standard deviation (σ), fitting coefficient (R 2 ), and the number of pixels (N) are given.
Uncertainty in glaciological mass balance measurements at individual point locations originate from errors in stake readings, snow/firn probing, and the sinking of the stakes [42,50]. Following [51], the uncertainty of the glaciological method raised with the total number of sampling sites, corresponding error of ablation measured in ice, and firn was estimated to be about 0.14 w.e. a −1 and 0.27 w.e. a −1 , respectively, and the error in accumulation measurements was 0.14 w.e. a −1 , based on independent geodetic mass balance. Uncertainty in a given point can be calculated based on the law of error propagation [42]. Thus, errors of ablation measurements in ice ( ) and firn ( ) are evaluated using 0.14/ and 0.27/ , respectively, where and are the total number of ablation stakes and snow pits, respectively; errors in accumulation measurements ( ) are Figure 6. Spatial distribution of surface elevation changes over stable non-glacier terrain (a-c) and corresponding frequency distributions of these changes (d-f). The median (µ), the standard deviation (σ), fitting coefficient (R 2 ), and the number of pixels (N) are given.
Uncertainty in glaciological mass balance measurements at individual point locations originate from errors in stake readings, snow/firn probing, and the sinking of the stakes [42,50]. Following [51], the uncertainty of the glaciological method raised with the total number of sampling sites, corresponding error of ablation measured in ice, and firn was estimated to be about 0.14 w.e. a −1 and 0.27 w.e. a −1 , respectively, and the error in accumulation measurements was 0.14 w.e. a −1 , based on independent geodetic mass balance.
Uncertainty in a given point can be calculated based on the law of error propagation [42]. Thus, errors of ablation measurements in ice (σ ice ) and firn (σ firn ) are evaluated using 0.14/ √ N ice and 0.27/ √ N firn , respectively, where N ice and N firn are the total number of ablation stakes and snow pits, respectively; errors in accumulation measurements (σ acc ) are determined based on 0.21/ √ N c , where N c is the number of snow pits. The total uncertainty in point annual mass balance (σ glac.p ) is calculated as: Here, we mainly considered the measured errors in the ablation area and the calculated error σ glac.p was about 0.03 m w.e. a −1 . Due to surface elevation change in the ablation area is the arithmetic sum of mass balance and the emergence velocities (σ vel ), the uncertainty in emergence velocity can be calculated applying a reciprocal density conversion: The resulting values are listed in Table 2 and discussed in Section 5.1.

Glacier Surface Elevation and Geodetic Mass Changes
As shown in Figure 7, TLS-derived DEM differencing gives detailed spatial distributions of glacier surface elevation changes for the latest three mass-balance years. The ablation area of the Muz Taw Glacier showed clearly surface lowering and negative geodetic mass changes, with a value of −1.33 ± 0.13 m w.e a −1 , −1.64 ± 0.09 m w.e a −1 , and −1.01 ± 0.11 m w.e a −1 , for the mass-balance years 2017-2018, 2018-2019, and 2019-2020, respectively. Spatial changing patterns were maintained as generally similar for all of the three investigated periods with the most negative surface lowering values ranged from −5.0 to −4.0 m a −1 near the glacier terminus and less pronounced thinning in the upper elevations. Measured glacier thinning was remarkably higher for the second period with 92.2% of the elevation change falling in the range −3.0 to −1.0 m, which was significantly more negative compared to the first and third mass-balance year (corresponding proportion was 79.6% and 60.2%, respectively). These temporal differences in observed glacier surface elevation and mass changes were consistent with summer (June-September) air temperature (11.4, 12.2, and 11.6 • C in summer 2018, 2019, and 2020, respectively) recorded by Jeminay (website: https://www.tianqi.com/, accessed on 1 March 2021), which showed a similar response to the observed climatic forcing compared to glaciers in other regions, demonstrating that glaciers in the Sawir Mountains were also sensitive to the warming climate.
The altitudinal distribution of glacier elevation changes is a combination of glacier surface mass balance and dynamic processes (emergence/submergence velocities). The altitudinal profiles of glacier surface elevation changes were depicted in Figure 8 and showed thinning continuously decreased with increasing of altitude, which was a commonly changing behavior and agreed well with the knowledge of existing in situ measurements [8]. In addition, reduced surface lowering occurred at the glacier terminus, which was consistent with the glaciological mass balance versus elevation of most reference glaciers reported by [8]. and 60.2%, respectively). These temporal differences in observed glacier surface elevati and mass changes were consistent with summer (June-September) air temperature (11 12.2, and 11.6 °C in summer 2018, 2019, and 2020, respectively) recorded by Jeminay (we site: https://www.tianqi.com/, accessed on 1 March 2021), which showed a similar sponse to the observed climatic forcing compared to glaciers in other regions, demonstr ing that glaciers in the Sawir Mountains were also sensitive to the warming climate. The altitudinal distribution of glacier elevation changes is a combination of glac surface mass balance and dynamic processes (emergence/submergence velocities). T altitudinal profiles of glacier surface elevation changes were depicted in Figure 8 and showed thinning continuously decreased with increasing of altitude, which was a commonly changing behavior and agreed well with the knowledge of existing in situ measurements [8]. In addition, reduced surface lowering occurred at the glacier terminus, which was consistent with the glaciological mass balance versus elevation of most reference glaciers reported by [8]. In order to identify aspect controls on the patterns of glacier surface elevation changes, we calculated the distribution by the aspect of the glacier area (in 2017) compared with the distribution of mean glacier surface elevation change. Results suggested that the mean glacier surface lowering on the eight aspects was in good agreement with the glacier mass loss of the three mass-balance years. The majority of the observed glacier area had aspects north through east, while more negative mean glacier surface changes were dominant for west, southwest, and south aspects with the exception of surface changes on the west aspect for the mass-balance year 2017-2018, confirming that aspect controlled glacier ablation, since it was possible to influence the glacier surface albedo and the TLS had potential for detailed investigation of the controls.

Surface Ablation and Emergence Velocity
In situ measurements of single-point mass balance were available by using the ablation stakes. Here, we selected the well-measured data in the mass balance year 2018-2019  In order to identify aspect controls on the patterns of glacier surface elevation changes, we calculated the distribution by the aspect of the glacier area (in 2017) compared with the distribution of mean glacier surface elevation change. Results suggested that the mean glacier surface lowering on the eight aspects was in good agreement with the glacier mass loss of the three mass-balance years. The majority of the observed glacier area had aspects north through east, while more negative mean glacier surface changes were dominant for west, southwest, and south aspects with the exception of surface changes on the west aspect for the mass-balance year 2017-2018, confirming that aspect controlled glacier ablation, since it was possible to influence the glacier surface albedo and the TLS had potential for detailed investigation of the controls.

Surface Ablation and Emergence Velocity
In situ measurements of single-point mass balance were available by using the ablation stakes. Here, we selected the well-measured data in the mass balance year 2018-2019 to reflect glacier surface ablation and emergence velocity. Glaciological measurements revealed that significant surface ablation occurred at the row of C (−3.06 m w.e. a −1 at the ablation stake of C3), and the ablation at the row of "A" was relatively less negative, this changing pattern was consistent with TLS-derived glacier surface elevation changes (Figures 7 and 8). Stakes in the row of "D-E" exhibited surface ablation between −2.50 w.e. a −1 and −1.50 w.e. a −1 , less negative mass balance occurred at the uppermost stakes (H1 and H2), which was near with the TLS-derived glacier mass balance at the same elevation bands.
The calculated emergence velocities ranged from 0.11 to 0.86 m a −1 , which slightly compensated for the surface ablation at the ablation area. The higher velocities happened at the ablation stakes of C1 and D1. The spatially distributed emergence velocity showed pronounced spatial variability, which decreased from the lower elevations to the upper glacierized regions (Figure 9). Spatial patterns in emergence velocities may be related to glacier surface slope; big slope or steep terrain was conducive to the downward migration of glacier mass and resulted in large horizontal velocities and small emergence velocities. Our comparisons between the emergence velocity and glacier surface slope at the ablation stakes exhibited an excellent correlation ( Figure 10). Note that the spatial distribution of measured sites was insufficient, since a very small number of well-measured ablation stakes can be used, which may affect the correlation coefficient. So, the numerical results need more sites to validate it. In addition, the left side around the ablation stakes of C1 and D1 had thicker glacier ice with thickness in the range of 70-90 m according to ground penetrating radar measurements in 2013 [16]. To our knowledge, glacier thickness makes the most important contribution to ice velocity, which may result in higher velocities at the left regions [52]. The movement of tributary glaciers at the right side can slightly increase the flow of mass; thus, the ablation stake of E3 was characterized with higher velocity. However, those were only quite qualitative descriptions, and we needed more observations to validate the interpretations. of glacier mass and resulted in large horizontal velocities and small emergence velocities Our comparisons between the emergence velocity and glacier surface slope at the ablation stakes exhibited an excellent correlation ( Figure 10). Note that the spatial distribution of measured sites was insufficient, since a very small number of well-measured ablation stakes can be used, which may affect the correlation coefficient. So, the numerical results need more sites to validate it. In addition, the left side around the ablation stakes of C1 and D1 had thicker glacier ice with thickness in the range of 70-90 m according to ground penetrating radar measurements in 2013 [16]. To our knowledge, glacier thickness makes the most important contribution to ice velocity, which may result in higher velocities at the left regions [52]. The movement of tributary glaciers at the right side can slightly increase the flow of mass; thus, the ablation stake of E3 was characterized with higher velocity. However, those were only quite qualitative descriptions, and we needed more observations to validate the interpretations.

Accuracy of DEM Differencing and Geodetic Mass Balance
The fixed scan positions used in this study can enhance the veys, and we selected windless atmospheric conditions to reduce

Accuracy of DEM Differencing and Geodetic Mass Balance
The fixed scan positions used in this study can enhance the stability of the TLS surveys, and we selected windless atmospheric conditions to reduce the vibration of the TLS device while scanning. Nevertheless, the TLS itself and glacier surface terrain can introduce some errors into the acquired point clouds. The systematic error of the TLS cannot quantify well, so the instrument manufacturers usually give simplified and constant values for the ranging accuracy and precision. For instance, the precision of Riegl VZ ® -6000 TLS reported has been reported to be about 10 mm [31], which completely satisfies the measurement accuracy in glaciology. Terrain-induced errors refer to glacier surface characteristics or geometry and the incident angle between the target and the scanner (e.g., slope, aspect, and size of the targets) [53,54]. The applied TLS can achieve a scanning range of 6 km provided that the target size exceeds the laser beam footprint, perpendicular incident angle, and good atmospheric visibility [31]. For the purpose of this study, the established scan positions were situated at different elevations and directions relative to the glacier surface. Dry weather conditions and overcast sky were also selected to reduce the influence of atmospheric refraction [31]. Note that the laser beam footprint had reached up to the top elevations of the Muz Taw Glacier under these premises, but some voids existed at the accumulation area. This study chose the best-surveyed data at the ablation zone to guarantee the good quality of the TLS-derived DEMs.
Surface elevation changes over stable non-glacier terrain are used to evaluate the accuracy of DEM differencing, which exhibited normal distribution ( Figure 6). Our statistics suggested that more than 80% and 90% of the total pixels with elevation changes fell in the confidence intervals of 2 × sigma and 3 × sigma, respectively, indicating that the accuracy of DEM differencing should be limited to the decimeter level. This may be related to fixed scan positions and the use of multi-temporal registration [38]. According to the point clouds over stable terrain, here, multi-temporal registration was used indispensably to correct the misalignment of two consecutive scan campaigns in the horizontal and vertical orientation via ICP algorithms due to the lack of in situ-measured ground control points [21,39].
On average, the uncertainty of TLS-derived geodetic mass balance was 0.11 m w.e. a −1 for the three investigated periods ( Table 2). The accuracy of our results was close to the studies that used high-precision 3D laser scanning, such as Hintereisferner in Austria, with a mean error of 0.07 m w.e. a −1 [47], four very small glaciers in Switzerland, with a mean error of 0.13 m w.e. a −1 [21], and Echaurren Norte Glacier in Chilean Andes, with a mean error of 0.09 m w.e. a −1 [55]. However, the uncertainty estimated in our study was smaller than the similar studies that used coarse DEMs (e.g., [56][57][58]), confirming the high precision of the TLS surveys. Note that this study only focused on glacier mass balance in the ablation area, and uncertainties in estimated accumulation had not been considered, which probably improved the accuracy of geodetic mass balance [59].

Spatiotemporal Variability of Glacier Surface Elevation and Mass Changes
It is widely recognized that glacier mass balance is very sensitive to air temperature [60]. The aforementioned high air temperature in summer 2019 was one of the key factors for clearly negative surface elevation and mass changes in the second mass-balance year. Other factors, such as the topographical conditions, can also influence the spatial patterns in TLS-derived surface elevation changes. Our results suggested that altitudinal distributions of glacier surface slope and elevation changes exhibited opposite tendencies. The lower elevations presented steep terrain with slopes in the range of 16-18 • , and clearly negative surface elevation changes were observed in those elevations (Figure 11a). A good correlation existed between the mean glacier surface slope and elevation changes over the three periods (Figure 11b). The steep terrain can increase the removal of large snow depositions and is detrimental to glacier accumulation. Note that the mean glacier surface slope at different elevation bands showed an increasing trend with continuous glacier ablation, which can contribute to accelerated glacier mass loss.

Influence of Emergence Velocity and Mass Balance on Surface Lowering
In order to estimate the influence of each component in the ablation zone quantitatively, the annual mass balance measured at the nine stakes in the mass-balance year 2018-2019 was compared to the emergence velocities ( Figure 12). TLS-derived glacier surface elevation changes were dominated by surface ablation (the average value of stake height changes was −2.18 m), and the mean contribution of emergence velocity was only limited to be 0.33 ± 0.11 m, which slightly compensated the surface ablation with a mean percentage of 14.9%. The findings suggested that dynamic thickening had small contributions to glacier surface evolution at the ablation area, and glacier surface lowering was mainly controlled by mass balance; such findings could happen at other glaciers in the Sawir

Influence of Emergence Velocity and Mass Balance on Surface Lowering
In order to estimate the influence of each component in the ablation zone quantitatively, the annual mass balance measured at the nine stakes in the mass-balance year 2018-2019 was compared to the emergence velocities ( Figure 12). TLS-derived glacier surface elevation changes were dominated by surface ablation (the average value of stake height changes was −2.18 m), and the mean contribution of emergence velocity was only limited to be 0.33 ± 0.11 m, which slightly compensated the surface ablation with a mean percentage of 14.9%. The findings suggested that dynamic thickening had small contributions to glacier surface evolution at the ablation area, and glacier surface lowering was mainly controlled by mass balance; such findings could happen at other glaciers in the Sawir Mountains. With the continuous glacier mass loss (concomitant with ice thinning) and the reduction of accumulation area ratio, emergence velocity will gradually decrease, which in turn results in rapid glacier surface lowering [61]. For instance, a decline in emergence (−0.11 m) velocity had been observed at Qiyi Glacier in the Qilian Mountains from the periods 1975-1985 to 1986-2002 [62], and similar phenomena were also observed at Lirung Glacier in the Nepal Himalaya since 2000 [29].

Limitations and Further Work
Although the applied TLS can provide a scanning range of 6 km and the scanned area exceeded 4.0 km 2 for the four years, our study only focused on a portion of the ablation area of (accounting for about 17.6% of the total glacier area) the Muz Taw Glacier, since some voids occurred at the upper elevations. This is mainly because of the limited scanning angle of the TLS. Auxiliary surveys from other techniques such as a cut-price unmanned aerial vehicle (UAV) can be required. Consequently, the combination of TLS and UAV had been done in many studies (e.g., [22,67,68]). A 30 m-high glacier monitoring tower had been placed at the fore of the glacier, which had a good visual angle relative to the entire glacier surface. Therefore, the TLS device can be mounted at the tower to increase the scanning range. For further work, we will integrate the TLS, UAV, and glacier monitoring tower to acquire the entire glacier surface terrain at the seasonal and monthly scales. In such a case, high-resolution and high-precision geodetic glacier mass balance can be yielded, which will bring insight into the mass balance processes and their mechanisms. Moreover, accurate quantification of the submergence/emergence velocities is important for geodetic mass balance estimates and simulation of the glacier dynamic process; therefore, we propose to strengthen the glaciological mass-balance measurements by using stakes and snow pits.

Stakes
Emergence velocity Surface ablation Surface elevation changes Figure 12. Emergence velocity (green), in situ measured surface lowing due to ablation (blue) and TLS-derived glacier surface elevation changes at each ablation stake (red).
Here, we compared our results to previous similar studies to discuss the dynamic characteristics of different regions and types of glaciers and thus to investigate the possible mechanism of dynamic thickening. The estimated annual emergence velocity (0.33 m a −1 ) at the ablation area of the Muz Taw Glacier, which was close to the corresponding values of Urumqi Glacier No. 1 (0.20-0.30 m a −1 ) in the Tien Shan [63] and Qiyi Glacier (0.32-0.43 m a −1 ) in the Qilian Mountains [62], but it was significantly smaller than that of Parlung No. 4 Glacier (3.7 m a −1 ) in the southeastern Tibetan Plateau [64], the debriscovered ablation area of Khumbu Glacier (5.7 m a −1 ) in the Nepal Himalaya [29], the ablation area of Hintereisferner (5.0 m a −1 ) in the Central Alps [65], Kesselwandferner (0-5.0 m a −1 ) in the Ötztal Alps [27], and Mer de Glace (4.0 m a −1 ) in the French Alps [26]. Glaciers in northwestern China are dominated by cold type, and they are characterized by relatively smaller mass accumulation and ablation compared with temperate glaciers in the southeastern Tibetan Plateau and Alps [66]. Particularly for the Muz Taw Glacier, very less annual precipitation recorded by Jimunai Meteorological Station directly determined limited glacier accumulation, and our field observations also confirmed this. Thus, the weak dynamic thickening can potentially be attributed to poor mass accumulation.

Limitations and Further Work
Although the applied TLS can provide a scanning range of 6 km and the scanned area exceeded 4.0 km 2 for the four years, our study only focused on a portion of the ablation area of (accounting for about 17.6% of the total glacier area) the Muz Taw Glacier, since some voids occurred at the upper elevations. This is mainly because of the limited scanning angle of the TLS. Auxiliary surveys from other techniques such as a cut-price unmanned aerial vehicle (UAV) can be required. Consequently, the combination of TLS and UAV had been done in many studies (e.g., [22,67,68]). A 30 m-high glacier monitoring tower had been placed at the fore of the glacier, which had a good visual angle relative to the entire glacier surface. Therefore, the TLS device can be mounted at the tower to increase the scanning range. For further work, we will integrate the TLS, UAV, and glacier monitoring tower to acquire the entire glacier surface terrain at the seasonal and monthly scales. In such a case, high-resolution and high-precision geodetic glacier mass balance can be yielded, which will bring insight into the mass balance processes and their mechanisms. Moreover, accurate quantification of the submergence/emergence velocities is important for geodetic mass balance estimates and simulation of the glacier dynamic process; therefore, we propose to strengthen the glaciological mass-balance measurements by using stakes and snow pits.

Conclusions
Glacier meltwater in the Sawir Mountains is an important freshwater resource for populations and hydro-economies in Jeminay and its surrounding counties; mass-balance measurements are fundamental to evaluate the impact of glacier changes on water resources. However, studies on detailed glacier mass changes are still vacant. This study presents high-precision evolution of annual surface elevation and geodetic mass changes in the ablation area of the Muz Taw Glacier (Sawir Mountains, China) over the latest three consecutive mass-balance years (2017-2020) based on multi-temporal terrestrial geodetic surveys. We mainly analyzed the spatiotemporal pattern of surface elevation and geodetic mass changes and the influences of emergence velocity on glacier surface lowering and discussed the possible mechanism of glacier dynamic thickening.
Our results revealed clearly surface lowering and negative geodetic mass changes, spatial changing patterns were maintained as generally similar for all of the three investigated periods with the most negative surface lowering (approximately −5.0 to −4.0 m a −1 ) around the glacier terminus and less pronounced thinning in the upper elevations, and the altitudinal increasing gradient was 0.69 m, 0.60 m, and 0.95 m per 100 m, for the mass-balance years 2017-2018, 2018-2019, and 2019-2020, respectively. The gradient of altitudinal elevation changes was commonly steep at the low elevations and gentle in the upper-elevation parts (the turning point occurred at the altitude of 3210 m a.s.l. for all of the three periods), and reduced surface lowering was observed at the glacier terminus. The majority of the more negative mean glacier surface changes were dominant for west, southwest, and south aspects, confirming that aspect had a major impact on glacier ablation. Calculated emergence velocities ranged from 0.11 to 0.86 m a −1 with pronounced spatial variability and higher velocities at the left and upper right sides, which were mainly controlled by surface slope, ice thickness, and the movement of tributary glaciers.
High air temperature in summer 2019 was one of the main factors for clearly negative surface elevation and mass changes in 2018-2019, and our foundlings suggested that glacier surface terrain can also influence surface elevation and mass changes. Emergence velocities made slight compensations to the surface ablation at the ablation area with a proportion of 14.9%, indicating that dynamic thickening had small contributions to the glacier surface evolution. The cold Muz Taw Glacier is located in a continental climate setting characterized by relatively smaller mass accumulation and ablation compared with temperature glaciers in the southeastern Tibetan Plateau and Alps, which can probably result in weak dynamic thickening. Our study only focused on a portion of the ablation area of the Muz Taw Glacier. In the follow-up work, we will integrate the TLS, UAV, and glacier monitoring tower to realize the whole glacier surface surveys at the seasonal and monthly scales, and submergence/emergence velocities should be carefully considered for geodetic mass balance estimates and simulation of glacier dynamic processes.