Glacier Surface Velocity Retrieval Using D-InSAR and Offset Tracking Techniques Applied to Ascending and Descending Passes of Sentinel-1 Data for Southern Ellesmere Ice Caps , Canadian Arctic

The Terrain Observation by Progressive Scans (TOPS) acquisition mode of the Sentinel-1 mission provides a wide coverage per acquisition with resolutions of 5 m in range and 20 m in azimuth, which makes this acquisition mode attractive for glacier velocity monitoring. Here, we retrieve surface velocities from the southern Ellesmere Island ice caps (Canadian Arctic) using both offset tracking and Differential Interferometric Synthetic Aperture Radar (D-InSAR) techniques and combining ascending and descending passes. We optimise the offset tracking technique by omitting the azimuth offsets. By doing so, we are able to improve the final resolution of the velocity product, as Sentinel-1 shows a lower resolution in the azimuth direction. Simultaneously, we avoid the undesired ionospheric effect manifested in the data as azimuth streaks. The D-InSAR technique shows its merits when applied to slow-moving areas, while offset tracking is more suitable for fast-moving areas. This research shows that the methods used here are complementary and the use of both to determine glacier velocities is better than only using one or the other. We observe glacier surface velocities of up to 1200 m year−1 for the fastest tidewater glaciers. The land-terminating glaciers show typical velocities between 12 and 33 m year−1, though with peaks up to 150 m year−1 in narrowing zones of the confining valleys.


Introduction
Remote sensing is becoming an increasingly important tool in climate research [1], with Synthetic Aperture Radar (SAR) being one of the techniques that has experienced a faster growth.The European Space Agency (ESA) has long supported various satellite Earth Observation missions with SAR sensors onboard such as the European Remote Sensing (ERS-1 and ERS-2) and ESA Environmental Satellite (ENVISAT) satellites.Continuing this long-term policy of providing continuous and consistent observational data, Sentinel-1A was launched in 2014, with its twin satellite, Sentinel-1b, launched in 2016.Both are equipped with a C-Band sensor [2].The Sentinel mission is comprised of a set of three different observational platforms intended to provide measurements that will contribute to better understanding the Earth System [3].The Sentinel mission is part of the Copernicus program (formerly known as Global Monitoring for Environment and Security, GMES), which comprises three components: space, in situ and services.The latter component is in charge of defining the Sentinel-1 observation scenario [4].The space component is intended to provide free of charge easily accessible data for registered users [5].The Sentinel-1 mission maintains continuity with key characteristics of former GMES missions (e.g., ERS and ENVISAT) including stable and accurate, well-calibrated products while at the same time enhancing other key mission aspects such as data reliability, revisit time and geographic coverage [6].A major improvement of the Sentinel-1 mission has been the reduction of the revisit time of the sensor, which has dropped from 35 days for ENVISAT to 12 days for Sentinel-1A, and six days when considering the full Sentinel-1 constellation.Furthermore, the coverage has increased from 80 km (for image mode in ERS and ENVISAT) to 250 km in the case of Sentinel-1 with its interferometric wide-swath (IW) mode (Table 1).The Terrain Observation by Progressive Scans (TOPS) acquisition mode of Sentinel-1 allows for a substantial increase in geographical coverage at the expense of the increasing difficulty associated with co-registering overlapping bursts.A small co-registration error can produce azimuth phase ramps.This problem is solved by applying a spectral diversity method over the overlapping areas of contiguous bursts [7,8].This type of SAR acquisition implies the antenna steering from aft to fore during burst acquisition leading to an azimuth variant Doppler centroid [9], which is corrected with a deramping operation of the Doppler spectrum.
There are several examples in the literature that apply the interferometric technique to TOPS acquisitions [10,11].The application of interferometry to ascending and descending passes is also well documented for glacier surface velocity estimations [12][13][14].There are plenty of articles dealing with glacier velocity retrieval using the offset tracking technique [15][16][17][18][19], including the particular case of Sentinel-1 TOPS acquisitions [20,21], the latter reference providing a comprehensive velocity field for the entire Greenland Ice Sheet.Furthermore, the combination of both techniques has also been covered with different approaches (e.g., merging of results or using least squares [22][23][24][25]) and also for entire ice sheets, such as the study by Mouginot et al. [26] for the Antarctic Ice Sheet.Finally, Hu et al. [27] give a holistic approach on the state of the art of SAR techniques for surface deformation retrieval.
The aim of the paper is to show the performance of the Sentinel-1 SAR Single Look Complex (SLC) Interferometric Wide-swath level-1 product (with a nominal resolution of 5 and 20 m in the across-and along-track directions, respectively) for retrieving glacier surface velocities using offset tracking and differential interferometry (D-InSAR) techniques.We estimate the surface velocity using differential interferometry where this technique proves to be more useful, namely, land-terminating glaciers and ice caps characterised by a low amount of deformation, which prevents decorrelation.The offset-tracking technique, in turn, proves to be more useful for fast-flowing glaciers.We apply the latter technique to Sentinel-1 imagery from ascending and descending passes [16], with the particularity of avoiding the use of azimuth offsets for such velocity estimates.The reason behind this is that Sentinel-1 TOPS IW acquisitions have their worse resolution in the azimuth direction.Additionally, we avoid the ionospheric effect strongly affecting this direction [28,29].By using this approach, we try to tailor the existing techniques to the particular characteristics of Sentinel-1, while avoiding the cumbersome corrections of ionospheric disturbances.Furthermore, the Sentinel-1 acquisition plan often includes ascending and descending passes from the same glaciated areas (Nunavut, Alaska, Himalayas, etc.), therefore making this approach suitable to be readily used in those cases.We selected as a test site the southern Ellesmere ice caps, Nunavut, in the Canadian Arctic, because it combines both land-terminating and tidewater glaciers, and presents wide glacier-free areas that allow a finer co-registration between acquisitions.

Study Area and Data Sources
Ellesmere Island is located in the Canadian High Arctic, neighbouring the western coast of Greenland (Figure 1).Ellesmere is part of the Northern Canadian Arctic Archipelago (NCAA) whose glacierized area in 2000 was around 104,000 km 2 [30] and contains one-third of the global volume of land ice outside Greenland and Antarctica [31].The climate characteristics of this location are low amounts of precipitation during the winter time (300-600 mm year −1 ) and increasing surface temperatures linked with a shift to a dominantly easterly positioned July vortex associated with an increased frequency of tropospheric ridging over the Canadian High Arctic [32,33].This implies that the surface mass balance (SMB) of the NCAA is governed by precipitation, in the form of snow, and increasing meltwater runoff accompanied by a shrinking capability of the snowpack and the firn layer to buffer mass loss through refreezing of percolating meltwater [34].Surface mass budgets in this region are currently showing an increasing negative trend [37].Between 2005 and 2009, CAA ice caps suffered the most negative trends in their observational records [38].These trends show that the SMB in the NCAA is moving from a mass loss predominantly dominated by calving during the 20th century to a contemporary mass loss driven by meltwater runoff [34,[38][39][40].
Studies focusing on the regional glacier dynamics revealed that surge-type glaciers, with a relatively long surge phase, are common in the Canadian High Arctic.Many of them are tidewater glaciers located in the high-precipitation coastal areas [39,41,42].Independent of this surge behaviour, glacier ice velocities during the summer are typically an order of magnitude larger than the mean annual velocities, and summer calving rates are 2-8 times larger than the annual average [41,43].Other studies point to a high interannual variability of the velocities of the large outlet glaciers [44].van Wychen et al. [39] showed that velocity changes associated with the termination or initiation of surges can rapidly change the rate of ice discharge to the oceans.Finally, many research papers indicate a clear increase in ice mass losses and associated contribution of Canadian Arctic glaciers to sea-level rise during recent years.This increase in mass loss is mostly caused by surface melt and runoff, and not by glacier dynamics [38][39][40][41]45].
We frame our study within this climatic and geographical context.Specifically, we focus on glaciers located in the Prince of Wales (POW) and the Manson Icefields, and some eastern glaciers of Sydkap Ice Cap (Figure 1).The glaciers in this area are both of land-terminating and tidewater types.The different nature and dynamics of these glaciers, together with the presence of wide ice-free areas (which allows performing a finer co-registration between acquisitions) makes this location optimal to test and apply different SAR surface velocity estimation methodologies.Furthermore, the low amount of precipitation in the area, together with the fact that all acquisitions correspond to northern hemisphere winter months, characterised by lower ice velocities, increases the chance of a successful application of various techniques.
The surface velocities on Ellesmere Island were obtained from a total of 14 Sentinel-1 SAR TOPS IW Level-1 Single Look Complex images taken during February and March 2016 (Table 2).These type of products provide a geo-referenced image (using accurate altitude and orbital information from the satellite), in slant-range geometry, processed in zero-Doppler.The image is normally composed of three sub-swaths with each sub-swath comprising normally nine consecutive bursts, which overlap in azimuth.Burst synchronisation is needed for interferometry and for accurate offset tracking [46].The data were retrieved in two different tracks, with six and nine scenes each, during ascending and descending passes, respectively.The incidence angle, θ, was in all cases close to 33 degrees, while the α angle between ascending and descending tracks was around 40 degrees.The time span between corresponding acquisitions was 12 days, which halves to six days when considered the ascending and descending passes.We used the horizontal transmit and horizontal receive (HH) channel, which preserves better the amplitude features and has also a higher signal-to-noise ratio as compared to the horizontal transmit and vertical receive (HV) channel, making the former better suited for glacier ice motion retrieval.The TOPS acquisition mode improves the performance of already existing SAR imaging algorithms such as ScanSAR mode [7,47].During a TOPS mode acquisition, the SAR antenna is steered in the azimuth direction from aft to fore with a constant rate.This mode of observation has several advantages; in particular, the measured targets are observed with the whole azimuth antenna pattern, which also reduces the scalloping effect [47].The main disadvantage is that it imposes a more restrictive approach for the co-registration procedure and for the interferometric processing.The peculiarities of the co-registration with the spectral diversity method have been addressed by Grandin [8], while the procedure for obtaining interferometry from TOPS mode images has been dealt with by various authors [8][9][10].
For the application of the differential interferometry to the SAR imagery, a fine digital elevation model (DEM) of the island was required.For this purpose, we used the freely accessible Canadian DEM (CDEM) designed by Natural Resources Canada (NRCan) with a resolution of 0.75 and 3 arc second in the south-north and east-west directions, respectively [48].The DEM was used to produce the differential interferometry results and for geocoding and co-registering the imagery employed for the intensity offset tracking technique.This is the best currently available DEM from the area and it dates back to the 1960s.We want to highlight here the degrading effects that this could have specially on the results of the D-InSAR technique [12].

GAMMA Remote Sensing Offset Tracking
We used GAMMA software [49] (GAMMA Remote Sensing AG, Gümlingen, Switzerland) for processing the SAR Sentinel-1 acquisitions.This software package is extensively used as a tool for retrieving ground deformations and surface displacements in many areas, including glacier-covered regions [18,50].We exclusively applied the intensity offset tracking algorithm.The speckle tracking technique has its own peculiarities when it is applied to TOPS products, in particular when applied to non-stationary areas such as glaciers or ice caps [20] (e.g., TOPS azimuth antenna scan entails phase deramping and azimuth common band filtering [9]).Orbital accurate ephemerides are becoming a common resource in new space missions, as is the case for the Sentinel-1 [21,51] mission, which allows for a seamless generation of wide regional surface displacement products.
Sentinel-1 TOPS mode images first need to be co-registered before any offset tracking or interferometric algorithms are applied on them.This fine co-registration procedure requires the use of a fine DEM within the same area of the acquisition image extent [49].The co-registration begins with obtaining the DEM in SAR coordinates followed by the matching algorithm and the spectral diversity method (which includes considering the interferometric phase) applied over the overlapping areas of the bursts.The co-registration requirements are quite stringent and a co-registration quality of 1/1000 of a pixel in the azimuth direction is required for the phase discontinuity between consecutive bursts edges to be less than three degrees [9].The co-registration performance can be improved when masking areas with expected azimuth displacements (e.g., glaciers or ice caps).This prevents fine co-registration methodologies to be applied to wide ice caps where no stable ground is found [20].This was a critical factor for choosing Ellesmere Island as a suitable area to apply a full co-registration procedure using the matching and the spectral diversity methods.When no interferometry is required and the only derived product is to be offset tracking displacements, the Sentinel-1 ephemerides are sufficient for comprehensive, region-wide results [21].
After a full co-registration is achieved, a deramping of the SLC images for correcting the azimuth phase ramp is required to apply oversampling in the offset tracking procedures [49].Once the above-mentioned steps are done, the offset tracking technique is the same as for normal stripmap mode scenes [15][16][17][18][19]. Surface displacements can be obtained in ground coordinates (e.g., slant rage and azimuth directions), which are finally geocoded using a lookup table derived from the use of the DEM and the image parameter file.
We used a matching window of 320 by 64 pixels (1200 by 1280 m) in range and azimuth directions, with an oversampling factor of two for improving the tracking results.The sampling steps were of 40 by eight pixels and the resolution of the final velocity map was 200 by 160 m in range and azimuth directions.The geocoding was completed using the Canadian Digital Elevation Model.We used a bicubic spline interpolation to generate the geocoded grid.This grid was transformed into a Tiff file using a functionality within GAMMA software.Finally, the determined velocities were manually checked for blunders and mis-matches.

GAMMA Remote Sensing D-InSAR
The application of interferometry and differential interferometry to TOPS mode data has been dealt with by various authors [8][9][10][11].The D-InSAR technique requires a full co-registration of the reference and slave scenes.For this purpose, we proceed as in the offset tracking case.Generation of wide area surface deformation products is possible thanks to highly accurate ephemerides products [51].However, caution should be taken when considering multiple frames along track, making sure to account for the variation of azimuth co-registration error in this direction [52].The use of stationary areas is paramount for achieving a successful co-registration.Therefore, we foresee limitations for fine co-registration methodologies in the absence of such areas.
When all of the above-mentioned conditions are met, interferometry is applied in the same way as conventional interferometry to stripmap imagery (i.e., phase filtering, masking and unwrapping).We used the minimum cost-flow algorithm for unwrapping the filtered phase.When calculating the differential interferogram by using a DEM, the derived differential interferometry should be carefully checked for the presence of phase jumps between subsequent bursts [9,49].For the computation of the differential interferometric product, we used sampling steps of 10 by 2 pixels (50 by 40 m) in the range and azimuth directions, respectively.The geocoding was completed using the Canadian Digital Elevation Model.We used a bicubic spline interpolation to generate the geocoded grid.This grid was transformed into a Tiff file using a functionality within GAMMA software.

Offset Tracking Case
In a study by other authors aimed at obtaining the full 3D displacement vector from ascending and descending passes using offset tracking [16], all of the offsets (both in range and azimuth) were used within a least-squares approach for calculating the velocity vectors.
We avoid using the azimuth offsets altogether, for two main reasons.First, the resolution of Sentinel-1 IW images in the azimuth direction is four times worse than in the range direction.Second, we avoid the ionospheric effect, which degrades substantially the performance of the offset tracking algorithm in the azimuth direction creating the so-called azimuth streaks [28,29].The ionospheric effect will have a stronger influence on Sentinel-1 datasets due to the shorter repeat cycle of the constellation (which implies less amount of glacier movement in the retrieved signal).As in the case of D-InSAR, we will obtain a 2D velocity vector by applying a transformation from a nonorthonormal basis to the corresponding geographical coordinate system [12].By doing so, we tailored the existing methodologies to the special characteristics of Sentinel-1 TOPS product.
The horizontal velocity vector v h can be expressed as a function of the ascending and descending across-track vectors (see Figure 1d): We identify the range displacements from the ascending and descending passes as the first and second terms on the right-hand side of Equation (1).At this stage, we used the intensity offset tracking range displacements obtained with GAMMA software and described in Section 3.1.The main problem in the transformation is that the basis (â, d) is nonorthonormal.Therefore, a relation between the projections of the horizontal velocity vector and the ad coordinates needs to be established.These relations are given in Equation ( 2) below, where the matrix form is also displayed: The first two expressions in Equation ( 2) are inferred manipulating the equations that describe the two-dimensional transformation between the ascending and descending coordinate systems [12].Next, we transform the projected velocity vector v h from the ad coordinate system to xy coordinates (Equation ( 3)): Finally, Equation ( 4) shows the application of the whole transformation from the nonorthonormal across-track vectors to the xy coordinate system: The α and β angles necessary to perform the transformation were obtained using the GAMMA software.The program provides a function to calculate the look-vector orientation producing an output grid with the SAR look-vector orientation as a function of position.This grid was geocoded back to the geographical coordinate system.Once the rasters of the look-vector direction from the ascending and descending images were obtained, we had all the necessary information to proceed with the computation of the velocity components.The raster algebra to transform the coordinate system ad to the geographical coordinate system xy was performed entirely using QuantumGIS software 2.18.2,Quantum GIS Development Team and an ad hoc code to automatize the processing with Python.The histograms and plots were completed using Octave (4.2.1, GNU project, Free Software Foundation).

D-InSAR Case
The retrieval of surface displacement from ascending and descending SAR pairs has been covered by several authors [12][13][14].In all cases, for retrieving a 3D full vector, an assumption of surface parallel flow was made.The downside of this approach is that ice sometimes does not flow parallel to the surface.However, when considering only the interferometric phase difference due to the horizontal displacement, one can derive the horizontal 2D velocity vector from ascending and descending passes.The important step in this procedure is the transformation from a nonorthonormal basis (the across-track ascending and descending vectors) to the geographical coordinate system [12].We proceeded in the same way as for the intensity offset tracking case.The D-InSAR across-track ascending and descending horizontal displacements were obtained using GAMMA software.The angles α and β were obtained in the same way as for the intensity offset tracking case.The raster algebra was done with QuantumGIS software.By using this approach, we do not readily get the full three-dimensional velocity vector, but the missing vertical component can be easily obtained afterwards using the surface gradient.Moreover, the vertical component is generally small.

Offset Tracking Velocity Mapping of Southern Ellesmere's Ice Caps
In general, our results show that the performance of the amplitude matching algorithm is good, and the method is able to resolve the glacier velocities in most areas.There is a large variability in the observed velocities, due to the different dynamic behaviour of the various glaciers.Some glaciers present a surge-type behaviour [42] or high seasonal and interannual variability of their velocities [43].We can broadly group our studied glacier into three groups: (1) fast flowing with winter velocities greater than 600 m year −1 ; (2) medium flowing with winter velocities within 50-200 m year −1 ; and (3) slow flowing with speeds less than 50 m year −1 .Trinity and Wykeham glaciers (Figure 2) belong to the fast flowing group with speeds up to 1200 m year −1 for Trinity and 600 m year −1 for Wykeham.To the medium flowing group belong glaciers such as Ekblaw Glacier (Figure 1c) or neighbouring glaciers such as Stygge or Cadogan (north and south of Ekblaw, respectively; not shown in figure), reaching winter velocities of 80 and 60 m year −1 , or the southern glacier of Manson Icefield, with speeds up to 75 m year −1 .This is currently the fastest glacier of the Manson Icefield because Mittie Glacier, which can have velocities greater than 1 000 m year −1 when in full surge [42], is currently in its quiescent phase [41].Most of the slow flowing glaciers are land-terminating glaciers located on the western side of the POW Icefield (Figure 3) and on Sydcap Ice Cap.However, some of the land-terminating glaciers in the western POW Icefield have larger velocities, as happens with unnamed West 3, 4 and 5 (Supplementary Materials Figure S1), with velocities up to 75, 150 and 100 m year −1 , respectively.These high velocities are due to the narrowing of the valleys confining the glacier flow.In the case of unnamed West 4, the velocity increase is more noteworthy because of its larger accumulation area.However, the feature tracking algorithm is not able to resolve the surface velocity in slowest-moving zones of the ice cap's accumulation areas (Figure S1) and in some zones of the Trinity glacier (Figure 2).In the case of the accumulation areas, this can be attributed to the absence of visible features due to the snow cover.In the case of the Trinity Glacier, the underlying reason is the strong surface velocity gradients that deform the surface features (e.g., at the junction between the main glacier trunk and its principal tributaries, at the corners where the glacier changes direction or at the glacier margins. The feature tracking algorithm is best suited for resolving areas of fast glacier movement.This is clear when comparing the feature tracking results with those from differential interferometry (Figures S1 and S2).However, there is also a limit on the amount of movement that the feature tracking algorithms are able to resolve.Nagler et al. [21], using an in-house developed algorithm, were able to resolve velocities up to 12 m d −1 when using Sentinel-1 IW data (the same sensor that we are using).With increasing image resolution (smaller pixel size), larger velocities could be resolved.Anyway, this restriction is not relevant to us; because 12 m d −1 is equivalent to 4380 m year −1 , a velocity well above the largest glacier velocities observed in the Canadian Arctic.
We turn now our attention to the performance of offset tracking using ascending and descending tracks without using the azimuth offsets.As expected, the areas of stable ground do not present the typical azimuth streaks that one would expect when using azimuth offsets [28,29].In the scenes that we used, the orientation of the ascending and descending across-track vectors was close to the horizontal (e.g., 12 • below the west-east direction for the ascending case).This makes the retrieval of surface velocities for glaciers oriented in the north-south direction more challenging.Nevertheless, glaciers oriented in this direction (e.g., most of the tributaries of Ekblaw Glacier) still present homogeneous surface velocities (Figure 1c).Furthermore, when analysing the surface velocity gradients of Ekblaw Glacier in the northern direction (Figure 1c), the better performance of our approach is noticeable, which avoids the use of azimuth offsets.We can observe how the northern component of the velocity increases and decreases depending on the direction of glacier flow.By contrast, when using the standard approach that uses azimuth offsets, we were unable to resolve the velocity gradients for this glacier.Please also note that the standard approach requires the correction of the ionospheric effect, which increases the complexity when considering Sentinel-1 TOPS acquisition mode.This is because each sub-swath presents its own ionospheric-induced azimuth pattern.When no information on the total electron content of the atmosphere is available, we need to have stable ground available in different sub-swaths, a condition not easily fulfilled in many glacierized areas.Furthermore, when using a modelled trend of the ionospheric disturbance for subtracting the ionospheric effect, residues that degrade the azimuth offset results always remain.

D-InSAR Velocity Mapping of Southern Ellesmere's Ice Caps
Glacier surface velocity results using ascending and descending passes show good performance on all land-terminating glaciers, such as the western glaciers of POW Icefield (Figures 3a and S2) and the eastern glaciers of Sydcap Ice Cap.Moreover, there are two tidewater glaciers that also show good results when applying this technique, namely the South Margin glacier (POW Icefield) and Mittie East and West arms (Manson Icefield), as well as their tributary glaciers.All of them present surface velocities below 50 m year −1 and smooth surface velocity fields.D-InSAR performs better than offset tracking in areas with low glacier velocities.Furthermore, it is able to resolve small surface velocity gradients.The low amount of precipitation in this region [32,33] allows for a successful application of D-InSAR in the region's ice mass accumulation areas.This is the case of POW Icefield, where we see a fine and continuous ice velocity field.It is precisely this fine resolution that allows us to discern the ice divide between the eastern and western flowing glaciers (Figure S2).We note that in some cases the Randolph Glacier Inventory (RGI) [36] ice divides do not match well with our surface velocity field estimates.This fact is noticeable e.g., in Taggart Lake glacier and in the unnamed West 4 and 5 glaciers (Figure S2).This fine resolution also allows us to recognize different glacier flows (tributaries) in some western land-terminating basins of the POW Icefield.Examples are Taggart Lake glacier, which shows three distinctive flows (one in its northern side and two in the southern one), Unnamed West 2 glacier and Unnamed West 4 glacier (Figure 3a).
Regarding surface velocities for individual glaciers, Unnamed West 1 reaches 24 m year −1 and Taggart Lake glacier reaches 33 m year −1 in their main branches.None of the three different ice

Comparison of Offset Tracking and D-InSAR Results
If we compare the results from offset tracking and differential interferometry (Figure 3), the more irregular velocity patterns produced by the offset tracking technique are noticeable (e.g., there is some degree of variability on stable ground).Additionally, we perceive a lower resolution of the technique, either due to the size of the matching window or to the quality of the amplitude features.This lower resolution prevents the offset tracking technique from producing the continuous and clearly defined glacier surface flow fields produced when the differential interferometry technique is used.However, the offset-tracking technique allows us to derive region-wide results with less constraints as compared with other methodologies (e.g., resilience to de-correlation or fast surface velocities).
Differential interferometry produces a high-quality, continuous and smooth velocity field.The main shortcomings of this technique include the inability of the minimum cost flow algorithm to resolve strong surface gradients and its deterioration with longer time spans between acquisitions (e.g., increased de-correlation).The former situation would be exemplified by the case of the glacier tongues of glaciers Unnamed West 3, Unnamed West 4 and Taggart Lake north, whose velocities cannot be retrieved with the D-InSAR algorithm (Figure 3).Offset tracking, on the contrary, is capable of capturing these type of velocity gradients, making this methodology especially useful to retrieve region-wide surface velocity fields.
The areas where the differential interferometry performs best are those with subtle transitions between parallel or adjacent surface flows or in accumulation areas with lower velocities.In such areas, we notice the difference between offset tracking and D-InSAR.While the former hardly resolves the velocity gradients (and, when it does, it produces a nonsmooth field; see the accumulation area of Taggart Lake glacier), D-InSAR results illustrate the capabilities of this technique in accumulation areas with low amount of precipitation, which helps to preserve the correlations between subsequent acquisitions.

Comparison with Previous Studies
We have produced regional velocity estimations for the Southern Ellesmere Island ice caps using intensity offset tracking and differential interferometry from Sentinel-1 acquisitions.What is special in our application of the intensity offset tracking technique is that we retrieved the surface velocity fields exclusively from range offsets using ascending and descending passes.We also applied differential interferometry for the same areas in order to establish a comparison on the performance of both methodologies in glaciers with different dynamic regimes [42,43].Sentinel-1 IW TOPS mode data was already tested with the offset tracking methodology by Nagler et al. [21] giving successful region-wide results for the ice velocity field of Greenland.The application of the D-InSAR technique to TOPS mode acquisitions from ice caps is still a challenge due to the need for stable ground to achieve a fine co-registration with the spectral diversity method, aimed to avoid undesired azimuth phase ramps in the resulting products [8][9][10].
The application of the intensity offset tracking approach to TOPS mode acquisitions was already dealt by Dall et al. [20], who acknowledged the need for azimuth deramping and azimuth common band filtering when applying speckle tracking.These two latter steps are not required in the feature tracking that we apply.The application of the intensity offset tracking technique to ascending and descending passes was covered by Fallourd et al. [16].However, their approach considered the use of all azimuth and range velocity components solving by least-squares, while our approach discards the azimuth offsets from the estimation to avoid undesired ionospheric effects and, simultaneously, to avoid Sentinel-1 shortcomings such as its lower azimuth resolution.We did this by transforming a nonorthonormal coordinate system (e.g., ascending and descending across-track vectors) to a geographical coordinate system [12].
The results of the application of the proposed method show that there is a noticeable improvement in the resolution of the components of the surface velocity together with an increase of the final precision.We avoid the ionospheric effects and therefore, considering the improved accuracy of the ephemeris, the remaining error sources are restricted to the co-registration, the matching procedure and the DEM related geocoding error.In addition, there are still some de-correlated areas in the accumulation zone devoid of results due to the lack of features for tracking.We would like to highlight that using azimuth offsets with Sentinel-1 TOPS data not only worsens the final velocity product but also increases the difficulty of correcting the azimuth streaks which present a different pattern for each sub-swath.
There are several articles devoted to the topic of interferometry using TOPS mode data [10,11].Furthermore, the use of ascending and descending passes for retrieving ice surface velocities is also well covered in the literature [12][13][14].The main obstacle that we faced when applying D-InSAR to TOPS mode data was to guarantee a successful fine co-registration using the spectral diversity method between contiguous bursts.The choice of Ellesmere Island as a test site was key for overcoming this limitation because of the wide areas of stable ground that are found next to the ice masses.
Our results show that differential interferometry is optimal for areas of low movement with no (little) surface change (e.g., accumulation areas with low amount of precipitation, and land terminating glaciers).The method produced smooth velocity fields for most land terminating glaciers of Ellesmere southern ice caps.On the other hand, the method failed to work in areas with strong velocity gradients and areas of fast glacier movement (e.g., the minimum cost flow (MCF) algorithm did not resolve a continuous velocity field [12]).The comparison between methodologies illustrates the difficulties of MCF to resolve the velocity field of a few tongues from the western glaciers of the POW Icefield.
Regarding the comparison with velocities presented in previous studies in the region, focusing on those based on scenes closest in time to ours, we note that the velocities are very similar, especially for the fastest glaciers.Our maximum winter (Feb-March 2016) velocities of 1200 m year −1 near the terminus of Trinity Glacier or 600 m year −1 for Wykeham Glacier are very similar to those of 1250 and 500 m year −1 (respectively) reported by van Wychen et al. [41] for winter 2014-2015 (from speckle tracking of RADARSAT-2 data), or those of 1200 and 650 m year −1 (respectively) given by Millan et al. [40] for winter 2015-2016 (from speckle tracking of Sentinel-1a data).Similarly, for Ekblaw Glacier, both our results and those of van Wychen et al. [41] show similar maximum velocities of ∼100 m year −1 .

Offset Tracking Case
The three main sources of errors in the offset tracking technique are: (1) the matching procedure (which is a function of the co-registration between images, template size, and the quality of the image features) [21]; (2) the ionospheric effect and its influence on the azimuth offsets in the form of azimuth streaks [28,29]; and (3) the geocoding error (with the high quality ephemerides of Sentinel-1 data, this error is reduced to the quality of the DEM used for the topographic correction [21]).
The first type of error is obviously unavoidable due to the intrinsic nature of the tracking algorithm.However, the ionospheric errors account for a big share of total velocity errors, providing a possible and effective way of reducing offset tracking errors.Sentinel-1 is equipped with a C-band sensor which is less affected by this disturbance than L-band sensors such as those on board of Advanced Land Observing Satellite-Phased Array type L-band Synthetic Aperture Radar (ALOS PALSAR) [29].L-band sensors, however, have higher penetration depth in snow and ice, which improves correlation between scenes and therefore also interferometric results, as has been analyzed by Rignot et al. [53].On the other hand, a Sentinel-1 repetition cycle of 12 days has the side effect of increasing the share of signal corresponding to ionospheric effects.Various authors have acknowledged the influence of the ionospheric effect on the offset tracking results.Nagler et al. [21] estimated this effect to account for ∼0.08 m day −1 displacement when considering the low velocity areas of Greenland.This estimate dropped to 0.02 m day −1 when averaging over several acquisitions.Other authors [26,29], when evaluating RADARSAT-1 offset tracking data from Antarctica, estimated this error to be between 0.016 and 0.042 m day −1 .
We performed an analysis of our results in Ellesmere Island, focusing on the performance of the algorithm on stable ground (e.g., on ice-free areas on the western side of POW Icefield) where both matching and ionospheric components of the error are present and can be quantified.We calculated the root mean square error (RMSE) of the normal intensity offset tracking results (e.g., including azimuth offsets and the related ionospheric effect) and also of our suggested methodology using exclusively range offsets from ascending and descending passes.The results show that the former approach casts an RMSE value of ∼0.051 m day −1 (the mean error being ∼14.2 m year −1 ), while the latter drops to 0.012 m day −1 (the mean error being ∼3.3 m year −1 ) (Figure 4).If we consider that both sources of error are independent, the share of error that could be attributed to ionospheric disturbances in the first case amounts to ∼0.05 m day −1 , a value falling within those of Nagler et al. [21] and Gray et al. [29] estimates.We emphasize that, when using only range offsets extracted from ascending and descending passes with Sentinel-1 data, errors represent a fourth of the intensity offset tracking technique error budget.Furthermore, our technique still improves the multiple acquisition averaging, by halving the error of the latter.

Differential Interferometric Case
The InSAR technique has several error sources.As for intensity offset tracking, the co-registration is a typical error source.For instance, the co-registration of the ascending and descending scenes using the DEM could have significant errors.These errors could be due to a low quality DEM, to shortcomings in the along-track timing or to inaccurate ephemeris [12,54].The largest share of the error corresponds to the interferometric baseline definition, whose influence can be observed as linearly varying errors.Sentinel-1 data is no longer affected by these two types of orbit-related errors because of its high quality orbital data [21,51].The phase unwrapping algorithm plays a role in a correct and continuous definition of surface velocity gradients.
In differential interferometry, a DEM is used to remove the topographic phase component from the interferogram.Therefore, a low quality DEM combined with long baseline acquisitions can also lead to worsened velocity estimates [14].
Most glaciers of this region have directions approximately perpendicular to the along-track satellite direction (especially all the western POW Icefield glaciers).This optimizes the application of interferometry to these glaciers.All the above mentioned elements contribute with different shares to the total error budget.This makes it difficult to give a proper error estimate.Nonetheless, Luckman et al. [55] gave an estimate of the error between 0.012 and 0.017 m day −1 for the specific case of ERS-1 satellite data applied to the Himalayas.We consider our performance to be better both because of Sentinel-1 improved orbital data and to the more gentle topography of Ellesmere Island and its positive impact on the DEM.

Conclusions
The proposed intensity offset tracking algorithm applied to ascending and descending passes, which disregards the azimuth offsets, has demonstrated its merits to resolve velocity gradients in the along-track direction, improving the resolution of the offset tracking and, at the same time, avoiding the ionospheric disturbances on the data.In this way, we have improved the accuracy of the derived surface velocity product, whose main source of error is reduced to the matching-related uncertainties.The D-InSAR interferometry has also shown a good performance in slowly moving areas, providing high resolution, smooth and continuous velocity fields for land-terminating glaciers.
In view of these results, we acknowledge the complementarity of both techniques [27,55].Offset tracking can be applied in areas of faster glacier movement and velocity gradients, with some limitations, while differential interferometry is suited to areas with low velocities and small velocity gradients.Moreover, offset tracking is less affected by de-correlation than differential interferometry.We see here an opportunity to develop a hybrid velocity product combining DInSAR and offset tracking results in regions where either one method or the other, or both, perform best, following the lines suggested by Joughin [22] and Liu [23].The recently deployed Sentinel-1B satellite will improve the interferometric capabilities of the Sentinel-1 constellation by reducing the de-correlation and the amount of movement between consecutive acquisitions [12].Furthermore, we see a possible application of the ascending and descending passes to speckle tracking.This approach would allow for obtaining surface velocity fields within featureless accumulation areas, while avoiding the ionospheric effects typical of Arctic and Antarctic areas [21].

Figure 1 .
Figure 1.(a) Main ice masses of Ellesmere Islands within the Arctic [35].The glacier outlines are from the Randolph Glacier Inventory (RGI) version 5.0 [36]; (b) zoom of red rectangle in panel a, showing our three main study areas (I, II and III) in Southern Ellesmere; (c) northern velocity component for Ekblaw Glacier obtained using intensity offset tracking.The zone shown corresponds to the red rectangle labelled I in panel b.The background images in this panel c were taken between 10 and 17 February 2016; (d) schematics of the vectors and angles involved in the non-orthonormal estimation of the ice flow velocity.All maps in this paper use the Universal Transverse Mercator (UTM) projection of the zone 17 north, and the reference ellipsoid is WGS84.

Figure 2 .
Figure 2. Glacier ice velocities for Trinity and Wykeham glaciers obtained using intensity offset tracking applied to ascending and descending passes.The zone shown corresponds to the red rectangle labelled II in panel b of Figure 1.Velocities calculated from images acquired between 22 February and 12 March 2016.Velocity map overlaid on top of Sentinel-1 Synthetic Aperture Radar (SAR) images acquired between 22 and 29 February 2016.

Figure 3 .
Figure 3. (a) Glacier ice velocities for the western glaciers of Prince of Wales Icefield obtained using Differential Interferometric Synthetic Aperture Radar (D-InSAR).The zone shown corresponds to the red rectangle labelled III in panel b of Figure 1; (b) Centreline velocities for longitudinal section A-A' in panels a, c; (c) ice velocities for the same zone of POW Icefield obtained using intensity offset tracking; (d) centreline velocities for longitudinal section B-B' in panels a, c.The velocity maps are overlaid on top of Sentinel-1 SAR images taken between 10 and 17 February 2016.

Figure 4 .
Figure 4. Glacier surface velocity errors from a traditional intensity offset tracking technique (average 14.7 m year −1 ) and from intensity offset tracking applying ascending and descending passes (average 3.3 m year −1 ).

Table 1 .
Sentinel-1 Beam modes with incidence angles, nominal resolution and scene size.

Table 2 .
Sentinel-1 acquisitions used in this study.

Acquisition Date Beam Mode and Number Orbit Number Slice Number Pass Direction Methodology
-InSAR velocity mapping for southern Ellesmere Ice Caps.