Surface Characteristics, Elevation Change, and Velocity of High-Arctic Valley Glacier from Repeated High-Resolution UAV Photogrammetry

: Unmanned Aerial Vehicles (UAVs) are being increasingly used in glaciology demonstrating their potential for the generation of high-resolution digital elevation models (DEMs) that can be further used for the evaluation of glacial processes in detail. Such investigations are especially important for the evaluation of surface changes of small valley glaciers, which are not well-represented in lower-resolution satellite-derived products. In this study, we performed two UAV surveys at the end of the ablation season in 2019 and 2021 on Waldemarbreen, a High-Arctic glacier in NW Svalbard. We derived the mean annual glacier surface velocity of 5.3 m. The estimated mean glacier surface elevation change from 2019 to 2021 was − 1.46 m a − 1 which corresponds to the geodetic mass balance (MB) of − 1.33 m w.e. a − 1 . The glaciological MB for the same period was − 1.61 m w.e. a − 1 . Our survey includes all Waldemarbreen and demonstrates the efﬁciency of high-resolution DEMs produced from UAV photogrammetry for the reconstruction of changes in glacier surface elevation and velocity. We suggest that glaciological and geodetic MB methods should be used complementary to each other.


Introduction
Estimation of elevation and mass changes, velocity, and extent of glaciers are the key properties for the understanding of their responses to climate change and future development. Only multitemporal studies give clues to the dynamics of seasonal and annual changes, their character, and their impact on glacier evolution, which is important for future models and predictions [1]. Glaciers around the globe including mainland Norway, Svalbard, and the Kaffiøyra region in NW Spitsbergen where the glacier of this study is located [2][3][4][5][6][7][8] have been constantly retreating (except surging events), thinning at fronts, and showing more negative mass balance (MB) in the twenty-first century [9]. Thus, monitoring glacier changes are crucial for the understanding of contemporary and future dynamics of the cryosphere and climate system. It has been found that smaller glaciers experience more negative surface MB in Svalbard [7]. Thereby, small glaciers are especially sensitive and prone to melting and even disappearance in this century [10] but their changes are less well resolved, frequently due to the lower resolution and noise of satellite-derived products [6].
Analysis of digital elevation models (DEMs) is one of the most used methods for glacier surface change detection [1,11,12]. Differences between two or more DEMs can be used to calculate geodetic MB. There are multiple techniques in DEM acquisition for

UAV Survey
UAV surveys were conducted on August 11-13, 2019, and on August 20, 2021. The surveys were executed with two DJI UAV platforms-Phantom 4 Pro v2.0 (2019) and Phantom 4 RTK (2021). The technical specifications of both platforms are reported in Table  1. In 2019, an automated mission was planned in the Drone Harmony application. In order to get a constant flight altitude, the terrain following function was used. As a source of terrain elevation values, the ArcticDEM [60] was incorporated in Drone Harmony and Field studies of the glaciological MB of Kaffiøyra region glaciers have been carried out for 25 years, each year since 1996, and provide some of the longest measurement series from Svalbard [4]. The mean annual MB of the Waldemarbreen was: −0.72 m w.e. in 1997-2015, and −0.85 m w.e. in 2011-2015 [4]. In 2019, the mean annual MB decreased even to −1.06 m w.e. [32]. Waldemarbreen is a valley glacier facing the southwest and it has one accumulation area (cirque) that has decreased significantly in recent decades. The accumulation area ratio during 1996-2009 has varied from even 0% to 48% [58]. The annual equilibrium line altitude (ELA) of Waldemarbreen has been increasing with an estimated trend of +7.4±2.6 per year from 1996-2015, reaching 525 m a.s.l. in 2015 [4].
The meteorological records of the Kaffiøyra plain show that the mean air temperature between 2013-2017 was −2.0 • C, and the summer temperature has been increasing by 1.2 • C from 1975 to 2017 [59]. Since 1909, the margin of Waldemarbreen has retreated by 818 m, which is equal to 20.1% of its length [4]. From 2000From -2009, the retreat of the glacier margin accelerated to 11 m a −1 [3]. Till 2019 the glacier was characterized by extensive icings near its front all year round but now they are lost completely during the summer. The existence of icings has been speculated to reveal the polythermal nature of the glacier [54], and it has been demonstrated by Karušs et al. [32] that Waldemarbreen is indeed a polythermal glacier with a relict basal temperate ice layer in its upper reaches. Despite being mainly Remote Sens. 2022, 14, 1029 4 of 20 cold, the glacier retains an englacial drainage system within cold ice mainly originated by a hydrofracturing mechanism, rather than a cut and closure mechanism [32].

UAV Survey
UAV surveys were conducted on 11-13 August 2019, and on 20 August 2021. The surveys were executed with two DJI UAV platforms-Phantom 4 Pro v2.0 (2019) and Phantom 4 RTK (2021). The technical specifications of both platforms are reported in Table 1. In 2019, an automated mission was planned in the Drone Harmony application. In order to get a constant flight altitude, the terrain following function was used. As a source of terrain elevation values, the ArcticDEM [60] was incorporated in Drone Harmony and the glacier extent KML file was used for the mission planning area. Since the terrain-following function was just introduced in Drone Harmony right before the survey date, there were many interruptions in software execution that resulted in a 3-day-long mission. Ground control was established with 49 GCPs (Table 2), sparsely scattered around the glacier extent. 25 of them were used as checkpoints (Figure 2a). GCP coordinates were acquired with the Emlid Reach RS2 GNSS receiver. The base station for the GNSS data correction was placed near the polar station of the Nicolaus Copernicus University, on a permanent reference point. The aerial survey was executed at 100 m elevation above the Arctic DEM surface. The true survey elevation varied in different parts of the glacier since its surface elevation was changed from the Arctic DEM acquisition time. From our estimates, the true elevation was in the order of 120-150 m above the glacier surface. Despite that, the ground sampling distance varied around 3.2-3.8 cm/px. The survey was planned with 85% forward and 85% side overlaps. The amount of overlap was chosen, based on previous experience and flight time optimization. More details on aerial and GNSS surveys are given by Karušs et al. [32].

GCPs and Flight Trajectory Corrections
During both surveys, Emlid Reach RS2 GNSS receivers were used for GCP coordinate acquisition. The final coordinates were obtained by using the post-processing kinematic technique in RTKlib software. Precision estimates and measurement quality for GCPs are reported in Table 2. In 2021, GCPs were used in combination with flight trajectory logs to provide the ground reference. The log files from Phantom 4 RTK were processed in Emlid Studio software with the reference data from the local base station. Final image locations were calculated by Aerotas XLS file, which includes the coordinate interpolation between time marks, as well as the position shift, based on camera sensor position. Flight trajectory quality was estimated by GNSS measurement status value (FIX or FLOAT). Usual uncertainty estimations, such as RMSE are not suitable in this situation, because image acquisition is instantaneous. Flight trajectories in 70% were calculated with In 2021, DJI GroundStation was used for automated flight planning and execution. The glacier extent from the 2019 survey was used as a mission area. The same Arctic DEM was used for the terrain following. Since Phantom 4 RTK was used for this mission, the amount of GCPs points was reduced to 19 ( Figure 2b). GCP coordinates were obtained with the Emlid Reach RS2 receiver. The base station was located on the same point as in the 2019 mission. Since Phantom 4 RTK has the same camera, as Phantom 4 Pro v2.0, the Ground Sampling Distance (GSD) values were almost the same. Side and front overlap values were set to 75% each. The glacier territory was divided into 4 equal parts in a way that each part could be covered with 2 batteries at maximum. In such a way, the whole mission was executed in one day.

GCPs and Flight Trajectory Corrections
During both surveys, Emlid Reach RS2 GNSS receivers were used for GCP coordinate acquisition. The final coordinates were obtained by using the post-processing kinematic technique in RTKlib software. Precision estimates and measurement quality for GCPs are reported in Table 2. In 2021, GCPs were used in combination with flight trajectory logs to provide the ground reference. The log files from Phantom 4 RTK were processed in Emlid Studio software with the reference data from the local base station. Final image locations were calculated by Aerotas XLS file, which includes the coordinate interpolation between time marks, as well as the position shift, based on camera sensor position. Flight trajectory quality was estimated by GNSS measurement status value (FIX or FLOAT). Usual uncertainty estimations, such as RMSE are not suitable in this situation, because image acquisition is instantaneous. Flight trajectories in 70% were calculated with FIX status and 30% FLOAT. Areas with FLOAT solution were located near the mountains and in the highest parts of the glacier (glacier cirque). These glacier parts usually have obscured the sky. Average PDOP values were estimated as good~2.0.

Orthophotos and Digital Surface Models from UAV Survey
Agisoft Metashape 1.7.4 was used for photogrammetric processing. Both surveys were processed on the same system-Intel Xeon E5-2640 v4CPU, 64 GB RAM, NVIDIA GeForce RTX 2080 Ti. The processing parameters were set according to the authors' previous experience [31,32,[34][35][36][37] and other research publications [17,43,61]. The processing flow is presented in Figure 3. The processing includes a standard procedure of UAV imagery in Agisoft Metashape. The key difference from other procedures, described in previous publications [32,[34][35][36][37] is in filtering disabling during the dense cloud generation. The purpose of this step is to reduce the decrease of points in the regions of poor imagery data, e.g., the upper part of the glacier that was covered with fresh snow. Instead, the confidence calculation was enabled, and points were filtered manually afterward, based on the confidence value. The points with confidence less than level 5 were removed. Some regions in the upper part of the glacier were excluded from this procedure in order to save the points in problematic areas. As a result, DEMs and orthomosaics were exported from Agisoft Metashape with the following resolution and accuracy: 2019 models-10.360 cm RMSE on control points, 27.287 cm RMSE on checkpoints, DEM resolution 12.5 cm/pix; 2021 models-2.417 cm RMSE on control points, 38.703 cm RMSE on checkpoints, DEM resolution 11.9 cm/pix. Both DEMs and orthomosaics were exported in WGS 84/UTM zone 33N coordinate system and imported in QGIS.
Produced orthomosaics were analyzed in QGIS and their co-location was estimated by comparing the overlapping areas. Since the aim of the survey was to map the glacier, the periglacial territory was mapped only in some parts. In co-location analysis, only these parts were compared. Since they represent the border areas of the survey, their analysis must be carried out with caution. Borders of the photogrammetric survey usually carry the biggest uncertainties, because of a smaller image count. In some parts, the image count was less than 4 images. Possible changes in extra-glacial terrain must be also considered. Fortunately, the Waldemarbreen glacier is adjoining steep mountain slopes, and they were partly visible on the orthomosaics. Both DEMs were analyzed by comparing elevation values in extra-glacial parts. As a result, both DEMs were co-located. Areas with bigger uncertainties were compared with model precision estimates. Precision estimates were of the same order as the uncertainties. Possible co-location errors on the glacier surface were analyzed by comparison of model precision values. Since both models were developed under the same coordinate reference system, using the same base reference point and post-processing methodology their uncertainty sources can be considered the same. regions in the upper part of the glacier were excluded from this procedure in order to save the points in problematic areas. As a result, DEMs and orthomosaics were exported from Agisoft Metashape with the following resolution and accuracy: 2019 models-10.360 cm RMSE on control points, 27.287 cm RMSE on checkpoints, DEM resolution 12.5 cm/pix; 2021 models-2.417 cm RMSE on control points, 38.703 cm RMSE on checkpoints, DEM resolution 11.9 cm/pix. Both DEMs and orthomosaics were exported in WGS 84 / UTM zone 33N coordinate system and imported in QGIS. Produced orthomosaics were analyzed in QGIS and their co-location was estimated by comparing the overlapping areas. Since the aim of the survey was to map the glacier, the periglacial territory was mapped only in some parts. In co-location analysis, only these parts were compared. Since they represent the border areas of the survey, their analysis must be carried out with caution. Borders of the photogrammetric survey usually carry the biggest uncertainties, because of a smaller image count. In some parts, the image count was less than 4 images. Possible changes in extra-glacial terrain must be also considered. Fortunately, the Waldemarbreen glacier is adjoining steep mountain slopes, and they were partly visible on the orthomosaics. Both DEMs were analyzed by comparing elevation values in extra-glacial parts. As a result, both DEMs were co-located. Areas with

Glacier Extent, Surface Elevation Change, and Mass Balance
QGIS raster calculator was used in the production of DEM of difference (DoD) subtracting DEM of 2019 from DEM of 2021. The output was clipped within the Waldemarbreen extent (2019) polygon, which was manually digitized from the generated orthomosaic. The accuracy of the glacier extent thus should not be higher than pixel size. Considering the high resolution of the orthomosaic (~3 cm), this error can be considered insignificant. However, a higher error, which cannot be precisely calculated, could be related to issues delineating the glacier outline in some places, where it was completely covered by debris. DoD was used to calculate the geodetic MB using an ice density of 910 kg m −3 for the conversion into water equivalent (w.e.). This value of ice density is commonly used in the MB calculations for the glaciers in the Kaffiøyra region [4]. A 1-m-resolution geodetic MB raster was smoothed using the focal statistics tool in ArcMap to generate the contour lines.
The glaciological method selected to evaluate the MB of Waldemarbreen consisted of direct measurements based on a network of ablation stakes, supported by observations of the snow cover in snow profiles (snowpits) [56,62]. Surface ablation was measured in 2019-2021 at the same 21 points. All the ablation poles were inserted to a depth of 10 m using a steam-driven Heucke ice drill. An ice density of 910 kg m −3 was used to convert ablation thickness to water equivalent. Where the snow was found on the glacier, snow density was measured using a steel cylinder and sample weight was determined using professional dynamometer scales and applied to the computations. It should be stressed that in the analyzed period, snow and firn were nearly not observed on the glacier at the end of the ablation season, thus, for the geodetic MB, only ice density was used.
The specific MB refers to an individual point balance on the glacier. The ablation converted to water equivalent for the individual stakes was interpolated to cover the whole surface area of the glacier. Mean values of ablation for the whole area of the glacier were determined using the method proposed by Østrem and Brugman [63]. The balance year adopted for this study lasted from October to September of the following year (stratigraphic system), including the entire accumulation and ablation seasons [64].

Surface Velocity from UAV Products
UAV SfM approach for glacier surface change and velocity detection bears a compromise of ultrahigh-resolution high accuracy data and local area scale. The potential of these data in the field of glacier velocimetry is still being investigated [1,44,47]. Many software packages for glacier velocimetry applications require satellite aerial data and are not opti-Remote Sens. 2022, 14, 1029 8 of 20 mized for high-resolution UAV inputs (e.g., CARST, autoRIFT, GIV) [65][66][67]. It is possible to use UAV-derived results in this software, but, for optimal results, they must be sufficiently downscaled. There were some successful applications of the ImCorr plugin [68] from the QGIS software package on UAV data and glacier velocity calculations from rock-covered glaciers [69,70]. The specifics of rock-covered glaciers are in a large number of features that can be tracked with image-based particle tracking algorithms. Waldemarbreen is covered with rocks only near lateral margins, which are close to the mountain ridges. However, there are other features presented on the glacier surface, such as crevasses, crevasse traces, and supraglacial channels. Since the ImCorr plugin showed promising results on UAV data from rock-covered glaciers, we decided to use it on Waldemarbreen data. ImCorr algorithm allows estimating the glacier surface velocity by raster image correlation between two surveys. This algorithm is based on particle tracking analysis (also known as "image-to-image cross-correlation", "feature tracking") [68]. It searches for patterns with similar textures within images by using a cross-correlation. Fleischer et al. [71] mentions that sub-pixel precision can be achieved with this algorithm. Displacement vector and distance can be measured between matched patterns. To estimate the glacier velocity, DEMs were resampled to 1 m resolution and hillshade rasters were produced. Hillshades were used as source files for computation. Multiple runs of the ImCorr algorithm were performed with different search chip sizes and reference chip sizes. The search chip size varied from 16 to 256 cells and the reference chip size from 16 to 128. The best results should be achieved if the search window size is big enough to represent only homogeneous displacement without local deformations [72]. As the surface of Waldemarbreen mostly consists of bare ice and only in some parts is covered with rocks, search patterns on such glacier will mainly consist of crevasses, crevasse traces, supraglacial streams, and some traces of primary stratification. Therefore, the search chip size must be big enough to include such features. Multiple experiments showed better performance with high search chip size-64 x 64 cells and reference chip size 128 x 128 cells. Another study mentions the same settings to be suitable for such investigations [70]. For optimization purposes, a grid spacing of 50 m was used. Correlation vectors were filtered using a moving average spatial filter. Displacement vectors with values in a direction different by more than 15 degrees were counted as unambiguous values and were removed. Additional visual analysis was carried out by examining orthomosaics and hillshades with the corresponding correlation vectors. A surface velocity map was produced by Inverse Distance Weighted interpolation of the obtained values in QGIS.

Glacier Surface Elevation Change and Geodetic Mass Balance
Our results show that the average glacier surface elevation change from 2019 to 2021 was −2.92 m, which equals −1.46 m a −1 . A maximum elevation change of −11 m was reached at a small spot near the southern glacier slope due to the retreat of a steep ice cliff, but overall maximum surface change values were between −7 and −8 m at isolated sites near the southwestern margin of the glacier. The glacier surface change was calculated to the geodetic MB using an ice density of 910 kg m −3 for the whole glacier as no snow was left on the glacier in 2020 due to extreme melting. The accuracy of the geodetic MB resulting from the accuracy of the DEMs equals ±0.391 m w.e.
The derived cumulative geodetic MB in 2019-2021 was −2.66 m w.e. The mean annual geodetic MB was −1.33 m w.e. a −1 . The most negative geodetic MB was reached near the southwestern glacier margin (Figure 4a). These more negative geodetic MB values correlate well with the glacier aspect ( Figure 5a) and slope (Figure 5b)-the southwestern marginal part of Waldemarbreen has the steepest slope except for the glacier cirque. This slope also faces the south indicating the substantial influence of the slope exposition on the MB. A wide central part of Waldemarbreen is characterized by relatively uniform geodetic MB. The area of the less negative geodetic MB is located in the middle part of the relatively flat glacier cirque, but this changes at the highest and steepest facing slope of the cirque. The

Glacier Surface Characteristics
The altitude of the glacier surface decreases from ~530 to ~135 m. The steepest part of the surface of Waldemarbreen is located on the northern exposition slope of the southern part of the glacier cirque (Figure 5b). Overall, the inclination of the glacier surface is small, and the mean slope angle is 8 degrees.
The glacier surface is cut only by a few permanent supraglacial streams (Figures 5a,  6). The two largest ones begin to develop only some 600 m from the glacier terminus from meltwater accumulating in a flat glacier surface area ( Figure 5). The third one collects meltwater from the higher part of the glacier, which, firstly, flows into the transverse crevasse and then follows the glacier surface gradient to the south. All these streams are meandering, but the pattern of meanders seems to undergo radical changes yearly. For example, one of the largest supraglacial streams was almost straight in 2021 but highly meandering in 2019 (Figure 5c,d). In this part of the glacier, the surface ablation reaches 2.5-3.0 m w.e. per year, thus preventing the further incision of supraglacial streams. We observed that the meandering of supraglacial streams sometimes can be connected to transverse crevasses, which locally divert surficial drainage even by 90° (Figure 5c,d). The ice velocity estimated from the redistribution of surface crevasses in this area reaches 6 m per year (Figure 6d). The meltwater from the glacier cirque seems to mainly penetrate the glacier through crevasses, thus preventing the formation of permanent supraglacial streams. This crevassed area of the Waldemarbreen is very well visible on the map of aspect (Figure 5a), and it has been attributed to the hydrofracturing and formation of temperate ice by Karušs et al. [32].
The southern part of the ablation area of Waldemarbreen is characterized by a steeper slope (Figure 5b), and it is reflected in the surface drainage pattern. The majority of supraglacial streams in this part of the glacier follows the glacier surface gradient and flows towards the south, although the mean aspect of the glacier is southwest (226°) (Figure 5a).

Glacier Surface Characteristics
The altitude of the glacier surface decreases from~530 to~135 m. The steepest part of the surface of Waldemarbreen is located on the northern exposition slope of the southern part of the glacier cirque ( Figure 5b). Overall, the inclination of the glacier surface is small, and the mean slope angle is 8 degrees.
The glacier surface is cut only by a few permanent supraglacial streams (Figures 5a and 6). The two largest ones begin to develop only some 600 m from the glacier terminus from meltwater accumulating in a flat glacier surface area ( Figure 5). The third one collects meltwater from the higher part of the glacier, which, firstly, flows into the transverse crevasse and then follows the glacier surface gradient to the south. All these streams are meandering, but the pattern of meanders seems to undergo radical changes yearly. For example, one of the largest supraglacial streams was almost straight in 2021 but highly meandering in 2019 (Figure 5c,d). In this part of the glacier, the surface ablation reaches 2.5-3.0 m w.e. per year, thus preventing the further incision of supraglacial streams. We observed that the meandering of supraglacial streams sometimes can be connected to transverse crevasses, which locally divert surficial drainage even by 90 • (Figure 5c,d). The ice velocity estimated from the redistribution of surface crevasses in this area reaches 6 m per year (Figure 6d).
The meltwater from the glacier cirque seems to mainly penetrate the glacier through crevasses, thus preventing the formation of permanent supraglacial streams. This crevassed area of the Waldemarbreen is very well visible on the map of aspect (Figure 5a), and it has been attributed to the hydrofracturing and formation of temperate ice by Karušs et al. [32].
The southern part of the ablation area of Waldemarbreen is characterized by a steeper slope (Figure 5b), and it is reflected in the surface drainage pattern. The majority of supraglacial streams in this part of the glacier follows the glacier surface gradient and flows towards the south, although the mean aspect of the glacier is southwest (226 • ) (Figure 5a).

Glacier Surface Velocity
The glacier has moved 10.6 m on average from August 2019 to August 2021. The mean annual velocity was 5.3 m. The maximum movement reached 28 m for the study period and 14 m per year (Figure 7). The mean flow direction was 247 degrees (southwest) that coincides with the main direction of the glacier, and it ranges from 180 to 300 degrees. However, the general pattern of flow velocity and direction does not completely follow the main direction of the glacier itself (SW) but rather agrees with the steepest slope gradient (Figure 5b), which has more SSW aspect (Figure 5a) in the central-southern part of the glacier. This is the same direction followed also by the supraglacial streams in this part of the glacier. The flow velocity increases directly at the lower end of the glacier cirque, where the most crevassed part of Waldemarbreen begins. This explains the existence of narrow but still open crevasses in this area because sufficient tensile stress is needed for these crevasses to develop and remain open. The flow velocity remains high up to the southern lateral margin but decreases in the northern part of the glacier and near its front.

Accuracy and Prospects of RTK UAV Surveys
Direct image georeferencing with the use of a UAV RTK system is still a novice approach in the field of photogrammetry. There are still a very small number of articles that discuss the usage of RTK systems and their performance to date [49][50][51]. Štroner et al. [50] stated that the standalone use of RTK UAVs without ground control points can introduce errors with drastic magnification in the vertical component of the model. Such usage can be performed with caution and only with pre-calibrated metric cameras. However, once the single GCP is introduced in the survey, the error is minimized. In their tests, Štroner et al. [50] were able to achieve RMSE of ~1 GSD. Our survey showed that the RTK system can achieve RMSE at the same level, but only on control points. RMSE of the checkpoints is still large (3 GSD in XY error and more than 10 GSD in Z error). If we remove 2 checkpoints with the highest RMSE values, we can get down to ~9 GSD (0.242 m). The largest errors were detected on GCPs that were placed in the upper reaches of the glacier with bright snow cover and therefore lacked the proper number of tie points between the images.
It should be noted that many comparisons of RTK UAV systems and non-RTK systems are done in rural or urban areas [49][50][51]. Glaciated areas have their specific attributes-lack of contrast, poor color space, repetitive patterns, etc. Such terrain characteristics are possible obstructions for the SfM algorithm processing and can enlarge RMSE values. On the other side, additional GCP introduction into the survey can control the varia-

Accuracy and Prospects of RTK UAV Surveys
Direct image georeferencing with the use of a UAV RTK system is still a novice approach in the field of photogrammetry. There are still a very small number of articles that discuss the usage of RTK systems and their performance to date [49][50][51]. Štroner et al. [50] stated that the standalone use of RTK UAVs without ground control points can introduce errors with drastic magnification in the vertical component of the model. Such usage can be performed with caution and only with pre-calibrated metric cameras. However, once the single GCP is introduced in the survey, the error is minimized. In their tests, Štroner et al. [50] were able to achieve RMSE of~1 GSD. Our survey showed that the RTK system can achieve RMSE at the same level, but only on control points. RMSE of the checkpoints is still large (3 GSD in XY error and more than 10 GSD in Z error). If we remove 2 checkpoints with the highest RMSE values, we can get down to~9 GSD (0.242 m). The largest errors were detected on GCPs that were placed in the upper reaches of the glacier with bright snow cover and therefore lacked the proper number of tie points between the images.
It should be noted that many comparisons of RTK UAV systems and non-RTK systems are done in rural or urban areas [49][50][51]. Glaciated areas have their specific attributes-lack of contrast, poor color space, repetitive patterns, etc. Such terrain characteristics are possible obstructions for the SfM algorithm processing and can enlarge RMSE values. On the other side, additional GCP introduction into the survey can control the variability of the errors. It is obvious that UAV RTK systems can drastically improve the overall framework of UAV fieldworks on the glaciers and more comparisons must be added in the future, especially in the glaciated areas.
The absence of a stable internet connection in polar regions makes RTK usage without a local base station unavailable. In this situation, it is possible to use DJI D-RTK 2 base station and make a connection between the UAV and the base station, however, such a setup requires a geodetic point near the survey area to conduct correct georeferencing. Another approach is to use post-processing kinematics and process the UAV data with other GNSS correction sources, such as third-party GNSS base stations. In this situation, there is no need for a connection between the UAV and the base station. The base station can be placed anywhere near the survey area or even RINEX log files from the NTRIP caster can be used instead. Taddia et al. [51] compared the RTK and PPK approaches in coastal area mapping and concluded that the PPK approach can perform at the same level as RTK if additional oblique imagery is being introduced in the survey. They [51] stated that double-grid flight plans with oblique imagery can deliver the same high accuracy results without GCPs. In our study, the Emlid Reach RS2 receiver was used in base station mode, and post-processing was done in Emlid Studio software. The main problem in this workflow is the inability to control the GNSS measurement solution status during the flight. FLOAT solution status may have bigger RMSE values and can be offset from the absolute position. In our survey, additional GCPs were introduced to compensate for this possibility. In comparison with the study by Taddia et al. [51], our survey area was enclosed by high mountains that obstruct the sky view. In such conditions, the possibility to get poor GNSS measurements is higher than in the coastal environment, for example. This factor is also being magnified by the fact that Svalbard is located in high latitudes, where GNSS signal quality is worse because of limited satellite coverage, poor trajectory geometry, and frequent ionospheric scintillation [73]. We conclude that the use of GCPs is a necessity to add additional control to the outcomes of the survey and must be considered obligatory in RTK surveys in difficult, especially polar, environments.

Compression of Direct Glaciological and Geodetic Mass Balance
Our results show that the mean annual geodetic MB was −1.33 m w.e. in 2019-2021 (Figure 4a), while the glaciological MB was −1.61 m w.e. a −1 (Figure 4b). Differences between these two methods are widely discussed in the scientific literature [5,[74][75][76][77], and there are many variables influencing the possible inconsistencies although a close agreement is usually reported. Nevertheless, it has been noted that each of these methods reveals different processes because the glaciological method measures only surface (point) balance contrary to the geodetic method where the englacial and subglacial melting is included as well [75]. Still, the direct glaciological MB probably remains the most widely used and precise method due to its consistent methodology [64,76], and has the longest time series recorded globally [4,5,8,76].
The geodetic MB method is increasingly used nowadays due to the availability of DEMs generated from satellite-derived data [39] and more recently due to the repeated UAV surveys [1,12,[44][45][46][47]. However, despite the increasing popularity of UAV surveys, several issues remain to precisely capture the glacier surface changes and MB. The main issue could be related to the variable survey times. For example, UAV surveys in this study were done in 2019 and 2021, skipping 2020 due to the encumbrance of the COVID-19 pandemic. This created a situation where the cumulative MB for 2019-2021 is well represented by both glaciological and geodetic methods, which results in comparable annual values, although the direct measurements reveal that the annual balance in 2020 and 2021 was substantially different (−2.28 m and 0.95 m w.e. a −1 , accordingly) due to the extreme melting in 2020 [57].
This finding emphasizes the role of the great interannual variability of MB of Svalbard glaciers, which can be detected only by annual data.
Discrepancies in both methods could also arise due to differences in measurement and survey times. Although the UAV surveys were done in August, the last direct measurements of ablation stakes were performed even later in September. It should also be noted that the geodetic method can be biased by snow events and missing or inaccurate information about snow, firn, and ice density as already pointed out by Fischer [75]. In our case, the geodetic MB could be slightly overestimated as we used the ice density of 910 kg m −3 for the whole glacier due to the extraordinary melt in 2020, when almost all snow melted on Waldemarbreen [57] but some remained in the uppermost reaches of the glacier in 2021.
Glaciological MB has been reported as unable to capture all the local spatial variability of glacier-wide changes due to a measurement network [77]. This could be supported also by this study, where some local MB perturbations are visible on a map of the geodetic MB but missing on the interpolated map of glaciological balance (Figure 4). The geodetic MB map seems to better characterize the variable changes of the glacier surface due to local variations of slope steepness, for example. This is well represented in the southernmost part of the marginal zone of the glacier, where the geodetic MB is able to reflect the sudden increase in the surface changes, which could be mainly related to the southern exposition of this glacier surface slope. This study shows that the geodetic MB calculated from surface changes in high-resolution DEMs (obtained from UAV surveys), can be complementary to the glaciological MB and both methods should be used together, if possible.
A clear trend of decreasing MB of Waldemarbreen is visible from this and previous studies. Glaciological MB was −0.72 m w.e. in 1996-2015, but already −0.85 m w.e. in 2011-2015 [4]. The trend of the average change of MB with time was estimated to be −0.040 m w.e. a −1 for 1996-2015 [4] increasing to −0.062 m w.e. a −1 for 1996-2020 [57], although it has been very variable between individual years [4], due to complicated local patterns of snow accumulation, melt and weather that has been also represented by large variations in ratios of the accumulation area [58]. Negative MB from 1996, which becomes more negative [4] correlates well with the measured increase in the summer temperature by 1.2 • C from 1975 to 2017 in the Kaffiøyra region [58,59] and indicates a clear connection between ongoing climate change and glacier melt in NW Svalbard.

The Pattern of Surface Velocity of Waldemarbreen
We found that the mean annual velocity of Waldemarbreen was 5.3 m during 2019-2021 reaching a maximum of 14 m per year ( Figure 7). As Waldemarbreen is a polythermal glacier with mainly cold ice and basal temperate ice layer only in the upper reaches [32], the majority of flow components should be related to the internal deformation, which in turn is governed by the gradients of glacier bed, ice thickness, and surface slope. In the case of Waldemarbreen, we observe that the general ice flow direction is related to the steepest surface slope gradient. It is very well manifested in the central-southern part of the glacier, where the flow direction is south-southwesterly and even southerly near the lateral ice margin (Figure 5b). This emphasizes the role of the glacier slope steepness in the governance of the local surface flow directions. However, the fastest velocity and the most rapid increase is related to the zone, where the temperate basal ice layer is located [32]. Thus, the increase of the glacier velocity and fastest speeds are hypothesized to be related to the lubrication of the bed and at least a partial component of basal sliding as well. This idea is supported by the well-developed englacial drainage network, hydrofracturing process, and presence of pressurized water beneath the Waldemarbreen [32,55]. Similar ice bed lubrication has been suggested, for example, by Kraaijenbrink et al. [44] for a debris-covered Himalayan glacier and by Melkonian et al. [78] for glaciers in Alaska. All these mentioned examples showed the peak velocities during the summer related to higher liquid water input from increased melt rates and rainwater. Our study does not allow evaluation of the seasonal characteristics of the surface velocity thus only advising this further research opportunity, which could be particularly interesting due to variations of the thermal structure of Waldemarbreen.
The area of the highest flow velocities is also related to the most crevassed part of the glacier. This suggests that the component of the internal ice deformation is significant as well in this site. We found that the pattern of glacier velocity is also related to differences in ice thickness and subglacial topography. These were investigated by Karušs et. al. [32] during the field campaign in 2019. Generally, the flow velocity begins to increase in the area, where the thickest ice is located, although this is the least inclined part of the glacier surface (Figure 5b), and then it follows the ice thickness contours (See Figure 5 in [32]). Higher flow velocities in the southern ablation area of the glacier coincide also with the longitudinal trough located under this part of the glacier and the generally thicker ice there emphasizing the main role of the subglacial topography on the general glacier velocity field.

UAV Surveys for Glacier Surface Velocity Estimation
In the previous studies, it was noted that the accuracy of glacier velocimetry, obtained from imagery data, is affected by many physical and optical parameters of particular glaciers and data [79,80]. For example, the inaccuracy of velocity measurements for the rock glaciers is affected by rough micro-topography, snow coverage, steep slopes, and lack of terrain contrast [79]. Jouvet et al. [80] concluded that the accuracy of ice surface displacement fields depends on the time duration between two data acquisitions. With longer time windows, surface displacement is bigger, resulting in smaller relative errors [80]. This statement is true for glacier areas with lower speeds as well. Regions with lower speeds can produce inaccurate results due to small displacement values.
The quality of co-registration of source data is another key factor for accuracy estimation. Several different approaches can be followed during the imagery or DEM coregistration before velocity field calculations. CARST (Cryosphere and Remote Sensing Toolkit) software uses predefined areas with an assumed absence of movement for the calibration of the aerial images [66], e.g., bedrock. This is a successful approach for satellite imagery processing because each data timeframe consists of one or several images. The uncertainty of the data is described mostly with the satellite sensor and atmospheric ambiguities. Aerial imagery produced from SfM photogrammetry techniques depends on internal model ambiguities that can consist of several error sources: (1) camera sensor and lens imperfections; (2) georeferencing accuracies; (3) human error during GCP placement on imagery; (4) reprojection error. The bigger errors are usually distributed around the borders of the survey areas, where the image count participated in the reconstruction is the lowest [81]. Waldemarbreen does not have any bedrock features except for those outside the glacier and even these not always can be considered stable features due to frequent mass movements on steep and unstable slopes. The bedrock is exposed around Waldemarbreen and is situated on the border of our UAV surveys. Due to this fact, the approach of calibration with a fixed area is not suitable because it can provide another source of errors. Our approach was a co-registration of DEMs based on their geolocation. ImCorr algorithm [68] is matching small sub-scenes in two input images. The images must have the same pixel grid size. To achieve it, we resampled both DEMs to equal resolutions and cropped them to the same grid size. The downside of the results of the ImCorr algorithm is the absence of original geolocation data. This problem is solved in other feature tracking packages, such as autoRIFT [65]. Processing algorithms of autoRIFT allow velocity calculations based on geographical units. However, this package is suitable for velocimetry based on single aerial imagery and does not allow us to evaluate unambiguities from SfM models.
The other issue of glacier velocity estimation from feature tracking was mentioned by Karimi et al. [82] who noted that intense ice melting can prevent feature tracking algorithms to find a coherent displacement pattern. Due to volumetric glacier changes, some of the surface features can propagate towards the opposite direction of the flow e.g., the projection of closed crevasses on aerial imagery can be shifted backward along with the flow because of ice thinning. This factor is a crucial consideration for feature tracking velocimetry of glaciers where crevasses and other supraglacial features are the main key points and should be addressed hereafter. A successful attempt in glacier velocity estimation has been made from terrestrial laser scanning point clouds [15]. With the use of the multiscale model-tomodel cloud comparison (M3C2) algorithm [83], Ulrich et al. [15] were able to compute velocities for each point pair. The benefits of cloud-to-cloud comparison compared with optical imagery feature tracking approach is in 3D vector calculations. The results of this method are presented as 3D vectors and allow to estimate velocities in all directions. Since photogrammetry-derived point clouds could be used in the same manner, they could be considered in future studies as an alternative to optical image feature tracking.

Conclusions
In this study, two repeated UAV surveys were conducted at the end of August in 2019 and 2021 on a small valley glacier, Waldemarbreen, NW Svalbard. Our UAV surveys include all areas of the glacier and demonstrate the efficiency of this method in the means of time, cost, resolution, and accuracy for the generation of DEMs and following reconstructions of the changes in glacier surface elevation, geodetic MB, and velocity. This further suggests the applicability of this technique for high-precision surveys of dynamics and changes of small valley glaciers in Svalbard and the Arctic.
A large number of GCPs (~50) was used in 2019 to ensure the accuracy of DEM and orthomosaic, while less than half of GCPs were used in 2021 due to the application of the RTK UAV platform (DJI Phantom 4 RTK), thus allowing the execution of all flights for one day. The accuracy of RTK flights was more than 4 times better than non-RTK flights, based on control point RMSE (0.0242 m in 2021 survey against 0.1036 m in 2019). We suggest that the use of GCPs is a necessity to add additional control to the outcomes of the survey and must be considered obligatory even in RTK UAV surveys, especially in polar environments.
Glacier surface elevation change, geodetic MB, and surface velocities were obtained from UAV data. The mean annual glacier surface elevation change in 2019-2021 was −1.46 m. It was converted to the geodetic MB of 1.33 m w.e. a −1 . The direct glaciological MB was calculated for the same period, and it was −1.61 m w.e. a −1 . The two methods were analyzed, and they show comparable mean MB, although the direct measurements reveal substantially different balances in 2020 and 2021 of −2.28 m and 0.95 m w.e. a −1 , accordingly, which were not captured from UAV surveys due to the missing survey in 2020. Such extreme melting in 2020 was observed for the first time since 1996. Our results emphasize once more the role of a great interannual variability of MB of Svalbard glaciers. However, the more and more negative MB of Waldemarbreen is clearly observed in the latest years indicating the strong connection between air temperature increase and glacier melt. The glacier extent between 2019 and 2021 decreased by 3%. It was found that the geodetic MB contrary to glaciological MB better characterize variable changes in the glacier surface due to local variations of slope steepness, aspect, and debris content. Not all local spatial variations of glacier surface changes are captured by the glaciological MB due to the measurement stake network. This study shows that the geodetic MB calculated from high-resolution DEMs significantly complements the direct glaciological method and vice versa, thus both should be used together, if possible.
The mean glacier surface displacement derived from UAV surveys was 5.3 m a −1 , reaching a maximum of 14 m a −1 in the central part of Waldemarbreen. The pattern of glacier velocity can be mainly related to the differences in the ice surface steepness, ice thickness, subglacial topography, and internal deformation. We emphasise also the possible impact of the lubrication of the bed and a component of basal sliding, which is supported by the presence of pressurized water beneath Waldemarbreen.