Multisensor validation of tidewater glacier flow fields derived from SAR intensity tracking

Following the general warming trend in Greenland, an increase in calving rates, retreat and ice flow has been observed at ocean-terminating outlet glaciers. These changes contribute substantially to the current mass loss of the Greenland Ice Sheet. In order to constrain models of ice dynamics as well as estimates of mass change, detailed knowledge of geometry and ice-flow are needed, in particular on the rapidly changing tongues of ocean-terminating outlet glaciers. In this study, we validate velocity estimates and spatial patterns close to the calving terminus of such an outlet derived from an iterative offset 5 tracking method based on SAR intensity data with a collection of three independent reference measurements of glacier flow. These reference data sets are comprised of measurements from differential GPS, a Terrestrial Radar Interferometer (TRI) and repeated UAV surveys. Our approach for the SAR-velocity processing aims at achieving a relatively fine grid spacing and a high temporal resolution in order to best resolve the steep velocity gradients in the terminus area and to exploit the 12 day repeat interval of the single-satellite Sentinel-1A sensor. Results from images of the medium-sized ocean terminating outlet 10 glacier Eqip Sermia acquired by Sentinel-1A and RADARSAT-2 exhibit a mean difference of 11.5% when compared to the corresponding GPS measurements. An areal comparison of our SAR velocity-fields with independently generated velocity maps from TRI and UAV showed a good agreement in magnitude and spatial patterns, with mean differences smaller than 0.7 md-1. In comparison with existing operational velocity products, our SAR-derived velocities showed an improved spatial velocity pattern near the margins and calving front. There 8% to 30% higher surface ice velocities are produced, which has 15 implications on ice fluxes and on mass budget estimates of similar sized outlet glaciers. Further, we showed that offset tracking from SAR intensity data at relatively low spatio-temporal sampling intervals is a valid method to derive glacier flow fields for fast-flowing glacier termini of outlet glaciers and, given the repeat period of 12 days of the Sentinel-1A sensor (6 days with Sentinel-1B), has the potential to be applied operationally in a quasi-continuous mode.

In this paper, we are investigating the limits of offset tracking methods and demonstrate that accurate estimates are possible even close to the calving front when choosing appropriate template sizes. As validation datasets, for the derived flow velocity information (magnitude and spatial pattern) we use field measurements based on three independent methods, namely differ-5 ential GPS (dGPS), a terrestrial radar interferometer (TRI), and high resolution imagery from an Unmanned Aerial Vehicle (UAV). The glacier studied is Eqip Sermia, a medium sized marine-terminating outlet glacier in the southwest of Greenland (cf. Fig. 1). We demonstrate that flow velocity estimates generated at a relatively fine ground sampling distance are more accurate close to the glacier's terminus compared to operational, ice-sheet wide ice velocity products which tend to underestimate the glacier's dynamics and thus also the calving flux.

Study Area
For this study, the medium-sized ocean terminating outlet glacier of Eqip Sermia (69 • 48' N, 50 • 13' W), located in Western Greenland, was observed by multiple sensors. Given the available historical geometry and flow velocity survey data, dating back to 1912 (e.g. de Quervain andMercanton, 1925), Eqip Sermia offers ideal preconditions. Eqip Sermia has a calving front roughly 3.5 km wide and 30-200 m high. The long-term flow speed at the terminus was stable for almost a century at about 15 3 md -1 (Bauer, 1968), followed by an acceleration towards the end of the 20th century. Between 2000 and 2005, Eqip Sermia accelerated by 30% as well, doubling the discharge (Rignot and Kanagaratnam, 2006;Kadded and Moreau, 2013;Lüthi et al., 2016). More detailed velocity fields from the recent decade (Joughin et al., 2008(Joughin et al., , 2010 indicate strong spatial variations in flow in the terminus area and a strong acceleration towards the calving front (Lüthi et al., 2016;Catania et al., 2018).

SAR Data
The increasing availability of freely and openly available SAR data over the ice sheets at high spatial and temporal resolutions allows circumvention of acquisition issues intrinsic to optical systems that have been available for several decades. Making use of the all-weather, day/night imaging capabilities of SAR sensors, cloud cover or illumination effects do not interfere with the acquisition schedule. The source data consists of a dual-pol (HH/HV) Sentinel-1A/B Single Look Complex (SLC) C-band 25 (5.405 GHz) time-series starting in October 2014, accessed through the Copernicus Open Access Hub (Torres et al., 2012). Out of these more than 200 acquisitions, 5 scenes complement ground measurement data acquired during a field campaign that took place in August 2016. These S-1A satellite products have a 12-day repeat cycle and all interferometric wide swath (IW-mode) SLC products used were acquired from the same ascending relative orbit (cf. Table 1). After this validation campaign, the S-1B satellite was commissioned, lowering the S-1 repeat interval to 6 days. In addition to the S-1A data sets, a total of 20 Figure 1. Eqip Sermia Glacier in West Greenland. The yellow area indicates the region evaluated in this study using spaceborne SAR images (S-1A and RADARSAT-2), with the red line showing the extent of the mosaic built up using airborne UAV and the blue polygon depicting the area with mean ground-based backscatter coherence >0.7 acquired from the terrestrial radar interferometer (TRI; its position is indicated with the blue triangle). The purple crosses mark the GPS tracker positions. In orange and brown the flowlines across the tongue are depicted (cf. Fig. 14 and 15). The green line depicts the central flowline with distance markers every 2500 m, starting at the glacier's terminus (cf. (Project CSA-SOAR-EU-16821). Based on this allowance, two SLC scenes were acquired using RS-2's Ultra-Fine wide mode operating also at the frequency of 5.405 GHz with a temporal baseline of 24 days (cf. Table 1). The detected HH polarized SAR images from both sensors were geometrically terrain corrected using Range-Doppler geocoding (Meier et al., 1993) based on the "GIMP" Digital Elevation Model (DEM). The Greenland Ice Mapping Project (GIMP) DEM has a grid spacing 5 of 30×30 m (Version 2.1; Howat et al., 2014), which was oversampled to 2.5×2.5 m. As we operated in the DEM geometry, no separate co-registration was performed. No tiepoints were employed during geocoding, as the geolocation accuracy was sufficient (Schubert et al., 2017). The HH polarization was chosen due to its higher signal-to-noise ratio.

SAR Intensity Tracking
With the increasing availability of space-borne interferometric data and using well-established methodologies, measurements 10 of ice-velocities at high spatio-temporal resolutions and accuracies of a few meters per year are performed on a regular basis for  (Gray et al., 1998;Joughin, 2002). As this study focuses on the glacier's flow dynamics at the fast-moving glacier terminus area, interferometric methods were not expected to be viable with a temporal baseline of 12 days. In addition to the temporal decorrelation, InSAR is limited to observing the line of sight. By evaluating both ascending 5 and descending tracks (Joughin et al., 1998) and making an assumption of surface-parallel flow of the surface ice one can retrieve a full 3D displacement using two different tracks. As only ascending geometries were available in this instance for Eqip Sermia, we used a speckle tracking approach to derive the glacier's movements. This method makes use of the backscattered speckle pattern within image patches from subsequent, co-registered image acquisitions to derive two-dimensional offset values by calculating the normalized cross-correlation between the image patches (Gray et al., 2001;Strozzi et al., 2002;Joughin, 10 2002). As this approach does not rely on phase information, using instead the detected SAR image, phase decorrelation caused by meteorological conditions or incoherent and/or rapid flow does not influence the velocity estimation and therefore allows retrievals at higher ice speeds and longer orbit repeat intervals (Gray et al., 2001;Strozzi et al., 2002). Nevertheless, strong changes in the amplitude of the backscattered signal (e.g. due to substantial changes in the presence of surface melt water) may result in a deterioration of the velocity estimates.  Figure 2. Amplitude Match Methodology for 12-day Sentinel-1A SAR products. The same methodology was used for RADARSAT-2 image pairs, but with an increased search window size of 261×261 pixels to accommodate the 24-day repeat orbit.

Outlier Detection and Process Iteration
To cull out bad matches, the two-step approach outlined in Fig. 3 was implemented, using a long-term flow velocity average product (V med ) and a correlation value threshold. For this V med , offsets in X-and Y-direction were calculated using a 101×101 pixel (252.5×252.5 m) template patch size and search window sizes based on the temporal baseline of each of the 5 256 image pairs available between 2014/10/11 and 2018/03/18 (cf. Table 2). Subsequently, a median flow velocity (V med ) and median flow angle map ( med ) were computed based on these 256 flow fields. The correlation signal-to-noise-ratio (SNR) was calculated by dividing the maximum height of the correlation peak (C max ) by the mean level of the correlation function (C mean ) for every pixel (Strozzi et al., 2002).
where ∠ dif f = |∠ med − ∠ map | and ∠ map representing the mapped flow angle from the image pair. For every pixel flagged as an outlier, the intensity tracking was reiterated with the reference patch size increased by 10 pixels in each dimension. The patch size of the search region was defined for every outlier pixel as a function of V med , limiting the search region to three times the expected displacement, thus decreasing the chance of missing the peak correlation. Following the outlier detection, iterative intensity tracking for outlier pixels was repeated for a maximum of five times or until no further outlier pixels were detected. The application of the iteration procedure drastically reduced the number of void pixels, while enabling reasonably small patch sizes and therefore meaningful results even close to the lateral glacier margin and in particular at the calving front.
Following the outlier detection and process iteration step, V map was downscaled to a pixel spacing of 100×100 m by applying 5 a median filter to account for the spatial correlation of adjacent pixels caused by overlapping template patches in the initial V map (cf. Sect. 2.2).

GPS Velocity Data
For validation purposes, seven low-cost single-frequency continuous GPS receivers (Wirz et al., 2013;Buchli et al., 2012) were deployed on the glacier using a photovoltaic system in combination with a battery (cf. Fig. 1). In addition to the receivers on the glacier, a base station was deployed on bedrock. The devices were installed on June 29, 2016: five were recovered by August 25, 2016, a sixth one was recovered in Summer 2017 and GPS solutions at 24 h intervals were calculated at the Geodesy and Geodynamics Lab of ETH Zurich (Wirz et al., 2013). To derive robust average velocities over the relevant image-pair periods, the 24 h GPS solutions were further processed using a three-step approach that begins with the culling of outliers, followed by temporally averaging both, the latitudinal and longitudinal positions as well as the resulting velocities, in the manner described 15 by Ahlstrøm et al. (2013). The velocities of the three sensors located in parts of the glacier with flow velocities higher than 1 md -1 (GPS sensors EG09, EG14 and EG19) were chosen for accuracy assessment, as lower flow velocities exhibit relatively low signal-to-noise ratios (SNR) when intensity tracking was applied to two immediately successive S-1A SAR acquisitions.

UAV survey
To acquire additional ground-truth data, an unmanned aerial vehicle (UAV) was flown over Eqip Sermia on three occasions in 20 August 2016 (cf. Table 3). The UAV surveys were carried out using a SenseFly eBee system, a light, fixed-wing UAV with a wingspan of 96 cm (senseFly SA, 2017). The images were acquired using a modified Sony Cyber-shot DSC-WX220 with 18.2 MP and a sensor size of 6.17×4.55 mm. For every image, the approximate 3D position as well as the UAV's orientation (i.e. roll, pitch, and yaw) were annotated based on information from the on-board GPS and inertial measurement unit (IMU) (senseFly SA, 2016SA, , 2017. A total of 6 flights per mission were carried out in parallel strips within 3 h, resulting in an average of more than 5 overlapping images (i.e. 70% longitudinal/lateral overlap per image), with reduced redundancy near the outer 5 limits of the acquisition areas. For each flight, the images acquired by the eBee UAV were processed by applying the structure from motion technique using Pix4Dmapper Pro (Eltner and Schneider, 2015), resulting in an optical mosaic (RGB) and a digital surface model (DSM) of the area with a ground sampling distance (GSD) of ∼0.17 m (cf. Table 3). The resulting datasets were georeferenced using 5 manually surveyed ground control points (GCP). As it was not possible to reach the northern lateral bedrock due to heavy crevassing, these GCPs were restricted to one side of the glacier, resulting in small residual differences 10 between the geolocation accuracy of some mosaics. To reduce differences between the mosaics, an affine transformation based on a total of 40 manually selected GCPs was performed with acquisition eBee1 (cf. Table 3) as a master reference and the other two acquisitions as slaves.
Due to differences in illumination between the acquisitions eBee1 and eBee3 (cf . Table 3), image matching using the different optical bands was not feasible. We therefore used the DSMs to derive shaded reliefs and input these to the image matching 15 algorithm, to minimize matching errors caused by illumination differences between the data sets. Flow velocity and direction were estimated by applying the approach described in Sect. 2.2. Mosaics eBee1 and eBee3 with a temporal baseline of 4 days (cf . Table 3) were chosen, as only those two included the tongue's full extent. Despite the difference in pixel spacing between the SAR images and the UAV's mosaics, the reference window size remained the same (101×101 pixels, ∼17×17 m), together with a larger search window size of 501×501 pixels (∼85×85 m) to accommodate the glacier's movement during the 4 days 20 between the acquisitions.

Terrestrial Radar Interferometer Data
In order to reference flow velocity and geometry information with high spatial and temporal resolution at the calving front, a terrestrial radar interferometer (TRI; Caduff et al., 2015) was set up on stable ground 5 km south of the calving front with an unobstructed view of the glacier. The TRI system used was a GAMMA Portable Radar Interferometer (GPRI; Werner et al., 2008), operating at Ku-band (17.2 GHz). The device operates as a real-aperture radar interferometer, having one transmitting and two receiving antennas. The glacier was scanned at 1 min intervals for 8 consecutive days. Occasional data gaps were caused by hard-drive issues. The resulting radar intensity and phase measurements were processed further by application of 5 a standard workflow to determine the displacements in line-of-sight (e.g. Caduff et al., 2015;Lüthi et al., 2016). According to Voytenko et al. (2015) the absolute velocity errors were <0.5 md -1 with averaging times of ∼1 hour even for distant points in a humid atmosphere, with a range resolution of 0.75 m and a linearly scaling azimuth resolution of 35 m at 5 km distance (Werner et al., 2008).

10
A typical example of a velocity field derived from a 12 day Sentinel-1A repeat acquisition in August 2016 is shown in Fig.   4. In the main upstream tributary trunk of the glacier flow speeds were in between 1 and 2.5 md -1 , decreasing to zero towards the lateral margins. Within the last 5 km towards the calving front, the flow in the center line strongly increases to maximum values reaching up to 7 md -1 but with a rapid decrease towards the lateral margins. Within the last few hundred meters of the calving front, velocities dropped again substantially due to image templates containing mixed information from the glacier and 5 the sea (or ice-mélange in winter), leading to erroneous underestimations.
Note that the velocity field in Fig. 4 is from one single 12 day pair and was only filtered for outliers (magnitude and direction) and downscaled to 100×100 m but no further smoothing was applied, which explains the data voids and somewhat noisy appearance.

Comparison of GPS data to satellite-derived velocities 10
Using the GPS based flow velocity measurements as a first ground control dataset, we were able to assess the accuracy of the offset tracking approach applied to satellite products. In a first comparison, the flow velocities derived from spaceborne SAR satellites (cf. Table 1, example in Fig. 4) were plotted against the mean GPS flow speeds of the corresponding periods, as shown in Fig. 5. A total of 13 SAR derived velocity estimates were compared with the GPS derived flow velocities, five periods for GPS tracker EG09 and four periods each for GPS trackers EG14 and EG19 (cf. Table 4). To improve the SNR of the intensity tracking estimates, the integration time of the velocity estimates based on S-1A was doubled to 24 days for the two GPS stations located in the upper part of the glacier (i.e. EG14 and EG19) where the velocities are relatively low (1-1.5 md -1 ).
In addition to the improvement in SNR and thus a decrease in pixels flagged as outliers (cf. Fig. 4), this allowed for direct comparisons with results from RS-2 which also had a 24 day orbital repeat. The GPS velocities were averaged over the same 5 periods as the satellite observation interval. Comparing the flow velocities, they showed generally good agreement (R 2 : 0.674), with a mean relative difference of 11.5% between the two methods and an RMSE of 0.202 md -1 . When comparing the average of both, the GPS measurements and SAR derived flow velocities, for the whole campaign duration (about two months), the mean relative differences were reduced to 9.1% for EG09, 8.4% for EG14, and 5.4% for EG19.

Comparison of satellite-derived to UAV-derived velocities 10
The UAV-derived velocity field (2016/08/21 and 2016/08/25, Table 3), shown in Fig. 6, covered the lower fast flowing 5 km of the glacier tongue and in general showed a very similar spatial pattern in flow speed to the SAR-derived velocities. The main differences were that the flow field was much smoother and that it was additionally able to significantly better resolve the acceleration towards the terminus where flow speeds of 12 md -1 are reached. These discrepancies can be attributed to the much higher spatial resolution and hence differences in patch size. This improved spatial resolution even allowed resolution of 15 discontinuities in flow speed near the calving front related to deep crevasses and rifts. The UAV-derived flow field ( Fig. 6) also confirms the decreased but non-zero flow along the orographic left margin on the tongue already indicated in the SAR data ( Fig. 4).
In order to quantitatively (pixel by pixel) compare the flow speeds from the SAR with the UAV, the resolution differences between the radar and eBee datasets needed to be resolved. Therefore, the flow velocity map based on the UAV data was downscaled to match the 100 m grid size of the SAR based result. A mask was then applied for the areas outside of the glacier tongue to exclude stationary (e.g. moraine) and incoherently moving areas (e.g. open water) from the statistics (cf. black outline 5 in Fig. 7). The mask was manually traced using a Sentinel-2 scene acquired on 2016/08/03 as a reference.
The comparison between the two methods showed generally good agreement in most parts of the glacier tongue with differences between ±1 md -1 . However, a discrepancy was observed near the calving front as well as at the side of the glacier, where differences exceeded ±2 md -1 . Due to these larger differences, the standard deviation was 1.915 md -1 with a mean difference of -0.689 and a median of -0.281, showing a shift towards higher flow velocities based on the UAV data. These differences were 10 expected, as the initial 101×101 pixels (i.e. 252.5×252.5 m) template size crossed the glacier's boundaries at the calving front as well as at the sides, resulting in miscorrelations and therefore incorrect flow speeds from the SAR data, in particular close to the calving front. When possible border regions were excluded using a 250 m buffer around the glacier mask, the statistical values improved marginally (cf. red line in Fig. 7), resulting in a standard deviation of 1.576 md -1 and mean and median values of -0.626 and -0.317 md -1 respectively.
One should note that the time-periods of data acquisition were not identical (4 days UAV vs. 12 days SAR), but temporal variations within these time-scales reached their minimum at upstream GPS EG09 and were only substatial within 500 m of the terminus in the continuous TRI data (personal communication Andrea Walter, 26 Oct 2018).

Comparison of UAV-derived to interferometrically measured velocities 5
The TRI derived velocity map of the lower tongue was not continuous due to shadows from topography with respect to the TRI line of sight. The general spatial patterns were very similar to the UAV data, confirming the spatial gradients towards the lateral margins as well as the strong rapid step-wise acceleration within the last 1 km towards the front. The flow at the calving front is now also very well resolved and there maximum speeds in line of sight of 14 md -1 are reached.
Again for a detailed quantitative comparison of velocity measurements from the TRI measured interferometrically with those 10 derived from the UAV's hillshaded DSMs, the UAV-derived velocity field was first projected into the line-of-sight direction relative to the TRI's position, resampling the data onto the same grid. The comparison of these projected estimates and the interferometrically derived TRI flow velocities (cf. Fig. 9) show a close correspondence of the data for all areas aside from the very front and the lateral edges of the UAV's acquired area. The frontal differences can be explained by the differences in the measurement techniques, measurement times and spatial resolutions used. The image matching algorithm used with the UAV 15 data relies on features visible in both the reference and search scenes. Due to calving events, features at the glacier front may no longer exist in subsequent acquisitions, causing phantom image correlations. As the TRI relies on direct interferometric measurements at 1 minute intervals, it is less susceptible to errors caused by calving. As mentioned above, the point density of the drone's DSM was reduced at the margins of the acquired area, resulting in error-prone image matching results in those areas. Note that the non-zero ice flow at the orographic left lateral margin (in contrast to the right margin) was again very well reproduced by both acquisition methods.
Inspecting the overlapping areas between the two sensors within the glacier's marginal boundary (black line in Fig. 9), eBee 5 showed a mean velocity of 1.982 md -1 vs. 2.62 md -1 for the TRI. The mean difference between the two datasets was 0.633 md -1 ; the median difference was -0.007 md -1 (cf. Fig. 9). The standard deviation was 2.9 md -1 . After applying a 250 m buffer around the glacier's margin (cf. red line in Fig 9), the statistical values changed to a standard deviation of 1.844 md -1 and mean and median differences of 0.372 and 0.129 md -1 respectively.

10
The results and comparisons to other data presented here show the feasibility but also the limitations of operational offset tracking using Sentinel-1 intensity data to estimate glacier flow dynamics at relatively high spatial (100×100 m) and temporal sample intervals (12 to 24 days). While the applicability of interferometric approaches has been demonstrated in the past, its applicability is limited by temporal decorrelation in cases of fast movements which is particularly the case on tongues of tidewater outlet glaciers. The S-1/GPS comparison highlights the validity of the feature tracking approach with a mean velocity 15 difference of 11.5% between the datasets, agreeing with the reported differences of 9.7% by Ahlstrøm et al. (2013). Further, in comparison to the UAV-derived velocities (which agrees with TRI), our approach is able to represent the spatial pattern of ice flow towards the fast-flowing calving front well. Specifically, the main acceleration in flow towards the front as well as gradients to the margins are well reproduced. However, some edge effects mostly at the calving front remain, due to overlaps of the templates used in the cross-correlation with non-glaciated areas, but this effect is limited to a narrow zone determined by 5 the chosen template dimension and is mostly filtered out when considering flow directions as well. The non-zero velocity at the orographic left side of the terminus is reproduced surprisingly well by our 12 day SAR estimates, indicating that the inability of deducing strong velocity gradients within template patches mentioned by Nagler et al. (2015) does not impact our velocity results substantially or is compensated enough by the chosen smaller patch sizes. The good representation of spatial velocity gradients implies that strain rate fields are also robust which is crucial for constraining ice flow or calving models or process 10 studies.
However, strong changes in backscatter occurring between two acquisitions (e.g. due to changes in temperature causing surface melt or precipitation) can cause a deterioration of the results, resulting in data voids after the outlier detection. This was especially the case for the acquisition of 2016/08/25 that was within a warm period with temperatures never falling below freezing for a week, while the corresponding acquisition on 2016/08/13 occurred during a period of pronounced diurnal variations of the temperature. Detection of small flow velocities (<2 md -1 ) using our offset tracking from SAR intensity data can be impaired by low SNR, resulting in unreliable and noisy velocity estimates. An increase of the time between acquisitions for slower parts of the glacier, for example doubling the period to 24 days, can alleviate these issues, at the cost of reduced temporal resolution.

Influence of sample interval 5
To improve the representation of the velocity gradients even close to the glacier's terminus, a short sample interval followed by a downsampling step to e.g. 100×100 m is beneficial, although at a higher computational cost. When choosing too large a sampling step size (e.g. 40×40 px, i.e. 100×100 m), large gradients in flow velocity (such as areas close to the calving margin) might not be resolved, whereas a fine sample interval increases the chance of calculating offsets right up to the calving front without the reference window overlapping into the fjord area (cf. Fig. A9). Given this paper's focus on validation of velocity estimates and spatial patterns using high-resolution reference data, the increased computational cost caused by a 2×2 pixel sample interval was acceptable. For processing at e.g. ice-sheet scale, choosing a coarser sample interval (e.g. 20×20 or 40×40 pixels) at the cost of some detail is advisable.

Uncertainties and sensitivities to patch size and acquisition period
As reported by Nagler et al. (2015) three main sources of error exist when using offset tracking for ice velocity estimation:

15
-Errors in the matching procedure, errors due to ionospheric disturbances, geocoding errors.
Errors in the matching procedure are not only influenced by the image pair's co-registration and the quality of the amplitude features, but also by the chosen size of the template (Nagler et al., 2015). This choice is not straightforward, as its optimality 20 depends on the presence and prominence of features within. A smaller patch size can produce better results in regions with strong velocity gradients, yet suffer from increased noisiness. In contrast, a patch size that is too large might cause a blurring of the velocities, resulting in a trend towards lower values. This issue is investigated along a center profile in Fig. 10, corresponding to spatial extents of about 150×150 m, 250×250 m, and 350×350 m. A larger template size substantially reduces the noise, but increases the area affected by edge effects such as the deceleration artifact at the calving front. 25 Furthermore, the temporal integration time of subsequent intensity trackings influences the results. While shorter integration times suffer from a higher noise level compared to products averaged over longer periods, they can be used to investigate shortterm changes in flow dynamics of a glacier. This is illustrated in the sensitivity analysis of Fig. 11 showing SAR velocity results along the center flowline for different temporal averaging window sizes. Temporal averaging over several 12 day acquisitions substantially smooths the data and in particular reduces the artefact of velocity reduction at the calving front. Regarding the 30 geocoding process, i.e. the transformation from slant range to map projection, the errors introduced are primarily caused by DEM inaccuracies, as the geolocation accuracy of the S-1A/S-1B products has been shown to be well within the 7 m absolute location accuracy requirements specified by the European Space Agency (Schubert et al., 2017;Miranda et al., 2018;Piantanida et al., 2018). The DEM used for the project was the Greenland Ice Mapping Project (GIMP) DEM generated from data nominally from around 2007 at an original resolution of 30×30 m. We oversampled to 2.5×2.5 m. Since then, the slope and shape of glaciers has changed due to rapid thinning over the last decade.  Rizzoli et al., 2017), the surface lowering rates were 6 ma -1 close to Eqip Sermia's calving front and 2 ma -1 at 17 km from the terminus, amounting to a maximum of ∼70 m of horizontal location error or less than one pixel in our product. Changes in slope and shape of the glacier need to be accounted for as well when comparing threedimensional flow magnitudes or flow velocities assuming surface parallel flow, as they can introduce biases in the estimated magnitude of surface velocities (Nagler et al., 2015). Using the calculated lowering rates stated above, this results in a surface 10 slope change of 0.15 degrees in 9 years. This value is almost identical to the one reported for Jakobshavn Glacier by Nagler et al. (2015), reporting a bias of 0.5% in the magnitude of surface velocity caused by surface lowering. Further, as Eqip Sermia is a fully grounded glacier, errors due to tidal influences can be neglected.
As we relied solely on the feature tracking methodology for the reasons explained in Sect. 2.2, errors are expected to be larger compared to those derived from interferometric-based approaches alone (Strozzi et al., 2002;Short and Gray, 2014). When comparing the derived mean error of 11.5% between our GPS measurements and the velocities derived from intensity tracking, our findings agree with the 9.7% difference reported by Ahlstrøm et al. (2013) using intensity tracking. Similar errors were found when comparing the SAR derived velocities to the high-resolution UAV data over the glacier tongue, with a mean and 5 median difference of 12.4% and 8.5%, respectively. The differing spatial resolutions of the ground truth data sets used (i.e. UAV, TRI, GPS) resulted in additional uncertainties introduced during the resampling of the data to match the spaceborne acquisitions.
An analysis of the offsets in the easting and northing directions calculated over stable, non-moving terrain north and south of Eqip Sermia showed stable results with a mean velocity of <0.01 md -1 and a standard deviation of <0.3 md -1 in both the easting 5 and northing directions over the 13 month period corresponding to the time span of the Greenland Ice Sheet CCI product (Nagler et al., 2015) (cf. Fig. 12).

Comparison to operational ice velocity products
Due to the large spatial coverage of currently available, operational ice flow velocity products such as products available from the Greenland Ice Sheet CCI (Nagler et al., 2015) and the National Snow & Ice Data Center's MEaSUREs Greenland Ice Sheet 10 Velocity Map (Joughin et al., 2010(Joughin et al., , 2015(Joughin et al., , updated 2018, these products are only available for specific glaciers and for specific times at a high temporal resolution and do not cover Eqip Sermia. For our observation period, the monthly MEaSUREs Ice Velocity products are available at a spatial resolution of 200×200 m, while the Greenland-wide ice velocity map from the Greenland Ice Sheet CCI is only available on a yearly basis with a grid spacing of 500×500 m. In order to avoid effects from differing time periods, we compiled time-averages over the Greenland Ice Sheet CCIs 13 months time period (2015/10/01 -15 2016/10/31) (Nagler et al., 2015) based on the velocities from all available 40 Sentinel-1 image pairs and 13 monthly MEa-SUREs ice velocity products (Joughin et al., 2015(Joughin et al., , updated 2018. Both operational products differ significantly from the flow velocities calculated in this study, both along the center flowline of the glacier (cf. Fig. 13) and along the two cross profiles (cf. Fig. 14 and 15). In addition to different processing methods, these differences are also likely a result of spatial smoothing in the operational products, and are most strongly pronounced towards the fast flowing frontal part of the glacier, where strong 20 velocity gradients occur (cf. Fig. 13). There, and in contrast to our SAR-and UAV-derived velocities, the operational products indicate a slight or substantial deceleration in the vicinity of the terminus (cf. Fig. 13), likely an artefact from boundary effects Smaller templates are better in capturing the velocity gradients occurring towards the glacier front, resulting in slightly higher flow velocities.
Furthermore, the area affected by overlaps with surrounding areas of the glacier is diminished, resulting in reliable values closer to the glacier front. Bigger template sizes tend to result in smoother results. and smoothing. The operational product appears to underestimate the velocities of the main tongue up to 3 km behind the front by about 10% to 20% as an effect that spans the whole width of the tongue (cf. Fig. 14). Our UAV-derived velocities confirm this underestimation over the tongue and are even slightly higher than our SAR-results which may be an effect of the 25 different acquisition period over 4 days in the early Arctic summer. Further upstream, the discrepancies between SAR and the operational products generally decreases, but near the centerline our SAR estimates at the location of the lowest GPS (EG09, Fig. 15) were still significantly higher. These findings are also valid when comparing the monthly MEaSUREs product from August 2016 with time-averages based on our Sentinel-1A image pairs (cf. Fig A2-A6). Despite the differences in flow velocities between our maps and the operational products, there is good agreement between the different products on direction along 30 the flowline (calculated from the X-and Y-Offsets; cf. Fig. A8).
The above differences, calculated for the period between 2015/10/01 and 2016/10/31 (and similarly for the monthly MEa-SUREs product), imply that using the available, operational glacier flow velocity data sets for estimation of ice discharge with e.g. a flux-gate approach will result in an underestimation of ice flux (between 7% in comparison to MEaSUREs and 28% when compared to Greenland Ice Sheet CCI across the tongue, cf. Fig. 1), as such fluxes are calculated using surface velocity observations to approximate horizontal, depth-averaged ice velocity (Osmanoglu et al., 2013). This underestimation cancels out if the focus is for example on changes in ice flux over time (Howat et al., 2011;Rignot et al., 2008). However, mass budget  Fig. 1) for different temporal integration periods. Shorter periods are more prone to noisy results, but offer higher temporal resolution, while an increase in temporal integration time smoothens the results. methods take the difference between the absolute discharge and the surface mass balance integrated over the upstream catchment and an underestimation in flow then systematically underestimates mass loss. Considering that the mass loss from the 5 mass budget calculation is only a fraction (a few 10%) of the total discharge at the terminus (Enderlin et al., 2014;Rignot and Kanagaratnam, 2006), this 7% to 28% underestimation in near terminus ice flux would substantially affect mass loss estimates.
Of course this issue is less pronounced if the flux gates are located in the slower flowing parts upstream, but then an extra estimation of mass changes downstream is still required. Given the parameter settings used to produce the operational products (i.e. template window size, sampling step size, ground sampling distance), this underestimation in ice flow near the terminus, 10 may likely apply also to other similar medium-sized outlet glaciers and hence have an impact on mass loss estimates of the whole Greenland Ice Sheet.
Our near-terminus flow-fields will also imply higher frontal strain rates (compared to the operational products) which affects observational constraints for models of flow dynamics of calving and our understanding of terminus dynamics (Nick et al., 2009Choi et al., 2018;Morlighem et al., 2017;King et al., 2018). 15

Conclusions
The in situ measurements using a terrestrial radar interferometer (TRI), an unmanned aerial vehicle (UAV), and continuous GPS measurements used in this study were acquired during summer 2016 on a medium sized, ocean-terminating outlet glacier in West Greenland. This data was successfully used for validation of a standard SAR-based velocity offset tracking approach with a specific focus on resolving the fast-flowing terminus area. We could show good agreement in magnitude and spatial pattern 5 between multiple independent ground measurements and the flow velocities derived from different spaceborne C-band SAR sensors. The validity of the velocities derived using an iterative intensity tracking approach was demonstrated also for areas close to the glacier's calving front. Accurate near terminus velocity fields and related strain-rates are crucial for investigating the calving process, for constraining flow models and also ultimately for assessing ice flux and mass loss.
While there is a good agreement between the different data sets, caution should be exercised close to the glacier's margin, where the detection of its movements can be influenced by the chosen template size. While the size of this overlapping region can be minimized by choosing smaller template sizes, this introduces noise into the resulting velocity estimates due to erroneous correlations. Similarly, the upper boundary of the template size is limited by increased blurring, especially in regions where high 5 velocity gradients occur. These limitations could be mitigated using an iterative matching approach with increasing template sizes. Filtering based on long-term directional information seemed successful in removing outliers, albeit at the cost of data voids, in particular towards the glacier margins.
Finally, we were able to demonstrate the feasibility of our offset tracking approach using spaceborne SAR intensity data to derive glacier flow velocities much closer to the calving front than standard operational products. Given the short orbitalrepeat and illumination independency of Sentinel's SAR sensor, our approach has excellent potential for quasi-continuous operational derivation of accurate flow speed near calving termini and thus to provide velocity time-series for the analysis of terminus dynamics and ice sheet mass changes.   Fig. 1) for different sampling step sizes derived from S1A scenes acquired between 2016/07/08 and 2016/08/25 (48 days), starting at the glacier's terminus.