Surface velocity of the Northeast Greenland Ice Stream (NEGIS): Assessment of interior velocities derived from satellite data by GPS

. The Northeast Greenland Ice Stream (NEGIS) extends around 600 km upstream from the coast to its onset near the ice divide in interior Greenland. Several maps of surface velocity and topography in the interior Greenland exist, but their accuracy is not well constrained by in situ observations. Here we present the results from a GPS mapping of surface velocity in an area located approximately 150 km from the ice divide near the East Greenland Ice-core Project (EastGRIP) deep drilling 15 site. A GPS strain net consisting of 63 poles was established and observed over the years 2015-2019. The strain net covers an area of 35 km by 40 km, including both shear margins. The ice flows with a uniform surface speed of approximately 55 m a -1 within a central flow band with longitudinal and transverse strain rates in the order of 10 -4 a -1 and increasing by an order of magnitude in the shear margins. We compare the GPS results to the Arctic Digital Elevation Model and a list of satellite-derived surface velocity products in order to evaluate these products. For each velocity product, we determine the bias and 20 precision of the velocity compared to the GPS observations, as well as the smoothing of the velocity products needed to obtain optimal precision. The best products have a bias and a precision of ~0.5 m a -1 . We combine the GPS results with satellite-derived products and show that organized patterns in flow and topography emerge in the NEGIS when the surface velocity exceeds approximately 55 m a -1 and are related to bedrock topography. Here, we present results from a geodetic surface program to characterize surface topography and ice flow of an interior section of NEGIS in an area near its onset in North Central Greenland, and to assess remote sensing products from this interior area of the Greenland Ice Sheet. The area is located approximately 150 km from the ice stream onset and centered around a reference 55 stake (75°38’ N, 35°58’ W) located 300 m from the East Greenland Ice-core Project (EastGRIP) deep drilling site. We compare our GPS derived heights and surface velocities with the ArcticDEM elevation model (Porter et al. 2018), as well as with 165 published and experimental remote sensing velocity products from the NASA MEaSUREs program, the ESA Climate Change Initiative, the PROMICE project, and three experimental products based on data from the ESA Sentinel-1, the DLR TerraSAR-X, and USGS Landsat satellites, in order to validate and assess these products (the complete list and references are given in 60 section 4 below). We use the GPS derived horizontal surface velocities and strain rates in combination with the remote sensing velocity products to characterize the ice stream flow, shear margins and structure of NEGIS near its onset. program, the PROMICE program, and experimental data products from MEaSUREs, DTU-Space and AWI. For each product, we calculate the bias, the standard deviation relative to the GPS derived surface velocities, and the spatial smoothing that minimizes the standard deviation. Our assessments show:


Introduction 25
The discharge from Greenland's marine terminating outlet glaciers has increased over the last decades and contributed to the increasing mass loss from the Greenland Ice Sheet (Mouginot et al., 2019;Mankoff et al. 2019;Imbie Team, 2019). During the same period, many outlet glaciers have accelerated and thinned in response to changes in atmospheric and oceanic forcings, thereby adding to the dynamic mass loss (Bevis et al. 2019;Khan et al., 2015). Further dynamic thinning and acceleration in ice flow at marine outlet glaciers can potentially propagate inland, and activate the vast high-elevation and slow moving interior 30 part of the ice sheet, thereby leading to additional mass loss (Mouginot et al., 2019).
2 Fast flowing ice streams drain a significant fraction of the ice from the Greenland Ice Sheet into marine outlet glaciers and they thereby connect the interior parts of the ice sheet with the margins. The fast flow involves basal sliding and friction at the bed and along the shear margins, but the understanding of the mechanisms controlling ice stream dynamics and their connection to the surrounding slow-moving ice is incomplete (Minchew et al., 2019;Minchew et al., 2018;van der Veen, 35 2018, Gillet-Chaulet et al. 2016). In the interior, in situ observations of surface movement are sparse and limited to a few locations (e.g. Hvidberg et al. 1997;Hvidberg et al. 2003), and satellite-derived observations of surface velocity and elevation change are limited by their temporal and spatial resolution and the lack of validation data (Joughin et al., 2017). A small surface thickening is observed since 1995 from satellite altimetry in the interior, but it is not clear whether it is due to increased precipitation or ice dynamical changes (Mottram et al. 2019). As a result, there is a significant uncertainty in the projections 40 of the future response of the interior areas of the Greenland Ice Sheet to changes at the marine outlet glaciers (Imbie Team, 2019; IPCC, 2019).
The North-East Greenland Ice Stream (NEGIS) drains a basin in northeast Greenland with an area of about 16% of the total area of the Greenland Ice Sheet into three main marine outlet glaciers, Nioghalvfjerdsfjorden glacier (NG), Zachariae Isstrøm (ZI), and Storstrømmen Glacier (SG) (Fig. 1). NEGIS extends around 600 km upstream of its outlet glaciers to its onset near 45 the ice divide in the interior of northern Greenland. The inferred mass loss from NEGIS has increased since 2003 (Mouginot et al. 2019). This is mainly due to a rapid retreat of ZI since it lost its floating tongue in 2003 and a slow retreat of NG (Khan et al. 2014;Mouginot et al. 2015), while SG has slowed down after its surge around 1980 (Mouginot et al. 2018). If the marginal loss continues and induces dynamical thinning and acceleration upstream along NEGIS, it could potentially activate the interior parts of NEGIS (Khan et al., 2014;Choi et al. 2017). The onset of NEGIS in the interior may be related to the 50 geothermal heat flux and subglacial drainage system in the area (Karlsson and Dahl-Jensen, 2015), but the sensitivity of the system to the ongoing marginal mass loss is not well known.
Here, we present results from a geodetic surface program to characterize surface topography and ice flow of an interior section of NEGIS in an area near its onset in North Central Greenland, and to assess remote sensing products from this interior area of the Greenland Ice Sheet. The area is located approximately 150 km from the ice stream onset and centered around a reference 55 stake (75°38' N, 35°58' W) located 300 m from the East Greenland Ice-core Project (EastGRIP) deep drilling site. We compare our GPS derived heights and surface velocities with the ArcticDEM elevation model (Porter et al. 2018), as well as with 165 published and experimental remote sensing velocity products from the NASA MEaSUREs program, the ESA Climate Change Initiative, the PROMICE project, and three experimental products based on data from the ESA Sentinel-1, the DLR TerraSAR-X, and USGS Landsat satellites, in order to validate and assess these products (the complete list and references are given in 60 section 4 below). We use the GPS derived horizontal surface velocities and strain rates in combination with the remote sensing velocity products to characterize the ice stream flow, shear margins and structure of NEGIS near its onset.
3 2 GPS data and methods

GPS stake network
The surface program includes a repeated in situ survey using the Global Positioning System (GPS) with a strain net consisting 65 of 63 stakes. Observations cover the years 2015-2019. The stake network was established in 2015 with 16 stakes, including a central reference stake (75°38' N, 35°58' W) located 300 m from the EastGRIP deep drilling site, and gradually expanded in 2016, 2017 and 2018 to include 63 stakes (Fig. 2). The growing network of stakes were measured by GPS every year from 2015 to 2019 and supplemented with additional temporary stakes that were measured only once. The layout of the stakes was designed to provide 1) transects of flow velocities along and transverse to the flow, and 2) longitudinal and transverse strain 70 rates in the center of the fast flow and at both shear margins. To fulfill these requirements, the stake network contains sets of stakes placed in a diamond shape centered around the mid-point of NEGIS and at both shear margins. The stake network extends 35 km along NEGIS and 40 km across NEGIS, thereby covering the entire 25 km width of NEGIS and extending across both shear margins into the slower moving regions outside the ice stream. The purpose of the additional stakes added in 2018 was to obtain detailed information of strain rates across a topographic surface undulation northwest of NEGIS (a 20-75 30 km dark/bright pattern perpendicular to NEGIS, Fig. 2a). All stake observations are included in this analysis.
The GPS observations were carried out with a Leica GX1230 GPS receiver with data acquisition lasting minimum 1 hour, and typically 2-4 hours. The GPS antenna was mounted on the top of each stake, and the height above the surface was measured manually. The stakes were 3.5 m long aluminum stakes, which were drilled approximately 2 m below the surface and extended when needed due to a continuous snow accumulation in the area (approximately 0.3 m of snow equivalent per year; Vallelonga 80 et al. 2014). All stakes established in 2015, 2016 and 2017 were extended during the observational period when the antenna heights decreased below 1 m above the surface. A few stakes were moved and/or replaced due to camp activities.
The GPS observations were post-field processed using the open source software package ESA/UPC GNSS-lab tool (gLAB) (Sanz Subirana et al. 2013;Ibáñez et al. 2018). We use the Center for Orbit Determination in Europe (CODE) final orbit/clock product, which includes Earth rotation parameters. We took the antenna phase center offset and variation into account. Receiver 85 clock parameters are modelled, and the atmosphere delay parameters are modelled using the CODE maps for the ionosphere and ESA's Niell Mapping function with Simple Nominal for the troposphere. We applied solid Earth tidal corrections using the IERS Convention's degree 2 tides displacement model (Sanz Subirana et al. 2013). Ocean tidal correction is not implemented in the gLAB processing tool, and for our interior site the associated error is estimated to be within 1 cm. The coordinates are computed in the IGS14 frame. We use the software in static mode and developed an automated protocol in 90 order to perform a systematic Precise Point Positioning (PPP) processing of the stake observations. The PPP approach can introduce systematic errors if the stake is moving (King, 2004). To optimize our processing protocol and evaluate timing estimates and position uncertainties, we observed the central reference stake at the EastGRIP site (red dot in Fig. 2) over extended periods each season and compared separate 1-hour static, 24-hour static, and kinematic solutions. We found that the 24-h static solution performed better than the average position of a 24-hour kinematic solution. With a maximum observed 95 4 surface speed of approximately 60 m a -1 , the uncertainty related to the static solution is estimated to be < 2 cm. We estimate the combined uncertainty of our GPS positions to be within 3 cm. We process the stake observations from each year, including multiple observations of some of the stakes within the annual field seasons.
The gLAB processing protocol was assessed by comparing processing results from two one-hour observations with PPP processing results from the GIPSY-OASIS version 6.4 software developed at the Jet Propulsion Laboratory (JPL). We use JPL 100 final orbit products, which include satellite orbits, satellite clock parameters and Earth orientation parameters. The orbit products take the satellite antenna phase center offsets into account. Receiver clock parameters are modeled, and the atmospheric delay parameters are modeled using the Vienna Mapping Function 1 with VMF1 grid nominal (http://vmf.geo.tuwien.ac.at/) (Kouba, 2008;Boehm et al, 2006). Corrections are applied to remove the solid Earth tide and ocean tidal loading. The amplitudes and phases of the main ocean tidal loading terms are calculated using the Automatic 105 Loading Provider (http://www.oso.chalmers.se/~loading) (Scherneck and Bos, 2002) applied to the FES2014 ocean tide model including correction for center of mass motion of the Earth due to the ocean tides. The site coordinates are computed in the IGS14 frame (Altamimi et al, 2016). All GIPSY-OASIS processing results were within < 1.7 cm of the gLAB processing results (Table 1) and within our estimated uncertainty of 3 cm.

GPS derived velocities and surface elevations 110
To derive the horizontal surface velocity, a linear fit was performed to the observed Northing and Easting positions, respectively (projected to the National Snow and Ice Data Center (NSIDC) Sea Ice Polar Stereographic North and referenced to WGS84 horizontal datum (EPSG:3413)) assuming a constant displacement rate. A small tilt of the stake can lead to uncertainties in the horizontal velocity. We take this into account by including an unknown horizontal shift in the position of stakes that were vertically extended, and we neglect any other changes in the tilt. For each stake, the shift was determined 115 independently from the other stakes, and the linear fit and the shift were determined simultaneously. The estimated shifts are in the range 0.05 m to 0.2 m, and often exceeds 0.10 m. The surface velocity was calculated by assuming that the flow is along the surface, thereby neglecting vertical movement. We estimate the uncertainty of the derived velocities due to the combined uncertainty of the GPS observations and the method to be in the order of 10 -2 m a -1 . As a horizontal reference position of the stakes, we use the estimated horizontal position of the stakes on 2017-01-01, assuming a constant horizontal displacement of 120 each stake over the observational period. We select a common reference for the network in order to consistently derive horizontal strain rates and assess surface elevations, but the reference date is not an accurate timestamp due to the different initial dates of the stakes.
We also estimate a mean GPS derived surface elevation of the stakes to be used below for the assessment of satellite based observations. In the interior areas of the Greenland Ice Sheet, climate driven variations in snow accumulation and firn 125 compaction lead to seasonal and interannual variations of the surface elevation that are not resolved by our annual GPS observations. We estimate the mean GPS derived surface elevation as the mean of the individual observations at each stake, neglecting trends over the observational period due changes in snow accumulation, firn processes or ice dynamical changes.

5
The resulting horizontal stake velocities are shown in Figure 2, and the reference positions and horizontal velocities are listed in supplementary Table S1. The GPS derived surface elevations and the magnitudes of the horizontal stake velocities, are 130 shown along three transects across NEGIS in Figure 3 and one transect along NEGIS in Figure 4. The stake velocities show that the surface speed is relatively constant at approximately 56.6 m a -1 along the centerline of NEGIS, and is above 55 m a -1 in central flow band wider than 10 km. The GPS derived surface elevations reveal 20 m deep topographic troughs at the shear margins. The direction of the fast flow at the center line is 33.5 from North.

GPS derived strain rates 135
After having derived horizontal surface velocities, we calculate horizontal strain rates, which are essential in understanding the ice flow pattern and internal stratigraphy of the ice stream and its surroundings. We calculate the horizontal principal strain rates in 32 different triangular sections within the GPS strain net. Each triangle is defined by a combination of three GPS stakes and assumes a linear velocity field within the triangle, i.e. constant strain rates within the triangles (Fig. 5). The principal strain rates are generally in the order of 10 -4 a -1 in a wider than 10 km central flow band along NEGIS, as well as in the slow moving 140 areas outside NEGIS. In the two shear margins, horizontal principal strain rates increase by an order of magnitude and reach a maximum in the northern shear margin of 3.8ꞏ10 -3 a -1 (horizontal extension) and -3.6ꞏ10 -3 a -1 (horizontal compression), and in the southern shear margin of 3.6ꞏ10 -3 a -1 (extension) and 4.3ꞏ10 -3 a -1 (compression). In both shear margins, the principal strain rates are orientated at an angle of approximately 35 relative to the direction of the flow, due to a combination of longitudinal extension, transverse compression and a high shear strain rate along the shear margins. The principal strain rates 145 are slightly higher in the triangles north of the central flow line of NEGIS than in those south of the central flow line, probably because the northern shear margin is wider than the southern shear margin and not captured as precisely by the GPS strain net as the southern shear margin. We estimate the uncertainty of the strain rates averaged over the triangles (>2 km) to be in the order of 10 -5 a -1 .
Along a transect across NEGIS, we calculate the horizontal strain rate tensor along the direction of the flow at three stakes. 150 The three stakes are located in the northern shear margin, in the center (EastGRIP site), and in the southern shear margin, respectively (Fig. 5). The strain rates along the direction of the flow are calculated as the mean of the rotated strain rate tensors in the four adjacent triangles, and they are 2 km horizontal averages, corresponding to the dimensions of the adjacent triangles. The normal strain rate components along the direction of the flow at these three stakes are plotted in Figure 3b. At the central flow line at the EastGRIP site, the normal strain rates are (0.90.2)ꞏ10 -4 a -1 in the longitudinal (along flow) direction, 155 and (0.90.5)ꞏ10 -4 a -1 in the transverse direction. The  2km-average longitudinal and transverse strain rates relative to the local flow direction in the northern shear margin are (0.81.9)ꞏ10 -3 a -1 and (0.80.4)ꞏ10 -3 a -1 , respectively. The horizontal shear strain rate is (2.10.9)ꞏ10 -3 a -1 . The 2 km-average longitudinal and transverse strain rates relative to the local flow direction in the southern shear margin are (0.81.6)ꞏ10 -3 a -1 and (1.30.5)ꞏ10 -3 a -1 , respectively, and the horizontal shear strain rate is (2.60.8)ꞏ10 -3 a -1 . The resolution of the GPS strain net is limited by the position of the stakes and may not capture the 6 peak strain rates at the shear margins, but the sharp transition in the southern shear margin stand out in the relative velocity pattern in Supplementary Figure S1 and in the strain rates along the transect across NEGIS (Fig. 3b).

Satellite-derived digital elevation model
The GPS derived surface elevations are used to validate the accuracy of the ArcticDEM release 7. The ArcticDEM is a digital 165 elevation model (DEM) based on stereo auto-correlation techniques on visual imagery from Worldview satellites (Porter et al., 2018). The resolution of the ArcticDEM is 2m with a bias of less than 5m (Noh and Howat, 2015). The timestamp of the ArcticDEM in the EastGRIP area is 2017 (estimated), and the vertical accuracy has not been verified (Porter et al., 2018), thus overlapped by the GPS observation period.

Satellite-derived surface velocity products 170
The GPS derived surface velocities are used to validate and assess the accuracy of several available ice velocity products derived from satellite data from the interior of the Greenland Ice Sheet. The ice velocity can be derived from space using data from synthetic aperture radar (SAR) or optical sensors. Optical feature tracking can provide velocities in very high resolution from coherent pairs of visual images. In the interior of the Greenland Ice Sheet, the surface is mostly featureless and SAR processing methods between pairs of SAR images are useful to derive surface velocities, either based on speckle tracking or 175 on phase displacements from interferometric synthetic aperture radar (InSAR).
We include several experimental ice velocity products in our assessment, as well as several one-year and multi-year products constructed from various remote sensing sources and methods. For each type of velocity product, we have also calculated a long-term average of all the velocity maps, and included this in our assessment. In total, we include 165 velocity products from the following sources: year v1). This map was derived from InSAR, SAR, and Landsat 8 optical imagery data using a combination of speckle-tracking, InSAR and optical feature tracking methods, supplemented with balance velocities near the ice divides where the flow speeds are <5 m a -1 . The data are provided with a resolution of 250 m. In the interior of the ice sheet, the estimated errors of this product are up to ~2 m a -1 , and reported to be <1 m a -1 in areas where InSAR is used. 185 2) NASA MEaSUREs Greenland Ice Sheet winter velocity maps (September-May) from InSAR data v2, 2000-2018 (Joughin et al., 2010(Joughin et al., , 2018) (MEaSUREs InSAR v2). These maps was derived entirely from data obtained by CSA RADARSAT-1, JAXA ALOS, and DLR TerraSAR-X/TanDEM-X (TSX/TDX) satellites, as well as from ESAs C-band SAR data from Copernicus Sentinel-1A/B. The maps were produced using an integrated set of SAR, speckle-tracking and interferometric algorithms (Joughin, 2002). The data are provided with a resolution of 200 m and the error is estimated to be <10 m a -1 . 7 3) NASA MEaSUREs Greenland Annual and Quaterly Ice Sheet velocity maps from SAR and Landsat v1, 2015-2018(Joughin et al., 2018 (MEaSURE SAR&Landsat v1). These maps are derived from SAR data obtained by DLR TerraSAR-X/TanDEM-X (TSX/TDX) and ESA Copernicus Sentinel-1A/B satellites, and from the USGS Landsat 8 optical imagery using a combination of speckle-tracking, InSAR and optical feature tracking methods (Joughin et al. 2018 (Nagler et al., 2015). 5) PROMICE Greenland velocity maps, 2016-2019 from SAR (Solgaard and Kusk, 2019; available from www.promice.dk).
These products are derived from ESA Copernicus Sentinel-1A/B SAR data using offset tracking  by 200 employing the operational interferometric post processing chain (IPP) (Kusk et al., 2018;Dall et al., 2015). Each map is a mosaic consisting of both 12-day and 6-day pairs within two Sentinel 1A cycles, and thus the temporal resolution of the product is 24 days. A new map is available every 12 days. The spatial resolution is 500 m, and the estimated error is 10-30 m a -1 . 6) DTU-Space experimental Sentinel-1A/B Greenland Ice Sheet velocity product, from InSAR (DTU-Space-S1). This product is derived from SAR data acquired by ESA Copernicus Sentinel-1A/B satellites in the period from 2019-01-01 to 2019-01-18 205 from two ascending and three descending tracks. Eight 6-day pairs and five 12-day pairs were processed using the in-housedeveloped Interferometric Post Processing chain (IPP) (Kusk et al., 2018). The spatial resolution is 50 m, and the estimated errors are <1 m a -1 (Andersen et al., 2020). 7) AWI experimental TerraSAR-X (TSX) Greenland velocity product, from InSAR (AWI-TSX). The velocity field was derived from SAR interferometry obtained by DLR TSX by combining data from ascending and descending satellite orbits 210 following well-established methods (e.g. Joughin et al., 1998). Three interferograms were formed from descending satellite data acquired between 2016-09-07 and 2016-10-01 and another three from ascending satellite data acquired between 2017-10-24 and 2018-01-03. All interferograms have a temporal baseline of 11 days with perpendicular baselines varying between 25 and 180 m. Due to the latter a certain topography induced phase difference is present in the interferograms, which was removed with the help of the global DLR TanDEM-X DEM with a 30 m grid resolution. The topography corrected interferograms were 215 unwrapped using GAMMA's minimum cost flow algorithm  and combined to 3D velocity maps assuming surface parallel ice flow. In order to set the relative velocity estimates to absolute values seed points were extracted from the MEaSUREs multi year v1 dataset, and adjacent velocity fields where patched together using the average value in their overlapping areas. The final product was gridded to 30 m spatial resolution. The AWI-TSX product is developed for this study. USGS Landsat 4, 5, 7, and 8 optical imagery using the auto-RIFT feature tracking processing chain described in Gardner et al. (2018). The final product was gridded to 120 m spatial resolution. The images suffers from x-and y-geolocation errors of 8 margin areas or to the median reference velocity in slow-moving areas. In interior Greenland, the MEaSUREs Greenland 225 Annual Ice Sheet Velocity Mosaic from SAR and Landsat version 1 velocity product is used as the reference velocity (Joughin et al., 2010). For the assessment here, we derived a multi-year velocity product from 1985 to 2018, averaged from the annual products. In our observed area, the data is mainly derived from Landsat 8, and we therefore also derived an additional 6-year average of the annual product from 2013 to 2018, covered by the Landsat 8 imagery. These two products were included in the assessment. 230

Comparison between GPS data and a satellite-derived digital elevation model
We compare GPS derived surface elevations with the surface elevation sampled from the ArcticDEM release 7 at the stake positions (WGS84 ellipsoidal heights) and find an agreement within 1m, except for one stake with a deviation of >1.5 m (Supplementary Figs. S2 and S3). Minor differences between the two datasets could be due to variable snow accumulation 235 through the years, leading to seasonal and interannual variability in surface elevation, which is captured differently by the two datasets due to their different timestamps. The outlier is located in the exceptionally deep and narrow trough in the southern shear margin (Fig 3a). Local topography effects at these stakes could possibly be due to interpolation or shadow effects in the Arctic DEM (Porter et al., 2018). The difference between the 63 GPS derived surface elevations and the ArcticDEM at the location of the GPS stakes is 0.48 m (mean) and 0.47 m (median) with a standard deviation of 0.53 m, confirming the low 240 uncertainty of the ArcticDEM (Noh and Howat, 2015;Porter et al., 2018).

Comparison between GPS data and surface velocity products derived from satellite data
The assessment consists of an inter-comparison between the GPS derived velocities of the stakes from the period 2015-2019 and the interpolated surface velocity at the location of the stakes from the satellite-derived velocity products. For each velocity product, we determine the accuracy (the bias) and precision (the standard deviation, i.e. the root mean square difference (RMS) 245 after removing the bias), between the GPS derived velocities and the satellite-derived velocities at the location of the stakes.
In addition to the direct inter-comparison between the GPS derived velocities and the satellite-derived velocities, we also investigate the variability of the satellite-derived velocity products. In order to do so, we perform a spatial smoothing of the satellite-derived velocity product with a running mean filter with a smoothing length, and we then vary the smoothing length in order to determine the optimum smoothing length that minimizes the standard deviation (RMS) between the GPS 250 observations and the velocity product. The results of the inter-comparison for the top ten products (sorted according to the standard deviation) are listed in Table 2 (with a complete overview of the results from all products in Supplementary Table   S2), and it is illustrated in Figure 6. 9 It is important to note that the different time stamps and temporal coverage ∆ of the observations are not taken into account in the inter-comparison. In satellite-derived velocity products with longer temporal coverage, possible temporal variability 255 and/or noise are smoothed, and there is a clear relationship between increasing temporal coverage and decreasing bias, standard deviation and optimum smoothing length. Similarly, spatial smoothing can remove noise. The improvement of the products with the temporal coverage ∆ is significant, with the bias decreasing approximately linearly with √∆ , as illustrated in Figure   6. Some long-term products were calculated as averages of short-term products, i.e. based on more observations, which would also help reduce the temporal variability and noise of these products compared to short-term products. The bias of all the 165 260 products are in the range ~0.3 -40 m a -1 , with a standard deviation in the range ~0.4 -22 m a -1 . The velocity products already include some smoothing as part of their production, but additional smoothing both temporally or spatially, for most products, reduced the standard deviation. After applying optimum spatial smoothing, the standard deviation is reduced to a range of ~0.4 -10 m a -1 . The optimum smoothing length is typically in the order of 500-3000m.
As part of the assessment, we use the whole set of satellite-derived surface velocity products to trace flow lines along NEGIS. 265 We use a starting point at the central reference stake near the EastGRIP site, which is located in the center of our observed area in a relatively narrow section of the NEGIS ice stream. We trace the flow lines upstream into the slower moving areas where flow converges into NEGIS and downstream into faster flow where the ice stream widens (Fig. 7). The flow lines are gradually displaced depending on their bias and fluctuate depending on their standard deviation.

Assessment of surface velocity products derived from satellite data
In the interior regions of the Greenland Ice Sheet, validation of satellite derived ice velocity and surface elevation products is generally limited due to lack of in-situ data. Our GPS stake network provides a unique dataset for validation in the interior accumulation area of the ice stream, and it represents a range of velocities and velocity gradients over one order of magnitude in the NEGIS ice stream, the shear margins and the surrounding slow moving areas. However, the assessment is restricted due 275 to the limited spatial extent of the GPS data, and our conclusions may not apply to margin areas with very fast flow, seasonal variability or high surface slopes.
In our comparison, the DTU-Space-S1 experimental product stands out among all the investigated products with its short temporal coverage (~10-20 days), low bias of 0.35 m a -1 and the low standard deviation of 0.55 m a -1 , which can be reduced to 0.53 m a -1 after optimum smoothing of 354 m. The AWI-TSX experimental product stands out because of its minimum 280 standard deviation of 0.39 m a -1 of all the investigated products, its high spatial resolution, which results in a very low optimum smoothing length of 10 m, i.e. no further smoothing is needed to reduce the noise, and the bias of the AWI-TSX product is 0.51 m a -1 , also among the lowest of the investigated products. Both these products are based entirely on InSAR processing methods. of 0.77 m a -1 and a standard deviation of 0.50 m a -1 . Since this product is already a 20-year average, the optimum smoothing length is only 200 m and only slightly reduces the standard deviation to 0.48 m a -1 . It is notable, that the bias of this product is similar to several other MEaSUREs products with shorter temporal coverage, while the standard deviation of this product is smaller than the other MEaSUREs products. If the interior of the ice sheet changes slowly over time, the differences between the temporal stamp of the GPS observations and of the multi-year velocity product covering 1995-2015 may become important. 290 However, MEaSUREs winter velocity map from 2008-2009, the MEaSUREs InSAR v2 product, performed very similar to the MEaSUREs multi-year v1 product, with a bias of 0.89 m a -1 , standard deviation of 0.46 m a -1 , and an optimum smoothing length of 51 m. The winter velocity product from 2008-2009 is based on InSAR and stands out with its low standard deviation and a relative short temporal coverage of nine months. The similar agreement between these products and the GPS derived velocities suggest that the velocity in the interior part of NEGIS has not changed significantly in the last decade. 295 The five products with a minimum bias are the MEaSUREs ITS_LIVE 6-year average product with a bias of 0.31 m a -1 , the MEaSUREs combined SAR&Landsat v1 1-year product for 2015 with a bias of 0.33 m a -1 , the DTU-Space-S1 18-day product from 2019 with a bias of 0.35 m a -1 , the ESA CCI annual velocity product from 2015-2016 with a bias of 0.43 m a -1 , and the 24-days PROMICE product from 2018-02 with a bias of 0.46 m a -1 . Of these, the DTU-Space-S1 product stands out as mentioned above. MEaSUREs ITS_LIVE product stands out with its long temporal coverage, low standard deviation, very 300 low optimal smoothing length, and because it is the only product in our study entirely based on optical feature tracking. The other products have a time stamp that overlap with the first one to two years of the GPS observation period, but their standard deviations are much higher due to the SAR processing techniques. The MEaSUREs combined SAR&Landsat v1 product has a standard deviation of 1.85 m a -1 , which reduces to 1.65 m a -1 after optimum smoothing over 1224 m. The ESA CCI product performs very similar with a standard deviation of 1.94 m a -1 , which reduces to 1.11 m a -1 after optimum smoothing length 305 over 1185 m. The PROMICE product with its very high temporal resolution of 24 days has a standard deviation of 5.39 m a -1 , which reduces to 2.6 m a -1 after optimum smoothing length over 2264 m. The optical product (MEaSUREs ITS_LIVE) can obtain a comparably high accuracy when averaged over long time intervals (several years), but the standard deviation is slightly higher than the radar based products. Mouginot et al. (2017) also derived 315 optical ice velocities from Landsat-8, and concluded similarly that the quality of the products derived from optical satellite sensors is comparable to data obtained with SAR sensors.

Inferred flowlines from satellite-derived products
Knowing the accurate flowlines of an ice sheet is useful for many applications, such as defining the outlines of drainage basins or to identify the source area for ice flowing through a specific survey site. For studies related to the internal stratigraphy and 320 ice properties, e.g. in ice cores or radar profiles, it is essential to know the upstream flow path in order to infer the deformation history of the internal layers. However, minor uncertainties and bias in satellite-derived velocity products can severely affect flow lines traced along the velocity field, as these uncertainties can displace the flowline and propagate along the flow line (Fig. 7). The flowlines for the products with minimum bias are marked (Fig. 7), showing that flowlines can only be reliably traced if the bias is small. These products have a low bias of 0.31 m a -1 to 0.43 m a -1 or approximately 1% of the surface speed. 325 This is particularly critical when the flow is strongly convergent or divergent. We notice that the back-trajectories diverge more than the forward trajectories, and we attribute this to the higher uncertainty of the upstream lower velocities compared to downstream. As a result, it may be better to use surface slopes instead of surface velocity products to trace flow trajectories in slow moving areas.

Estimated errors of satellite-derived strain rates 330
Strain rates are derived from the satellite-derived products as derivatives of the velocity fields and have therefore higher errors.
Our assessment provides an estimate of the strain rates error depending on the resolution from the standard deviation of the velocity product. For velocity products with a standard deviation in the order of 0.5 m a -1 , the strain rate uncertainty is in the order of 10 -3 a -1 on a 500 m grid, but it could be improved by smoothing the velocity product using the optimum smoothing length from the assessment. We compare the GPS derived strain rates with strain rates from MEaSUREs multi-year v1 velocity 335 product in transects across the NEGIS (Fig. 3) and along NEGIS (Fig. 4). We calculate the satellite-derived strain rate tensor directly from the gridded 250 m resolution velocity product without further smoothing according to the optimum smoothing length (included in Table 2), and rotate according to the local flow direction in order to determine the strain rates along the flow. While the GPS derived strain rates are limited in resolution and do not exactly capture the maximum strain rate at the southern shear margin, they do capture the enhanced strain rates in the shear margins. The fluctuations of the satellite-derived 340 strain rates are less than 10 -3 a -1 , thus confirming our estimated uncertainty above. The satellite-derived strain rates capture the high resolution strain rate peaks of approximately 3-4ꞏ10 -3 a -1 in the shear margins (Fig. 3) and approximate 2ꞏ10 -3 a -1 along NEGIS (Fig. 4).

Structure and flow of NEGIS
The assessment of satellite-derived velocity and height products inform of the accuracy and limitations of the satellite-derived 345 products in the interior regions of the Greenland Ice Sheet. In our study area the best products have a bias and precision of less than approximately 0.5 m a -1 , i.e. about 5% of the smallest observed GPS derived velocities of around 10 m a -1 in the slowmoving areas north of NEGIS. Knowing the limitations of the satellite-derived products, we are now able to combine the GPS derived velocities and strain rates with satellite-derived data to characterize spatial patterns in surface structure and ice flow in the interior part of the NEGIS ice stream. We discuss here the observed patterns. 350 The flow and surface topography across the NEGIS ice stream reveal a distinct 25 km wide fast-flowing ice stream near the EastGRIP site, which is sharply marked at both sides in speed, strain rates and surface geometry (Fig. 3). The cross sections show a central 10 km wide section with an almost uniform speed of 55 m a -1 , and well-defined shear margins at both sides with a width of about 5 km separating the ice stream from the surrounding slow-moving ice (Fig. 3). The velocities are above 20 m a -1 on the southern side where a broad flow field is merging with NEGIS, and approximately 10 m a -1 on the northern 355 side (Fig. 2 and 3). The strain rates are at a level of approximately 10 -4 a -1 within the ice stream and in the surrounding slow moving areas outside the ice stream. In the shear margins, they increase by an order of magnitude to a maximum value of approximately 10 -3 a -1 . The remarkably uniform velocities and low strain rates in the fast flowing central band of NEGIS with narrow shear zones at the margins with enhanced strain rates is characteristic of ice stream flow (e.g. Minchew et al., 2018).
In our study area in NEGIS, Holschuh et al. (2019) proposed that thermal softening of ice is present in the shear margins, 360 despite the relatively low strain rates. The surface topography reveals a 30-40 m deep lowering coinciding with the fast flow within NEGIS with well-defined deep troughs marking the shear margins. These deep shear margin troughs form due to a combination of enhanced longitudinal stretching and shear as the ice flow enters the fast flowing ice stream from both sides at an angle of ~15, accelerates and turns, and an enhanced firn densification in the shear margins due to the enhanced horizontal deformation (Riverman et al., 2019a). 365 The location of the shear margins cannot be clearly linked to the bedrock topography in the area (Christianson et al., 2014;Franke et al., 2020). Christianson et al. (2014) proposed that the shear margins of NEGIS are controlled by a self-stabilizing mechanism related to gradients in the subglacial hydropotential due to the surface troughs that restrict widening of the ice stream, and the internal stratigraphy suggests that the shear margins have been relatively stable during the Holocene (Keisling et al., 2014). Detailed maps of bedrock topography in the area reveal subglacial landforms proposed to be related to basal 370 erosion due to the fast flow (Franke et al., 2020), and elongated bedforms aligned with the flow (Franke et al., 2020;Riverman et al., 2019b). These elongated bedforms are seen in the transects across NEGIS as 100-300 m undulations in bedrock topography (Fig. 3), and they appear here to be related to the location of the shear margins. The southern very well-defined shear margin trough is consistently located above a local bedrock low in the three cross sectional profiles spanning a 5 km distance along NEGIS (Fig. 3). The northern broad shear margin trough is located over a wide bedrock valley, and the shear 375 margin trough narrows from a wide double-trough to a single trough in the three cross sectional profiles over a 5 km distance along NEGIS, as the bedrock valley over the same distance narrows (Fig. 3). Thus, our observations support that these bedforms are related to the shear margins (Franke et al., 2020;Riverman et al., 2019b), but further studies are needed to fully understand the conditions at the shear margins.
The ArcticDEM surface topography of the NEGIS ice stream shows that an organized spatial pattern of wavy undulations 380 develops perpendicular to the NEGIS flow in the area around EastGRIP (Fig. 8). The undulations develop within the fast flowing central flow band of NEGIS in a 25 km section along NEGIS where the surface velocity remain at a level of 13 approximately 55 m a -1 . Upstream from this section, the flow accelerates over tens of km with an acceleration of approximately 10 m a -1 over 10 km, i.e. longitudinal strain rates of ~10 -3 a -1 . The undulating patterns start forming as ice velocity exceeds a threshold velocity of approximately 55 m a -1 , and as the ice flows over a 200-300 m bedrock transition to a bedrock plateau of 385 approximately 200 m elevation and widens (Franke et al., 2020), suggesting that the undulations are related to the bedrock topography (Fig. 4). The undulations in surface slope are connected to undulations in the longitudinal and transverse surface strain rates, and to some degree related to undulations in bedrock topography (Fig. 4). Similar organized undulating patterns in driving stress was previously reported in Antarctic and Greenland ice sheets in fast flowing areas (Sergienko and Hindmarsh, 2013;Sergienko et al., 2014) and related to patterns in basal stress located in areas with significant sliding. These previous 390 studies attributed the patterns to be the result of instabilities related to subglacial water beneath a sliding glacier, and our results support that bedrock topography plays a role in relation to these undulations.

Conclusion
We We compare our GPS derived heights and surface velocities with the ArcticDEM height model (Porter et al., 2018), as well as 400 published and experimental remote sensing velocity products in order to validate and assess these products. We include surface velocity products from the MEaSUREs program, the ESA CCI program, the PROMICE program, and experimental data products from MEaSUREs, DTU-Space and AWI. For each product, we calculate the bias, the standard deviation relative to the GPS derived surface velocities, and the spatial smoothing that minimizes the standard deviation. Our assessments show: -Among the top five surface velocity products with lowest standard deviation compared to the GPS derived surface velocities, three are entirely based on InSAR (DTU-Space, 2019;AWI-TSX, 2016MEaSUREs winter velocity by InSAR v2, 2008-2009, and two are combined products averaged over a multi-year period (MEaSUREs multi-410 year product, 1995(MEaSUREs multi-410 year product, -2015MEaSUREs multi-year SAR andLandsat, 2014-2018).
-SAR based surface velocity products from ESA CCI, PROMICE and MEaSUREs can obtain comparable precision compared to the GPS derived surface velocities if they are averaged over longer time periods (years) and smoothed spatially, and they generally obtain a low bias.
products if it is averaged over long periods (several years), but the bias is slightly higher.
Overall, the assessments show that for interior velocity estimates, the InSAR based products stand out with higher resolution in time and space and low errors. For all products, longer observation time improves the products in these interior areas where surface velocity has not changed significantly over the last decade.
The assessment inform of the accuracy of the satellite-derived velocities, and thereby allow us to evaluate the use of these 420 products for investigations of flow patterns in the interior regions of the Greenland Ice Sheet. We show that satellite-derived strain rates can capture high-resolution spatial signals at the shear margins and within the fast flowing part of NEGIS, despite the high uncertainty in the order of 10 -3 a -1 . We show further that the strain rate peaks along NEGIS are part of a regular undulating pattern forming in surface slope and strain rates when the surface velocity exceeds approximately 55 m a -1 , and we argue that the formation of these undulations appears to be related to bedrock topography. 425 We derived flowlines from the satellite-derived velocity products and showed that even a minor bias in these products can severely affect the path of the flowlines, in particular in slow-moving areas. We conclude that reliable flowlines can only be derived from satellite-derived velocities with a low bias compared to the surface speed, and that surface slopes may produce more realistic flowlines than satellite-derived velocities in slow-moving areas.
The study demonstrates that it is important to know the limitations of the satellite-derived products. We conclude that available 430 satellite-derived products are sufficiently accurate to allow a detailed analysis of the ice flow in the interior part of NEGIS, which can contribute to understand the flow near its onset in interior North Greenland, and ultimately to improve projections of its future response to mass loss at the margins.
Data availability: The GPS derived positions and surface velocities presented in this paper are included in the supplementary 435 Table S1. The assessment results for the 165 velocity products are included in the supplementary Table S2.   Table 2. The assessment results for the ten velocity products with the smallest standard deviation (RMS). The products are sorted with increasing RMS. Notice that the AWI-TSX, the DTU-Space-S1, and the MEaSUREs ITS_LIVE products are also among the ten products with the smallest bias. The velocity bias is determined for both the x-and y-direction (see supplementary table S2). The bias listed here is the length of the velocity bias vector, i.e. the average rate of change in distance between poles moving with the satellite based velocity field compared to poles moving with the GPS velocities. The complete list of assessment results can be found in Supplementary Table S2.    2b). For each cross section, the three plots are: Left: Surface velocity (blue) and surface elevation (red). Middle: surface strain rates relative to the local flow direction, along flow (blue full) and transverse to flow (blue dotted) with positive for stretching and negative for compression, and surface elevation (red). Right: Surface elevation (red) and bedrock topography (blue). GPS observations are shown as black circles/squares at three stakes marked in Fig. 5. Surface elevation is from the ArcticDEM (Porter et al., 2018), surface velocity and strain rate profiles are derived from MEaSUREs multiyear Greenland Ice sheet velocity v1 product (Joughin et al., 2016(Joughin et al., , 2017, and the bedrock topography is from BedMachine v3 (Morlighem et al., 2017a(Morlighem et al., , 2017b. The vertical gray lines in b) panels indicate the position of the central stake near the EastGRIP site.  (Porter et al., 2018), MEaSUREs multi-year Greenland Ice Sheet Velocity v1 product (Joughin et al., 2016(Joughin et al., , 2017, and BedMachine v3 (Morlighem et al., 2017a(Morlighem et al., , 2017b. GPS derived surface elevations, velocities and strain rates are shown in red. The vertical gray line indicates the position of the central stake near the EastGRIP site. Figure 5. Horizontal principal strain rates for 32 different triangular sections within the GPS array. The principal strain rates are plotted at the centroids of each triangle (black circles) with black lines indicating positive strain rates (extension) and green lines indicating negative strain rates (compression). The GPS-derived surface velocities (red arrows) are plotted at the stakes (blue dots). At three GPS stakes (yellow dots), strain rates along the direction of the flow are calculated. Figure 6. Results of the assessment of a list of 165 satellite-derived surface velocity products: a) The mean bias of the velocity product compared to the GPS derived velocities at the location of the GPS poles, b) the standard deviation of the velocity products relative to GPS derived velocities (RMS), and c) the optimal smoothing length of the velocity product that minimizes the standard deviation ( ). All results are shown as a function of ∆ , the timespan of the velocity product. Notice that some products have been averaged over time to provide results with longer temporal coverage. The gray lines suggest a linear dependency of the bias, RMS and on temporal coverage ∆ .   (Porter et al., 2018), c) longitudinal strain rate in a -1 along the direction of the flow, and d) transverse strain rate in a -1 along the transverse direction to the flow, both from MEaSUREs multi-year Greenland Ice Sheet Velocity v1 product (Joughin et al., 2016(Joughin et al., , 2017.