Summer Mass Balance and Surface Velocity Derived by Unmanned Aerial Vehicle on Debris-Covered Region of Baishui River Glacier No. 1, Yulong Snow Mountain

: Glacier retreat is a common phenomenon in the Qinghai-Tibetan Plateau (QTP) with global warming during the past several decades, except for several mountains, such as the glaciers in the Karakoram and the western Kunlun Mountains.

. Photograph of glacier terminal. Its position was measured. In pictures, (a) denotes the glacier terminus of Baishui River Glacier No. 1; (b) denotes the glacier mass avalanching from the glacier front; and (c) denotes a ground control point at the end of the glacier.
In recent years, Unmanned Aerial Vehicles (UAVs) have been used to understand the dynamic features of glacier surface by repeated UAV photogrammetry, especially for debris-covered glaciers [9,21,25]. For example, the heterogeneous patterns of mass loss, surface velocity, and surface temperature on the debris-covered Lirung Glacier in Central Nepal were assessed by high-resolution images from an UAV, which revealed the important catalytic role for supraglacial lakes and ice cliffs [9,20,26,27]. In addition, an assessment of image network geometry and strength without Ground Control Points (GCPs) was also carried out, which suggests that the glacier images without GCPs have a good ability to describe the glacier surface features when providing a larger of feature points or matches points among the UAV photogrammetry images [25]. Though several studies on the UAV monitoring vegetation have been conducted to understand the ecosystem environment in the QTP [28,29], there has been no special study on the glacier melting and surface velocity in the region. The glaciers in the southeast QTP are significantly controlled by the Indian and southeast monsoons. These glaciers have stronger glacier ablation and accumulation, high velocity, large crevasses, and debris. However, the dynamic features and mass loss mechanisms of the debris-covered and crevasses-developed glaciers are still unclear. The features and mechanisms play an important role in understanding the relationship between glacier and climate.
During the past several decades, the in situ records of glacier melting and retreat have been conducted in the central glacier surface and terminus of BRG1 in the southeast QTP, respectively. The in situ data of debris-covered glacier are lacking, so we did not know much about the glacier surface features of mass loss and move. These surface features are very important to assess the glacier mass balance, understand the distribution of crevasses and cliffs, and reveal the relationship between glacier melting and local climate change. Herein, we carried out the two tests using an UAV to survey BRG1 on 20 May 2018 and 22 September 2018, respectively. Using UAV monitoring, the Orthomosaic and Digital Elevation Model (DEM) were obtained with high spatial resolutions in the debris-covered glacier during the summer seasons of 2018. The Orthomosaics and DEMs of a debris-covered glacier during the summer filled the data gap and were very important to assess the glacier mass loss and surface velocity in the area.

Study Area
BRG1 is the largest glacier in the YSM, located in the southeast TP ( Figure 2). The total glacier area disappeared by 64.02% in the YSM during the period 1957-2017 and decreased to 4.48 km 2 in In recent years, Unmanned Aerial Vehicles (UAVs) have been used to understand the dynamic features of glacier surface by repeated UAV photogrammetry, especially for debris-covered glaciers [9,21,25]. For example, the heterogeneous patterns of mass loss, surface velocity, and surface temperature on the debris-covered Lirung Glacier in Central Nepal were assessed by high-resolution images from an UAV, which revealed the important catalytic role for supraglacial lakes and ice cliffs [9,20,26,27]. In addition, an assessment of image network geometry and strength without Ground Control Points (GCPs) was also carried out, which suggests that the glacier images without GCPs have a good ability to describe the glacier surface features when providing a larger of feature points or matches points among the UAV photogrammetry images [25]. Though several studies on the UAV monitoring vegetation have been conducted to understand the ecosystem environment in the QTP [28,29], there has been no special study on the glacier melting and surface velocity in the region. The glaciers in the southeast QTP are significantly controlled by the Indian and southeast monsoons. These glaciers have stronger glacier ablation and accumulation, high velocity, large crevasses, and debris. However, the dynamic features and mass loss mechanisms of the debris-covered and crevasses-developed glaciers are still unclear. The features and mechanisms play an important role in understanding the relationship between glacier and climate.
During the past several decades, the in situ records of glacier melting and retreat have been conducted in the central glacier surface and terminus of BRG1 in the southeast QTP, respectively. The in situ data of debris-covered glacier are lacking, so we did not know much about the glacier surface features of mass loss and move. These surface features are very important to assess the glacier mass balance, understand the distribution of crevasses and cliffs, and reveal the relationship between glacier melting and local climate change. Herein, we carried out the two tests using an UAV to survey BRG1 on 20 May 2018 and 22 September 2018, respectively. Using UAV monitoring, the Orthomosaic and Digital Elevation Model (DEM) were obtained with high spatial resolutions in the debris-covered glacier during the summer seasons of 2018. The Orthomosaics and DEMs of a debris-covered glacier during the summer filled the data gap and were very important to assess the glacier mass loss and surface velocity in the area.

Study Area
BRG1 is the largest glacier in the YSM, located in the southeast TP ( Figure 2). The total glacier area disappeared by 64.02% in the YSM during the period 1957-2017 and decreased to 4.48 km 2 in 2017 [22]. The elevation of the glaciers ranged from 5361 m a.s.l. to 4395 m a.s.l. The local climate is controlled by the Indian and southeast monsoons over the summer, while it was controlled by the south branch of the Westerly Belt during the winter and spring [24]. The precipitation is mainly concentrated in the warming seasons, i.e., June to September, which accounts for the 80%~90% of the total of annual precipitation [30]. The average air temperature is 5.08 • C in the monsoon season [31]. The glacier length decreased to 1.90 km from 2.65 km in 1957, and the glacier area decreased to 1.32 km 2 in 2017 [22]. In addition, the area of debris-covered glacier accounts for 12-13% of the total area, and the debris is mainly concentrated on the glacier surface below 4600 m a.s.l., i.e., the region of repeated UAV surveyed in the study. Note that the debris-covered glacier region had previously been a field-data-free zone because of its inaccessibility and the coverage of crevasses and cliffs to observe in situ in the glacier tongue ( Figure 1). 2017 [22]. The elevation of the glaciers ranged from 5361 m a.s.l. to 4395 m a.s.l. The local climate is controlled by the Indian and southeast monsoons over the summer, while it was controlled by the south branch of the Westerly Belt during the winter and spring [24]. The precipitation is mainly concentrated in the warming seasons, i.e., June to September, which accounts for the 80%~90% of the total of annual precipitation [30]. The average air temperature is 5.08 °C in the monsoon season [31]. The glacier length decreased to 1.90 km from 2.65 km in 1957, and the glacier area decreased to 1.32 km 2 in 2017 [22]. In addition, the area of debris-covered glacier accounts for 12-13% of the total area, and the debris is mainly concentrated on the glacier surface below 4600 m a.s.l., i.e., the region of repeated UAV surveyed in the study. Note that the debris-covered glacier region had previously been a field-data-free zone because of its inaccessibility and the coverage of crevasses and cliffs to observe in situ in the glacier tongue ( Figure 1).

Unmanned Aerial Vehicle
Considering the cost-effectiveness and aircraft stability, we used the Phantom 4 Pro (i.e., four-rotor drone) to survey the glacier surface. The UAV was produced by DJ-Innovations (DJI) and equipped with a CMOS image sensor with size of 1-inch and effective pixels of 20-megapixel, i.e., FC6310 [32]. The maximum flight time of the DJI Phantom 4 Pro was ~30 minutes in the normal temperature environment with a low-elevation region, while it ranged from 15 minutes to 30 minutes in the alpine glacier region. The maximum relative height was 500 m, while the maximum control distance was ~7 km. The drone was equipped with a Global Positioning System (GPS) and Global Navigation Satellite System (GLONASS), which provided accuracy position information. In addition, the controller and display equipment were required in order to make the normal drone flight. The drone hovered and took a photo in the air with a horizontal accuracy of ±0.3 m and vertical accuracy of ±0.1 m [32]. We found that the best time for UAV flight measurement was between 9 a.m.~12 a.m local time by analyzing many times of observations and experiments. Lower windspeed and less clouds occurred in the period, which benefitted the glacier survey using UAV flight. During the glacier survey, each flight time was set to 15-20 minutes in this region due to the

Unmanned Aerial Vehicle
Considering the cost-effectiveness and aircraft stability, we used the Phantom 4 Pro (i.e., four-rotor drone) to survey the glacier surface. The UAV was produced by DJ-Innovations (DJI) and equipped with a CMOS image sensor with size of 1-inch and effective pixels of 20-megapixel, i.e., FC6310 [32]. The maximum flight time of the DJI Phantom 4 Pro was~30 minutes in the normal temperature environment with a low-elevation region, while it ranged from 15 minutes to 30 minutes in the alpine glacier region. The maximum relative height was 500 m, while the maximum control distance was 7 km. The drone was equipped with a Global Positioning System (GPS) and Global Navigation Satellite System (GLONASS), which provided accuracy position information. In addition, the controller and display equipment were required in order to make the normal drone flight. The drone hovered and took a photo in the air with a horizontal accuracy of ±0.3 m and vertical accuracy of ±0.1 m [32]. We found that the best time for UAV flight measurement was between 9 a.m.~12 a.m local time by analyzing many times of observations and experiments. Lower windspeed and less clouds occurred in the period, which benefitted the glacier survey using UAV flight. During the glacier survey, each flight time was set to 15-20 minutes in this region due to the lower air temperature.

Flight Description
Two surveys were conducted on the debris-covered area of the BRG1 on 20 May 2018 and 22 September 2018 (Figure 1c), approximately corresponding to the beginning and end points of the ablation season, respectively. In the first survey, three flights were operated between 9:45 a.m. and 11:00 a.m., while six flights were finished between 9:30 a.m. and 12:00 a.m in the second survey ( Table 1). The surface elevation of the monitored glaciers ranged from approximately 4400 m a.s.l. to 4700 m a.s.l. The flight height of UAV measurement was set to~150 m above the glacier surface. However, the altitude of flight path had to be adjusted for each flight in order to keep a constant height above the glacier surface, especially in fights 3-6 (detailed information of the flights is shown in Figure 2 and Table 1). During the UAV flight measurement, the flight and photo collection were operated using the execution software of DJI GS Pro (https://www.dji.com/). For both campaigns, each image was acquired after focusing on the glacier surface, with a resolution of 5472 pixels × 3078 pixels. To meet the condition of image mosaic in space, the image the overlap was set to 80% in both the lateral and longitudinal directions with respect to the UAV flight path.

Ground Control Points
It was difficult to arrive at the debris-covered glacier surface in the glacier front, as travel paths were built on the south lateral moraine because the YSM was developed into a mature glacier resort. During the flight, 18 Ground Control Points (GCPs) were distributed on the south lateral moraine and bare rock at the end of the glacier. Approximately 60% of the GCPs was used for Orthomosaic and DEM generation, while the remaining 40% was used to validate the results of the generated models. That is, 6 of the 16 GCPs and 7 of the 18 GCPs were used as Ground Validation Points (GVPs) in the two surveys, respectively ( Table 1). The treatment of these GCPs was the same in terms of the process of generating Orthomosaic and DEM images because these points were fixed during the study period. The output Orthomosaics and DEMs, with the spatial resolutions of 0.05 m and 0.09 m, respectively, were resampled to 0.10 m in order to reduce the effects of any remaining motion blur as well as JPEG artifacts for further processing. In the procedure, the techniques of dealing with the georeferenced UAV photogrammetry were similar to those described by Immerzeel et al. [9], Kraaijenbrink et al. [26], and Rossini et al. [33].
In addition, the GCP and GVP coordinates were collected using differential global navigation satellite system (dGNSS) geodetic receivers (two CHCNAV X9 smart antennas and one CHCNAV controller HCE300). Then, the coordinates were corrected during post-processing [34]. Two identical smart antennas were used simultaneously: A base station and a rover. The base station was installed on the south lateral moraine of BRG1 and occupied this position for a cumulative time of over 16 h. The rover was used to measure the inner corner points of each GCP and GVP, which were marked using red spray paint to ensure visibility. The dGNSS used during the surveys had a horizontal precision of 2.5 mm + 0.5 ppm and a vertical precision of 5 mm + 0.5 ppm according to CHCNAV X9 specification Remote Sens. 2020, 12, 3280 6 of 15 data. The coordinate of each point was measured by taking a measurement every second with at least four satellites in the field of view of the receiver according to the proposed protocol [35]. During the record processing, at first, a measurement point was called as a dynamic point due to the coordinate changing in real-time with more than four satellites. After, the point changed into a floating point after several minutes and recorded when it changed into a fixed point. It was recorded after several minutes (around more 6 min) and an error threshold of 0.05 m was set during the measurement process. In other words, these precision errors of GCPs and GVPs were very small compared to the geodetic accuracy of the base station. Hence, they were negligible [9].

Orthomosaic and DEM Generation
Georeferenced photos collected during each UAV survey were assessed to obtain the vertical photogrammetry photos. The selected images must be cloudless and high-transmittance images, i.e., with sufficient quality and overlap from multiple positions and/or angles. These images were input and dealt with in the software Pix4Dmapper in order to generate the products of Orthomosaic and DEM. The data processing is shown in Figure 3. In addition, the process of images dealt with and output for the Orthomosaic and DEM was finished using Pix4Dmapper, which was fully described on the official homepage (https://www.pix4d.com/product/pix4dmapper-photogrammetry-software).
and GVPs were very small compared to the geodetic accuracy of the base station. Hence, they were negligible [9].

Orthomosaic and DEM Generation
Georeferenced photos collected during each UAV survey were assessed to obtain the vertical photogrammetry photos. The selected images must be cloudless and high-transmittance images, i.e., with sufficient quality and overlap from multiple positions and/or angles. These images were input and dealt with in the software Pix4Dmapper in order to generate the products of Orthomosaic and DEM. The data processing is shown in Figure 3. In addition, the process of images dealt with and output for the Orthomosaic and DEM was finished using Pix4Dmapper, which was fully described on the official homepage (https://www.pix4d.com/product/pix4dmapper-photogrammetry-software).
When the georeferenced photos were chosen and input, the GCPs coordinates were also needed in the professional software. These images were then adjusted and aligned using an image feature recognition algorithm, which automatically detects and matches unique image feature characteristics (i.e., tie points) in the overlap region among these images [36]. The generation of Orthomosaic and DEM covered the flight region, which needed to be validated using the GVP coordinates. Subsequently, the Orthomosaic of the debris-covered glacier was manually interpreted and clipped from the generated Orthomosaic, including the non-glacier region. In addition, the DEM image of a debris-covered glacier was also generated. The data process was conducted, and the data processing progress is consistent with other studies [9,25,26,33,37], especially the structure-from-motion applied to create Orthophtos and DEMs from the images collected during each UAV survey [38,39]. It is worth noting that the density of 900 kg/m 3 of glacier was used for converting changes in glacier surface elevation to changes in glacier mass balance.

COSI-Corr for Glacier Velocity
COSI-Corr (Co-registration of Optically Sensed Images and Correlation) is a software tool developed to co-register pairs of optical remotely sensed images (aerials or satellites), which is used to look for seismic ground deformations, glacier flows, horizontal changes between multitemporal images, and so on [40,41]. The main tasks performed by COSI-Corr are precise orthorectification, When the georeferenced photos were chosen and input, the GCPs coordinates were also needed in the professional software. These images were then adjusted and aligned using an image feature recognition algorithm, which automatically detects and matches unique image feature characteristics (i.e., tie points) in the overlap region among these images [36]. The generation of Orthomosaic and DEM covered the flight region, which needed to be validated using the GVP coordinates. Subsequently, the Orthomosaic of the debris-covered glacier was manually interpreted and clipped from the generated Orthomosaic, including the non-glacier region. In addition, the DEM image of a debris-covered glacier was also generated. The data process was conducted, and the data processing progress is consistent with other studies [9,25,26,33,37], especially the structure-from-motion applied to create Orthophtos and DEMs from the images collected during each UAV survey [38,39]. It is worth noting that the density of 900 kg/m 3 of glacier was used for converting changes in glacier surface elevation to changes in glacier mass balance.

COSI-Corr for Glacier Velocity
COSI-Corr (Co-registration of Optically Sensed Images and Correlation) is a software tool developed to co-register pairs of optical remotely sensed images (aerials or satellites), which is used to look for seismic ground deformations, glacier flows, horizontal changes between multitemporal images, and so on [40,41]. The main tasks performed by COSI-Corr are precise orthorectification, image co-registration, and correlation, and it is specially designed to retrieve the sub-pixel displacements between optical images. This useful tool has been widely used and developed [42]. In recent years, the COSI-Corr has been widely used to compare the displacement of glacier surfaces using multitemporal images [16,17,26,33,43], including aerial, Landsat, and MODIS images.
To assess the displacement of debris-covered glacier surveyed repeatedly, the two Orthomosaics of glacier surfaces were clipped to the same spatial range with a resolution of 0.10 m. First, the post-event and pre-event images were input to the COSI-Corr tool in order to compare the displacement image. The calculating step in the X and Y directions was set as five in the correlation procession, respectively. The output image then included three bands with a resolution of 0.50 m, respectively, i.e., the East/West band, North/South band, and SNR (Signal-to-Noise Ratio) band. Note that the value of SNR denotes the possibility of glacier displacement, ranging from 0 to 1 [44]. The closer the SNR was to 1, the larger the displacement of the glacier surface. On the contrary, the ground surface had almost displacement when the SNR was lower than the threshold value. As shown in Figure 4, the number of SNRs at different values was counted. Second, an SNR threshold of 0.956 was chosen to deal with some abnormal values using the function "Discard/Replace Image Values" in Cosi-Corr, such as the missing value (NaN). Filtered bands of the East/West band and North/South band were also needed to perform the Non-Local Means Filter (NLM), as shown in the workflow process ( Figure 3). Third, the NLM algorithm filtered data by stepping through the dataset one value at a time. At each datapoint, a small region of surrounding values (typically 5×5 or 7×7; here, 7×7, was adapted) was compared with other non-local values. This approach can remove the repetitive and redundant information in the image (e.g., additive white Gaussian noise) and has the effect of preserving true signals and features [45,46]. Finally, to obtain a relatively reliable result, the standard deviation was applied three times to constrain the boundaries of the glacier surface velocity (i.e., the "Stacking" process in Figure 3), of which the detailed technological process was shown in Figure 3. image co-registration, and correlation, and it is specially designed to retrieve the sub-pixel displacements between optical images. This useful tool has been widely used and developed [42]. In recent years, the COSI-Corr has been widely used to compare the displacement of glacier surfaces using multitemporal images [16,17,26,33,43], including aerial, Landsat, and MODIS images.
To assess the displacement of debris-covered glacier surveyed repeatedly, the two Orthomosaics of glacier surfaces were clipped to the same spatial range with a resolution of 0.10 m. First, the post-event and pre-event images were input to the COSI-Corr tool in order to compare the displacement image. The calculating step in the X and Y directions was set as five in the correlation procession, respectively. The output image then included three bands with a resolution of 0.50 m, respectively, i.e., the East/West band, North/South band, and SNR (Signal-to-Noise Ratio) band. Note that the value of SNR denotes the possibility of glacier displacement, ranging from 0 to 1 [44]. The closer the SNR was to 1, the larger the displacement of the glacier surface. On the contrary, the ground surface had almost displacement when the SNR was lower than the threshold value. As shown in Figure 4, the number of SNRs at different values was counted. Second, an SNR threshold of 0.956 was chosen to deal with some abnormal values using the function "Discard/Replace Image Values" in Cosi-Corr, such as the missing value (NaN). Filtered bands of the East/West band and North/South band were also needed to perform the Non-Local Means Filter (NLM), as shown in the workflow process ( Figure 3). Third, the NLM algorithm filtered data by stepping through the dataset one value at a time. At each datapoint, a small region of surrounding values (typically 5×5 or 7×7; here, 7×7, was adapted) was compared with other non-local values. This approach can remove the repetitive and redundant information in the image (e.g., additive white Gaussian noise) and has the effect of preserving true signals and features [45,46]. Finally, to obtain a relatively reliable result, the standard deviation was applied three times to constrain the boundaries of the glacier surface velocity (i.e., the "Stacking" process in Figure 3), of which the detailed technological process was shown in Figure 3.

Accuracy of DEMs and Orthomosaics
The Phantom 4 Pro coordinates every image it takes using its GPS/GLONASS module, which has been carried out on th alpine glacier of Bernina Massif in the Raethian Alps and has presented high accuracy in DEMs and Orthophotos [33]. To deeply assess the accuracy of generating Orthomoasaics and DEMs, six GVPs and seven GVPs were used to check the imagery accuracy in May and September 2018, respectively. The errors of northing and easting were lower than 11 cm in both surveys, while the vertical error was lower than 12 cm ( Table 2). The total Root Mean Square Error (RMSE) of GVPs in the XY directions (i.e., horizontal direction) was lower than 14 cm in both

Accuracy of DEMs and Orthomosaics
The Phantom 4 Pro coordinates every image it takes using its GPS/GLONASS module, which has been carried out on th alpine glacier of Bernina Massif in the Raethian Alps and has presented high Remote Sens. 2020, 12, 3280 8 of 15 accuracy in DEMs and Orthophotos [33]. To deeply assess the accuracy of generating Orthomoasaics and DEMs, six GVPs and seven GVPs were used to check the imagery accuracy in May and September 2018, respectively. The errors of northing and easting were lower than 11 cm in both surveys, while the vertical error was lower than 12 cm ( Table 2). The total Root Mean Square Error (RMSE) of GVPs in the XY directions (i.e., horizontal direction) was lower than 14 cm in both surveys, and the total RMSE was lower 18 cm. Therefore, the generated DEMs and Ortomosaics showed good accuracy, which met the conditions of computing the images (e.g., glacier surface velocity and elevation changes). In addition, the credibility of results from comparing images in different times is described in the discussion section.

Orthomosaic and Digital Elevation Model
Based on the UAV images of drone flights on 20 May and 22 September 2018, we derived detailed Orthomosaics and DEMs for the debris-covered glacier surface with a resolution of 0.10 m, respectively ( Figure 5). The first monitoring was implemented on 20 May 2018, and the debris-covered area of 0.08 km 2 in front of the glacier was monitored (Figure 5b). The UAV flew three sorties due to the lack of special experience, including flight 1, flight 2, and flight 3 (Figure 5a). The first monitored DEM of glacier surfaces ranged from 4573 m a.s.l. to 4426 m a.s.l. (Figure 5c surveys, and the total RMSE was lower 18 cm. Therefore, the generated DEMs and Ortomosaics showed good accuracy, which met the conditions of computing the images (e.g., glacier surface velocity and elevation changes). In addition, the credibility of results from comparing images in different times is described in the discussion section.

Orthomosaic and Digital Elevation Model
Based on the UAV images of drone flights on 20 May and 22 September 2018, we derived detailed Orthomosaics and DEMs for the debris-covered glacier surface with a resolution of 0.10 m, respectively ( Figure 5). The first monitoring was implemented on 20 May 2018, and the debris-covered area of 0.08 km 2 in front of the glacier was monitored (Figure 5b). The UAV flew three sorties due to the lack of special experience, including flight 1, flight 2, and flight 3 (Figure 5a). The first monitored DEM of glacier surfaces ranged from 4573 m a.s.l. to 4426 m a.s.l. (Figure 5c

Mass Balance of Debris-Covered Glacier
The first monitored glacier front (Figure 5b) was covered by thicker debris. In order to accurately assess the elevation change of the glacier surface, the same regions in the twice-monitored glacier surface were compared (Figure 6a). During the period from 20 May to 22 September, the elevation of the glacier front surface dropped, on average, by 6.58 m ± 3.70 m (arithmetic mean ± standard deviation). The maximum decrease in glacier surface elevation was 32.51 m, while the maximum elevation increase was 9.68 m. We also found that some stronger areas of glacier surface elevation decrease happened at the crevasse edge and in ice cliff regions. The elevation change was used to calculate the mass balance of the glacier front using an ice density of 900 kg/m 3 . The results

Mass Balance of Debris-Covered Glacier
The first monitored glacier front (Figure 5b) was covered by thicker debris. In order to accurately assess the elevation change of the glacier surface, the same regions in the twice-monitored glacier surface were compared (Figure 6a). During the period from 20 May to 22 September, the elevation of the Remote Sens. 2020, 12, 3280 9 of 15 glacier front surface dropped, on average, by 6.58 m ± 3.70 m (arithmetic mean ± standard deviation). The maximum decrease in glacier surface elevation was 32.51 m, while the maximum elevation increase was 9.68 m. We also found that some stronger areas of glacier surface elevation decrease happened at the crevasse edge and in ice cliff regions. The elevation change was used to calculate the mass balance of the glacier front using an ice density of 900 kg/m 3 . The results showed that the mean mass balance of the glacier front was −5.92 m w.e. ± 3.33 m w.e. (arithmetic mean ± standard deviation) (Figure 6b), in which the maximum value of mass balance was 8.77 m w.e., while the minimum value was −29.26 m w.e.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 15 showed that the mean mass balance of the glacier front was −5.92 m w.e. ± 3.33 m w.e. (arithmetic mean ± standard deviation) (Figure 6b), in which the maximum value of mass balance was 8.77 m w.e., while the minimum value was −29.26 m w.e.

Summer Surface Velocity
During the summer, i.e., from 20 May to 22 September 2018, we produced two high-resolution Orthomosaics of the overlap region with the area of 0.08 km 2 in the glacier front. The displacement of the glacier surface in the same region was computed by the Cosi-corr tool in ENVI. As shown in Figure 7a, we found that the glacier mainly moved from west to east, with a mean glacier surface displacement of 18.30 m ± 6.27 m (arithmetic mean ± standard deviation) during the study period. The maximum displacement of the debris-covered surface was 37.04 m, while the minimum displacement was 0.81 m. The glacier displacement in the overlap region showed a significant difference in space. We found that the distribution of the relatively small or large displacements was dispersed in space. Notably, the relatively stronger displacements did not occur in the terminus of the glacier front. For this point, we inferred that the discrepancy in surface displacement was mainly controlled by the local ice surface slope, ice thickness, ice density, and temperature [47].
Here, 125 days between the two UAV monitoring tests in the summer were used to compute the glacier surface velocity per day. The mean daily velocity of the glacier debris-covered surface was 0.14 m/d ± 0.05 m/d (arithmetic mean ± standard deviation) during the summer. The maximum velocity of the glacier surface was 0.30 m/d, while the minimum velocity was 0.01 m/d. In addition, the melt velocity of the glacier debris-covered surface was also assessed (Figure 7b). The mass balance per day ranged from −0.23 m w.e./d to 0.07 m w.e./d, with the mean mass balance of −0.05 m w.e./d ± 0.03 m w.e./d (arithmetic mean ± standard deviation).

Figure 7.
Glacier surface displacement of the debris-covered glacier derived from the two surveys in the glacier front. In plots, (a) denotes the surface displacement distance of the debris-covered glacier

Summer Surface Velocity
During the summer, i.e., from 20 May to 22 September 2018, we produced two high-resolution Orthomosaics of the overlap region with the area of 0.08 km 2 in the glacier front. The displacement of the glacier surface in the same region was computed by the Cosi-corr tool in ENVI. As shown in Figure 7a, we found that the glacier mainly moved from west to east, with a mean glacier surface displacement of 18.30 m ± 6.27 m (arithmetic mean ± standard deviation) during the study period. The maximum displacement of the debris-covered surface was 37.04 m, while the minimum displacement was 0.81 m. The glacier displacement in the overlap region showed a significant difference in space. We found that the distribution of the relatively small or large displacements was dispersed in space. Notably, the relatively stronger displacements did not occur in the terminus of the glacier front. For this point, we inferred that the discrepancy in surface displacement was mainly controlled by the local ice surface slope, ice thickness, ice density, and temperature [47].
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 15 showed that the mean mass balance of the glacier front was −5.92 m w.e. ± 3.33 m w.e. (arithmetic mean ± standard deviation) (Figure 6b), in which the maximum value of mass balance was 8.77 m w.e., while the minimum value was −29.26 m w.e.

Summer Surface Velocity
During the summer, i.e., from 20 May to 22 September 2018, we produced two high-resolution Orthomosaics of the overlap region with the area of 0.08 km 2 in the glacier front. The displacement of the glacier surface in the same region was computed by the Cosi-corr tool in ENVI. As shown in Figure 7a, we found that the glacier mainly moved from west to east, with a mean glacier surface displacement of 18.30 m ± 6.27 m (arithmetic mean ± standard deviation) during the study period. The maximum displacement of the debris-covered surface was 37.04 m, while the minimum displacement was 0.81 m. The glacier displacement in the overlap region showed a significant difference in space. We found that the distribution of the relatively small or large displacements was dispersed in space. Notably, the relatively stronger displacements did not occur in the terminus of the glacier front. For this point, we inferred that the discrepancy in surface displacement was mainly controlled by the local ice surface slope, ice thickness, ice density, and temperature [47].
Here, 125 days between the two UAV monitoring tests in the summer were used to compute the glacier surface velocity per day. The mean daily velocity of the glacier debris-covered surface was 0.14 m/d ± 0.05 m/d (arithmetic mean ± standard deviation) during the summer. The maximum velocity of the glacier surface was 0.30 m/d, while the minimum velocity was 0.01 m/d. In addition, the melt velocity of the glacier debris-covered surface was also assessed (Figure 7b). The mass balance per day ranged from −0.23 m w.e./d to 0.07 m w.e./d, with the mean mass balance of −0.05 m w.e./d ± 0.03 m w.e./d (arithmetic mean ± standard deviation).  Here, 125 days between the two UAV monitoring tests in the summer were used to compute the glacier surface velocity per day. The mean daily velocity of the glacier debris-covered surface was 0.14 m/d ± 0.05 m/d (arithmetic mean ± standard deviation) during the summer. The maximum velocity of the glacier surface was 0.30 m/d, while the minimum velocity was 0.01 m/d. In addition, the melt velocity of the glacier debris-covered surface was also assessed (Figure 7b). The mass balance per day ranged from −0.23 m w.e./d to 0.07 m w.e./d, with the mean mass balance of −0.05 m w.e./d ± 0.03 m w.e./d (arithmetic mean ± standard deviation).

Accuracy of DEMs and Orthomosaics
Using GCPs can improve the accuracy of generating DEM and Orthomosaics from aerial mapping surveys. The ground control points can be measured by some methods, i.e., traditional surveying methods, LiDAR, an existing map, and even Google Earth [48][49][50]. Many researchers think that adding more ground control points takes more time but creates more accurate drone maps, but does it? To test how many GCPs are really required for accurate results, the Nevada Department of Transportation (DOT) set out a series of tests. A DJI Phantom 4 Pro drone and Pix4D mapper software were used during the experiment and a covered area of~0.14 km 2 (i.e., 34 acres) was measured. The results showed that 11-13 GCPs would likely be the most efficient number for geometry. The Northing and Easting RMSEs were 1.52 cm (i.e., 0.05 feet) when the number of GCPs ranged from 18 to 4 GCPs. The vertical RMSE was about 1.52 cm (i.e., 0.05 feet) from 18 to 14 GCPs and consistently around 3.05 cm (i.e., 0.1 feet) from 10 to 4 GCPs. In addition, reducing the number of GCPs by four to a total of seven located on the outside of the study region led to a maximum vertical difference of 13.41 cm (i.e., 0.44 feet) and vertical RMSE of 7.31 cm (i.e., 0.24 feet). A more detailed description of this assessment process was presented in the official document [32]. Herein, the number of GCPs in our study met the condition of generating the high-resolution images. The total RMSE of GCPs was lower than 20 cm, which was slightly higher than DOT experiment because the topography of the glacierized region was more complicated.

Changes in the Elevation and Mass Balance
It was widely accepted that the elevation change of the glacier surface is translated into the change of glacier mass balance to better understand the glacier melting. Here, we justly discussed the mass balance with respect to understanding glacier surface change. The annual mass balance of the BRG1 ranged from −1.94 to 2.26 m w.e. during the period 1952-2017 according to the glaciology measurement [22,51]. Based on the observation at 4700 m a.s.l. from 2008 to 2013, the mass balance of the BRG1 ranged from −3.64 m w.e. to 3.10 m w.e. from June to September, especially the mass balances of −3.53 m w.e. and −3.64 m w.e. in the melt seasons of 2008-2009 and 2009-2010, respectively [51]. The mean mass balance in situ at 4700 m w.e. was around −3.22 m w.e. in the melting seasons from 2008 to 2013, and the maximum and minimum mass loss were −3.64 m w.e. and −2.63 m w.e., respectively. In addition, the seasonal fluctuation of the glacier mass balance gradient, which changed with elevation, was lower than that of the annual change. For example, the mean mass balance gradient from July to September 2013 was 0.72 m w.e./(100m) from 4600 m a.s.l. to 4900 m a.s.l, which was significantly lower than 1.45 m w.e./(100m) from 2012 to 2013. In this paper, the mean glacier surface elevation of our repeatedly monitored glacier area was around 4551 m a.s.l. (Figure 5c). We deduced that the mean mass balance of the melting season ranged from −4.71 m w.e. to −3.70 m w.e. on the debris-covered glacier using the mass balance and mass balance gradient of the melting season during the period 2008-2013. Correspondingly, the mean mass balance in the melting season was around −4.29 m w.e. during that timespan.
In addition, the debris-covered glaciers in the Himalaya, being similar to glacier in YSM in terms of glacier surface feature, may have a spatially averaged rate of surface height change that is similar to those observed on bare-ice glaciers, despite the protected effects of thick debris. Those surfaces of debris-covered glaciers usually develop many ice cliffs. However, the mechanisms controlling the formation and survival of cliffs remain largely unknown [52]. An assessment of understanding the distribution and characteristics of these surface features on the debris-covered glacier was lunched in the Lantang Glacier located in the Himalaya, using UAV to obtain high-resolution Orthomosaics and elevation models [37]. The results indicated that the cliffs play an important role in the mass loss of debris-covered glaciers, as they melt comparatively faster than others. The mass loss from ice cliffs has hardly been considered using the traditional glacier mass balance method. In other words, the mean mass balance of −4.29 m w.e., calculated during the period of the melt seasons at the front of BRG1 in YSM from 2008 to 2013, was lower than that of the actual melt. That is, the actual mass balance of melt season should have been lower than −4.29 m w.e. Therefore, the computed result of mass balance (−5.92 m w.e. ± 3.33 m w.e.) from UAV surveying in terms of the melt season was credible and close to the actual situation. After all, the air temperature was warming from 2013 to 2018 [22].

Glacier Surface Velocity
To understand the glacier surface velocity, a campaign monitoring the displacement of 16 mass balance stakes distributed the regions from 4800 m a.s.l. to 4600 m a.s.l. of the BRG1 was carried out during the period of July-October 2016 [53]. The Trimble GeoXT of manual dGPS was used in the campaign and its accuracy was within 0.10 m. The results showed that the monthly glacier surface velocity was at 2.28 m/month in terms of the entire glacier. The mean daily velocity of 16 stakes was 0.11 m/day during the melt season, while the mean daily velocity of the debris-covered glacier was 0.14 m/day ± 0.05 m/day using UAV repeat photogrammetry. Similar studies also found that the surface velocity at the end of the glacier was faster than that of the glacier surface with higher elevation [54][55][56]. For example, the observation in the Yanzigou glacier in the Gongga Mountains, located in the northeast region of the YSM, showed that the change of glacier surface velocity on the ablation area was more significant than that at the end of the glacier [56]. Besides, we also found that the glacier surface with the highest velocity was not the same as the glacier region with the most loss of glacier mass balance (Figure 7b). Overall, the surface velocity of the debris-covered glacier computed using the high-resolution images from UAV photogrammetry showed high accuracy and was very close to the actual velocity.

Advantages and Disadvantages of the UAV Monitoring
In the past decades, optical Very High Resolution (VHR) and Synthetic Aperture Radar (SAR) imagery have been rapidly produced and have represented effective tools in the glacier monitoring [57]. Although these spaceborne platforms had short revisiting times (i.e., ranging from one day to a few days), the time and spatial resolutions of satellite imagery still did not match the means of UAVs. Compared with spaceborne platforms, the advantages of UAV monitoring are significant. For example, the UAV photogrammetry was applied for to survey, monitor, and obtain characteristics of glacier surface elevation, surface velocity, supra-glacial lakes, crevasses, and unstable cliffs [57,58]. The mass loss from ice avalanches in glacier fronts had been ignored, and the phenomenon was not revealed until the application of UAV [37]. Meanwhile, the disadvantages of UAV monitoring also cannot be ignored [57]. For example, a fixed-wing platform needs a proper takeoff and landing sites, which is sometimes difficult to locate in rugged mountainous terrain. The length of UAV flight has been relatively short. The longest flight duration of a multirotor drone cannot currently exceed one hour. In addition, power consumption is faster in the glacier region than that in the plain region due to the low-temperature environment. The area of UAV monitoring is also limited. To make sure the high-resolution and high-overlap rate for images, the UAV has to fly at a special mid-and low-altitude over the glacier surface. Thus, future development of UAV monitoring should focus on these problems. Increased battery capacity is important to address short aircraft endurance. However, it is more important to focus on the bridge function of UAV monitoring between the ground observation and satellite technology of glacier monitoring.

Conclusions
UAVs are an important tool with respect to the study glacier melt and dynamic in ablation zones. Here, two UAV experiments were conducted to measure the glacier topography for the debris-covered glacier of the BRG1 on the YSM. Using UAV flight photogrammetry, we obtained the high-resolution Orthomosaic and DEM images with a resolution of 0.10 m. The images provided detailed information of the unobservable glacier surface. By comparing the different time images, we computed the dynamic features of debris-covered glacier on the ablation zones in the summer. During the period from 20 May to 22 September in 2018, the elevation changes of the debris-covered glacier ranged from −32.51 m to 9.68 m, while the mean mass balance decreased by 6.58 m ± 3.70 m. Correspondingly, the mean mass balance of the glacier front was −5.92 m w.e. ± 3.33 m w.e. In addition, during the summer, we found that debris-covered glacier surface displacement ranged from 37.04 m to 0.81 m, with a mean displacement of 18.30 m ± 6.27 m. The mean daily velocity was 0.14 m/d ± 0.05 m/d. The results provide the new insights to understand glacier dynamic features for the debris-covered glacier, which reveal the accuracy mass balance and glacier displacement and also provide validation for glacier modeling and remote-sensing assessments.