Dynamic Monitoring of Laohugou Glacier No. 12 with a Drone, West Qilian Mountains, West China

: Laohugou glacier No. 12 (LHG12), located in the northeast of the Qinghai–Tibet Plateau, is the largest valley glacier in the Qilian mountains. Since 1957, LHG12 has shrunk signiﬁcantly. Due to the limitations of in situ observations, simulations and investigations of LHG12 have higher levels of uncertainty. In this study, consumer-level, low-altitude microdrones were used to conduct repeated photogrammetry at the lower part of LHG12, and a digital orthophoto map (DOM) and a digital surface model (DSM) with a resolution at the centimeter scale were generated, from 2017 to 2021. The dynamic parameters of the glacier were detected by artiﬁcial and automatic extraction methods. Using a combination of GNSS and drone-based data, the dynamic process of LHG12 was analyzed. The results show that the terminus of LHG12 has retreated by 194.35 m in total and by 19.44 m a − 1 on average during 2008–2021. The differential ablation leading to terminus retreat distance markedly increased during the study period. In 2019–2021, the maximum annual surface velocity was 6.50 cm day − 1 , and during ablation season, the maximum surface velocity was 13.59 cm day − 1 , 52.17% higher than it is annually. The surface parameters, motion, and mass balance characteristics of the glacier had signiﬁcant differences between the west and east branches. The movement in the west branch is faster than it is in the east branch. Because of the extrusion of the two ice ﬂows, there is a region with a faster surface velocity at the ablation area. The ice thickness of LHG12 is decreasing due to intensiﬁed ablation, leading to a deceleration in the surface velocity. In large glaciers, this phenomenon is more obvious than it is in small glaciers in the Qilian mountains.


Introduction
As one of the components of Earth's climatic system, the cryosphere is the most sensitive to climate change and has a significant impact on ecosystems, the environment, and sustainable social and economic development [1][2][3]. As the most important component of the cryosphere, glaciers and ice sheets store about 70% of the freshwater resources in the world [4]. Glacier changes are a key point in glaciology research and also represent an important factor in the climate, water resources, and sea level changes [5][6][7][8][9]. Glaciers are the most important freshwater resource in the arid regions of Northwest China [10]; at the same time, glaciers are not only an important driving force but also record and act as indicators of climate change [11,12]. Due to global warming, glaciers around the world have shrunk significantly, and this is irreversible [13,14]. Over the past few decades, the area, length, and mass balance (MB) of global mountain glaciers have been reduced consistently, something that has been detected by both satellite and ground observations; however, inter-annual and regional differences can be observed [15].
The Qilian mountains are located in the northeast region of the Qinghai-Tibet Plateau and are the natural boundary of the Qinghai and the Gansu provinces. They provide a source of fresh water via the inland river, especially the Hexi corridor. They are also an important ecological barrier of the northwest region of China. In recent years, the glaciers in the Qilian Mountains have shown massive volume losses, with the contribution rate of glacier melt runoff being over the critical point [16,17]. In the Qilian mountains, there are 2683 glaciers that are 1597.81 km 2 in area and 84.48 km 3 in volume [18]. Previous research used ground survey or remote sensing methods to describe how glaciers in the Qilian mountains have shrunk significantly; however, because of the limitations of satellite remote sensing data, there have been some differences in terms of the study period and results, and many studies have focused on changes in glacier area [19][20][21][22][23]. Investigations of glacier volume or surface elevation changes acquired using methods, such as volume-area scaling or DSM difference estimations, still had high levels of uncertainty [24][25][26][27].
Satellite remote sensing data with global coverage can effectively monitor glacial terrain. However, temporal and ground resolutions, the cost of high-resolution data with large scales, and the influence of cloud or snow cover limit the application of satellite remote sensing in glaciology research. Drone-based monitoring can easily overcome those limitations. The potential for using drones in glaciology is huge, and drones could play a great supporting role in field-based and satellite remote sensing-based glaciology observations [28][29][30][31][32][33]. Glaciers located in high mountainous areas are always under harsh conditions, such as having a hostile climate or poor approachability, and this is a huge limitation in field observations along with inadequate labor and funding, creating great difficulties in the dynamic research on layers. Observations of Laohugou Glacier No. 12 (LHG12) are more comprehensive than others in the Qilian Mountains. Observational parameters, such as mass balance and surface velocity, among others, were investigated by glaciological methods, but because of low accessibility in the ice serac region, the mid-low part of LHG12, observation data are always interrupted [34,35]. Due to the limitations of in situ observation, the investigation of LHG12 in the Qilian mountains with a thermomechanically coupled ice flow model experienced errors at the terminus region [36]. In this study, consumer-level, low-altitude micro drones were used to conduct repeated photogrammetry at the lower part of LHG12, and a digital orthophoto map (DOM) and a digital surface model (DSM) with centimeter-scale resolution were generated. LHG12 s dynamic process was analyzed; we found that there have been significant differences in the dynamic processes between the east and west branches, and the movement in the west branch is faster than it is in the east branch, even below the confluence area. In the previous monitoring and simulation studies, the differences in dynamic processes between the east and west branch were ignored. The aim of this paper is that provide accurate data support for glacier dynamic simulations for further research.

Study Area
LHG12 (glacier code: 5Y448D0012) is located in the northeast region of the Qinghai-Tibetan Plateau (Figure 1a). It is 9.7 km in length, 20.37 km 2 in area, and is the largest valley glacier in the Qilian mountains [37]. LHG12 has a typical continental climate that is characterized and influenced by westerlies all year round. It is also located at 4200 m a.s.l., has a daily mean air temperature in the summer (from 1 June to 30 September) above 0 • C, and precipitation that mainly occurs from May to September. The average annual temperatures of the air and glacial surface are −2.6 • C and −6.4 • C, respectively, with 317 mm w.e. of annual precipitation. The maximum monthly precipitation was 142 mm w.e. in July, accounting for 36% of the annual precipitation in 2019 [38]. Due to increases in air temperature, especially in autumn and winter, the ice temperature and climate sensitivity of LHG12 have increased, and LHG12 has shrunk significantly since 1957 [39]. Groundpenetrating radar (GPR) data show that the bed topography of the east branch is generally more rugged than that of west branch, while the surface slopes of both tributaries are Remote Sens. 2022, 14, 3315 3 of 17 gentle; in the upper confluence area, the ice thicknesses of the east branch and west branch were 122 m and 157 m, respectively. Furthermore, due to transverse compression and convergence from the two ice flows, the ice depth is 162 m at the center of the confluence area [40].
July, accounting for 36% of the annual precipitation in 2019 [38]. Due to increases in air temperature, especially in autumn and winter, the ice temperature and climate sensitivity of LHG12 have increased, and LHG12 has shrunk significantly since 1957 [39]. Groundpenetrating radar (GPR) data show that the bed topography of the east branch is generally more rugged than that of west branch, while the surface slopes of both tributaries are gentle; in the upper confluence area, the ice thicknesses of the east branch and west branch were 122 m and 157 m, respectively. Furthermore, due to transverse compression and convergence from the two ice flows, the ice depth is 162 m at the center of the confluence area [40].    [41] and Mavic 2 Pro (M2P) [42] produced by the DJI Innovation Company were used as the image platform, and Pix4D Capture (Android edition) was used as the route-planning software. Drone observations of LHG12 started in 2017. At first, we used P4P to test the flight control, time, and altitude, and the battery capability had also been tested. To increase efficiency and reduce discharge from the lithium-ion batteries in high-altitude and low-temperature alpine areas and ensure flight safety, each flight time was not more than 13 min. At the terminus region, the effects of the ground control points (GCPs) with different numbers and distributions were tested to determine the accuracy of the drone-based photogrammetry, and an observation scheme suitable for drone-based photogrammetry at the mid-lower regions of LHG12 was determined [43]. In this study, five GCPs were deployed along the glacier's mid-line (Figure 1c), and except for in 2017, each drone survey included 12 routes and captured about 1500 images. A total 6 drone surveys were conducted. The overlap was set to 75% in Pix4D Capture, and the wind speed was no more than a moderate breeze at every site. An overview of the flights is shown in Table 1. Images covered about 2 km 2 ( Figure 1c). The drone flew according to the planned routes, and captured the images automatically; meanwhile, the GCPs deployed on the surface of LHG12 in the survey region were measured by the Globe Navigation Satellite System (GNSS) receivers. The measurement accuracies of SOUTH S86 and SOUTH INNO7 GNSS receivers produced by the SOUTH [44] set in real-time kinematic (RTK) mode were ±3 cm vertically and ± 1.5 cm horizontally. GNSS was used with the RTK mode to measure the terminus position of LHG12 in the period from 2008-2013 [37] and 15 sticks (along the mid-ice flow, 6 sticks were set at the west branch surface, and 9 sticks were set at the east branch surface) that has been placed in the ice were also measured as surface velocity markers. The ground sample distance (GSD) of the drone-based photogrammetry was dependent on flight height (H) as well as the sensor parameters and could be estimated by Functions (1) and (2) [45]. Some lens manufacturers provide the focal length (F 35 ) as the 35 mm equivalent. It is not the 35 mm equivalent, but instead the real focal length that should be used in Pix4Dmapper. In order to find the real focal length, some computations are needed. In the case of a 4:3 ratio system, the formula for the real focal length FR is: F 35 is the focal length that corresponds to the 35 mm equivalent; FR is real focal length (mm); SW is the real sensor width. The relationship between H (m) and GSD (cm/pixel) is: where imW is the image width (pixel).

Image Processing
Five GCPs and all of the aerial photos were imported to the Pix4D Mapper software, and each GCP was marked in at least three photos. Finally, the DSM and DOM with a 5 cm GSD were acquired. The projection of the DSM and DOM was UTM 47N, and the ellipsoid of projection was WGS84. We used QGIS [46] software to process the DOM and DEM, and the glacier extent and terminus position were extracted by the DOM acquired in 2018, 2019, and 2021. For 2017, the DOM was only extracted in the terminus position. Slope, aspect, and roughness were extracted by DSM acquired in 2018, 2019, and 2021. Images were clipped by the same range layer in the survey region.

Glacier Surface Displacement Extraction
The ImGRAFT software package was used to extract the LHG12 surface displacement from July 2018 to September 2021. ImGRAFT is an open-source package that encompasses georeferencing, georectification, and terrestrial feature tracking [47]. In order to determine displacement, we used the templatematch function with the standard normalized crosscorrelation (NCC) algorithm in ImGRAFT. NCC was outperformed in areas with poor visual contrast and in areas with thin clouds or changing snow conditions from one image to the next. Different illumination conditions in two images can lead to mismatches [48]. During the flight missions, the illumination conditions inevitably changed, leading to an inconsistent DOM gray scale ( Figure 2c). There were also huge displacement extraction outliers ( Figure 2a). To avoid this issue, hill-shade images generated by DSM were used as input. Those hill-shade images highlighted LHG12 s surface texture while avoiding different illumination conditions at the same time. However, there were still some outliers in the displacement-extracted results ( Figure 2b). In those areas, the glacier surface mostly consisted of bare ice and lower roughness, leading to more outliers. Therefore, the search window size had to be big enough to include surface texture features [49]. Considering the computing capability, the images were re-sampled to a 1 m resolution, and the search window size was set to 200, while the relative distance on the ice surface was 200 m.

Meteorologic Data
In 2009, an automatic meteorological station (AWS) was deployed at 4200 m a.s.l. beyond the research station which is 1.5 km from terminus of LHG12. The AWS had air temperature (1.5 m from ground) and precipitation (T-200B) sensors connected to a data logger (CR-1000), and air temperature and precipitation data have been recorded continuously since August 2009.

Uncertainty Estimation
In order to estimate of the glacier surface displacement extraction uncertainty, 100 ground validation points (GVPs) were selected (Figure 2b), displacement was extracted by ImGRAFT, and the visual feature track was contrasted. Image pairs with a certain time

Meteorologic Data
In 2009, an automatic meteorological station (AWS) was deployed at 4200 m a.s.l. beyond the research station which is 1.5 km from terminus of LHG12. The AWS had air temperature (1.5 m from ground) and precipitation (T-200B) sensors connected to a data logger (CR-1000), and air temperature and precipitation data have been recorded continuously since August 2009.

Uncertainty Estimation
In order to estimate of the glacier surface displacement extraction uncertainty, 100 ground validation points (GVPs) were selected (Figure 2b), displacement was extracted by Im-GRAFT, and the visual feature track was contrasted. Image pairs with a certain time interval contained information of glacial surface movement, the top of the ice serac was selected as the GVP, the coordinates of the same GVP on both of the image pairs were measured, and the displacement was calculated. The correlation between ImGRAFT and the visual feature track was relevant: R 2 was 0.96 (Figure 2e), and the median was slightly positive (Figure 2d). The root mean squared error (RMSE) was used to describe the displacement extraction error, and the RMSE was ±1.69 m. Using the combined error of the GNSS-RTK measurement system, the glacier surface displacement extraction error (Error SV ) could be calculated by Function (3): E Image is the glacier surface displacement extraction error and is equal to the RMSE; E GNSS is the error of the GNSS-RTK measurement system and we determined that Error SV was ±1.69 m in this study.

Changes in Terminus & Area
The terminus of LHG12 was briefly stable in the 1980s, and after that, it began to retreat significantly (Figure 3a) [37].  Table 2). In ablation season, the ice surface was eroded by strong fluviation, a combination of surface moraine distribution and ice flow, which caused evident differential ablation at the mid-lower part of LHG12. At the terminus of LHG12, some regions with weak ablation remained for one to two years, turning the edge into an irregular shape similar to a "bulge", and if this part of the ice ablated completely or separated from glacier within the next one or two years, it caused the terminus to retreat by a markedly increased distance. From 2008 to 2021, rapid terminus retreat appeared two times, in 2010 and 2015, respectively (Figure 3b), an interval of approximately five to six years. A drone image from 2021 showed that there was a new "bulge" caused by differential ablation, which might cause "rapid retreat" to occur again in the next one to two years.

Surface Morphology
The main orientation of LHG12 was northwest, and due to differential ablation, the aspect had variation in the ice serac region (Figure 4a). The effective range of the glacial surface roughness was zero to four. In the monitoring area, the ice serac region was the roughest region, and mid-upper parts were smooth. The roughness of the east branch was higher than that of the west branch, especially near the shear line. The surface roughness of the monitoring area increased slightly in 2018-2021 (Figure 4b). The effective range of the slope was 0-45 • . In 2018-2021, the average surface slope of LHG12 was 13.86-17.02 • . The glacier's slope increased significantly, and the rate of slope increase was higher in June-July than in July-September, due to strong ablation in the monitoring area (Figure 4c).

Surface Velocities
The main flow direction of LHG12 was towards the north. The confluence of the east and west ice flow had a huge influence on the ice flow direction. Above the confluence area, the azimuth of east ice flow was about 313 • , the azimuth of west ice flow was about 50 • , and the confluence of the two ice flows formed a shear line along the valley. Near the shear line, the flow direction of the west branch has turned 33.99 degrees to the northwest, and the flow direction of the east branch has turned 54.8 degrees to the north ( Figure 5).

Surface Velocities
The main flow direction of LHG12 was towards the north. The confluence of the east and west ice flow had a huge influence on the ice flow direction. Above the confluence area, the azimuth of east ice flow was about 313°, the azimuth of west ice flow was about 50°, and the confluence of the two ice flows formed a shear line along the valley. Near the shear line, the flow direction of the west branch has turned 33.99 degrees to the northwest, and the flow direction of the east branch has turned 54.8 degrees to the north ( Figure 5). The glacier's surface velocity was obviously influenced by extrusion of the two ice flows ( Figure 6). From the surface velocity extracted by three drone-based measurements in 2021, it can be seen that: the maximum annual surface velocity was 6.50 cm day −1 (Figure The glacier's surface velocity was obviously influenced by extrusion of the two ice flows ( Figure 6). From the surface velocity extracted by three drone-based measurements in 2021, it can be seen that: the maximum annual surface velocity was 6.50 cm day −1 (Figure 6b); the surface velocity accelerated during the ablation season, and the maximum was 13.59 cm day −1 (Figure 6a); the surface velocity in the summer was doubled compared to the annual, and decelerated at the end of ablation season, with a maximum surface velocity of 9.60 cm day −1 . The surface velocity that was higher in the confluence area than in other areas (Figure 6c). The lower portions of the east and west branches and the middle of the confluence region were the areas with faster movement. On both sides of the shearing line, the surface velocity was different.
At the glacier surface, six cross sections (CS) perpendicular to the main flow direction and two longitude sections (LS) along the main flow directions of the east and west branches were selected, and the surface velocities of each period and altitude on 6 September 2021 were extracted according to those parameters (Figure 7). Because of factors, such as topography, there were different characteristics of surface velocity included in those parameters, but a higher surface velocity than the annual average during ablation season was a common feature of all of them. At CS A, near the terminus, the ice was thin, the surface velocity was slow, and no obvious differences were observed; at CS B, the surface velocity was slow at the edges of LHG12 and faster along the west side; at CS C, located below the confluence region, the ice flow direction changed according to the terrain. The surface velocity was higher in the middle of CS C and lower along both sides of CS C. At the east side of the glacier, the ice was squeezed by mountain, the elevation of the ice surface was increased, and the surface velocity was significantly lower; at CS D, located in the mid-upper confluence region, the west branch was squeezed by the mountain, the flow direction changed, the elevation of the western ice surface was increased, and the surface velocity was lower in the west and higher in middle and east in the summer; at CS W-E, located in the lower part of the west branch, the western elevation of the ice surface was higher than the eastern branch, but the trend in the surface velocity was the opposite; at CS E-E, located in the lower east branch, the eastern elevation of the ice surface was higher than the western elevation, and the surface velocity had the same trend. The elevation trend was gentle at LS W but rolling at LS E, and the annual surface velocities at LS W and LS E increased with altitude and fluctuated in the summer. The trends in the seasonal surface velocity and air temperature of LHG12 were more consistent (Figure 8). Movement accelerated significantly during the ablation season: it was higher than the annual speed but slowed down during the accumulation season. Because of the huge difficulty of making basal sliding observations, we could not analyze the relationship between melt water and basal sliding for LHG12 directly.
Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 18 6b); the surface velocity accelerated during the ablation season, and the maximum was 13.59 cm day −1 (Figure 6a); the surface velocity in the summer was doubled compared to the annual, and decelerated at the end of ablation season, with a maximum surface velocity of 9.60 cm day −1 . The surface velocity that was higher in the confluence area than in other areas (Figure 6c). The lower portions of the east and west branches and the middle of the confluence region were the areas with faster movement. On both sides of the shearing line, the surface velocity was different.  trend. The elevation trend was gentle at LS W but rolling at LS E, and the annual surface velocities at LS W and LS E increased with altitude and fluctuated in the summer. The trends in the seasonal surface velocity and air temperature of LHG12 were more consistent (Figure 8). Movement accelerated significantly during the ablation season: it was higher than the annual speed but slowed down during the accumulation season. Because of the huge difficulty of making basal sliding observations, we could not analyze the relationship between melt water and basal sliding for LHG12 directly.   (Figure 9). The surface velocity of the east branch decreased, but it did not slow down as obviously in the west branch [50].  1960-1961, 2008-2009, 2010-2011, and 2011-2013, respectively (Figure 9). The surface velocity of the east branch decreased, but it did not slow down as obviously in the west branch [50]. measurement points were mostly distributed along the eastern branch of LHG12. The maximum surface velocities of the east branch at 4750-4800 m a.s.l. were 36 m a −1 , 32.4 m a −1 , 26.4 m a −1 , and 26.6 m a −1 , in 2008-2009, 1960-1961, 2010-2011, and 2011-2013, respec tively; the west branch observation points on the upper part of confluence region were located at about 4650 m a.s.l. The annual surface velocities were closer to or greater than those of the east branch; 35.6 m a −1 , 32.6 m a −1 , 34.7 m a −1 , and 39.0 m a −1 , in 1960-1961 2008-2009, 2010-2011, and 2011-2013, respectively (Figure 9). The surface velocity of the east branch decreased, but it did not slow down as obviously in the west branch [50].

Changes of Surface Elevation
In 2018-2019, ice surface elevation of the survey area decreased by an average of 1.22 m. In 2019-2021, it by an average decreased 3.28 m, decreasing by about 1.65 m per yea ( Figure 10). The characteristics of surface elevation difference results were quite differen between east and west branches. In the monitoring area, the area of western branch wa larger than that of the eastern branch. The regions in which the elevation increased were   Figure 10). The characteristics of surface elevation difference results were quite different between east and west branches. In the monitoring area, the area of western branch was larger than that of the eastern branch. The regions in which the elevation increased were located on the upper part of the confluence area of the western branch and the lower part of the confluence area of eastern branch. In these areas, the morphology of the eastern branch was more complex. After the east and west branches merged with each other, the eastern branch narrowed rapidly, especially in the lower part of confluence area, and the surface elevation of eastern branch increased significantly under the force of the bedrock [40]. The cross profile of the glacier is narrowed by the bedrock. The ice flow needs to get over the bedrock; meanwhile, the ice that comes from upstream has not ablated completely, which led to decelerate the surface velocity and increase the surface elevation.

Changes of Surface Elevation
Remote Sens. 2022, 14, x FOR PEER REVIEW 13 of 18 located on the upper part of the confluence area of the western branch and the lower part of the confluence area of eastern branch. In these areas, the morphology of the eastern branch was more complex. After the east and west branches merged with each other, the eastern branch narrowed rapidly, especially in the lower part of confluence area, and the surface elevation of eastern branch increased significantly under the force of the bedrock [40]. The cross profile of the glacier is narrowed by the bedrock. The ice flow needs to get over the bedrock; meanwhile, the ice that comes from upstream has not ablated completely, which led to decelerate the surface velocity and increase the surface elevation.

Climate Background
AWS data from Jan 2010 to Dec 2020 were analyzed in this study ( Figure 11). We calculated the monthly negative and positive cumulative, minimum, maximum, and average air temperature and divided the solid and liquid precipitation by the air temperature (3.5 °C) [51]. In 2010-2021, the positive cumulative temperature increased by 0.02 °C month −1 , and the negative cumulative temperature decreased by 0.14 °C month −1 . This is

Climate Background
AWS data from Jan 2010 to Dec 2020 were analyzed in this study ( Figure 11). We calculated the monthly negative and positive cumulative, minimum, maximum, and average air temperature and divided the solid and liquid precipitation by the air temperature (3.5 • C) [51]. In 2010-2021, the positive cumulative temperature increased by 0.02 • C month −1 , and the negative cumulative temperature decreased by 0.14 • C month −1 . This is higher than the positive cumulative temperature, as the amplitude, minimum, maximum, and average air temperature increased by 0.011 • C month −1 , 0.004 • C month −1 , and 0.005

Discussion
Increased air temperature, especially in the autumn and winter, could cause an increase in the borehole ice temperature and make the glacier more sensitive to climate change [39]. The AWS data show that in 2010-2021, the negative cumulative temperature had decreased and the positive cumulative temperature had increased. The variation amplitude in the negative cumulative temperature was obviously higher than it was in the positive cumulative temperature, which means that there was a serious decrease in the "cold storage" of LHG12. In 2010, 2016, and 2017, the positive cumulative temperature was higher, the negative cumulative temperature was lower, the liquid precipitation percentage was higher, and LHG12 experienced strong surface melt. Differential ablation on the surface of LHG12 may be exacerbated by supraglacial runoff and streams. In those years, there was strong ablation for a long time; therefore, a "bulge" gradually formed at the terminus of LHG12 over the next 1-2 years. This may be the main reason for the rapid retreat observed in the terminus in every few years.
In the Qilian Mountains, glaciers gradually become smaller, and the glacier surface velocity becomes slower from east to west (Figure 12). In recent decades, due to the climate warming, the glaciers in the Qilian Mountains have experienced huge mass loss, reduced ice thickness, and a reduced surface velocity. The surface velocity reduction trend

Discussion
Increased air temperature, especially in the autumn and winter, could cause an increase in the borehole ice temperature and make the glacier more sensitive to climate change [39]. The AWS data show that in 2010-2021, the negative cumulative temperature had decreased and the positive cumulative temperature had increased. The variation amplitude in the negative cumulative temperature was obviously higher than it was in the positive cumulative temperature, which means that there was a serious decrease in the "cold storage" of LHG12. In 2010, 2016, and 2017, the positive cumulative temperature was higher, the negative cumulative temperature was lower, the liquid precipitation percentage was higher, and LHG12 experienced strong surface melt. Differential ablation on the surface of LHG12 may be exacerbated by supraglacial runoff and streams. In those years, there was strong ablation for a long time; therefore, a "bulge" gradually formed at the terminus of LHG12 over the next 1-2 years. This may be the main reason for the rapid retreat observed in the terminus in every few years.
In the Qilian Mountains, glaciers gradually become smaller, and the glacier surface velocity becomes slower from east to west ( Figure 12). In recent decades, due to the climate warming, the glaciers in the Qilian Mountains have experienced huge mass loss, reduced ice thickness, and a reduced surface velocity. The surface velocity reduction trend observed in the large glaciers is higher than the trend observed in the small glaciers in the Qilian Mountains [52]. The dynamics of LHG12 were more complex than those of other single cirque or v ley glaciers. Below the confluence region and between the two ice flows, there was a shea ing surface similar to a "deep ditch" on the surface of LHG12 (Figure 6d). On both sid of the shearing line, the glacier surface parameters, motion, and mass balance characte istics represented significant differences. In previous simulations, due to the limitatio of observation conditions, the ground resolution of the data or the remote sensing image and other factors, the dynamic differences in the glacier on both sides of the shearing li were ignored. This may be the main reason for the underestimation of ice surface velo ties in simulation results near the glacier terminus [36]. In further numerical simulatio of glaciers, the effects of the shearing surface need to be considered seriously.
In addition, there were some limitations in the current work. First, although the DS and DOM generated by drones had a high resolution at the centimeter scale, the area aerial photogrammetry could not cover the entire glacier's surface. In previous cases drone-based glacier monitoring, there were only a few glaciers that were smaller and th had a more accessible route covering the entire area [49,57]. The images acquired drones for those glaciers have a larger scale or are less accessible, and only the partia covered regions, such as the terminus region [32,58]. Second, because of the complex d namics of LHG12, the MB could not be calculated from the DSM difference. Even hig precision measurements of elevation change cannot resolve the spatial patterns of ma balance across individual glaciers [59]. The DSM difference is the combination of chang in the surface MB, the englacial and basal components, and mass redistribution as a resu of ice flow [60] and represents elevation change measurements integrating ice motion r ther than the local signals of a specific MB [61]. As such, we need to use continuity equ tions to calculate the MB in further studies.

Conclusions
With the combination of drone-based photogrammetry and GNSS survey, DOM an DSM images with a centimeter-level ground resolution were acquired, and detailed d namic characteristics of LHG12 were obtained. Compared to other single cirque or vall glaciers, the dynamic characteristics of LHG12 are more complex. The main conclusio are as follows: The dynamics of LHG12 were more complex than those of other single cirque or valley glaciers. Below the confluence region and between the two ice flows, there was a shearing surface similar to a "deep ditch" on the surface of LHG12 (Figure 6d). On both sides of the shearing line, the glacier surface parameters, motion, and mass balance characteristics represented significant differences. In previous simulations, due to the limitations of observation conditions, the ground resolution of the data or the remote sensing images, and other factors, the dynamic differences in the glacier on both sides of the shearing line were ignored. This may be the main reason for the underestimation of ice surface velocities in simulation results near the glacier terminus [36]. In further numerical simulations of glaciers, the effects of the shearing surface need to be considered seriously.
In addition, there were some limitations in the current work. First, although the DSM and DOM generated by drones had a high resolution at the centimeter scale, the area of aerial photogrammetry could not cover the entire glacier's surface. In previous cases of drone-based glacier monitoring, there were only a few glaciers that were smaller and that had a more accessible route covering the entire area [49,57]. The images acquired by drones for those glaciers have a larger scale or are less accessible, and only the partially covered regions, such as the terminus region [32,58]. Second, because of the complex dynamics of LHG12, the MB could not be calculated from the DSM difference. Even high-precision measurements of elevation change cannot resolve the spatial patterns of mass balance across individual glaciers [59]. The DSM difference is the combination of changes in the surface MB, the englacial and basal components, and mass redistribution as a result of ice flow [60] and represents elevation change measurements integrating ice motion rather than the local signals of a specific MB [61]. As such, we need to use continuity equations to calculate the MB in further studies.

Conclusions
With the combination of drone-based photogrammetry and GNSS survey, DOM and DSM images with a centimeter-level ground resolution were acquired, and detailed dynamic characteristics of LHG12 were obtained. Compared to other single cirque or valley glaciers, the dynamic characteristics of LHG12 are more complex. The main conclusions are as follows: (1) The terminus of LHG12 retreated by 194.35 m in total and by 19.44 m a −1 on average during 2008-2021. Differential ablation may have led to changes in the morphology of the terminus edge. This caused a marked increase in the terminus' retreat distance. (2) The dynamic characteristics of LHG12 were analyzed via drone-based photogrammetry, and differences in the surface velocity features between the west and east branches were found. The movement in the west branch is faster than it is in the east branch.
Because of the extrusion of the two ice flows, there is a region with a faster surface velocity at the ablation area. (3) In 2019-2021, the maximum annual surface velocity was 6.50 cm day −1 , and during ablation season, the maximum surface velocity was 13.59 cm day −1 , 52.17% higher than it is annually. The surface velocity in the west branch was faster than it was in the east branch. (4) With the intensification of LHG12 ablation in recent decades, the ice thickness is decreasing, leading to the surface velocity decrease. This phenomenon is more obvious in large glaciers than it is in small glaciers in the Qilian mountains.