HIGH-RESOLUTION SURFACE AND BED TOPOGRAPHY MAPPING OF RUSSELL GLACIER (SW GREENLAND) USING UAV AND GPR

This study presents the detailed survey of the northern marginal part of Russell Glacier, SW Greenland using the combination of unmanned aerial vehicle (UAV) photogrammetry and low-frequency ground penetrating radar (GPR) measurements. Obtained digital elevation model (DEM) and ice thickness data from GPR data allowed the generation of high precision subglacial topography model. We report uncertainties arising from GPR, GPS, and DEM suggesting sufficient accuracy for the reconstruction of glacier bed topography. GPR data and generated subglacial topography model does not reveal any possible Nye channel that could be incised into the bedrock, however, we were able to detect englacial tunnel that runs approximately parallel to the ice margin and possibly is a remnant of a tunnel that provided passage for ice-dammed lake waters during the latest jökulhlaups (2007, 2008). We also observe a radar-transparent layer up to 20 m from the glacier surface suggesting the boundary of cold/temperate ice or piezometric surface. The latter one is preferred due to the warm climatic conditions which are supposed to warm up possible winter cold wave. * Corresponding author


INTRODUCTION AND STUDY AREA
The region of the Russell and adjoining outlet glaciers of the Greenland Ice Sheet is one of the most studied areas in Greenland due to its accessibility and various scientifically intriguing features, but previous studies focus on the gathering of data along large catchment areas of glaciers (e.g. Morlighem et al., 2013;Lindbäck et al., 2014).
To obtain a precise model of subglacial topography, particular attention must be paid to the generation of precise DEM of the glacier surface. There are many ways, how the elevation of the glacier surface can be obtained, for example, interpolation of DEM from GPS points or DEMs derived from LiDAR data or photogrammetry. There are profits and downfalls for each methodpoor resolution of the data (GPS point grid and global DEM models), high acquisition costs (LiDAR data acquisition), unsuitable weather conditions (photogrammetry and LiDAR data acquisition). The remote polar locations and inaccessibility of electric power can limit the usage of the newest approaches. For example, Lamsters et al. (2016) demonstrated that an accurate map of the subglacial topography can be generated using only GPS measurements of glacier surface combined with ice thickness data from GPR, but it could be complicated if the ice surface is very articulated. The favour in this study was given to UAV photogrammetry, as a method that brings highresolution topography (HiRT) and takes minimum effort and time. UAV photogrammetry has become an accurate tool in glaciology for the generation of high-precision DEMs (e.g. Ryan et al., 2015;Bhardwaj et al., 2016;Ely et al., 2017;Rossini et al., 2018;Ewertowski et al., 2019;Jouvet et al., 2019;Lamsters et al., 2019Lamsters et al., , 2020. Here, we present the study where UAV photogrammetry products of glacier surface were generated and used in combination with very detailed ice thickness measurements by GPR. This approach is shown to be time and cost-effective for usage in Polar Regions where very detailed information about the ice surface and subglacial topography in a small area of interest is needed (e.g. Finlayson et al., 2018;Karušs et al., 2019).
Our primary focus is on the subglacial topography of study area and possible subglacial/englacial tunnels, which are proposed to drain the ice-dammed lake (e.g. Russel et al., 2011) but, as the digital elevation model (DEM) is necessary for the generation of subglacial topography model, we particularly asses the acquisition and accuracy of photogrammetrical data.
Our study area is located on the northern marginal part of the Russell glacier, SW Greenland ( Fig. 1), which is an outlet glacier draining the Greenland Ice Sheet (GIS). This marginal part of the glacier is particularly dark and rich in cryoconite material (e.g. Kalińska-Nartiša et al., 2017a;Stivrins et al., 2018) that is a major contributor to the darkening of the west GIS as well. Proglacial streams flowing from the Russell Glacier create the Sandflugtdalen valley-sandur, which is a large source for wind-transported dust (Kalińska-Nartiša et al., 2017b). The deposition of sandy and silty sediments in sandur is promoted by periodical jökulhlaups (Česnulevičius, Šeirienė, 2009;Russell, 1989Russell, , 2007, with the latest ones occurred in 2007 and 2008 (Russell et al., 2011). These jökulhlaups originate from the ice-dammed lake located near our study area ( Fig. 1) and drain through a ~1 km long englacial/subglacial tunnel (see Fig. 1 in Russell et al. (2011). Despite its small size, the study area comprises the marginal part of glacier and adjacent bedrock slope that is overlain by moraine ridges rising at least 50 m above ice margin. The surface height of the UAV study area changes from 395 to 516 m. The area of GPR measurements is smaller than the area of aerial survey (Fig. 1).

Figure 1. The location of the study area. © Sentinel Image
Copyright 2020, ESA The adjacent glacier surface has even more articulated relief of ice bumps and lows, crevasses and supraglacial river canyons that is not suitable for GPR measurements in the field. The compressional ice flow regime is responsible for visible stacked ice thrust planes and shear bands. Measuring the deformation of basal ice at the margin of Russell Glacier, Chandler et al. (2005) found that basal ice was actively deforming and recorded high rates of shear resulting from high effective stress due to strong longitudinal compression. He detected also stick-slip motion indicating sliding at the bed or shearing within discrete layers of debris-rich ice. It was reported also by Knight (1992) that strain rates increase markedly toward the margin, and compressive strain becomes increasingly important. These notions must be considered when interpreting GPR data regarding the internal structure of the glacier.
The previous mapping campaigns of the thickness and bed topography of the Russell Glacier catchment area has been made mostly using Operation IceBridge data and some other ground-based and airborne radar surveys (Morlighem et al., 2013;Lindbäck et al., 2014). These previous mappings include the whole ablation area stretching 80-90 km inland, and due to this the spatial resolution of ice thickness and bed topography DEMs was 250-500 m that is not comparable to the aim and small but very detailed mapping conducted in this study.
For UAV surveys, the studied area provides an excellent opportunity to test the DEM accuracy because of the large differences in height amplitude in a relatively small area and various surface characteristics including colour differences of white ice, dark ice, dark sediments and water surfaces.

METHODS
The expedition to the Russell Glacier, SW Greenland, was conducted at the end of July 2016. Aerial survey with small and low-cost unmanned aerial vehicle (UAV) and groundpenetrating radar (GPR) measurements of the study area were performed on July 27 with the main task to generate digital elevation model (DEM), orthomosaic, subglacial topography and ice thickness maps. GPS receivers were used for the measurements of the GPR profile points and UAV ground control points.

Aerial survey
In this study, DJI Phantom 3 Advanced quadcopter was used for image acquisition. In order to acquire pictures with optimal overlap between flights, automatic mission planning was done with Pix4Dcapture application. The flights were planned as parallel top-down missions and the launch point was on the highest point of the survey territory. It is known, that DJI Phantom 3 Advanced is calculating its altitude based on the barometric sensor relative to the initial launch point. Survey height was set at 80 meters, but, considering, that glacier surface was relatively articulated, and survey territory was spreading downwards the slope, the altitude of the UAV was deviating between 80 and 120 m. Ground sampling distance was also variable -3.46 -5.18 cm/pix. Due to the limitations of UAV battery capacity, four flights were needed to survey 0.45 km 2 territory. The speed of UAV was set to 9-10 m/s. The total amount of taken images was 453. Because of the shortage of batteries, the final flight was performed the next day.

GPR survey
To determine the subglacial topography of the study area, 2.8 km of GPR profile lines in total were recorded on July 27, 2016. Steep pronounced linear topographic features of glacier surface prevented GPR data collection parallel to glacier margin as a result altogether 10 GPR profiles were aligned perpendicular to the Russell glacier margin. The first profile was installed approximately 30 m of steep glacier margin that borders the ice-dammed lake, and each subsequent profile was set 50 m further along glacier margin. Each GPR profile was split into 50 m long sections. Altitude and coordinates of the start and end point of each GPR profile section were determined by using centimetre accuracy GPS system.
Geophysical measurements were performed with the Zond 12-e GPR system ( Fig. 2A). To obtain a clear reflection of the glacier bed as deep as possible, an hand held antenna with a low centre frequency (38 MHz) was used. During the data acquisition time window of 2000 ns was used, which allowed us to detect the reflection of a glacier bed up to 160 m beneath the ice surface. The GPR data were processed and analysed with Prism 2.5 software. During the processing of GPR data time-dependent signal gain function, background removal filter, and Ormsby band-pass filter with a low frequency cut off at 10 MHz and a high frequency cut off at 76 MHz was used. The GPR signal propagation speed and its corresponding dielectric permittivity (ε) value were determined by using englacial reflections represented as hyperbolas (Moore et al., 1999;Bradford, Harper, 2005). Altogether 9 hyperbolas were inspected. It was calculated that ε = 3.73 +/-0.42 at 99% confidence level. The calculated value of ε fits well in the interval of values of ε for temperate glaciers (3.2-4.6) reported by other researchers (Bradford, Harper, 2005;Blindow et al., 2010;Martín-Español et al., 2013). We used the propagation of uncertainty calculations as described by Navarro and Eisen (2009), to determine the error in calculated depth values. Calculated depth uncertainty near glacier margin, is approximately ± 1.14 m, while at the highest part it is ± 4.59 m.
During the analysis of the obtained GPR measurements, the two-way travel time for the basal reflection was determined along with each GPR profile. The distance between the individual readings varied depending on the complexity of the basal topography.

GNNS survey
The coordinates of start and end points of each GPR profiles were determined by GPS system Magellan Promark 3 that is composed of two GPS receivers. Working with this GPS system allowed us to take measurements and post-process results without GSM or radio support. GPS data postprocessing was carried out in the GNSS Solutions software package. Additional corrections from Kellyville (Kangerlussuaq) GPS station were downloaded and local base station coordinates were calculated after processing of the recorded observations against Kellyville GPS station. Local base station coordinates were calculated with the following precision: 0.131 m horizontal confidence interval, 0.092 m vertical confidence interval. All survey points were processed against the observations from the local base station (stop-n-go data acquisition). Mean confidence intervals for processed points were 0.349 m (st. dev. 0.322 m) horizontal confidence, 0.460 m (st. dev. 0.345 m) vertical confidence.

Development and validation of models
Photogrammetric processing was carried out in Agisoft Metashape 1.6. software. Computer configuration was as follows: Intel Xeon E5-2640@2.40 GHz, 64 GB RAM, NVIDIA RTI 2080 Ti GPU. Processing workflow and parameters were set according to the official Agisoft guidelines (Agisoft, 2020). The alignment accuracy was set to high; the key point limit was set to 100000, and the tie point limit was set to 0. The reconstruction parameters were as follows: dense point cloud qualityhigh, depth filtering modeaggressive.
Initial image alignment was performed with image reference preselection using UAV camera positions data. After alignment GCPs were used for more precise georeferencing. GCPs were divided into two groups -11 control points for model generation and 5 checkpoints for model validation (Fig. 3A). All necessary parameters for reference points (such as point vertical and horizontal accuracy) were provided to the algorithm via photogrammetry settings. After camera alignment, the sparse point cloud was filtered, using gradual selection. The input criteria for gradual selection were "Reconstruction uncertainty" with level <10 and "Projection accuracy" with level >10. RMSE reprojection error after model generation was 0.32 m. Sparse cloud resulted in ~4 000 000 points, dense point cloud in ~71 000 000 points. Orthomosaic resolution -4.24 cm/pix, DEM resolution 8.47 cm/pix. The RMS error of control points as given Agisoft Metashape was 23 cm. Model validation by checkpoints yielded the mean model deviation of 0.28 m ( Fig  3B). As a result, orthomosaic and DEM were produced and exported.
After model creation, additional steps were carried out for model precision estimation, following James et al. (2017) guide for photogrammetric uncertainty estimation by calculation of the precision maps. Since the algorithm provided in the original work is developed for Agisoft PhotoScan 1.2, and taking into account updated version of it on the developer's website for Agisoft PhotoScan 1.3, in order to work with it in Agisoft Metashape 1.6, the algorithm was refactored to comply with the last implementation of Agisoft Python API. No additional changes that could affect algorithm logic were made. Precision map calculation was summarised from 1000 iterations of ~200,000 points large sparse cloud "bundle" adjustments. The resulting sparse point clouds were processed by sfmgeoref software tool (James et al., 2012) into a point cloud that accompanies point with mean deviation values from all three axes. The point cloud was then interpolated into a triangulated network (TIN) in QGIS. In the final result, the only deviation from Z-axis is used, since it was the most important parameter for the precision of the final 3D model. The mean point variance throughout the model was 0.291 m with a standard deviation 0.093 m (Fig 3C). Precision estimates for GPR survey territory were as following: 0.186 m altitude variance with 0.043 m standard deviation (Fig 3D).
The strong deviation was observed near model edges that can be explained with a lesser amount of cameras used for the model production, as well as the lack of GCPs in these areas. GPR survey territory DEM has a lesser point deviations since it is situated in the centre of the complete photogrammetric model. In this case, these regions were situated in the middle of flight trajectories. Closer examination of these regions showed, that deviating points had lesser number of cameras used for point extraction than the points nearby that had smaller deviation (~3-4 cameras against 16-17 for points with lesser deviation). The source of these points was in two flights that were performed in two different days. Since optical conditions were different between two days and the time picked for UAV mission performance was different, it resulted in various shadows appearing on the survey territory that were inconsistent between flights. Extracted points in this area represented optical artifacts from both days that were not present throughout the survey.
To generate the final maps, the orthomosaic and DEM were converted to the UTM zone 22N projected coordinate system in ESRI ArcMap 10.6.1., where the models of the ice thickness and subglacial topography were created as well using kriging (ordinary) interpolation method. We used points with ice thickness obtained from GPR data and then determined the values of the ice surface elevation from the created DEM for the same points. The values of the measured ice thickness points were subtracted from the values of the ice surface elevation obtaining points with the altitude of the glacier bed, which were later used for the interpolation. 3D visualisation of the subglacial topography and ice surface was performed in the ESRI ArcScene application.
Main steps of data gathering and processing are summarized in block diagram below (Fig. 4).

RESULTS AND DISCUSSION
Altogether, obtained B-scans (Fig. 5) were of high quality and it was possible to trace not only reflections from the glacier bed but also to distinguish the boundary between the zone of intense and sparse GPR signal scattering (upper 15 to 20 m). The latter is a so-called radar-transparent layer, which is usually interpreted as cold ice (e.g. Gusmeroli et al., 2010;Reinardy et al., 2019). The alternative hypothesis would be that the boundary between the sparse-scatter zone and zone of intense scattering is related to the piezometric surface as proposed by Jania et al. (1996), Murray et al. (2000), Bredford and Harper (2005). If we assume that the upper part of the glacier consists of cold ice, it should be related to the winter cold wave that penetrates through thin ice in winter. We see this explanation as unlikely due to the overall warming of the western part of the GIS, where the melt extent and surface runoff has increased in the last decades (van As et al., 2012). Considerable surface melting can be observed in the field as well, suggesting that the possible effect of winter cold wave must disappear quickly due to warm summer temperatures. We prefer the interpretation of the observed GPR signal scattering layer as evidence for the piezometric water surface. This explanation agrees very well with the borehole data in the Russell glacier which have revealed the presence of pressurized subglacial water (Aaltonen et al. 2010). Our GPR data confirm that ice margin is definitely warm-based and that the surveyed part of the glacier is temperate (possible exception could be only the upper 15-20 m of ice where a radar-transparent layer exists, although we prefer the explanation of piezometric surface). Reflections visible as separate hyperbolas (Fig. 5) were interpreted as reflections from the englacial tunnel according to numerous studies (e.g. Murray et al., 2000;Baelum, Benn, 2011). It was possible to trace one distinct hyperbola through the almost all recorded B-scans except the latest one. The second hyperbola adjacent to the first one was visible in several B-scans. We interpret them as an englacial tunnel (Fig. 6) that possibly is a remnant of a larger tunnel that provided passage for ice-dammed lake water during the latest jökulhlaups that occurred in 2007 and 2008 (Russell et al. 2011). The possible remnant of this tunnel is visible outside our study area at the lower elevation (Fig. 2B). We do not detect any other hyperbolas that could be related to the possible drainage system in this particular part of the glacier. If we were able to record these well-visible hyperbolas in GPR data, this englacial tunnel must have been survived the internal deformation of moving ice suggesting a very slow ice movement near the glacier margin.
In some cases in the up-glacier direction, B-scans were more blurred and tracing of bed reflection was challenging. We assume that radar signal was altered by high water content in ice and complex/debris rich ice structure as well, formed due to compressional ice flow that was confirmed by Chandler et al. (2005).
We do not observe a distinct bedrock depression close to the ice margin that could be attributed to the possible Nye channel. Instead, the elevation of the subglacial slope gradually decreases in the direction of the glacier near the ice margin and increases further up-glacier (Figs 6 and 7). In the southern part of the study area, a wide subglacial depression is observed that stretches mainly parallel to the ice flow direction and its morphometry does not suggests meltwater erosional origin (Figs  6 and 7). Russell et al. (2011) suggested that there could be bedrock topography control on conduit location and that jökulhlaups flow through an incised bedrock-walled Nye channel. Our observations do not prove the aforementioned assumption of Russell et al. (2011) about the possible Nye channel but we have shown the existence of an englacial tunnel that could be a remnant of a tunnel, which was active in jökulhlaups.
Figure 7. 3D model of glacier surface (overlain by orthomosaic) and subglacial topography

CONCLUSIONS
Obtained results demonstrate that it is possible to collect highquality GPR data on the glacier even if the ice surface is uneven. A combination of DEM obtained from UAV photogrammetry and ice thickness data from GPR data allowed to create a high precision model of subglacial topography. We show that such an application can be used if the ice surface is complex. The usage of UAV in the polar environment allowed the generation of DEM and orthomosaic with high accuracy, fully sufficient for GPR data which include the depth uncertainty from ± 1.14 m near glacier margin and ± 4.59 m at higher elevations. The mean altitude deviation of DEM was 0.28 m as validated by checkpoints, and mean point variance throughout the model was 0.291 m, although the precision estimates exactly for a smaller GPR survey territory were 0.186 m of altitude variance with 0.043 m standard deviation. The precision of DEM in our case arises mainly form the GPS errors of GCPs due to usage of the particular GPS system.
GPR data suggest a warm-based ice margin with the radartransparent layer in the upper 15-20 m of the ice surface. We interpret the boundary between the upper sparse scatter zone and the lower zone of intense scattering as a piezometric surface. An alternative explanation would be the presence of cold ice due to the penetration of winter cold wave but this is unlikely due to the warm climatic conditions.
Our main finding is the englacial tunnel that runs approximately parallel to the ice margin and possibly is a remnant of a larger tunnel that provided passage for ice-dammed lake waters during the latest jökulhlaups that occurred in 2007 and 2008 (Russell et al. 2011). GPR data and generated map of subglacial topography does not show any possible Nye channel that could be incised into the bedrock.