Greenland ice velocity maps from the PROMICE project

,


Introduction
The Greenland Ice Sheet (GrIS) is a major contributor to sea-level rise, and approximately half of this contribution is due to ice dynamics (Shepherd et al., 2019). Thus, in order to constrain the on-going mass loss of the GrIS it is important to obtain ice-sheet wide observations of ice-flow velocities. High temporal and spatial resolution will further allow us to distinguish 20 between annual or sub-annual variations and long-term trends, aiding in improving our understanding of the processes behind the observed changes. This is especially important because the flow of glaciers and ice caps vary on a range of timescales in response to the seasonal cycles, climate change, or internal variability (e.g. Moon et al., 2020;Joughin et al., 2018;. In-situ measurements of ice-flow velocities (e.g. from GPS) are relatively sparse on the GrIS and are often of short acquired 6 or 12 days apart (see Sect. 3.1). To achieve a consistent coverage (see Sect. 6), each mosaic is based on velocity measurements from all possible 6-and 12-day pairs using data from Sentinel-1A and 1B within two consecutive Sentinel-1A orbit cycles (i.e. 24 days). A new map is produced for every Sentinel-1A cycle i.e. a new mosaic every 12 days. A given mosaic thus overlaps by 12 days with the previous and subsequent maps. The dataset is expanded continuously, and we aim to provide a new mosaic within 10 days of the last acquisition. However, during the winter campaigns where more data is acquired this 5 lag may be larger. The velocity provided at every grid point in the PROMICE ice velocity product is the weighted average of all velocity measurements available at that grid point within that 24-day period (see Sect. 4.3), and should be considered an average estimate of velocity over the 24-day period during which the radar images were acquired (see further discussion in Sect. 6). The start and end times of this period are given in the time_bnds variable in the PROMICE NetCDF product (see below). Figure 1 shows samples of the time series at different times during the year 2020.  Table 1. A quick look image for each mosaic is provided along with the dataset.

15
In the following, we present the characteristics of the Sentinel-1 data and introduce the input data that we use to generate the PROMICE ice velocity product.

Sentinel-1 SAR Data Characteristics
SAR sensors are well suited for polar observations because data collection is not impacted by the polar night or cloud cover.
The Sentinel-1 constellation currently consists of two satellites, Sentinel-1A and Sentinel-1B, equipped with identical C-band 20 (5.4 GHz) SAR sensors. Over the GrIS, the Sentinel-1 SAR mainly employs the Interferometric Wide Swath (IW) mode (De Zan and Guarnieri, 2006) allowing for generation of radar images with a resolution of approximately 3 m on ground in the slant range (line-of-sight direction) and 22 m in the azimuth (flight-path direction). The pixel spacing of the product is 2.3 m in slant range, and 14.1 m in azimuth.
The near-polar orbit has a repeat cycle of 175 orbits, corresponding to 12 days, with the two satellite orbits phased 6 days 25 apart. With the current observation schedule, the entire margin of Greenland is imaged every 12 days by both satellites (Fig. 2a).
Furthermore, the entire ice sheet is mapped from several additional tracks every winter from December to February, allowing Image pairs can be formed between acquisitions from the same satellite, i.e. S1A-S1A or S1B-S1B, with temporal baselines which are a multiple of 12 days, i.e. 12-days, 24-days etc. In addition, image pairs can also be formed from acquisitions obtained from two different satellites, i.e. S1A-S1B or S1B-S1A, with a temporal baseline which is an odd multiple of 6-days, i.e. 6-days, 18-days, etc. Although in principle the radar instruments on board the Sentinel-1A and -1B satellites are identical, there are some subtle differences, which should be taken into account when selecting the data pairs for processing, as detailed 5 in Sect. 10 and 5.2.

Input Data
The data used for generating the PROMICE ice velocity product are single-look complex (SLC) IW radar images (with annotation), supplied by the Copernicus Open Access Hub (https://scihub.copernicus.eu/). SLC images are focused SAR images referenced to the radar acquis ition geometry and have the highest resolution of the available IW product types (3 m × 22 m). 10 They are supplied as slices with a footprint of approximately 250 km × 250 km. Owing to the peculiarities of the IW mode, each slice is subdivided into three range swaths, named IW1-IW3, which are acquired in an interleaved (bursted) fashion; these swaths are stored in separate files. The SLC images are supplemented by restituted orbit files, available a few hours after acquisition, and precise orbit files, available 21 days after the data acquisition at https://scihub.copernicus.eu/gnss/. Since the PROMICE ice velocity product is typically generated before the latter become available, we use the restituted orbit files as we 15 have found that the difference between the restituted orbit files and the precice orbit files is insignificant for our data products (see Sect. 5.1).
In order to geocode measurements made in radar geometry, a digital elevation model (DEM) is used. We employ the Greenland Ice Mapping Project (GIMP) DEM based on the ASTER and SPOT 5 DEMs and AVHRR photoclinometry (Howat et al., 2014(Howat et al., , 2015, downsampled to 500 m spacing to match the resolution of the ice velocity product. 20

Methods
The data processing is carried out using the IPP processor, developed and maintained by DTU Space (Kusk et al., 2018).
Despite the name, the processor also performs offset-tracking for displacement measurements, which is the functionality used to generate the PROMICE product. It is a highly automated processing chain requiring little user intervention. The processing is described in the following section, and a high-level overview of the processing flow is shown in Fig. 3a. 25 To support the PROMICE product generation, a database with all available SLC products over Greenland is maintained. This is updated daily by searching the Copernicus Open Access hub. An automated system downloads all new SLC data to a central storage location. Product generation is initiated by an operator selecting a Sentinel-1A reference orbit cycle number.
All SLCs from Sentinel-1A in that and the following cycle (24 days of data) are first selected for processing. Additionally, all SLCs from Sentinel-1B acquired within the same 24-day timespan are selected for processing. Then, all possible SLC pairings 30 with a 6-or 12-day baseline are calculated, and the offset-tracking processing described in Sect. 4.1 is automatically carried out for each pair. When all pairs required for a product have been processed, the geocoding and error estimation described in Sect. 4.2 and 5 is performed for each pair, followed by fusion and mosaicking of all the pairs, as described in Sect. 4.3.

Offset-tracking
The offset-tracking procedure employed is similar to the one described in Strozzi et al. (2002) and estimates local shifts between two SLCs in radar geometry using normalized cross-correlation of intensity image patches.It is illustrated in Fig. 3b.

5
The SLC with the earliest acquisition time is used as the reference image. Prior to the processing, calibration constants in range and azimuth timing are applied to correct for the different geolocation biases observed in Sentinel-1A and Sentinel-1B SLCs (see also Sect.5.2). For 12-day pairs, these constants are identical for both SLCs and have no effect.
An output grid is defined in the reference image geometry, with a spacing of 40 pixels in the slant-range direction and 10 pixels in the along-track (azimuth) direction. These spacings correspond to approximately 150 m×150 m on the ground.

10
For intensity cross-correlation methods, the SLCs need not be coregistered and resampled to sub-pixel accuracy prior to the processing, as they rely on relatively large windows (several tens of pixels in each dimension). Instead, for the regular grid points in the reference SLC, we calculate the expected position of the corresponding grid points (assuming no motion) in the second SLC, based solely on SLC timing information, orbital state vectors, and the DEM. The grid points are selected on integer pixel positions in the reference SLC, but the corresponding grid points will generally not coincide with integer pixel locations 15 in the second SLC. To avoid resampling the second SLC, we round the corresponding grid points to their nearest integer pixel locations and save the fractional shifts, which are then added back to the offset measurement after cross-correlation.
At each grid point, surrounding image patches of 256×64 complex pixels (slant range×azimuth) are extracted in both SLCs.
This patch size has been chosen to maximize the coverage over different flow regimes and coherence levels, and means that the product has an effective spatial resolution on the order of 800-900m. As shown in (Boncori et al., 2018), an adaptive window 20 size approach similar to the one described in (Joughin, 2002) could provide for a locally finer spatial resolution when the data allows it, but this is currently not implemented in the IPP processor Each patch is deramped (Miranda, 2017), upsampled by a factor of two (in both range and azimuth) using FFT interpolation, and the intensity (magnitude squared of the complex pixel values) is derived. A normalized cross-correlation of the two upsampled real-valued patches is carried out, resulting in a correlation surface with values between 0 and 1. The integer shift between the two patches is then estimated by locating 25 the peak of the correlation surface, and a 9x9 neighbourhood surrounding the peak is upsampled by a factor of 4, again using FFT interpolation. Then the fractional shift is retrieved by fitting a parabola to the peak and its two surrounding pixels in each dimension, correcting finally for oversampling factors and accounting for the fractional shift initially estimated for the second SLC. A signal-to-noise ratio (SNR) for the peak estimate is calculated by dividing the correlation value of the peak with the mean of the surrounding pixels in the correlation surface (de Lange et al., 2007). The estimated 2-D shift, the peak normalized 30 cross correlation value (NCC), and the SNR are all saved for further processing.
The procedure described above will yield a shift estimate even if the two images are completely uncorrelated, so a culling of the estimated shifts is carried out. First, pixels with an NCC<0.05 or SNR<7 are set as invalid. Then, a further culling, is carried out on the on the range and azimuth shifts, using a normalized median test, as described in (Westerweel and Scarano, 2005).
For each measurement, U 0 , in a 5×5 neighbourhood, the median, U m , of the 24 surrounding measurements (U 1 , U 2 , . . . , U 25 ) is calculated (excluding U 0 ), and for each measurement in the neighbourhood, a residual, R i = |U i − U m | is calculated. The median, R m , of (R 1 , R 2 , . . . , R 24 ) is then calculated, and used to normalize the residual of U 0 so that R 0 = |U 0 −U m |/(R m + ), where is a minimum normalization level that accounts for cross-correlation noise. We use = 0.1 pixel, as suggested in (Westerweel and Scarano, 2005), and cull the measurement,U 0 , if R 0 exceeds a threshold of 5 for either of the range or azimuth 5 shifts. This value was found by experiments to remove most clearly visible outliers, without removing valid measurements.
Lower values removed more outliers, but had an adverse effect on measurement coverage. After the culling, small unconnected segments of pixels (<25 pixels) are removed, as these were found to often contain erroneous values. Some outliers may remain, especially in areas subject to surface melt, as the associated strong radar backscatter can create false correlation peaks. An additional culling based on time series statistics (see Sect. 4.3) is carried out on the final mosaicked product to further suppress 10 these outliers.
To aid in error estimation, the local standard deviations of the two shift maps (range and azimuth shifts) are estimated in a sliding 5×5 window, ignoring pixels with invalid measurements, and the shift maps are finally averaged by a 5×5 window.
The shift maps (in units of SLC pixels) and associated standard deviations are stored along with the SLC parameters and orbit information to be used in the subsequent processing.

Geocoding and horizontal velocity projection
The geocoding takes as input the shift maps and associated standard deviation maps output by the offset-tracking and the DEM.
Using the DEM and orbit information, the maps are interpolated to the output grid in map projection (see Sect. 2). The shifts and standard deviations are converted to velocity by multiplying with the SLC pixel spacing and dividing by the temporal baseline. At this stage, the velocities and standard deviations, even though provided on a georeferenced grid, are still measured 20 in the radar range/azimuth geometry. With a single pair providing only two velocity measurements, it is not directly possible to estimate three-dimensional flow. Instead, we assume surface parallel flow (SPF) and estimate the flow as described in the T be the three-dimensional velocity vector in map geometry, and v SAR = [v r , v a ] T be the velocity vector in radar geometry, with v r the range (line-of-sight) velocity and v a the azimuth (along-track) velocity. With the SPF assumption, the vertical velocity component becomes (Joughin et al., 1998): where ( ∂z ∂x , ∂z ∂y ) is the surface gradient derived from the DEM. The effective resolution of the velocity maps is on the order of the correlation window size, which corresponds to approximately 800 × 900 m on ground (see Sect. 4.1), so the resolution of the surface gradient map should approximately match this. The DEM is downsampled to the pixel spacing of the PROMICE product (500 × 500 m), and the gradient is derived using second order differences, which means the gradients are derived using 30 samples approximately 1000 m apart. The relation between the horizontal and the radar velocity can be written as: where angles φ and θ describe the orientation of the line-of-sight (LoS) vector pointing from the pixel under consideration to the sensor, with the horizontal angle φ measured counter-clockwise from the y-axis of the map projection and the elevation angle θ measured from the local horizontal plane to the LoS vector. The horizontal velocity components (and the associated standard deviation maps) can then be found by inversion of Eq. 2. Projection scaling factors are not applied to the velocities, so these represent physical velocities along the projection axes.

Fusion
The fusion step describes the process of combining and mosaicking the geocoded offset-tracking results on to a Greenlandwide grid. For every pixel on the output grid, we do a weighted averaging of the N valid velocity measurements from all pairs covering the pixel, using as weights the inverse of the measurement variances: wherev is the fused (x or y) velocity, v n is the (x or y) velocity measurement from pair n, and σ n its associated standard deviation. The estimated standard deviation of the pixel is then:

Culling
After all measurements have been fused and mosaicked, temporal culling is carried out to remove further outliers. This relies 15 on comparison of the measured value with an average value of all available measurements, based at the time of writing on more than 4 years of data. For each pixel, we reject the measurement, if: is the average velocity, v is a velocity constant preventing erroneous culling in areas with very low velocities, and k thr is a constant factor, setting the threshold for culling. A low value 20 of k thr will remove more outliers, but may also remove valid measurements in areas with strong seasonal variation, such as glaciers with significant speed up during the melt season. This effect is showcased for ice velocity along the flowline from Hagen Brae in north Greenland (Fig. 4a). The slow flowing outlet glacier experiences periods of speed up during summer, where velocities near the terminus increases more than 200%. At the same time surface melt inhibits processing parts of the data resulting in spikes in the ice velocity as evident in Fig. 4b. Figure 4c and d show how values of k thr =3 and 1 cull the data. 25 For k thr =3, the (real) summer speed up near the front is conserved, while the majority of spikes further inland are removed.
Applying a stricter value, k thr =1 removes not only outliers but also the real signal due to summer speed-up. In this case, 8 % of the pixels are culled while 4 % of the pixels are culled when applying k thr =3. For the PROMICE ice velocity product we apply, v = 20 m/yr, whereas k thr = 3. This choice of threshold is a balance between removing as much noise as possible without removing actual signal in the mosaics encompassing a wide range of ice dynamics. Figure 5 provides an example of the effect of applying k thr = 3 to a map from summer 2018 when surface melt influences the data quality. The unculled and culled maps are displayed in Fig. 5a and b. The location of the culled data points is shown in red in Fig. 5c. Note for example the removal of the noisy areas in the west and east Greenland ablation zone ( Fig. 5c, d, and 5 e) and locations influenced by ionospheric stripes in the slow moving interior. On average, ∼ 2 % of the pixels in a mosaic are culled using this procedure. More pixels are culled in summer mosaics than in winter mosaics (see Fig. 10a).

Error Sources and Estimation
In this section, we describe the error sources affecting the PROMICE ice velocity product in more detail. The error sources can be divided into five main groups: 10 1. Slowly varying errors, such as those caused by orbit errors (Sect. 5.1) or other timing biases in the products (Sect. 5.2).
2. Temporal decorrelation caused by changes in the radar backscatter between observations, see Sect. 5.3.
3. Ionospheric errors, resulting in localized, but spatially correlated errors in the measured azimuth shift, see Sect. 5.4.
4. Aliasing errors caused by the need to acquire two observations from which we infer displacement, and then velocity.
Any extreme velocity highs or lows will be smoothed out.
15 5. Errors due to tidal motion of floating glacier tongues.

Orbit Errors
Errors in the Sentinel-1 orbital state vectors provided by ESA will result in an apparent shift between the two SLCs in a pair, translating directly into biases on the velocity measurement. For Sentinel-1 data, absolute orbital errors are on the order of 5cm RMS when using the precise orbit product available after 21 days (Peter et al., 2017). The restituted orbits typically used in 20 the PROMICE product generation are available shortly after acquisition, and have a nominal accuracy of 10 cm RMS. This corresponds to 8.6 m/yr RMS for a measurement using a 6-day pair and 4.3 m/yr for a 12-day pair. To assess the difference between using restituted and precise orbits, we processed 18 different Sentinel-1 12-day pairs (9 S1A/S1A and 9 S1B/S1B pairs) acquired consecutively over an area in southwest Greenland where much of the scene in the IW1 swath consists of bedrock. The reason for using only 12-day pairs is to exclude the effect of S1A-S1B biases, which are treated instead in 25 Sect. 5.2. All pairs were processed twice, using either the precise orbits or the restituted orbits, with all other parameters identical. The processing carried out consisted of offset-tracking and geocoding (see Sect. 4.1 and 4.2). Averaging for each processed pair the measured range and azimuth velocities over the bedrock area, which can be assumed stationary, gives an estimated average residual velocity error for that pair. The mean of the 18 residual range and azimuth velocity estimates and associated standard deviations are listed in Table 2. We note that the range velocity bias (-0.5 m/yr) and azimuth velocity 30 bias (0.5 m/yr) do not differ between the precise and the restituted orbit files. The standard deviation in the range direction is 2.7 m/yr for the precise orbit files and 2.8 m/yr for the restituted orbit files, while the standard deviation in the azimuth direction is 3.3 m/yr for both orbit types. Overall, the error statistics for the two orbit types are almost completely identical, and the use of restituted versus precise orbits has an insignificant impact on the accuracy of the final velocity products.

5
With the commissioning of Sentinel-1B in late 2016 it became possible to generate ice velocity products with a 6-day temporal baseline by combining Sentinel-1A and Sentinel-1B data (see also 3.2). Although the SAR instruments are in theory identical, our analysis of the initially generated 6-day ice velocity products revealed velocity biases not present in 12-day (same satellite) pairs. To quantify this, an experiment was carried out using 37 Sentinel-1 pairs (9 6-day A/B pairs, 10 6-day B/A pairs, and 18 12-day A/A and B/B pairs grouped together). The 18 12-day pairs are identical to the pairs used in 5.1, thus, all pairs were 10 acquired over an area in southwest Greenland where much of the scene in the IW1 swath consists of bedrock. Averaging the measured range and azimuth velocity over the bedrock area, gives an estimated average residual velocity error for that pair. Figure 6a shows a scatterplot of the residual range velocity (V r ) versus residual azimuth velocity (V a ) for the 37 pairs. The corresponding mean and standard deviation values are listed in Table 3.
As expected, the 12-day statistics are identical to those observed in the orbit type comparison (see 5.1 and Table 2), with 15 a bias magnitude below 0.5 m/yr. For the 6-day pairs (middle part of Table 3), the bias magnitudes are significantly larger and change sign, depending on whether the first SLC in the pair is acquired from Sentinel-1A or Sentinel-1B. The standard deviations of the 6-day bias estimates are approximately two times those of the 12-day estimates, which is expected, as the velocity measurements are based on measurements of shifts between the images, which are then divided by the temporal baseline to arrive at velocities. The average bias magnitudes for the 6-day pairs are 8.8 m/yr in range and 28.8 m/yr in azimuth. 20 With the 6-day baseline, this corresponds to bias magnitudes on the measured shifts of 0.15 m in range and 0.48 m in azimuth.
These values are consistent with results obtained in detailed analysis of Sentinel-1 SLC product geolocation using corner reflectors, see (Schubert et al., 2017) and (Gisinger et al., 2020). The latter reports average shifts between Sentinel-1A and Sentinel-1B of 0.16 m in range and 0.40 m in azimuth, corresponding to velocities of 9.7 m/yr and 24.4 m/yr, respectively, for measurements using 6-day pairs, but also suggests that there may be a swath dependence of these shifts. Our analysis above 25 concerns only the IW1 swath, so for now, we use the constants from (Gisinger et al., 2020) mentioned above to calibrate the PROMICE product. The calibration is implemented as an adjustment to the timing annotation for the SLC products prior to the offset-tracking. Applying these calibration constants to the test dataset described above results in a significantly reduced bias on the 6-day measurements, as shown on Fig. 6b and in the bottom part of Table 3. The calibrated 6-day range velocities now have a mean bias magnitude of 0.8 m/yr, and the azimuth velocities a mean bias of 6.7 m/yr. In the final PROMICE ice 30 velocity product, the weighted averaging of 12-day pairs and both A/B and B/A 6-day pairs will tend to reduce the impact of any residual biases (see also 4.3).

Temporal Decorrelation
Temporal decorrelation is caused by changes in radar backscatter between acquisitions that reduce the correlation between the image patches which are cross-correlated in the offset-tracking procedure, leading to noisy or even missing measurements. The surface of the interior of the ice sheet is relatively homogeneous, with no large scale features, and the velocity measurement relies on preservation of the speckle pattern (coherence) between observations (Gray et al., 1998). Speckle is a property of 5 radar images, caused by variations in the sub-resolution structure of the imaged scene, resulting in large pixel-to-pixel intensity fluctuations in otherwise homogeneous areas. If the ice-flow is spatially uniform and the sensor track does not deviate excessively for the two acquistions (the latter is generally not a problem for Sentinel-1), the speckle pattern can be tracked between acquisitions using the cross-correlation procedure described in 4.1. Precipitation, surface melt, and steep spatial gradients in ice flow can all reduce the coherence and thus the ability to measure ice velocity in such areas. Often in the interior, the 10 signal-to-noise ratio is low, but since five velocity maps from each track are averaged to produce the PROMICE product, the noise can be reduced. In extended homogeneous areas of low coherence, the velocity measurements can become patchy, since many unreliable measurements will be discarded by the culling procedures described in Sect. 4.1 and 4.4.
On outlet glaciers, the rapid ice flow and associated deformation tends to destroy the coherence except for short temporal baselines. Here, the ability to measure displacement relies instead on the presence of larger scale features, such as crevasses,

15
which can be still be tracked between images with the cross-correlation procedure described in Sect. 4.1, even if there is no coherence.
Models that express the shift errors as function of coherence do exist (De Zan, 2014), but the coherence cannot directly be used to estimate errors or discard measurements, since velocity measurements can often still be made in non-coherent, fast moving areas, as mentioned above. In the PROMICE ice velocity product, the velocity error estimate is based instead on the 20 local standard deviation of the tracked shifts (see Sect. 4.1 and 5.6).

Ionospheric Errors
Ionospheric propagation errors arise due to spatial fluctuations (scintillations) in the ionosphere total electron content within the SAR synthetic aperture length (i.e. km-scale variations) (Gray et al., 2000). This is especially a problem in the near-polar regions. For a given image pixel, these fluctuations cause an azimuth variation in the raw signal phase, which is not accounted 25 for by the SAR focusing, resulting in an azimuth shift of the focused pixel. The varying propagation naturally also causes a shift in the range direction, but these shifts are much smaller (typically on the centimeter-level) than those observed in the azimuth direction (comparable to the azimuth pixel size, i.e. several meters (Mattar and Gray, 2002)). The shifts vary along the scene according to the ionosphere conditions along the satellite flight path, often present in only parts of the scene. Also the observed shifts are strongly correlated in the range direction, appearing as linear or slightly curved "streaks" superposed on the azimuth 30 shift map. In the PROMICE velocity mosaics, such streaks are readily identifiable by the human eye, appearing as roughly east-west oriented stripes of varying intensity. The velocity errors caused by ionosphere can exceed 200 m/yr, impacting 6-day pairs twice as much as 12-day pairs, since the shifts caused by the ionopshere do not depend on the temporal baseline. An example of the impact of the ionosphere can be seen in Fig. 7 showing a single-pair 6-day velocity measurement and a 12day velocity measurement, the former exhibiting significant ionospheric streaks. Both pairs share a common SLC (acquired 2016-10-11), suggesting in this case that the ionosphere effects can be attributed to the other SLC of the 6-day pair.
Methods for reducing the impact of the ionosphere on ice velocity measurements typically rely on the dispersive nature of the ionosphere delay, and have been applied to L-band interferometric ice velocity measurements (Liao et al., 2018). A method for 5 correcting azimuth shift measurements in Sentinel-1 data has been proposed in (Gomba, 2018), but has not been demonstrated for Sentinel-1 ice velocity measurements. Another method to reduce the impact of ionospheric effects is to exploit the fact that in some regions, measurements from both ascending and descending tracks are available, and in this case, ice velocities can be derived from only the range offsets -which are much less sensitive to ionospheric effects -and the SPF assumption (see Sect.4.2). In the standard S1 acquisition scenario (Fig. 2a), only two ascending tracks are acquired (the long track along the west coast of Greenland, and the track covering the northeast margin of the ice sheet up to the northernmost point), so the method is not always applicable everywhere. During the winter campaigns ( (Fig. 2b), this method will be applicable over a much larger part of the icesheet. Work is undergoing to include this method in a future update of the PROMICE product.
The mitigation of ionospheric effects in the PROMICE ice velocity product relies on culling and averaging. Pixels with large ionospheric errors, if present in regions with generally low velocities, will be removed by the temporal culling procedure 15 described in Sect. 4.3. In areas where multiple velocity observations are available, the weighted averaging in the fusion (see Sect. 4.3) will tend to reduce, but not completely remove, the ionospheric effects. We estimate that the ionospheric effects can cause a velocity error of up to 300 m/yr and will mainly affect the v y -component, which is roughly aligned with the azimuth direction due to the near-polar orbit of the Sentinel-1 satellites.

20
A few outlet glaciers in Greenland (Petermann and 79 Fjord being the most prominent) are characterized by having a floating tongue subject to tidal motion. The tidal motion introduces a vertical shift, with the sensitivity of this shift to the tidal signal increasing from 0 near the grounding line to 1 on the fully floating part of the tongue, a transition zone which is typically 5-10 km wide (Padman et al., 2018). This vertical shift will affect the ice velocity estimate, as the difference in the shift between the two radar acquisitions is projected on to the radar line-of-sight and interpreted as motion in the slant range direction. In 25 (Reeh et al., 2000), tidal-induced shifts of approximately +/-0.5 m were observed using GPS receivers placed on the floating part of 79 Fjord glacier. This could, for a 6-day pair, lead to errors in the ice velocity estimate of more than 50 m/yr, although the averaging of several acquisitions in the PROMICE product will tend to reduce this error. The effect is not modeled in the PROMICE product, so care should be taken when using the product on floating glaciers.

30
The error estimates provided with the PROMICE ice velocity product are derived from the local standard deviation of the underlying shift maps generated by the offset-tracking (see Sect. 4.1 and 4.3). As such, they do not account for slowly varying errors, such as those described in Sect. 5.1 and 5.2, and only to a limited extent for the impact of ionospheric errors, as these are locally correlated on the scale of the window size used to estimate the local standard deviations. Although this is not a complete error characterization, it was shown in Boncori et al. (2018) to provide the correct order of magnitude for the errors. Examples of relative error estimates accompanying the two ice velocity maps from Fig. 7 are shown in Fig. 8. The strong ionospheric streaks evident in the 6-day pair on Fig. 7 are seen to be reflected in the corresponding error estimate, although the magnitude is underestimated. In the central and lower left part of the maps, errors are seen to be generally higher on the 12-day pair, but 5 in a more diffuse pattern, even though the 12-day pair is less sensitive to a given shift error, due to the longer baseline. In this case, it is the higher temporal decorrelation of the 12-day pair that causes an increased noise level, which is also reflected in the error estimate. A Greenland-wide view of the error estimate for the PROMICE product is given in Fig. 9 for the same mosaics displayed in Fig. 1.
The slowly varying errors (Sect. 5.1 and 5.2) could potentially be corrected by calibrating the measured velocities using 10 ground control points (GCPs), either on stable terrain or in areas where the ice flow is known to vary little. In practice this is difficult to do in an automated system, as the calibration has to be carried out on the individual pairs, where the ionospheric and, to some extent, the temporal decorrelation errors associated with offset-tracking are often much larger than the slowly varying errors. If GCPs are unwittingly selected in areas affected by e.g. ionosphere, the GCP calibration could actually have a detrimental impact. A large number of GCPs, well distributed in the image, would be required to reduce the statistical noise, 15 but this can often not be achieved within the limited spatial coverage of a single pair. For this reason, the PROMICE ice velocity product is not calibrated using GCPs.

Properties of the PROMICE Ice Velocity Product
The PROMICE ice velocity product is designed as a compromise between good spatial coverage, high temporal resolution, and low noise. Other combinations of 6 day and 12 day pairs are possible resulting in a different temporal resolution and spatial 20 coverage. We explore other possibilities for products and compare them to the PROMICE ice velocity product with respect to coverage and noise. These products are time series of mosaics consisting of (see Table 4): 1. All 6 day pairs (no 12 day pairs) within 2 Sentinel-1A cycles (6dOnly).
Data coverage for each mosaic in each time series is displayed in Fig. 10a. We defined the coverage as the fraction of grid points that contains ice velocity data on the ice sheet in a given mosaic. We have included a time series in the analysis called All-pairsNoCull, which includes the same data as the PROMICE product, but has not undergone the culling procedure described in Sect. 4.4.The percentage of points that are culled in each mosaic of the PROMICE product using that method 30 is also shown in Fig. 10a. If all grid points on the ice sheet contains data then coverage is 1. All time series have close to full coverage during peak winter, where a campaign ensures full IW coverage of the ice sheet over a number of cycles. The coverage of the PROMICE product drops to ∼ 0.7 outside the campaigns with a low during summer months. Figure 10b shows the mean of the reported standard deviation for each mosaic in the time series. In general, when coverage is low the noise goes up and vice versa. The 6dOnly_1cycle is the time series with highest reported error estimate. All four time series have a large low-coverage area in the ice-sheet interior, where SAR data in IW mode is rarely acquired as is evident from Fig. 2. The same explanation is true for the smaller triangular areas in the Melville Bay area and northern Greenland as well as the Scoresbysund area (locations 1, 2, and 3 in Fig. 11a). However, the large area with low coverage along the southeast ice sheet margin as well as an area in southern Greenland, one north of Rink Glacier in west Greenland, and one 15 in the Melville Bay area all have routine SAR IW acquisitions every 6 days (locations 4, 5, 6, and 7 in Fig. 11 a). This will be discussed in the following.
The coverage and quality of each mosaic depend both on the SAR data coverage, the number of acquisitions going into each mosaic as well as on how the properties of the ice-sheet surface have changed between acquisitions (Sect. 5, Sect. 4.3 and Fig. 2a). The PROMICE product has the best coverage of the time series in Table 4 (excluding All-pairsNoCull) ( Fig. 10 and 11) 20 as it includes all the pairs contained in both 12dOnly and 6dOnly. Figure 10a shows that most often the 6dOnly has better coverage than the 12dOnly and for some extended periods of time, it is comparable to the PROMICE product. However, both 12dOnly and 6dOnly time series have periods with significant drops in coverage, while the PROMICE product still performs well. The difference in coverage is caused by differences in data acquisition and ice-sheet surface properties.
Not all tracks have both 12 day and 6 day coverage and often tracks in the interior are only covered by 12d pairs. This is 25 revealed by the lower coverage of the interior for 12dOnly compared to 6dOnly in Fig. 11. 6dOnly has better coverage along the ice sheet margins compared to 12dOnly, because coherence is more likely to be preserved for shorter temporal baselines as discussed earlier in Sect. 3.1. The PROMICE product mosaics thus have better coverage than both 6dOnly and 12dOnly, because they each have coverage where the other does not. However, even using the short 6 day temporal baseline, some areas consistently have low coherence and therefore rarely have ice velocity coverage. The largest of these areas is the southeast ice 30 sheet margin while the smaller areas include an area in southern Greenland, one north of Rink Glacier in west Greenland, and one in the Melville Bay area (Fig. 11). The areas are apparent in all time series, but are most pronounced in the 12dOnly series as a longer temporal baseline increases the probability of changes to the surface properties due to precipitation and/or surface melt. The areas discussed here largely coincide with the regions identified as high accumulation percolation areas (HAPA) by Vandecrux et al. (2019) studying firn properties. HAPAs are areas on the ice sheet characterized by frequent precipitation events and surface melt-water that percolates into the firn -both processes leading to loss of coherence. Figure 10 also shows the effect of performing the culling described in Sect. 4.4 on the time series as a whole. The All-pairsNoCull is the PROMICE Product without culling. The effect on the coverage is minor (Fig. 10a), however Fig. 10b shows that the average noise level is significantly reduced.

5
In the PROMICE ice velocity product each mosaic includes all possible 6-and 12-day pairs within two consecutive Sentinel-1A cycles and the timestamps supplied with the product lists the timespan of the product (first and last date) as well as the midpoint time as specified in Table 1. This information is true for the mosaic, but not for a given grid point. This is due to: -The SAR data (Fig. 2) is not acquired simultaneously over the GrIS as described in Sect. 3.1.
-Different areas on the ice sheet are covered by a varying number of tracks/varying amount of data acquired at different 10 times (also Fig. 2).
-Although data is acquired, the processor is unable to detect displacement for some pixels or larger areas due to loss of coherence or the processing of an image pair fails for various reasons and is therefore not included in the final mosaic.
Another point to keep in mind is that the mosaic is a weighted average of the processed pairs. This means that although the product has a high temporal resolution, the time series will be smoothed and likely miss short lived (real) peaks in velocity for 15 instance during summer.
The analysis from this section shows, that it is possible to provide a Greenland-wide ice velocity product with a higher temporal resolution than the PROMICE product (the 6dOnly_1cycle product), but also that this comes with the price of reduced spatial coverage and higher uncertainty. Creating a product spanning more than two Sentinel-1A cycles will have opposite effects. The two Sentinel-1A cycles choice for the PROMICE product is therefore a compromise between having reasonably 20 high temporal resolution and good coverage and reducing noise.

Validation
We validate the PROMICE ice velocity product against in-situ GPS measurements from the PROMICE automatic weather stations (AWS) (van As et al., 2011;Fausto and van As, 2019) and perform an analysis over stable ground. Only a limited number of GPS measurements are available since the data must overlap in time with the period of the PROMICE ice velocity 25 product and have a temporal resolution comparable to or higher than the PROMICE ice velocity product. Furthermore, the measurements are biased toward the slow moving parts of the ice sheet ablation zone. Locations are displayed in Fig. 12.
PROMICE AWSs measure a range of surface mass-balance components in the ablation zone of the GrIS. The stations are per design located in slow moving areas with an average flow generally lower than 100 m/yr (Fig. 12). The position of the AWS is measured every hour using a single frequency GPS receiver and a small ceramic patch active antenna. We use the 30 freely available hourly positions ((Fausto and van As, 2019)) to calculate velocities using a workflow similar to that described in GIScci-Consortium (2018): Daily positions of the GPS stations are calculated as a mean of the hourly positions for each day.
The velocity components are estimated using a weighted linear regression for each of the 24-day time spans of the velocity mosaics using the daily positions. The weights are inversely proportional to the number of hourly measurements going into the estimate of a daily position in order to account for gaps in the data.
Scatter plots of the satellite derived PROMICE ice velocity product (magnitude, v x and v y components) vs. PROMICE GPS 5 derived ice velocities are displayed in Fig. 13. The standard deviation of the difference between the GPS measurements and the satellite derived velocity (from here on referred to as the standard deviation) is calculated along with the mean difference (bias) between GPS and satellite velocity (see first line in Table 5). These values reflect not only the uncertainty of the satellite product but also that of the GPS derived velocity. The expected error of the satellite product is estimated to be 10-30 m/yr for individual pairs (GIScci-Consortium, 2013). The PROMICE product lies well within these bounds with a standard deviation 10 and bias of 20 m/yr and -3 m/yr for the v x -component and 27 m/yr and -2 m/yr for the v y -component, respectively. The larger standard deviation for the v y component is expected: Due to the general north-south orientation of the satellite tracks (Fig. 2), the v y component is aligned roughly parallel to the azimuth direction (satellite flight path), and Sentinel-1 IW SLC images (see Sect. 3.1) have a much lower resolution in azimuth than in the range (line-of-)sight direction). We note that due to the roughly east/west orientation of most Greenland glaciers the velocity range of the y-component in our validation is notably smaller than 15 that of the x-component.
We perform a similar analysis for the PROMICE product for the pixels over stable ground, where no movement is expected.
All pixels on ice-free terrain from all mosaics (each spanning 24 days) in the time series are included, which totals to more than 142 · 10 6 pixels. The resulting values of the standard deviation and bias are 8 m/yr and 0.1 m/yr for the v x -component and 12 m/yr and -0.6 m/yr for the v y -component, respectively (see Table 6). The values of the standard deviation are less than half 20 the values of the validation against GPS measurements, while the biases are significantly lower (in absolute value) and thus closer to zero. In this analysis, the standard deviation and bias are also largest for the v y -component as discussed above. shows, but similar to the analysis carried out on stable ground in Sect. 5.2 (see Fig. 6). This difference is mainly due to two things: 1) In the accumulation zone, changes at the ice sheet surface are most often due to snow fall and redistribution by wind.

30
In contrast, in the ablation zone, where all the PROMICE GPS stations are located, the surface properties are influenced by several factors e.g. melt, high accumulation rates, and rain. This influences the coherence of image pairs and thereby increases uncertainty in the velocity product in these areas.
2) The uncertainty of the GPS measurements reported in Hvidberg et al. (2020) is lower compared to the PROMICE GPS observations. This is due to both the longer temporal baseline between measurements as well as the data acquisition time of 2-4 hours per data point. The velocities derived from the PROMICE GPS 35 observations may therefore carry a non negligible part of the uncertainty in the validation. None of the other products in the comparison by Hvidberg et al. (2020) have a similar high temporal resolution as the PROMICE ice velocity product. However, a 3-year average of the PROMICE ice velocity product was also included in the analysis and Hvidberg et al. (2020) found that it had a standard deviation (0.7 m/yr) similar to other offset/feature tracking products covering longer timespans like the annual maps from ESA CCI (Greenland Ice Sheet velocity maps from Sentinel-1 (Nagler et al., 2015)) and MEaSUREs Greenland average of the PROMICE product is lower than for the annual maps, most likely because it includes more data. This is also a conclusion drawn by Hvidberg et al. (2020).
The 12dOnly and 6dOnly time series introduced in Sect. 6 have similar standard deviations as the PROMICE product, when 10 compared to the velocity derived from the PROMICE GPSs, whereas the higher temporal resolution product, 6dOnly_1cycle, has a significantly higher standard deviation. The PROMICE ice velocity product has the lowest standard deviation of the four. Using only 6 day pairs it is also possible to define a Greenland-wide product with a temporal resolution of 12 days, -the 6dOnly_1cycle product. It has the clear advantage of resolving the dynamics of the outlet glaciers even better, although this comes with the price of increased noise due to both the shorter temporal baseline and the geolocation bias as well as reduced 15 coverage of each mosaic ( Fig. 10a and b, Fig. 11 and Table 5). For studies concerned with changes in fast flow the increased temporal resolution may outweigh the downsides.
The uncertainty reported in the ice velocity product is lower than the values we found during our validation against GPS.
The average standard deviation found in Hvidberg et al. (2020) or the stable ground analysis is more comparable. The origin of some errors is such that the algorithm is unable to account for them. This is especially true for the spatially correlated errors 20 caused by ionospheric scintillations (Sect. 5.4), which are not fully estimated by the error estimation algorithm (Sect. 5.6). A second issue is the distribution of the PROMICE AWSs biased towards the slow flowing parts of the ablation zone as well as the uncertainty on the velocity estimates from these data.
For a time series of mosaics like the PROMICE ice velocity product, errors will vary both spatially and temporally due to the sources described in Sect. 5 as well as to variations in data coverage (Sect. 4.3 and Fig. 2a). A validation dataset which is 25 not biased towards slow flowing areas in the ablation zone but is representative of a larger range of flow regimes and surface conditions would help assess whether the reported product errors capture this correctly. The analysis above, however, shows that the size of the product errors are as expected.
8 Living Data: Updates and Improvements PROMICE will continue to distribute and update the PROMICE Ice Velocity product based on the Sentinel-1 data collected 30 and released by ESA and the Copernicus programme. We aim to deliver a clean and homogenous data product and offer the possibility of user-interaction and addressing issues with the data product. Associated with PROMICE, we have a user-contributable dynamic web-based data archive (GitHub), which list known data quality issues [https://github.com/GEUS-PROMICE/Sentinel-1_Greenland_Ice_Velocity]. On the GitHub page, we also offer the opportunity for data users to add and document new issues. Documenting dataset issues is often simpler than correcting them and future dataset versions will implement fixes to any verified issues as soon as they are done. All fixed issues will be tagged as closed and remain visible for new users.
We encourage users who are working with Sentinel-1 and the PROMICE ice-velocity data to search the issue database and 5 see if there are any known data issues relevant to their needs. We find it likely that there are issues unknown to us in the existing data and new issues may be found in the future data collection pipeline. We will do our best to improve the dataset with user-based help through the GitHub page.

Summary and Outlook
We have presented the PROMICE ice velocity product -a time series of GrIS wide velocity mosaics (September 2016 to 10 present) based on Sentinel-1 SAR data. The product has a 500 m spatial-and 24 day temporal resolution and is produced in an operational setup using the IPP processor. A new mosaic is produced every 12 days and is made available within 10 days of the last included acquisition. During the winter campaigns, this lag is larger due to the amount of data to be processed. Validation against PROMICE AWS GPS data show that the standard deviation of the difference between the ice velocity product and the GPS data is 20 m/yr and 27 m/yr for v x and v y -component, respectively. This is within the expected uncertainty range of 15 10-30 m/yr (GIScci-Consortium, 2013). However, we expect the actual values pertaining to the PROMICE ice velocity product to be lower as the PROMICE AWS GPS data carry a non negligible part of the uncertainty. This is indicated by the analysis carried out for pixels over stable ground, which showed a standard deviation of 8 m/yr and 12 m/yr for v x and v y -component, respectively as well as by the study of Hvidberg et al. (2020). Better spatially distributed validation data with low uncertainty would help assessing whether the processor captures the spatially and temporally varying uncertainty field correctly. 20 Ice velocities are retrieved by applying intensity offset-tracking to Sentinel-1 images acquired 6 and 12 days apart. The resulting velocity maps from all image pairs acquired during a 24-day period are temporally averaged and mosaicked to produce a consistent coverage. The processing chain is described in detail from the input data to the final outlier removal. We discuss the various error sources, which include biases and smoothly varying errors due to orbit and timing errors, noise-like errors due to changes in radar backscatter between radar acquisitions, and errors due to ionospheric scintillations. The error 25 estimation approach is also described.
We show how the product coverage varies temporally and spatially in response to variations in SAR data acquisitions and seasonal changes in surface properties. The southeast GrIS margin has good Sentinel-1 SAR data coverage, but often has gaps in the mosaics due to changes in the surface properties caused by surface melt and high precipitation rates, which hinder velocity retrieval. Other areas, like the small triangular area in the Melville Bay area, have low coverage in the mosaics simply 30 due to lack of SAR data.
The PROMICE Ice Velocity product will continue to update as long as the Sentinel-1 satellites are in operation. We will continue to make improvements to the product, and these updates will be posted at https://github.com/GEUS-PROMICE/Sentinel-1_Greenland_Ice_Velocity. Users are encouraged to add and document product issues or suggest improvements.
The PROMICE ice velocity product presented here was originally intended primarily to calculate ice discharge through marine-terminating glaciers of the GrIS as done in Mankoff et al. (2020). The PROMICE ice velocity product is thus less 5 suited for studying very short-lived changes in the velocity structure, as observed in-situ by e.g. Bartholomew et al. (2012) and Ahlstrøm et al. (2013) or through higher frequency acquisitions/non-mosaic products of satellite imagery as done by e.g. Sundal et al. (2013) and Davison et al. (2020). By not mosaicking all the individual image pairs like we do for the PROMICE ice velocity product, a much higher temporal resolution over a limited region is possible. Yet, the spatially comprehensive and temporally consistent nature of the PROMICE ice velocity product makes it attractive also for longer term large-scale 10 monitoring of the GrIS velocity structure and glacier dynamics as done by Vijay et al. (2019) and .
Author contributions. AS and AK designed and produced the PROMICE ice velocity product. AK, JPMB and JD developed the processing software. RSF and KDM set up the data-curation framework. AS and AK prepared the manuscript with contributions from all co-authors.
Competing interests. Authors declare that they have no conflict of interest.