High ‐ Resolution Monitoring of Glacier Mass Balance and Dynamics with Unmanned Aerial Vehicles on the Ningchan No. 1 Glacier in the Qilian Mountains, China

: Glaciers located in the Qilian Mountains are rapidly retreating and thinning due to climate change. The current understanding of small glacier mass balance changes under a changing climate is limited by the scarcity of in situ measurements in both time and space as well as the resolution of remote sensing products. Unmanned aerial vehicles (UAVs) provide an unparalleled opportunity to track the spatiotemporal variations in glacier extent at a high resolution and the changing glacier morphological features related to glacial dynamics. Five measurements were performed on the Ningchan No. 1 (NC01) glacier in the Qilian Mountains between 18 August 2017 and 13 August 2020. The glacier changes displayed in the digital orthophoto maps (DOMs) and digital surface models (DSMs) show a 7.4 ± 0.1 m a − 1 retreat of the terminus of NC01, a mass balance of − 1.22 ± 0.1 m w.e. a − 1 from 2017 to 2020, and a maximum surface velocity of 3.2 ± 0.47 m from 18 August 2017 to 26 August 2018, which clearly show consistency with stake measurements. The surface elevation change was influenced by the combined effects of air temperature, altitude, slope, and surface ve ‐ locity. This research demonstrates that UAV photogrammetry can greatly improve the temporal and spatial resolution of glaciological research.


Introduction
The retreat of mountain glaciers around the world is clearly a sign of global climate change [1]. Glaciers provide water for irrigation, power generation, sanitation, and religious activities for large populations and are important parts of rivers in High Mountain Asia (HMA) [2][3][4]. A continuous negative mass balance results in the decrease in the amount of ice in these mountain ranges [5][6][7][8]. Based on current climate projections, glacier mass loss will accelerate and meltwater runoff will increase in the coming decades [9,10]. Therefore, monitoring glacier extent, mass balance, and surface velocity are critical to understanding how climate perturbations affect glaciers [11]. Studies have shown that due to the large number of small glaciers, the impact of their changes is much greater than expected [12]. A deep understanding of small-scale glacier changes can provide a reference for future glacier models and predictions.
Field and geodetic methods can both be used to measure glacier extent, mass balance, and surface velocity, and each method has its own advantages and disadvantages [13]. Field-based glacier monitoring is hampered by glacier inaccessibility, harsh weather, and high costs associated with data acquisition [14]. Satellite remote sensing-based glacier monitoring is limited by the spatial and temporal resolution and availability of appropriate data [15], especially in small glaciers. The coarse resolution of most satellite data cannot be applied to calculate the length, area, surface elevation change, and velocity of a glacier [16,17]. Therefore, higher-resolution data are increasingly needed, whether they are related to the glacier dynamics or derived from models of glacier responses to climate change. Recently, cutting-edge technologies such as Terrestrial laser scanning (TLS) and UAVs have been widely used in geodetic mass balance estimates [18][19][20][21]. TLS and UAVs have higher precision than satellite images, implying that useful glacier changes can be derived in detail from images acquired over only a few days, even for a small glacier with a slow rate of change. This is a significant advantage because changes in the glacier surface over short time periods have been difficult to obtain [22,23]. Monthly fluctuations in glacier changes can now be observed from time series of UAV images, something that was previously only possible using radar images and only for fast-flowing in large glaciers [24,25]. Compared to TLS, UAVs provide measurements over large areas and can avoid the impact of terrain obscuration with ease and efficiency. In addition, the economic cost of UAVs is much lower.
At present, UAVs are being used by an increasing number of studies in natural science fields [23,26,27]. They are economically and technically well suited for surveying glaciers and have been proven to be a very useful glaciology tool [17,[28][29][30][31][32][33][34][35][36][37][38]. High-resolution DOMs and DSMs provided by UAVs have made it possible to describe glacier extent, mass balance, and velocity in detail [34], which fills the gap between field and geodetic methods.
Changes in the glaciers in the Qilian Mountains indicate a sensitive response to climate variation [39][40][41][42][43]. In recent years, glaciers in the Qilian Mountains have been quickly retreating, which may have affected the water resources in the Hexi Corridor [40,44]. Glaciers in the Qilian Mountains are small, with an average area of 0.6 km 2 , and 85% of the glacier areas are less than 1 km 2 [45]. Traditional remote sensing images are difficult to capture the dynamic changes of small glaciers. For this reason, high-resolution data are needed to determine how climate factors affect glacier retreat and to model future glacier mass changes in the Qilian Mountains. In this study, five DOMs and DSMs of the whole Ningchan No.1 (NC01) glacier in the Qilian Mountains were produced after UAV survey flights of the glacier from 2017 to 2020. The main objectives of this study were (1) to obtain detailed data on glacial area change, mass balance, and surface velocity through UAVs and validate them through independent stake measurements during the same period; (2) to combine multiple factors (air temperature, altitude, slope, surface velocity, etc.) and analyze the causes of dynamic changes in glacier surface elevation.

Study Area
The NC01 glacier is a mountain glacier located on the north slope of the eastern Qilian Mountains (Figure 1a). It had an area of 0.77 km 2 and an altitude range of 4200-4640 m a.s.l. in 1972 [45]. The annual average temperature is approximately −6.0 °C, and the mean annual precipitation (MAP) is >800 mm based on the automatic weather station on the NC01 lateral moraine (4450 m a.s.l.) [41]. In recent years, NC01 has experienced rapid and accelerated shrinkage [40,46]. The mean annual mass balance was approximately −0.7 and −0.9 m w.e. in 1972-2010 and 2010-2015, respectively [41,46]. The mean equilibrium line altitude increased from 4500 m in 1972 to 4680 m in 2010-2015 [41,45]. In addition, the area of the NC01 glacier had decreased to 0.39 ± 0.04 km 2 in 2014 [41]. There are several small ice streams and a deep canyons on the exposed ice in ablation season (Figure 1d). NC01 has relatively little debris covering the ice surface, but there are scattered boulders (Figure 1f,g).

Unmanned Aerial Vehicle Surveys
The Phantom 4 Pro quadcopter system equipped with a 20 megapixel fixed-lens camera from Dajiang Company was used in this research. It is a rotorcraft that can fly for approximately 20 min with each battery charge.
Five UAV surveys of the entire NC01 glacier region were performed on 18 August 2017, 18 July 2018, 26 August 2018, 18 August 2019, and 13 August 2020 (they are referred to below as 18 August 2017, 18 July 2018, 26 August 2018, 18 August 2019, and 13 August 2020, respectively). Each survey consisted of for flights to cover NC01 and were conducted between 10:00 and 12:00 Beijing time to maximize the flight stability and image quality, as the wind speeds usually increased during the day, especially after 14:00. The flights were planned using the DJI-GS-Pro flight planning software to maintain 66% lateral overlap and 88% longitudinal overlap between the images. And the average flight height was about 150 m.

Orthophotos and Digital Surface Models Generated from UAVs
Images (n ≈ 2300, Figure 2a) collected by each UAV survey were processed into DOMs and DSMs of the glacier and the surrounding environment using a structure-frommotion photogrammetry (SfM) workflow [47]. Further details on the SfM process can be found in Immerzeel et al. [29], Kraaijenbrink et al. [30], Wigmore and Mark [33], Rossini et al. [34], Verhoeven [48], and Kraaijenbrink et al. [49]. The camera positions and the orientation of each image and the image overlap were initialized using the GPS coordinates recorded by the UAV. Due to the influence of reconstruction error, there was some interference in the sparse point cloud. All points with a reprojection error greater than 0.5 pixels were removed [17]. The dense point cloud was used to construct a DSM with a resolution of 0.1 m ( Figure 2c). In addition, a DOM with a resolution of 0.1 m was created using elevation information (Figure 2b).

Glacier Extent, Mass Balance and Surface Velocity Derived from UAVs
Before any differential analysis, the multitemporal orthophoto images and the corresponding DSMs needed to be comatched to remove the horizontal or vertical offsets [50]. A common method that was proposed by Nuth and Kääb [51] was used here. This method adjusts the relationship between the DSM elevation differences and the aspect in the regions outside of glaciers and further adjusts the correlation. Based on this approach, in each pair of DSMs, one DSM was used as a reference ("master") and the other DSM was used as a "slave". The "slave" DSM moves iteratively along the x-, y-and z-axes in pixel fractions to minimize the standard deviation of elevation differences from the "master" DSM [52,53]. Shifts along the x-and y-axes are the also performed in the corresponding orthophoto images.
Even after the above correction, there are still elevation-dependent vertical biases in the two DSMs. This bias can be corrected based on the relationship between the elevation b difference and the maximum curvature over the stable terrain off of the glaciers [6,54]. After the adjustments, most of the off-glacier errors were approximately 0 m with a standard deviation of 0.12 m (Figure 3), which is an acceptable error [29].
The glacier extents in the five surveys were digitized using manual methods and adjusted images in ESRI ArcMap 10.2. The accuracy of digitizing of the glacier extent can be strongly influenced by the spatial resolution of the orthophotos [55]. We estimated the uncertainty with a pixel resolution of 0.1 m for the orthophotos. Approximately evenly spaced transects of clearly distinguishable surface features (70 to 80 points), i.e., primarily large boulders and small supraglacial streams, were visually selected on the orthophotos [29]. Following vector creation, a velocity surface was interpolated for the glacier using the kriging interpolation tool with a 0.47 m root-mean-square (RMS). The surface elevation changes were calculated using DEM differencing. A density of 850 ± 60 kg m −3 was employed to convert ice to actual mass balance (cf. Huss, 2013).

Mass Balance and Surface Velocity from the Glaciological Method
There have been 22 stakes in NC01 since September 2011 (Figure 1b), which were used to measure the mass balance and surface velocity of NC01. This is called the "glaciological method" [56]. Here, they were used to validate the accuracy of the UAV results. To measure the glacier mass balance and velocity using the stakes, a pair of 5800 real-time kinematic (RTK) differential GPS (dGPS) units manufactured by Trimble Company were used. The stakes have occasionally been lost, buried under snow, or broken by ice movement, but our observation network allowed for measurements with a minimum of 15 points. The glacier mass balance for NC01 was measured following the glaciological methodology described as ∑ , where si and bi represent the area between the two contours and the corresponding mass balance, respectively, and S is the total area of NC01 [41]. Surface velocities were measured by determining the locations of the stakes using dGPS. To compare the UAV data with the data from the stakes, we chose the time period that was consistent with the UAV surveys. The precision of the horizontal GPS measurements is approximately 0.1 m in glacierized areas [46]. Due to the readings and the tilt of the stakes, the mass balance at each stake has a 0.1 m uncertainty [17,57].

Glacier Terminus Retreat from 2014 to 2020
According to UAV surveys, 5 DOMs and DSMs were obtained. After registration, the glacier boundaries were digitized by manual interpretation. The glacier extent on 19 August 2014 derived by Cao et al. [41] was also used to obtain detailed information about changes in the glacier area and terminus. The glacier terminal retreated 46.9 ± 0.1 m from 2014 to 2020, and it retreated by 22.1 ± 0.1 m from 2017 to 2020, averaging 7.4 ± 0.1 m a −1 ( Table 1). The retreat rates from 2017 to 2018 (7.7 ± 0.1 m a −1 ) and 2018 to 2019 (8.4 ± 0.1 m a −1 ) were faster than those from 2019 to 2020 (6.0 ± 0.1 m a −1 ) (Figure 4).

Glacier Surface Elevation Change and Velocity
Though DEM differencing was generated by the UAV data, changes in the surface elevation were highly variable across NC01, including areas of both loss and gain ( Figure  5). The surface elevation change was −4.30 ± 0.12 m in the period of 18 August 2017-26 August 2020, with a mean annual mass balance of −1.22 ± 0.10 m w.e. The rate of glacier thinning varied from region to region on the surface of the glacier. On the whole, the rate of glacier thinning was more dramatic in the lower part of the glacier and slightly slower in the upper part of the glacier (Figure 6a), which was mainly controlled by temperature; the lower the altitude was, the higher the positive accumulated temperature, and the stronger the ablation.    Figure 7). The mean displacement that occurred in the 39-day interval accounted for 67.6% of the total annual displacement. In general, the surface velocity along the glacier flowline increased gradually from the upstream to the central region, reached its maximum, and then decreased gradually to the terminal edge [58,59]. According to the velocities we obtained from 18 July 2018-26 August 2018 and 18 August 2017-26 August 2018, the velocity distribution pattern of NC01 is generally consistent with the general pattern (Figure 6c).

Comparison of the Results of the UAV and Glaciological Method
Strictly speaking, the mass balance obtained by the glaciological method and the geodetic method are different. The former measures the surface mass balance [60], whereas the latter measures glacier elevation changes and mainly include surface mass balance and vertical velocity [58]. Vertical velocity refers to the upward or downward flow of the 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 [58]. The emergence and submergence velocity will have a certain influence on the dynamic process [61]. Dehecq et al. [62] pointed out that with continuous glacier loss and the decrease of accumulation area ratio (AAR), emergence velocity will gradually decrease. For example, this phenomenon has been observed in both the Lirung Glacier [63] in the Nepal Himalayas and the Qiyi Glacier in Qilian Mountains [64].
By comparing the stake ablation data and the surface elevation changes from the DSM differencing results, we found that in most stake positions, the two results corresponded very well, with an absolute difference of less than 0.3 m w.e. (Figure 8a). Due to the small size (less than 0.5 km 2 ), the slower surface velocity (less than 3.2 m a −1 ), and the ELA exceeding the maximum height of the glacier [41], the emergence velocity of the NC01 glacier should be very small. It is difficult to determine this by the difference between the glacier mass balance and the surface elevation change, because this value itself is within the error range. Therefore, the influence of emergence velocity on the results has been neglected. There is a considerable difference between only one stake ablation value and the DSM difference, mainly because that stake is located in an area with a relatively large slope. During the observation of this stake, the ablation was slightly overestimated.
The absolute difference between the velocity after interpolation through the feature points and the GPS measurements in stake positions is within 1 m. Compared to the precision of the surface elevation change, the precision of the velocity is slightly lower, which may be because it was obtained indirectly through interpolation, where some uncertainty was generated. Similarly, one stake displacement value does not correspond to the interpolation results of the feature points (Figure 8b). The main reason is that the feature points in the upper part of the glacier are relatively sparse, which leads to the interpolation results underestimating the velocity in some areas. The findings reveal that UAVs have a very high accuracy in glacier mass balance research and a relatively high accuracy in velocity research.

Factors Influencing Glacier Changes
Pan et al. [40] and Cao et al. [41,43,46] documented that glacier shrinkage in this study area can probably be attributed to the increase in air temperature. According to the datasets recorded at the Menyuan meteorological station (2924 m asl) located 20 km away from the NC01 glacier, the summer (JJA) air temperatures in 2018, 2019, and 2020 were 14.0 °C, 12.4 °C, and 12.2 °C, respectively, which were significantly higher than the mean values from 1972-2010 [46]. This partly explains the dramatic glacier terminus retreat and negative mass balance from 2017 to 2019.
As the air temperature is the main factor controlling glacier melting and the air temperature is related to the altitude, there is an obvious correlation between the mean glacier surface elevation change and altitude (Figure 6a), with the exception of a location that is approximately 400-800 m from the summit. This may be because the slope here is less steep and more radiation is received; hence, melting is stronger. In addition, the same pattern is observed for surface elevation changes during the period of 18 August 2017-26 August 2018 (Figure 6b). However, according to the characteristics of the 39-day interval (18 August 2018-26 August 2018) melt, melting was strong in the entire glacier area, which was related to the altitude but was not obvious. It may be that the summer temperature is so high (14.7 °C at Menyuan station) that the entire ice surface was in a state of intense melting, which also led to a large amount of melting during the aforementioned year (18 August 2017-26 August 2018) (Figure 6b).
The velocity along the flowline increased gradually from the upstream region to the central region and reached its maximum value at a location 800 m from the summit. However, a second peak appeared in the lower part of the glacier, about 1000 m from the summit. Cao et al. [65] also found that the maximum glacier velocity occurred in the lower location. There are two possible reasons for this occurrence: the slope is steeper (Figure  6c), and the glaciers are thicker [41] at this location. According to Glen's flow law, the surface velocity is positively related to slope and thickness [58], which may be the reason for the high velocity in this area.

Comparison with Previous Results
A comparison with previous results showed that the retreat rate at the terminus of the glacier showed late acceleration ( Table 2). The rate of terminal retreat has gradually increased from 4.0 m a −1 in 1972-1995 to 7.4 ± 0.1 m a −1 at present. Moreover, the mass balance of the glacier also became more negative, changing from −0.7 m w.e. a −1 in 1972-2010 to −1.22 ± 0.10 m w.e. a −1 at present. As climate change continues to accelerate, glaciers will continue to retreat at an accelerating rate. The same phenomenon occurs for many glaciers, such as in the Tien Shan, the Nyainqentanglha Mountains, and the Himalayas [66][67][68][69]. However, glacier velocity showed a decreasing trend; whether measured using stakes or interpolation from distinguishable surface features, the maximum speed of the glacier from 2010 to 2012 was 4.0 m a −1 [65], while from 2017 to 2018, the maximum speed of the glacier was only 3.2 ± 0.47 m a −1 . Dehecq et al. [62] compared the velocity to the mass balance of many glaciers in High Mountain Asia and found that in most regions, glacier velocity tended to slow and suggested that changes in velocity can be explained by changes in gravity driving stress, which are largely controlled by changes in ice thickness. Therefore, the slowing of glacier NC01 may also be due to continuous negative balance.

Conclusions
UAV surveys were conducted to reconstruct a glacier surface with accuracies higher than those previously reported using satellite images, allowing for small changes to be measured. Through five surveys over the NC01 glacier, we obtained detailed glacier characteristics for the past 3 years, including the retreat of the glacier terminal edge, the change in the glacier surface elevation, and the glacier velocity. The results show that the glacier terminal retreat was 7.4 ± 0.1 m a −1 , the annual mass balance was −1.22 ± 0.1 m w.e. from 2017 to 2020, and the maximum velocity of the glacier was 3.2 ± 0.47 m from 2017 to 2018. These data were validated through independent stakes measurements, which clearly showed the consistency of the results. Glacial thinning in the 39 days of the ablation season accounted for 47.8% of the thinning in the annual interval. The displacement occurring in the 39 days of the ablation season accounted for 67.6% of the total displacement. This implies that glacier melting and velocity in the summer are very intense. In addition, a comparison with previous results revealed that the retreat of the terminal edge and thinning of the glacier showed a trend of late acceleration, while the glacier velocity slowed as the glacier thinned.