1 Introduction

Glacier mass-balance provides significant information in understanding climate-glacier processes and its consequences such as water resource estimation and global sea level rise (IPCC, 2013). Numerous studies have been performed on glacier area, length, thickness (volume) and mass balance by using remotely sensed observations (Li KM et al., 2011 ; Gardelle et al., 2013 ; Shangguan et al., 2015 ) and in-situ measurements (Zemp et al., 2015 ). Mass-balance can be derived by using in-situ direct glaciological measurements and repeated geodetic surveys. The direct glaciological method monitors mass balance using ablation stakes and snow pits and net accumulation is measured by digging pits in the accumulated areas and stakes are drilled into ablation of the glacier surface (Cogley et al., 2011 ). Although the glaciological method provides quantitative results at high temporal resolution (Zemp et al., 2015 ), it is usually limited to individual glaciers. Direct measurements are time consuming and labor intensive. Thus, at present, only 41 glaciers have been studied for ongoing mass-balance series with more than 30 years of observations (WGMS, 2017).

Compared with traditional glaciological measurements, the geodetic method provides glacier volume changes by repeated mapping from ground, air or spaceborne surveys and subsequent differences of multi-temporal digital elevation models (DEMs) of glacier surface topography (Zemp et al., 2013 ). The geodetic mass-balance extends the glaciological sample both in space and time, which has become a widely used method to calculate glacier mass balance and evaluate the contribution of glaciers to runoff and periodically validate and calibrate long-term glaciological mass-balance series (Cogley, 2009; Cogley et al., 2011 ; Zemp et al., 2013 ). Whereas, the limited spatiotemporal resolution of existing DEMs usually influences the calculation of accurate geodetic mass balance. A new generation of terrestrial lidar techniques has offered dense point clouds and high resolution and precision DEMs of glacier surface (Fischer et al., 2016 ).

In China, glacier mass-balance monitoring mainly focus on Tien Shan, Qilian Shan and the Tibetan Plateau (Li ZQ et al., 2011 ; Wang et al., 2017 ; Zhu et al., 2017 ). Even though glaciers in Sawir Mountains have undergone an acceleration of mass loss and terminal retreat with the consequence of floods, these glaciers have received little attention at present (Huai et al., 2015 ). The present study aims to provide a detailed and precise mass balance and front variation monitoring of Muz Taw Glacier using an ultra-long-range terrestrial laser scanner (TLS). We then use Muz Taw Glacier as a case study to discuss the applicability of TLS systems to monitor the mass balance of mountain glaciers. Thus, the present study provides valuable data for poorly understood mass balance and front variation of mountain glaciers.

2 Study area

The Sawir Mountains are located south of the Ertix River (Wang, 1988). According to the first Glacier Inventory of China (GIC), there were a total of 21 glaciers with a total area of 16.84 km² in Chinese Sawir Mountains (Shi, 2005). Due to accelerated shrinkage, glacier area reduced by 42.74% at a rate of 0.14 km²/a from 1959 to 2013 (Huai et al., 2015 ). With the growth of numerous water resource management problems (e.g., flood protection, water supply, and operation of hydroelectric facilities) in this region, timely and correct assessment of glaciers in the Sawir Mountains is crucial to fill the gap, and it has an extremely important significance for Jeminay's economic development and people's livelihood.

Muz Taw Glacier (47.06°N, 85.56°E; glacier ID 5A259C0001) is a northeast-orientated valley glacier with an area of 4.27 km² and a length of 3.7 km according to the first GIC (Shi, 2005), located at the northern side of the Sawir Mountains (Figure 1a). The elevation from the glacier terminus to the highest point ranged from 3,090 m to 3,800 m a.s.l.. The snow line and ice volume of this glacier were about 3,310 m a.s.l. and 0.28 km³, respectively. The average ice thickness was 66 m. The mass-balance monitoring program was implemented beginning in 2014 (Figure 1b). From field observations, there was no debris on the glacier surface. Muz Taw Glacier has been dramatically shrinking due to extensive melting over the past several decades.

3 Data and methods 3.1 Riegl VZ^®-6000 ultra-long-range TLS

Traditional TLS systems use class 1 laser beams with wavelengths around 1,500 nm, with low reflectance over snow and ice cover, and the possible scanning distance is limited to a few hundred meters (Deems et al., 2013 ), which restricts the application of TLS surveys in the cryosphere (Rabatel et al., 2008 ). The Riegl VZ^®-6000 ultra-long-range TLS operates with class 3B laser beams (wavelengths around 1,064 nm), and the system is exceptionally well-suited for measuring snowy and icy terrain. Faster surveys (up to 222,000 measurements/s) are possible even at long range (>6,000 m) with unprecedented high accuracy and precision (RIEGL Laser MeasurementSystem, 2014a).

The Riegl VZ^®-6000 uses near-infrared laser signals to capture repeated dense point clouds of the glacier surface terrain, based on the principle of time-of-flight measurement with echo digitization and online waveform processing (RIEGL Laser Measurement Systems, 2013). The laser pulse transmitter emits a beam of laser pulse, which reaches the target object and then reflected back to the laser receiver, and the pulse is used to calculate the time difference of the pulse (RIEGL Laser Measurement Systems, 2014a).

3.2 Field measurements

On June 2, 2015 and July 25, 2016, Riegl VZ^®-6000 TLS measurements of Muz Taw Glacier were performed with two scanning stations (Figure 1c), coincident with same day of direct glaciological measurements. Each scanning station was fixed with a GPS-leveling point in order to prohibit ground motion. The instrument was mounted on a tripod placed on stable bedrock for each scan (Figure 2). Three-dimensional (3 D) coordinates of the two scanning stations were obtained using the real-time kinematic (RTK) global positioning system. The coordinated system of RTK surveys was the World Geodetic System 1984 (WGS84).

For TLS work, both data quality and acquisition time should be considered, thus surveying parameters are of significant importance. To avoid range ambiguity and to capture point clouds of the glacier surface to the maximum extent, coarse scanning was first performed with vertical and horizontal angles range of 60°–120° from zenith and 0°–360°, respectively, and corresponding laser pulse repetition rate was set as 50 kHz. Based on coarse scanning, laser pulse repetition rate was selected as 30 kHz for the implementation of fine scanning of the glacier surface terrain region for two scanning stations. We should note that the overlap percentage of two scans was more than 30% to satisfy the basic requirements for multi-station adjustment (MSA) (Mukupa et al., 2017 ). Surveying parameters are listed in Table 1.

3.3 Point cloud data processing

Here, we use RiSCAN PRO^® v1.81 software to perform the point cloud data processing (RIEGL Laser Measurement Systems, 2014b), which includes:

1 Direct georeferencing. An additional sensor (e.g., navigation sensors, total station) is used to transform the laser scanner's system (Mukupa et al., 2017 ). In this study, GPS-RTK surveying points of each scanning station were used to transform the Scanner's Own Coordinate System (SOCS) into a Global Coordinate System (GLCS). The theory can be illustrated as (Lichti et al., 2005 )

$\overrightarrow {{r_g}} = \overrightarrow {{r_0}} + R\left(k \right)\overrightarrow {{r_s}} $

(1)

$R\left(k \right) = \left[ {\begin{array}{*{20}{c}}{{\rm{cos}}\;k}&{{\rm{sin}}\;k}&0\\[3pt]{{\rm{sin}}\;k}&{\cos \;k}&0\\[3pt]0&0&1\end{array}} \right]$

(2)

where, $\overrightarrow {{r_g}} $ is the vector of the object in GLCS, $\overrightarrow{r_s} $ is the vector of the same object in SOCS, $\overrightarrow {{r_0}} $ is the vector of SOCS origin in GLCS and k is the azimuth between scanning station and backsight target.

2 Point cloud registration. The position of each scan was fixed after direct georeferencing. Whereas, due to the influence of orientation, the point clouds of the overlapped areas cannot coincide completely. The iterative closest point (ICP) algorithm was used to complete multi-station adjustment (MSA) (Besl and McKay, 1992; Zhang, 1994). The orientation of each scan was iteratively modified to calculate the best fit among the points patches based on least-squares minimization of residuals until point cloud layers of each scan coincided (Besl and McKay, 1992).

3 Point cloud vacuation and filtering. We used an octree algorithm to achieve point cloud data vacuation. In particular, the algorithm was applied to the registered scans to generate points with equal spacing (Schnabel and Klein, 2006; Perroy et al., 2010 ). Then a terrain filter was used to remove noise and non-ground points (RIEGL Laser Measurement Systems, 2014b).

In the support of ArcGIS10.1 software, TLS-derived DEMs were generated using natural neighbor interpolation with a resampled resolution of 1 m.

3.4 Glaciological mass balance measurements

The glaciological mass balance measurements of Muz Taw Glacier were performed in 2014. About 20 ablation stakes were drilled into the glacier surface evenly at different elevation bands using a stream drill (Figure 1c). Generally, a spatial distribution of single-point ablation or accumulation, snow density (if there is snow cover) were determined by stakes and snow pits at intervals of a month from the beginning of May to early September each year. Therefore, the measured programs include the stake vertical height over the glacier surface, thickness of superimposed ice, thickness and density of each snow/firn layer at individual stakes. By using in-situ measured densities of snow and ice, stake and snow pit changing values were converted into mass (water equivalent) to calculate single point mass balance (the summation of snow/firn, superimposed ice and glacier ice mass balance). Then point observations can be extrapolated to glacier-wide mass balance contour-line and isoline methods (Xie and Liu, 2010).

3.5 Geodetic mass balance 3.5.1 Calculation of geodetic mass balance

The geodetic mass balance was calculated by the subtraction of two high-resolution DEMs. The volume change ΔV was determined by summing the elevation change Δh_i at the individual pixel r of different time periods (Zemp et al., 2013 ).

$\Delta V = {r^2}\mathop \sum \limits_{i = 1}^N \Delta {h_i}$

(3)

where N is the number of pixels covering the glacier at the maximum extent, and r is the pixel size (1m×1m).

To estimate glacier mass balance (B_geo), observed geodetic volume change must be converted to specific mass balance (m w.e.) by using the following Equations (4) and (5):

${B_{{\rm{geo}}}} = \frac{{\Delta V}}{S}\cdot\frac{\rho }{{{\rho _{{\rm{water}}}}}}$

(4)

${\rm{S}} = \frac{{{S_{t1}} + {S_{t2}}}}{2}$

(5)

where ρ presents the average density of volume changes of 850±60 kg/m³ as recommended by Huss (2013), S is the mean glacier area over the PoR, the two survey dates t₁ and t₂, assuming a linear change (Thibert et al., 2008 ; Zemp et al., 2013 ).

3.5.2 Uncertainty analysis

Uncertainties in the geodetic mass balance include systematic or random errors (Zemp et al., 2013 ), those include: (1) point cloud data acquisition errors, including terrain and atmospheric environment (wind and moisture); (2) data processing errors and DEM creation, i.e., registration (multi-station adjustment), point cloud filtering and interpolation (smoothing terrain information) (Gabbud et al., 2015 ; Hartzell et al., 2015 ; Fischer et al., 2016 ); and (3) spatial bias of the two DEMs and sampling density of snow and ice (Nuth and kääb, 2011; Zemp et al., 2013 ; Thomson et al., 2017 ).

Point clouds of the glacier surface should be captured in sunny, windless weather to remove the influence of wind and moisture. The octree algorithm built the topological relationship of scattered points to reduce the point cloud, in order to preserve terrain information as much as possible. To eliminate the spatial bias of the two DEMs, scanning point clouds of 2015 were set as the reference; a relative registration selects the ground control points of the non-glaciated stable terrain to match the point clouds of the two periods.

Uncertainties of the average TLS-derived glacier surface elevation changes for individual glaciers were quantified based on geostatistical analysis methods (Rolstad et al., 2009 ), which can be approximated by

$\sigma _{\overline {\Delta h{\rm{TLS}}} }^2 = \sigma _{\Delta h{\rm{TLS}}}^2 \cdot \frac{{{S_{{\rm{cor}}}}}}{{5S}}$

(6)

where S is the glacierized area of the subtracted DEMs, which is same as the S defined by Equation (5) at the annual scale, S_cor is the area over which there are spatial correction errors in the subtraction of the two DEMs, and we conservatively assume S_cor = S (Fischer et al., 2016 ). The σ_ΔhTLS is the standard deviation for TLS-derived bedrock (non-glaciated ground) elevation changes.

The uncertainty σ_BTLS in the homogenized mass balance results from errors in the geodetic mass change (Huss et al., 2009 ):

${\sigma _{{\rm{BTLS}}}} = \pm \sqrt {{{\left({\overline {\Delta h{\rm{TLS}}} \cdot {\sigma _\rho }} \right)}^2} + {{\left({\rho \cdot {\sigma _{\overline {\Delta h{\rm{TLS}}} }}} \right)}^2}} $

(7)

where σ_ρ is the uncertainty of density. $\overline {\Delta h{\rm{TLS}}} $ is the mean of TLS-derived geodetic elevation changes and the related uncertainty relies on the accuracy of the DEMs. Derived values are given in Table 2.

4 Results and discussion 4.1 Glacier front variation (length changes)

Glacier front variation data are a key element for evaluating the regional representativeness of the few glaciological measurements at spatiotemporal scales (Zemp et al., 2015 ). Using Landsat TM/ETM+ images, previous studies show that the glacier terminus had an observed retreat of 269.19 m from 1977 to 2013, at a retreat rate of 7.28 m/a (Huai et al., 2015 ). However, the low spatial resolution images usually limited the accuracy of glacier front variations. As in-situ measurements of glacier terminus were poor, in this work, based on the visual interpretation of the TLS-derived high-resolution point clouds, we extracted six stable landmarks and calculated their corresponding distance to Muz Taw Glacier terminus, and those landmarks distributed evenly at the front of the glacier (Figure 3). Observed glacier front variations ranged from 8.4 to 10.6 m, with an average value of 9.3 m from June 2, 2015 to July 25, 2016 (Table 3), and the annual retreat rate was 8.1 m, which demonstrated an accelerated terminus retreat rate in recent years.

4.2 Glacier surface elevation changes and geodetic mass balance

The map of TLS-derived geodetic elevation changes revealed a serious ice loss at the ablation area (Figure 4) and the glacier tongue showed more negative ice loss than the upper part. The mean elevation change was −4.55 m for the ablation area between June 2, 2015 and July 25, 2016. Elevation changes measured along the main flow line allow more detailed insights into the characteristics of glacier behavior. The elevation differences showed well changing rule, that was to say, the ice thinning value slackened when the length increased (the distance from glacier terminus) (Figure 5), which agrees well with field observations. Figure 4 illustrates the distribution of individual ablation stakes (black crosses) and TLS-derived geodetic elevation changes (black digitals) were calculated at the corresponding stakes, which closely match direct glaciological elevation changes (stake height over glacier surface) from in-situ measurements. Considering the uncertainties of the two measuring methods, glacier ice movement and the influence of terrain factors, differences in surface elevation changes of the point measurements were neglected for this work, indicating that TLS-derived glacier surface elevation changes are reliable.

For the conversion of observed glacier elevation changes into the mass budget a density of ice/firn/ snow needs to be assumed. Here, a recommended density assumption of 850±60 kg/m³ is used and the density was provided by Huss (2013) using an empirical firn densification model with idealized surface mass balance forcing. Using Equations (3)–(5), the mass balance of the ablation area was determined, and we selected 375 stable points in the non-glaciated ground to evaluate the uncertainty of TLS-derived elevation changes and geodetic mass balance (Equations (6), (7)). The mass balance uncertainty of the observed period was ±0.378 m w.e.. and the calculated balance was −3.864±0.378 m w.e.. Compared to the mass balance, the uncertainty was very small, which indicates a high accuracy of TLS-derived DEMs.

4.3 Applicability and recommendation of the TLS system to monitor the mass balance

Due to high temporal-spatial resolution and ultra-long-range of Riegl VZ^®-6000, previous studies have shown that this technique can provide accurate mass balance data of very small glaciers in the Swiss Alps (Fischer et al., 2016 ) and surface melting values of Alpine glaciers at the seasonal and hourly scales (Gabbud et al., 2015 ). What's more, TLS can work without contact with objects and in a rugged environment, thus allowing measurements of dangerous or inaccessible areas of the glacier surface. Although we had inserted more than 20 ablation stakes on the Muz Taw Glacier surface, they were distributed evenly at different elevation bands. Some areas, such as the terrain of firn basin is precipitous, and the surface topography of the tributary glacier (right branch) is also dangerous (Figure 1c). Therefore, it is not easy for us to insert ablation stakes and dig snow pits in those inaccessible areas. However, glacial terrain and insufficient spatial distribution of stakes have significant influence on the calculation of glacier-wide mass balance. TLS technology makes it possible for us to monitor those inaccessible regions, so we recommend the application of Riegl VZ^®-6000 TLS for future mass balance monitoring of Muz Taw Glacier. To better monitor the mass balance, several recommendations need to be followed to improve the results of glacier measurements: (1) good weather (dry and windless atmosphere) is of great importance for field surveys; (2) stable and visual angles of the scanning stations are necessary between surveys, which could minimize registration errors among the scans of the same time; (3) adequate and visible areas of non-glaciated stable terrain are required to guarantee a good quality of the relative registration of different temporal point clouds; (4) TLS measurements should be implemented at the end of ablation season to decrease the influence of fresh snow cover; (5) the overlap percentage of each scanning station of the same period should be more than 30% to achieve an ideal mosaic of multi-station adjustment; and (6) monitoring of glacier density and glaciological mass balance at Muz Taw Glacier should continue. Long-term mass-balance series measurements are needed, especially for important regions where there is a shortage of in-situ monitoring data. The TLS-derived geodetic results can validate the reliability of glaciological balance and complement the glaciological method.

5 Conclusion and outlook

The Riegl VZ^®-6000 ultra-long-range TLS system is well suited for surveying snow- and ice-covered terrain and has a ranging accuracy better than 10 mm. The system can capture surface elevation changes and the mass of mountain glaciers at high spatiotemporal resolution, so the quality of TLS-derived geodetic net mass balance is very good. In the present paper, we have shown the applicability of an ultra-long-range terrestrial laser scanner to monitor mass balance of Muz Taw Glacier.

The TLS-derived high spatial resolution point clouds determined the front variation (glacier retreat) of Muz Taw Glacier of 9.3 m from June 2, 2015 to July 25, 2016. The geodetic elevation change was 4.55 m at the ablation area and the calculated balance was −3.864±0.378 m w.e.. Compared with in-situ measurements, the glaciological elevation change of individual stakes and the TLS-derived geodetic elevation change of corresponding points matched closely. Mean uncertainty in the TLS-derived mass balance was ±0.378 m w.e.. The uncertainty is attributed to data acquisition errors, data-processing and DEM creation errors. Thus, in further studies, we should take into account recommendations mentioned in Section 4.3 to generate high quality data and minimize the uncertainty of TLS-derived results.

Acknowledgments:

This research was supported by the National Natural Science Foundation of China (41601076, 41471058 and 91425303), the "Light of West China" program for Talent Introduction of Chinese Academy. We would like to thank the Tien Shan Glaciological Station for field surveys.