Edinburgh Research Explorer Intercomparison and validation of SAR-based ice velocity measurement techniques within the Greenland Ice Sheet CCI project

: Ice velocity is one of the products associated with the Ice Sheets Essential Climate Variable. This paper describes the intercomparison and validation of ice-velocity measurements carried out by several international research groups within the European Space Agency Greenland Ice Sheet Climate Change Initiative project, based on space-borne Synthetic Aperture Radar (SAR) data. The goal of this activity was to survey the best SAR-based measurement and error characterization approaches currently in practice. To this end, four experiments were carried out, related to different processing techniques and scenarios, namely differential SAR interferometry, multi aperture SAR interferometry and offset-tracking of incoherent as well as of partially-coherent data. For each task, participants were provided with common datasets covering areas located on the Greenland ice-sheet margin and asked to provide mean velocity maps, formulated, in particular concerning procedures which can reduce the impact of analyst decisions, and which are often found to be the cause of sub-optimal algorithm performance.

achievable spatial resolution of the measurements in each image dimension is in the order of several radar resolution cells, which in turn are between 2 m and 20 m in each image dimension for SAR acquisition modes relevant for ice-sheet studies. From the processing point of view, DInSAR requires image co-registration and 2D phase unwrapping, which are time-consuming and sometimes delicate steps for non-stationary and/or partially coherent scenes, particularly for imagery acquired in Terrain Observation by Progressive Scans (TOPS) acquisition mode ( [8,9]).
Offset-tracking methods measure both the LoS and the azimuth components of the displacement between two SAR acquisitions by estimating the mis-registration of corresponding points on ground due to motion. This is achieved by tracking incoherent amplitude or intensity features (incoherent offset-tracking) and/or coherent SAR speckle (coherent offset-tracking) in several ways. Amplitude or intensity cross-correlation (ACC and ICC in the following) [10], as well as translation estimation in the Fourier domain via phase correlation [11], are able to exploit both the presence of features and of correlated speckle, typically yielding results on the featureless, but possibly coherent, ice-sheet interior as well as on the crevassed, but often incoherent, glacier tongues. Improvements in terms of accuracy and spatial resolution can in some cases be obtained using coherence maximization [12] or complex cross-correlation (CCC) [13,14], although these methods are only applicable in practice to areas/datasets in which a high level of phase coherence is retained. The achievable accuracy and spatial resolution of all offset-tracking techniques is typically an order of magnitude worse than DInSAR. From the processing point of view, although these methods are computationally intensive, they are also straightforward to parallelize and can be automated to a high degree, since image coregistration and most importantly phase unwrapping are not required.
Multi Aperture Interferometry [15] measures the azimuth displacement component between two SAR acquisitions, and can be proven to be an alternative formulation of the azimuth spectraldiversity technique originally proposed for image co-registration [16]. The same measurement principle, also referred to as split-bandwidth interferometry [14], can also be applied in the range dimension, where it was first proposed for absolute phase estimation [17]. To our best knowledge, however, the range variant of this method has never been applied to measure ice-motion. Splitbandwidth techniques exploit the fact that image mis-registrations in the azimuth (or LoS) dimension are necessarily coupled with phase variations, which are proportional to the registration offsets and to the central Doppler (or range) frequency of the processed data bandwidth. Such phase variations can thus be estimated by forming interferograms between upper and lower azimuth (or range) spectral sub-bands, from which displacement are then computed by differencing and scaling. From the performance point of view, these methods achieve a higher spatial resolution and accuracy compared to amplitude or intensity offset-tracking for a given level of coherence ( [14,18]), and are more robust and computationally efficient than complex cross-correlation or coherence maximization. Although split-bandwidth techniques share the same phase-coherence requirements as DInSAR, and also require image co-registration, it is only for large displacements (comparable to the size of the SAR resolution cell) that phase unwrapping is required, so that this processing step can often be neglected.
Finally, for all of the above-mentioned techniques, it can be necessary to estimate and remove slowly-varying image-wide error trends due to satellite orbit and radar system uncertainties, and to the relative nature of the measurements in the DInSAR case [19,20]. This is typically done by estimating the parameters of suitable error models based on points of known position and velocity, referred to in the following as Ground Control Points (GCPs).

Experiment Description
The ice velocity round-robin algorithm inter-comparison experiment was organized into four different tasks, namely DInSAR (Task 1), MAI (Task 2), incoherent offset-tracking (Task 3), and partially-coherent offset-tracking (Task 4). Invitations to participate in one or more tasks were issued by the Greenland Ice Sheet CCI consortium to academic, research and business organizations, through announcements on the project website and on the CRYOLIST e-mailing list for snow and ice related research, as well as through personal communications. For each task participants were provided with the same input SAR and auxiliary datasets (Tables 1 and 2) and were asked to apply their preferred implementation of one of the three processing algorithms described in Section 2.1. A set of deliverables was requested for each task (Table 2) and participants were requested to fill in a feedback form providing information on their processing algorithm and parameters.    [21] were also provided. Both SAR image pairs are highly coherent due to the small spatio-temporal baselines (Figure 2a,b), except for an unphysical band of lower coherence for the I2 pair, due to missing lines in the raw data. Participants were asked to generate a LoS velocity map using one of two applicable DInSAR techniques, namely DEME and DD, and to provide quality parameters if applicable as well as the position of any GCPs used in the processing.
The Multi Aperture Interferometry dataset (Task 2) consisted of the previously mentioned I2 pair in Table 1. Users were provided with the Level-0 (raw) data products and with the same ancillary data as for Task 1. Participants were asked to generate an azimuth velocity map and supply any applicable quality parameters and GCP positions.
For the offset-tracking experiments (Tasks 3 and 4) image pairs I3 and I4 in Table 1 were provided. The I3 dataset consisted of an Environmental Satellite Advanced Synthetic Aperture Radar (ENVISAT ASAR) C-band image pair spanning two standard frames (about 200 km × 100 km) around the Upernavik Isstrøm glaciers in northwestern Greenland (Figure 1c). The Synthetic Aperture Radar raw data were made available as well as Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) precise state vectors and the GIMP DEM. Due to the C-band radar frequency and the relatively long temporal baseline, the I3 pair has very low coherence values on the ice sheet ( Figure 2c). The I4 dataset consisted of an Advanced Land Observing Satellite 1 Phased Array L-band Synthetic Aperture Radar (ALOS PALSAR) Fine Beam Single image pair, spanning 5 standard frames and covering an area of about 300 km × 70 km, featuring several large tidal outlet glaciers, such as Rink Isbrae, Store Gletscher and Sermeq Avannarleq (Figure 1d). The Synthetic Aperture Radar raw data and the GIMP DEM were provided to the participants. Due to the longer radar wavelength, the I4 pair shows several patches of relatively high coherence values, both on the bedrock outcrops and on slowly-moving areas on the ice sheet ( Figure 2d). Coherence is, however, also lost for this pair closer to the ice margin, and in particular on the tongues of the above-mentioned glaciers. For both offset-tracking tasks, participants were asked to generate LoS and azimuth velocity maps using their offset-tracking method of choice, and to provide any applicable quality parameters and the position of any GCPs used during the processing.  Table 1.

Validation Strategy
The SAR-based measurements provided by the participants were validated against in situ GPS measurements available from different sources. GPS coordinates and Cartesian velocity components are provided in Table S1, whereas their locations are shown in Figures 1 and 2. Neither the GPS  measurements nor their locations were made available to the participants beforehand. For Tasks 1 and 2, velocities were derived from the positions of 10 stakes drilled into the ice surface, which were surveyed at least twice between 1992 and 1995 [22]. Since the last observations were in late June 1995, and thus between 6 and 8 months before the acquisition of the I1 and I2 SAR datasets (Table 1), the measurements were linearly extrapolated to a reference date (1 February 1996), to account for the post-surge, long-term adjustment of the ice dynamics, as discussed in [22].
To compare the GPS measurements, recorded in a local ENU Cartesian coordinate system (ve, vn, vu), with the LoS and azimuth velocity components provided by the SAR-based methods, the GPS velocities were projected onto the LoS and azimuth directions according to the following transformation: ( ) = ( cos cos sin cos sin sin − cos 0 ) ( ) (1) in which vr represents the SAR LoS velocity component (positive for motion towards the radar), va the SAR azimuth velocity component (positive for motion in the satellite flight path direction), φ the angle between the East direction and the ground projection of the SAR LoS vector, and θ the elevation angle from the horizontal plane to the SAR LoS vector. The values of these angles at scene center are listed in Table 1.

Results
In the following subsections, the results received for each task are described and intercompared. Contributions are anonymized and participants are referred to through group numbers, which are randomly assigned for each experiment.

Measurements
For this task, results were provided by 6 groups, affiliated (in alphabetical order) to the Geological Survey of Denmark and Greenland (GEUS), the Italian National Institute for Geophysics and Volcanology (INGV) and the Technical University of Denmark (DTU). Some institutions provided more than one result, obtained with different processing algorithms. Image focusing was carried out by all groups with the software of [25], whereas the interferometric processing chains, which were used were based on different combinations of the software modules developed by [25] and on DTU's IPP processor [26]. The LoS mean-velocity maps generated by the participants are shown in Figure 3.
Group 1 applied a DEME approach to the I1 dataset (Table 1), whereas groups 2, 3, and 5 applied the same approach to I2. Groups 4 and 6 processed both I1 and I2 using a DD algorithm. The main processing parameters are detailed in Table 3.
Co-registration was carried out by all groups according to a 2-step procedure, initiated by geometric calculations exploiting the available precise state-vectors and subsequently refined by fitting a mis-registration model to the results of an intensity image patch cross-correlation carried out on a coarse grid. Algorithms differed in the number of cross-correlation windows, their size, and in the order of the mis-registration polynomial model (Table 3).
Interferometric phase and coherence were estimated with similar methods, the main differences lying in the averaging factors (number of looks) and in the use of adaptive phase filtering [27], which was carried out by groups 1, 2, 3, and 6.
Concerning phase unwrapping group 1 used a branch-cut algorithm [28], whereas all other groups used different implementations of the Minimum Cost Flow (MCF) algorithm [29]. Weights for the latter were based on coherence for groups 2, 3 and 6 and on both coherence and intensity edgedetection for groups 4 and 5. All groups used a coherence threshold to generate a phase unwrapping mask, with thresholds selected either manually (group 2) or as a function of the number of interferogram looks (all other groups).
Absolute phase estimation and orbital refinement were carried out by all groups using GCPs manually selected on bedrock outcrops. The number and spatial distribution of the GCPs, however, varied significantly ( Figure 4). In particular, group 2 used a single GCP to estimate only a constant phase calibration parameter, whereas all other groups estimated four parameters describing a linear baseline error model [19].
Quality parameters were provided by all groups (Figure 4). For groups 1 and 3 to 6, they consisted of maps of predicted LoS error standard deviation, estimated according to the framework of [30], whereas group 2 provided a coherence map.  . Task 1 quality parameters provided by groups 1 to 6 (a-f). For group 2 (b) the quality parameter is interferometric coherence, whereas for all other groups it is the predicted LoS velocity error standard deviation in m/y. Triangles represent the GPS stations used for validation. Crosses represent GCPs used for orbital refinement and absolute phase estimation.

Intercomparison
The received results were resampled to a common grid with a posting of 3 arcsec in latitude and 12 arcsec in longitude, corresponding to about 90 m × 90 m on ground. For comparison purposes, the differences of each LoS velocity map with respect to that provided by group 6 were computed ( Figure 5). The measurements provided by group 1 (Figure 3a) show a more limited coverage compared to other groups, as well as significantly higher velocity magnitudes on ice and on bedrock outcrops, whereas the results provided by the remaining groups show local differences, in particular in the bottom-left and bottom-right corners of the image, as well as slightly varying image-wide trends (Figure 5b-e).
The differences observed for group 1 are most likely due to the more conservative branch-cut phase unwrapping algorithm, concerning coverage (i.e., the percentage of valid measurements), and to the selection of several GCPs in non-stationary areas ( Figure 4a) as well as in areas affected by phase unwrapping errors, which cause a spatially varying bias in the measurements.
Concerning the results obtained by groups 2 to 6, some discontinuities are seen in the difference maps, e.g., in the lower-right of Figure 5b,d,e, which correspond to regions almost entirely delimited by areas of low-coherence in Figure 2a,b. These regions are prone to phase unwrapping errors, particularly if network approaches such as MCF are used [31], which could explain the observed differences. Results generated with a DEME approach applied to data pair I2 (Table 1), namely by groups 2, 3 and 5 (Figure 5b-d), show some common difference patterns with respect to the DD results based on the I1 and I2 pairs, obtained by groups 4 and 6. This could be due to the height estimation procedures, based respectively on an external DEM in the DEME case and on the joint inversion of height and displacement from two SAR data pairs in the DD case, and/or to the different sensitivity to GCP height errors, which affect only the DEME results [32].
Concerning the quality measurements, those provided by groups 3 to 6 are similar, except for the fact that the error predictions associated with the DD methods (Figure 4d,f) are slightly lower than those associated with the DEME ones ( Figure 4c,e), due to the insensitivity to GCP height errors. Interestingly, the error predictions for group 1 are an order of magnitude higher than other groups. This is probably due to the fact that the phase unwrapping error modelling of Mohr et al. ( [30]) identifies several GCPs to be in areas which are likely to be affected by single cycle phase unwrapping errors. Figure 5. LoS velocity differences between groups 1 to 5 (a-e) with respect to group 6.

Validation
Measurement accuracy was assessed on exposed bedrock, identified using the ice-cover mask of Citterio et al. ( [33]), and on ice-covered surfaces by comparison with the GPS positions of 10 stakes drilled into the ice surface [22]. The LoS-projected GPS velocities were up to 48 m/y in magnitude. Summary statistics are reported in Table 4, whereas the differences for each stake are detailed in (Table S2). Table 4 confirms that the results of group 1 are affected by a significant image-wide error trend, whereas those of all other groups are similar and show a level of agreement with the GPS measurements, which is typically of 1 m/y or better in terms of median and Median Absolute Deviation (MAD). Interestingly, the group which achieves the lowest deviation with respect to the GPSsurveyed stakes, namely group 2, shows a slightly higher variability on bedrock areas compared to groups 3 to 6. This is probably due to the use of a single GCP for orbital refinement (Figure 4b), which leads to larger errors further away from this point, where some of the bedrock areas are located.
Finally, it is interesting to analyze to what extent the quality parameters provided by the participants predict the validation results. For groups 3 to 6, indeed, most observed errors fall within two predicted error standard deviations ( Figure 4 and Table S2), whereas for group 1, although the predicted error standard deviations are higher by an order of magnitude, they still consistently underestimate the observed errors. This is due to the fact that the error prediction algorithm used by group 1 can account for errors in GCP velocity only if an appropriate uncertainty value is specified by the analyst. Furthermore, only single-cycle phase unwrapping errors are assumed to occur [30]. Concerning group 2, which provided a coherence map as a quality measurement, the values of this parameter are uniformly high and thus insensitive to the observed errors (Table S2).

Measurements
For this task, the I2 dataset in Table 1 was processed by 4 groups, applying different processing chains listed in alphabetical order: DTU [34], Massachussets Institute of Technology [15], University of Edinburgh [18], and the University of Seoul [35]. Image focusing was carried out with the software of Werner et al. ( [25]) and that of Rosen et al. ( [36]). The resulting azimuth mean-velocity maps are shown in Figure 6.
Multi-aperture processing was carried out by all groups by creating two azimuth sub-apertures for each acquisition, with frequency separations between 0.4 and 0.5 of the SAR Pulse Repetition Frequency (Table 5). Different approaches were followed to generate the sub-apertures. Group 1 focused each image by processing the full azimuth bandwidth, and subsequently carried out spectral whitening combined with azimuth common-band filtering and range spectral-shift filtering. Coregistration was performed according to the two-step procedure described in Section 3.1.1. Finally, the two azimuth sub-apertures were formed by azimuth band-pass filtering. In contrast, groups 2 to 4 carried out the focusing step twice for each raw dataset, selecting different azimuth center frequencies and processing bandwidths to jointly carry out azimuth common-band filtering and sub-aperture formation. The corresponding sub-aperture images for each acquisition were then co-registered following a similar two-step procedure as for group 1.
All groups performed low-pass filtering and decimation (multi-looking), followed by a local mean filtering in the case of groups 2 and 3 and by a Goldstein adaptive filter [27] in the case of groups 1 and 4. Group 4 carried out multi-looking and adaptive filtering on each sub-band interferogram, prior to MAI (difference) interferogram formation, whereas group 1 applied both filtering steps to the final MAI interferogram alone. An additional processing step was carried out by group 4, to reduce the coregistration errors induced by ionospheric propagation also known as "azimuth streaks" [37]. The algorithm which was applied is similar to the one proposed in Raucoules et al. ( [38]) for offset-tracking and consists in masking out areas with significant motion based on a user-defined threshold, rotating the image by a user-defined angle so that most ionospheric streaks are parallel to the slant-range direction (under the assumption that most streaks are almost straight lines), and finally isolating the streaks from residual motion and other error sources with a suitably designed spatial filter.
Phase unwrapping was carried out only by group 4, using an MCF approach [29]. MAI phase was then converted to azimuth displacement after performing an orbital refinement in the case of group 4, and applying different scaling factors, based either on the nominal antenna length (group 1) or on the sub-band frequency separation and the beam footprint (ground) velocity (groups 1,2,3).
Quality parameters were provided by all groups except group 1 and consisted of predicted error standard deviations (Figure 7), based on local moving-window estimates.

Intercomparison
The received results were resampled to a common grid with a posting of 3 arcsec in latitude and 12 arcsec in longitude, corresponding to about 90 m × 90 m on ground, and the difference of each azimuth velocity map with respect to that provided by group 4 was computed ( Figure 8). The difference maps are dominated by ionospheric streaks, since group 4 is the only group, which implemented a dedicated filtering step for their estimation and removal (Table 5). This is apparent also from Figure 6d, in which the results of group 4 show significantly lesser streak artefacts. Other differences are related to slowly-varying trends in both the range and azimuth dimensions and to the levels of spatially uncorrelated noise. Relative biases among the measurements can be due to different co-registration error models and parameter estimations, and to the orbital refinement phase calibration step in the case of group 4. Uncorrelated noise levels are instead influenced by the respective filtering approaches and parameters. In particular, groups 1 and 4 applied a higher degree of spatial filtering compared to groups 2 and 3. Concerning groups 1 and 4, although the same filtering algorithms were applied with similar parameters (Table 5), the results of group 4 contains a higher level of short-wavelength noise ( Figure S1d). This could be due to the application of the phase filtering steps on the sub-band interferograms in the case of group 4, and on the differential MAI interferogram in the case of group 1, as discussed further in Section 4.2.
Finally, the results of group 1 contain an error due to the handling of missing-lines by the focusing software (red line crossing GPS stake 6 in Figure 8a), which is also visible in the coherence map generated from the same focused image pair (Figure 2b).
Concerning the error predictions in Figure 7, these are significantly lower for groups 2 and 3, compared to group 4, even though the local variations would actually seem lower for the latter ( Figure S1). We speculate that this could be due to the use of smaller estimation windows for the computation of the local standard deviations on which the error estimates are based, although information on the window sizes was not provided by all participants and thus this hypothesis cannot be confirmed.

Validation
Measurement accuracy was assessed on exposed bedrock, identified using the ice-cover mask of Citterio et al. ( [33]), and on ice-covered surfaces, by comparison with the same GPS data used for the validation of Task 1 [22]. The azimuth-projected GPS velocities were within 112 m/y in magnitude (Table S3). Summary statistics are reported in Table 6, whereas the differences for each validation location are detailed in Table S3. Table 6 shows that group 4 achieved the highest level of agreement with the GPS measurements, as well as the lowest velocity variations on bedrock. Based on Figure 8, the main reason for this is the reduction of ionospheric streaks errors, which is also achieved to a lesser extent by the spatial filtering approaches applied by the remaining groups.
Concerning the quality parameters, the validation results of group 4 were essentially predicted by the error estimates ( Figure 7), with all but one observed error lying between two predicted error standard deviations (Table S3). Groups 2 and 3 instead systematically underestimated the observed errors, since these are mainly due to ionospheric streaks, which are not accounted for by locallyestimated standard deviations. No quality parameter was provided by group 1.

Measurements
For this task the ENVISAT ASAR I3 dataset in Table 1 was processed by 6 groups, applying different processing chains, listed in alphabetical order: DTU [39], Environmental Earth Observation Information Technology (ENVEO, [3]), GAMMA-RS [25], German Research Centre for Geosciences (GFZ-Potsdam, [40]), GEUS [25], Swansea University [41]. The raw SAR data were focused with the software of [25,36,41]. LoS and azimuth mean-velocity maps were provided by all participants and are shown in Figures 9 and 10 respectively.
The core processing algorithm was ICC for all groups, although several differences can be pointed out, as summarized in Table 7. Prior to offset-tracking, three groups carried out image coregistration, following the two-step approach described in Section 3.1.1, whereas the remaining groups skipped this step. Group 2 also carried out multi-looking and high-pass spatial filtering of the intensity images to enhance features. Concerning offset computation, groups 1, 4, and 6 chose the same oversampling factor, cross-correlation window size and a similar spacing between adjacent windows, leading to comparable spatial resolutions of the measurements. Group 2 chose window sizes corresponding to a larger ground coverage, whereas groups 3 and 5 applied larger oversampling factors prior to the cross-correlations, and chose smaller window sizes respectively in both image dimensions (group 3) and in the range dimension only (group 5). Regarding outlier culling, all groups applied thresholds to one or more quality parameters, namely cross-correlation Signal-to-Noise-Ratio (SNR), as defined in [41], Normalized Cross-Correlation (NCC) coefficient, and maximum expected velocity magnitude. Group 2 additionally carried out a culling based on the local variability of flow magnitude and direction. The culled offsets were smoothed (filtered) by groups 1 and 6, which applied a moving average window, and by group 3, which used anisotropic diffusion filtering [42].
Estimation of image-wide error trends (calibration) was carried out by groups 1 and 6 using GCPs selected on bedrock, respectively by visual inspection (group 1) and based on the ice mask provided in Citterio et al. ( [33]) (group 6). Group 4 implicitly performed a calibration step by fitting the coefficients of a 2D polynomial to high-SNR range and azimuth offsets computed on a coarse grid, and then co-registering the SAR images according to this model prior to offset-tracking. The positions of the GCPs used by each group are shown in Figure 11.
Quality parameters were not provided by all groups and were quite heterogeneous. They consisted of per-pixel LoS and azimuth error standard deviation estimates for groups 1 and 6, NCC coefficient amplitude for groups 3 and 5, and cross-correlation SNR [41] for group 4 ( Figure 11). The predicted error standard deviations provided by groups 1 and 6 were generated by applying the framework of [30] and modelling measurement errors through the local variance, scaled to account for the partial overlap of adjacent cross-correlation windows [13].

Intercomparison
The coverage of the measurements is the first apparent difference between the results. Due to the low interferometric coherence (Figure 2c), measurements were expected to be feasible only on feature-rich areas, namely on bedrock outcrops and close to the ice-sheet margin. Indeed group 5 provided results also towards the interior of the ice-sheet, although these are clearly noisy (Figures 9e  and 10e) and correspond to very low NCC values (Figure 11e). Groups 2 and 3 (Figures 9b,c and 10b,c) obtained almost complementary coverages, confined to glacier tongues in the case of the former and to rock and slow-moving ice-sheet areas for the latter. For group 2, this could be due to tight settings for the additional culling step based on local flow direction variability, which can be expected to be more sensitive to noise in areas where the velocity magnitudes are low. In contrast, the limited amount of good matches found by group 3 on glacier tongues seems due to the higher NCC threshold chosen for outlier culling, combined with the fact that the NCC values are reduced by the small crosscorrelation window sizes (Figure 11c). Finally, groups 1, 6 and 4 all provided results on rocky surfaces and close to the ice-sheet margin as expected, although with a decreasing measurement coverage. Based on the insets in Figures 9 and 10, group 1 appears to have used a less conservative hole-filling approach compared to both groups 4 and 6, although details of these procedures are not available to confirm this. Another lacking information is the culling procedure implemented by group 4, which is based on SNR thresholding and on a second unclearly described criterion, which is likely however to be the cause of the reduced measurement coverage compared to groups 1 and 6, since the latter used more conservative SNR thresholds.
To compare the location and spatial-correlation of the differences between groups, the received results were resampled to a common grid with a posting of 9 arcsec in latitude and 27 arcsec in longitude, corresponding to about 250 m × 250 m on ground. In Figures 12 and 13 the differences between the LoS and azimuth velocities compared to group 6 are shown for the area in the insets in Figures 9 and 10.  Differences are consistently higher for the azimuth measurements, probably due to the fact that for this dataset the highest velocity gradients, and hence challenges for the measurement algorithms, occur in this direction. Several spatially correlated differences are found in shear zones and more in general in areas adjacent to data voids. The fact that this holds also for participants, such as groups 1 and 6, whose measurements have the same spatial resolution and level of filtering, and also similar outlier culling thresholds ( Table 7), suggests that this is mostly due to different hole-filling procedures. Spatially uncorrelated differences, such as those of group 5 compared to group 6, are instead more likely to be due to lower outlier culling thresholds and a lesser degree of spatial filtering.
Concerning the quality parameters (Figure 11), the predicted error standard deviations provided by groups 1 and 6 are similar, and mostly reflect the local spatial variability of the velocity measurements, which are higher for the azimuth component compared to the LoS one. The parameters provided by groups 3, 4, and 5 characterize the quality of the cross-correlation peaks, which are influenced by the strength (signal-to-clutter-ratio) of the intensity features and/or by the interferometric coherence.

Validation
Measurement accuracy was assessed on exposed bedrock, identified using the ice-mask of [33], and on ice-covered surfaces, by comparison with the GPS measurements. Summary statistics are reported in Table 8, whereas the differences for each GPS station are detailed in Table S4.
For all groups, the velocity errors with respect to the GPS measurements are larger in the azimuth component compared to the LoS one, except for group 5, whose statistics are significantly affected by outliers. In contrast, the LoS and azimuth variations on bedrock are comparable and, except for group 2, they are an order of magnitude less than the errors observed at the GPS locations. Thus, the worse agreement of the SAR azimuth measurements cannot be due to ionospheric propagation errors, i.e., the azimuth streaks discussed in Section 3.2, which would equally affect measurements on bedrock and at the GPS locations. A possible explanation could be that the limited (kilometric) spatial resolution of the SAR measurements (Table 7) has a greater impact on the accuracy of the azimuth component, since the velocity gradients are higher in this direction compared to the LoS.
Concerning the differences discussed in Section 3.2.1, a velocity profile is shown in Figure 14, which crosses the location of 4 GPS stations along the dashed line in the insets of Figures 9 and 10. The negative differences between group 1 and most other groups in this area appear due to errors in the measurements of group 1, caused by excessive hole-filling as previously discussed. The higher levels of spatially uncorrelated noise of groups 4 and 5 compared to other groups are also apparent in Figure 14, as well as the unexpectedly high azimuth variability of group 2, which was also observed on bedrock as mentioned above, and the reasons of which are unclear.
Concerning the agreement between validation results and quality parameters, for groups 1 and 6, almost all observed errors fall within two predicted error standard deviations (Table S4). This is consistent with the previously mentioned hypothesis that the main error source for this dataset is the spatial variability of the underlying velocity field, particularly in the azimuth direction. The error predictions do in fact account for spatial variability, although just for the one between neighboring measurements, rather than for that within the coverage of the cross-correlation windows. For the groups which provided maps of NCC coefficient amplitude or cross-correlation SNR as quality parameters, there is not a clear relation between these and the observed errors (Table S4).
The LoS and azimuth mean-velocity maps provided by the participants are shown in Figures 15  and 16 respectively. Since group 7 provided the horizontal velocity components along the Cartesian axes of a National Snow and Ice Data Center (NSIDC) Polar Stereographic North projection (EPSG 3413), the LoS and azimuth components were derived for this group by using the surface parallel flow assumption and the GIMP DEM [21]. Based on a comparison with optical imagery, significant relative geocoding discrepancies were observed among groups, ranging between +/− 300 m on ground, so that constant 2D shifts were estimated and applied to inter-compare and validate the results.  8 (a-h). Positive values correspond to motion towards the radar along the LoS vector (from ground to satellite), the ground projection of which at scene center is shown by the arrow in (a). The insets show the GPS stations used for validation.
The participants applied several different offset-tracking techniques, of which the main processing steps and parameters are summarized in Table 9. Groups 1, 2, 4, 5, 6, and 8 all used ICC as a core calculation method, whereas Group 3 relied on ACC and Group 7 applied a combination of CCC and ACC.
Quality parameters, shown in Figure 17, consisted of per-pixel LoS and azimuth error standard deviation estimates for groups 1, 6 and 8, the amplitude of the NCC coefficient for group 5, and the cross-correlation SNR [41] for group 4. Groups 2 and 3 did not provide quality estimates, whereas LoS and azimuth quality parameters could not be derived for group 7, based on the information which was made available. The predicted error standard deviations provided by groups 1, 6 and 8 were generated by applying the framework of [30], adapted for offset-tracking as described in Section 3.2.1.

Intercomparison
The coverage of the measurements differs in slow-moving areas towards the interior of the icesheet, in fast-moving areas close to the glacier termini and in terms of the overall azimuth extent.
Towards the interior of the ice-sheet some results show significantly more data voids than others (e.g., Figure 15b,d,f). Several processing steps and parameter trade-offs can contribute to this. For instance, the conservative coverage of group 2 (Figures 15b and 16b) could be influenced by the highpass pre-filtering of the image intensities, since the ice-sheet interior is poor of intensity features, and/or by outlier culling based on the local variability of the flow direction, which was not applied by any other group. The limited coverage of group 4 compared to group 8, despite the choice of a lower cross-correlation SNR threshold (Table 9), is most likely due to the application of an additional culling procedure, which was however unclearly described. Interestingly, group 6, which used the same SNR threshold as group 8, also achieved a poorer coverage compared to the latter. This could be due to the coarser azimuth posting of the cross-correlations (150 m vs. 60 m), which reduces the chances of finding good cross-correlation peaks. The most comprehensive coverage was achieved by group 7 (Figures 15g and 16g) and is most likely due to the combination of a relatively dense posting for the cross-correlation computations and to the use of variable cross-correlation window sizes (Table 9).
In fast moving areas the coverage differences are apparent from the insets in Figures 15 and 16. Different strategies were applied by the participants to accommodate the full motion range. Group 2 used a large (2.5-3 km) search window size, which would be sufficient to track a 10 km/y motion. Groups 3, 7 and 8 used a priori information on the velocity magnitude to either increase the search window size (group 3), to selectively increase the cross-correlation window size and apply multilooking (group 7) or to select the position of the image patches to cross-correlate (group 8). Figure 17. Task 4 quality parameters. Groups 1 slant-range (a) and azimuth error standard deviation (b); group 4 cross-correlation signal to noise ratio (c); group 5 normalized cross-correlation amplitude (d); group 6 slant-range (e) and azimuth error standard deviation (f); group 8 slant-range (g) and azimuth error standard deviation (h). Triangles represent the GPS stations used for validation. Crosses represent GCPs used for velocity calibration. Groups 1, 4 and 6, which did not implement any of such strategies, achieved a limited coverage close to the termini of the outlet glaciers (panels a, d and f in Figures 15 and 16).
Finally, coverage differences can be noticed in the overall azimuth length of the ice velocity maps in Figures 15 and 16. Besides group 4 (panel d), which unintentionally processed only part of the provided dataset, these differences are due to the processed azimuth bandwidth at the image borders, which is typically a user-configurable parameter during the image focusing step.
Concerning the properties of the velocity measurements, the LoS maps differ primarily due to image-wide trends in the east-west direction. To highlight such differences all results were resampled to a common grid with a posting of 4.5 arcsec in latitude and 13.5 arcsec in longitude, corresponding to about 125 m × 125 m on ground, and differenced with respect to the results of group 8. Figure 18 shows the results for an area of interest around Store Gletscher and Sermeq Avannarleq (STO4 and AVA1 GPS stations respectively). Smaller scale differences are also apparent, particularly in the Rink Isbrae area (Figure 19), and could be caused by geocoding differences, particularly in the north-south direction, which were not corrected by the rigid 2D shift applied in the pre-processing procedure.
More subtle differences in the LoS measurements can be inferred from Figure 20, in which the velocities derived by all groups are shown relative to a common reference point and are wrapped so that each color cycle represents a 200 m/y variation. For instance, the results from group 2 (Figure 20b) show some anomalous bands parallel to the azimuth direction, whereas the measurements of group 7 show a greater fringe detail (spatial resolution) particularly in the glacier shear zones.   The azimuth velocities also differ primarily due to East-West trends ( Figure 21). Furthermore, the results of group 2 show a difference pattern with respect to other groups, which is correlated with the underlying azimuth velocity values (panel b in Figures 21 and 22). The measurements of all groups are affected by ionospheric streaks [37], as seen in Figure 23, although for group 4 these artifacts are reduced due to a dedicated filtering step (Table 9). Finally, differences are observed concerning the spatial resolution of the measurements, which can be inferred from the clarity of the fringes on Store Gletscher and on Sermeq Avannarleq in Figure 23. Concerning the quality parameters (Figure 17), the predicted error standard deviations provided by groups 1 and 6 are very similar, whereas those provided by group 8 are slightly lower, although the reason for this is not apparent based on the processing parameters. Since for this track the GCPs are located within 50 km of any measurement, the error predictions mainly reflect the local variability of the measurements. In contrast, the indices provided by groups 4, and 5 (Figure 17c,d), characterize the quality of the cross-correlation peaks, and are thus influenced by the signal-to-clutter ratio of amplitude/intensity features and by the level of phase coherence.

Validation
Measurement accuracy was assessed on exposed bedrock, identified using the ice-mask of [33], and on ice-covered surfaces, by comparison with GPS measurements. Summary statistics are reported in Table 10, whereas the differences for each GPS station are detailed in Table S5.
On bedrock outcrops, assumed to be stationary, median LoS velocities smaller than 1 m/y are measured by virtually all groups, whereas most groups show a negative bias in the azimuth component. Statistics on bedrock could not be computed reliably for group 5, since most of the variations found on bedrock fall within the quantization step of the measurements determined by the oversampling factor in Table 9, which corresponds to 3 m/y and 4.6 m/y respectively in the azimuth and LoS velocity components. For the remaining groups median absolute deviations lie between 1 and 4 m/y in LoS and 3 to 8 m/y in azimuth. A greater azimuth variability is expected due to the presence of azimuth streaks caused by ionospheric propagation [37]. Using the simplified approach of Raucoules et al. ( [38]), the contribution of this error source is found to be approximately within +/− 50 m/y ( Figure S2). For group 4 this error contribution is lower (Figure 23d), since an ionospheric filtering procedure was applied (Table 9).
On ice surfaces, in the area of the STO4 station, all groups underestimated the LoS-projected GPS measurement by an amount varying between 7.7 m/y (group 7) to 40 m/y (group 1) (Table S5). The profile in Figure 24c, taken from points β to β' in Figure 21, shows that the differences between groups are spatially correlated and not confined only to the vicinity of the STO4 station. In particular, around the peak of the profile, at about 20 km from point β, groups 5, 7, and 8 are in good agreement with each other, whereas the results from other groups differ significantly. This is very likely due to a higher azimuth resolution of the measurements in the shear zone of the glacier, which can be inferred from Figure 23e,g,h. In turn this could be caused by the use of smaller azimuth crosscorrelation window sizes (Table 9). We speculate that azimuth spatial resolution is also the main cause of the LoS differences seen at STO4. Concerning the azimuth measurements, ionospheric streaks explain most of the errors for groups 3, 5, 7, and 8 (Tables S5 and S6). For groups 1 and 6 an additional error contribution is provided by long-wavelength errors (Figure 21a,f), whereas group 2 ( Figure 21b) shows spatially correlated differences with respect to all other groups, which seem correlated with the underlying azimuth velocities, and could thus depend on the use of a different scale-factor in the conversion from pixel offsets to velocity values.

Station
Component Group  G1  G2  G3  G4  G5  G6  G7  G8  SAR-GPS  N a  LoS  2  3  3  1  3  1  3  2  Azimuth  2  3  3  1  3  1  3  In the area of the AVA1 station, concerning the LoS velocity component, group 7 shows quite a different profile compared to most other groups ( Figure 24a) and also the highest disagreement with the GPS measurement (Table S5). Based on a visual inspection this is due to a correlation of group 7′s results with topography, most likely due to the pre-processing step used to convert Cartesian velocities to LoS and azimuth components. Groups 1, 4, and 6, show LoS gradients which differ from the other groups (Figure 18a,d,f) and seem incorrect, based on the profile through GPS station AVA1 in Figure 24a, in which point α', located on a bedrock outcrop, has negative velocity values only for these groups. Such image-wide error trends are due to different strategies in the selection of the GCPs used for orbital refinement. Groups 1 and 6 manually selected points on bedrock outcrops, assuming these to have a zero motion, whereas group 4 selected points with a high cross-correlation SNR, although some of these do not fall in stationary areas (Figure 17c). Groups 3, 7 and 8 selected zeromotion GCPs automatically based on an ice-cover mask, whereas group 2 did not carry out orbital refinement. Concerning the azimuth measurements, ionospheric propagation errors (Table S6) account for the differences observed for group 5, and only partially for group 2 (Table S5), since the latter also shows an error pattern which is correlated with the velocity field ( Figure 21b). Groups 1, 3, 4, 6, and to a lesser degree groups 7 and 8, underestimate the GPS azimuth velocity at AVA1, even more if the ionospheric errors are accounted for, as seen in Figure 24b and Table S5. As for the STO4 case, the smaller errors of groups 5, 7, and 8 could be due to the higher azimuth resolution of the measurements.
In the area of the RNK2 station, concerning the LoS velocity component, all groups except group 4 overestimate the GPS measurement (Table S5 and Figure 24c). The reason for this is not clear. A speculation is that this could be due to orbital refinement inaccuracies due to a less favourable distribution of rock surfaces and thus of GCPs, compared to the STO4 and AVA1 validation sites. From Figure 19, it is seen that the results of groups 1, 2, 7, and 8 differ only locally, whereas groups 3, 4 and to a lesser extent also 5 show differences which are spatially correlated with the underlying LoS velocity pattern ( Figure 15). This could be due to relative geocoding shifts, which were not corrected by the constant shift model used in the pre-processing step. Groups 4, 5 and 6 also show differences with respect to image-wide East-West trends, due to the coregistration step and/or the GCP-based calibration of the offsets. Concerning the azimuth differences, group 2 ( Figure 22b) shows a bias in the measurements towards lower values, which is also visible in the profile in Figure 24f, whereas group 4 shows spatially correlated differences with respect to other groups, which could be due to relative geocoding shifts and/or to the coregistration step. Concerning the agreement between validation results and quality parameters, the observations are too few to draw reliable conclusions (Table S5). Two shortcomings of the approaches which were used can, however, be pointed out. Firstly, none of the quality parameters account for ionospheric errors, particularly those affecting the azimuth velocities, although this is the dominant error source on bedrock outcrops, where the GCPs used for calibration are located. Secondly, quality parameters, which reflect the strength of the cross-correlation peaks, e.g., cross-correlation SNR provided by group 4 (Figure 17), fail to characterize spatially-correlated errors, due for instance to the GCP-based calibration step (Figure 18d).

Discussion
The experiments discussed in Section 3 provide a comprehensive picture of the properties of the different SAR-based ice velocity measurements techniques. At the same time, it was unexpected for the results provided by the participants to differ often significantly in terms of accuracy and/or coverage. This is due to the variety of processing algorithms applied and to the influence of analyst choices on some critical processing steps.

Differential Synthetic Aperture Radar Interferometry
Concerning the DInSAR task (Section 3.1), most groups achieved LoS MADs better than 1 m/y (Table 4) and spatial resolutions of 45 m × 45 m on ground, using two-or four-pass approaches, applied to highly coherent ERS Tandem data. The most significant differences among groups consisted of image/wide trends due to the manual selection of GCPs for orbital refinement and absolute phase estimation, as well as to large-scale phase unwrapping errors ( Figure 5).
Errors in the GCP selection process could be reduced by automated procedures, such as those implemented by some of the tasks 3 and 4 participants, in which an ice-cover mask was used to avoid points in non-stationary areas. Methods to ensure a homogeneous spatial distribution of the GCPs throughout the scene have also been proposed in literature [44]. Moreover, for several missions which became operational in the last decade, the orbit control and overall geolocation accuracy are significantly improved compared to the sensors considered in this study (e.g., [45][46][47]), which allows to simplify GCP-based calibration approaches by reducing the number of error-model parameters, and thus the estimation errors.
In contrast, no general automated method is available for detecting phase unwrapping errors in single interferograms, although synergies with range split-bandwidth or offset tracking techniques could be exploited for sensors with sufficiently high range resolutions [14]. Since phase unwrapping errors are due to local violations of the assumption that interferometric phase gradients are smaller than +/−π radians, their occurrence can be greatly reduced by considering only acquisitions with small temporal and spatial baselines, thereby limiting phase noise due to temporal decorrelation and sensitivity to topography, and by applying these methods towards the interior of the ice-sheet, to limit the phase contribution of velocity gradients.

Multi Aperture Interferometry
Concerning the MAI task (Section 3.2), the best MAD accuracy achieved with respect to the GPS surveys was of 7 m/y (Table 6), and azimuth velocities maps with a spatial resolution of about 90 m × 90 m were generated. The algorithmic differences among groups, which most significantly affected performance, were the handling of ionospheric propagation errors [10] and the reduction of spatially uncorrelated noise due to loss of interferometric coherence.
The mitigation of ionospheric streak errors was successfully carried out by one of the participants, although the procedure which was used cannot be considered of general applicability, since moving areas were excluded from the estimation, an empirically assigned threshold on the MAI phase. For systems with a sufficiently high range bandwidth, and datasets/areas meeting the abovementioned conditions to reduce the occurrence of phase unwrapping errors, general ionosphericstreak correction approaches could be investigated based on the split-spectrum method [48][49][50], in which motion and ionospheric phase contributions are separated exploiting the frequency dependence of the latter.
Regarding the filtering of spatially uncorrelated noise, a difference among participants concerned the application of spatial filtering methods to each sub-band interferogram as opposed to their application to the final differential MAI interferogram. The superiority of the former approach was established for the case in which multi-looking is the only filtering method applied, based on analyses carried out on real datasets [35], as well as on theoretical derivations supported by numerical simulations [51]. However, the intercomparison between groups 1 and 4 in Section 3.2.2 suggests that these conclusions should be re-examined for the common case in which adaptive filtering methods, such as [27], are used. Unlike multi-looking, these techniques can in fact be expected to perform better on the differential MAI interferogram, in which fringes common to both sub-band interferograms (e.g., due to LoS motion) have been cancelled out.
The Multi Aperture Interferometry quality parameters provided by the participants were based on the local variability of the measurements. While this can provide a reasonable characterization for spatially uncorrelated error sources, such as decorrelation, other spatially-correlated error sources such as the above mentioned ionospheric streak errors and image-wide co-registration errors cannot be accounted for based on a single acquisition pair.

Offset Tracking
The results of both offset-tracking experiments (Sections 3.3 and 3.4) were the most heterogeneous. The highest levels of agreement with the GPS measurements yielded LoS and azimuth MADs of 20 m/y and 40 m/y respectively for the ENVISAT ASAR dataset (task 3), and of 10 m/y in LoS and 20 m/y in azimuth for the ALOS PALSAR dataset (task 4). The spatial resolutions can be considered at best to approach the cross-correlation window spacing and at worst twice the size of the correlation window size [52], which on average amounted to values between 200 m × 170 m and 2.7 km × 1.7 km (ground range × azimuth) for task 3, and between 160 m × 120 m and 640 m × 1100 m for task 4. For both processing tasks the limited spatial-resolution of the measurements was one of the major error sources for areas with high velocity gradients, in agreement with the SAR vs. GPS comparisons reported in [23]. Furthermore, sub-optimal outlier culling and hole-filling procedures were also relevant in particular for the void-rich task 3 dataset, whereas GCP calibration errors as well as ionospheric streaks were noticeable for task 4, due to the larger spatial extent of the measurements and the lower (L-band) radar frequency. Finally, both offset-tracking tasks intentionally provided a challenge regarding the measurement of very high velocities (i.e., >2 km/y) close to the glacier termini.
Concerning the choice of the cross-correlation window sizes, which typically determine also the spacing at which the cross-correlations are computed, several groups chose large values, assigned based on analyst experience, and constant for the processing of the whole dataset. Although this approach typically succeeds in providing good-quality cross-correlations in ice-sheet regions characterized by different levels of coherence and intensity features, the optimal window sizes depend on surface and radar properties, which can vary significantly. Within the validation area of the task 3 dataset (Section 3.2.3), the smallest window sizes which provide good quality measurements lie between 1.0 km and 1.5 km in both ground-range and azimuth, whereas for the task 4 dataset (Section 3.4.2), good cross-correlation peaks at the GPS locations are found already for window sizes between 300 m and 500 m in each dimension. Since these values cannot be known beforehand, a good compromise to preserve spatial resolution would seem to be the hierarchal approach of Joughin ( [13]), in which cross-correlations are attempted for each measurement with increasing window sizes until quality criteria are satisfied. Concerning the latter, the use of an NCC threshold, as proposed by Joughin ([13]), could be complemented by a second threshold on crosscorrelation SNR [41], to account also for off-peak cross-correlation values in the quality characterization.
Measurements in fast-moving areas were provided by groups which applied processing approaches capable of accommodating the full range of motion. This was done by adjusting either the search window size or its position, or both, possibly guided by a priori velocity information. Increasing the search window size has typically a significant computational cost and can lead to a larger percentage of false-detections due to the presence of competing cross-correlation peaks. In contrast the positioning of search windows based on a priori information, such as a seasonal velocity map, depends on the velocity differences compared to the latter, as well as on the pixel spacing and temporal baseline of the SAR acquisitions. This study does not provide clear evidence in favor of one approach or the other; however, the results of task 4 (insets in Figures 15 and 16) in particular show that in general one of these approaches should be implemented to derive measurements close to the glacier termini. Furthermore, since interferometric coherence is lost in fast moving areas, processing steps which enhance features with respect to radar speckle, such as the low-or high-pass filtering implemented by some of the participants, can be expected to improve the quality of the crosscorrelation peaks in these areas. However, clearer indications on this can be provided by dedicated studies, such as [41].
Regarding outlier culling, several groups obtained good validation results on static bedrock with approaches based on the thresholding of NCC magnitude and cross-correlation SNR, using the same thresholds for the datasets in tasks 3 and 4. Further culling steps based on the analysis of local velocity magnitude and/or flow-direction variations proved to be effective, but also very sensitive to analyst choices. Concerning the hole-filling procedures, the results of task 3 suggest that conservative approaches are preferable, to avoid unphysical interpolations in areas with strong velocity gradients such as shear margins.
Spatially-correlated measurement biases were observed, related to the selection of GCPs and to ionospheric propagation errors in the case of the azimuth measurements. Concerning the latter, the same considerations discussed above for MAI hold also for offset-tracking. However, a general ionospheric-streak correction procedure based on the split-spectrum principle, as discussed in Section 4.2 for MAI, would only be applicable to azimuth offset-tracking measurements in coherent areas, where interferometry and phase unwrapping can also be carried out. Concerning GCPs, the considerations discussed in Section 4.1 for DInSAR hold also for offset-tracking.
Finally, concerning the quality parameters associated with the offset-tracking measurements, error estimates were either provided in the form of intermediate quality parameters, such as phase coherence, NCC magnitude or cross-correlation SNR, or in the form of 1σ uncertainties, accounting for the local variability of the measurements and in some cases for the position of the GCPs used for image-wide velocity calibration. The latter provided in several cases the correct order of magnitude for the observed errors, although it can be considered only a partial error characterization, since uncertainties due to the limited spatial resolution in both LoS and azimuth, to image-wide coregistration errors and to ionospheric propagation remain unaccounted for. In contrast intermediate parameters, although useful for outlier culling, provided a very poor indication of end-to-end measurement errors, besides being of difficult interpretation for an end-user unfamiliar with the processing techniques.

Conclusions
In the current context in which SAR-based ice-velocity products are being made available by an increasing number of independent providers, it would be desirable from an end-user point of view, for these products to share similar state-of-the-art properties in terms of measurement spatial resolution, accuracy, coverage, and error characterization. However, the experiments described in this paper show that single-track LoS and azimuth ice-velocity measurements, which represent the building block of higher-level products distributed to the end-users, can have quite heterogeneous properties when generated by different processing chains and analysts.
The contribution of this study is two-fold. Firstly, a total of 24 SAR-based ice-velocity processing results are described, generated by 13 international research groups, and a total of 24 GPS measurements used for validation are provided. The results cover all SAR processing techniques currently used to carry out ice-motion measurements, namely DInSAR, MAI (azimuth split-bandwidth) and offset-tracking. Both the SAR-based results and the GPS data provide a unique benchmark to assess the performance of processing chains implementing these measurement techniques.
Secondly, the algorithmic steps and/or parameters responsible for performance variability are analyzed in detail for each technique. Interestingly, the causes for sub-optimal performance are most often related to the application of procedures, which are strongly influenced by analyst choices, e.g., concerning processing parameters. Useful indications to automate some key processing steps, while at the same time yielding improved performances, are obtained by analyzing the available processing results, as discussed in Section 4 and detailed in the inter-comparison and validation sections of each experiment (Section 3).
Finally, although this study was carried out on archive data from ESA and ESA Third Party Missions, the considerations presented in Section 4 concerning the processing algorithms are not sensor-specific and are thus equally applicable to current (e.g., Sentinel-1a,b) and forthcoming (e.g., NASA-ISRO) SAR missions. In particular, since the latter by design support the systematic application of interferometric techniques, this will allow the synergies between different measurement algorithms to be exploited to a much higher degree, thus contributing to define more robust and analyst-independent procedures for several critical processing steps, besides of course noticeably improving the spatial resolution and the accuracy of the final ice-velocity products.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Figure S1: Task 2 wrapped azimuth velocity maps, Figure S2: Task 4 estimated azimuth velocity ionospheric propagation errors, Table S1: GPS data overview,  Funding: This research as well as the APC were funded by European Space Agency contract 4000104815/11/I-NB. The SAR data was obtained through CAT-1P proposal 11339.