High-resolution ice thickness and bed topography of a land-terminating section of the Greenland Ice Sheet

. We present ice thickness and bed topography maps with a high spatial resolution (250–500 m) of a land-terminating section of the Greenland Ice Sheet derived from ground-based and airborne radar surveys. The data have a total area of ∼ 12 000 km 2 and cover the whole ablation area of the outlet glaciers of Isunnguata Sermia, Russell, Leverett, Ørkendalen and Isorlersuup up to the long-term mass balance equilibrium line altitude at ∼ 1600 m above sea level. The bed topography shows highly variable subglacial trough systems, and the trough of Isunnguata Sermia Glacier is overdeepened and reaches an elevation of ∼ 500 m below sea level. The ice surface is smooth and only reﬂects the bedrock topography in a subtle way, resulting in a highly variable ice thickness. The southern part of our study area consists of higher bed elevations compared to the northern part. The compiled data sets of ground-based and airborne radar surveys cover one of the most studied regions of the Greenland Ice Sheet and can be valuable for detailed studies of ice sheet dynamics and hydrology. The combined data set is freely available at doi:10.1594/pangaea.830314.


Introduction
The first radar measurements on the Greenland Ice Sheet were collected in the 1960s (Evans, 1963;Waite and Schmidt, 1962). Since then, various campaigns have measured the elevation of the ice-covered bedrock (Bogorodsky et al., 1985;Evans and Robin, 1966;Letreguilly et al., 1991). The first compilation of bed elevation data over the whole Greenland Ice Sheet by Bamber et al. (2001) consisted of 5 km gridded maps of ice thickness and bed topography. Numerous surveys have increased the data density or filled in the gaps in the data in this grid. An updated bed map, with six different data sources, was recently published as a 1 km grid (Bamber et al., 2013a).
The size of many outlet glaciers in Greenland is small and bed elevation data sets with a higher resolution than is currently available are required for modelling ice sheet dynamics. On a regional scale, high-resolution digital elevation models (DEMs) of the bed also allow subglacial hydrological pathways and drainage basins to be determined with confidence (e.g. Wingham et al., 2006;Wright et al.,2008) and subglacial landforms and landscapes to be studied in detail (e.g. Bamber et al., 2013b;Bingham and Siegert, 2009;King et al., 2009).
Recent high-resolution measurements of ice thickness have focused on mapping the fast-flowing marineterminating glaciers (e.g. Plummer et al., 2008;Raney, 2009) that drain the majority of the Greenland Ice Sheet, while the Published by Copernicus Publications.
typically slower, land-terminating glaciers have received less attention. Land-terminating glaciers and their catchments, however, provide ideal study areas for investigating the response of ice sheet dynamics to atmospheric forcing, as they are isolated from marine influences such as calving and submarine melt. Furthermore, under a warming climate tidewater outlet glaciers are expected to retreat inland, causing a larger portion of the ice sheet to be land-terminating in the future (Sole et al., 2008). A higher-resolution map of ice thickness and bed topography (< 1 km grid) of a land-terminating region of the Greenland Ice Sheet is therefore timely.
Here, we present ice thickness and bed topography DEMs of a land-terminating section of the Greenland Ice Sheet, based on a compilation of ground-based and airborne radar surveys. The DEMs have a total area of ∼ 12 000 km 2 at a resolution of 250-500 m. Our combined data set of groundbased and airborne radar surveys is available for integration into databases of bed elevation on Greenland (e.g. Bamber et al., 2013a) and can be valuable for detailed studies of ice sheet dynamics and hydrology. Our collected data will be used in a project that aims to improve the current understanding of hydrogeological processes associated with continental-scale glaciations, including the presence of permafrost and the advance/retreat of ice sheets.

Study area
The study area is located in West Greenland and includes the informally named Isunnguata Sermia, Russell, Leverett, Ørkendalen, and Isorlersuup glaciers and their catchment areas. The radar survey extends a farther 100 km south of the Isorlersuup Glacier and 90 km inland to approximately the 21year mean mass balance equilibrium line altitude (ELA) at ∼ 1600 m a.s.l. (van de Wal et al., 2012). The glaciated area is one of the most studied regions of the Greenland Ice Sheet with studies of mass balance (e.g. van de Wal et al., 2012), dynamics (e.g. van de Wal et al., 2008;Bartholomew et al., 2011;Palmer et al., 2011;Sole et al., 2013), and supraglacial lakes (Doyle et al., 2013;Fitzpatrick et al., 2014). Recently, two studies have published DEMs of the Isunnguata Sermia and Russell glaciers (Jezek et al., 2013;Morlighem et al., 2013), based on the IceBridge data set (Leuschen and Allen, 2010), also used in this study (see Sect. 2.2). These studies cover the northern part of our study area, to an extent of ∼ 25 % of our maps. In comparison, our data cover the whole ablation area and contain, in addition to the IceBridge data set, two previously unpublished data sets.

Data and methods
The data set in this study was compiled from three different sources (Fig. 1). We collected ground-based radar surveys during spring 2010 and spring 2011 which we combined with two airborne radar data sets, collected by the Technical University of Denmark in 2003 (previously unpublished)  and by the NASA IceBridge project between 2010 and 2012 (Leuschen and Allen, 2010). In the following sections, we describe the methods used to acquire, assimilate and interpolate each data set into the final product. The system parameters of each data set are summarized in Table 1.

Ground-based radar surveys
During spring (April-May) in 2010 and 2011 some 1500 km of common-offset radar profiles were collected with two ground-based impulse radar systems, hereafter referred to as the UU (Uppsala University) data set. Each radar system consisted of resistively loaded half-wavelength dipole antennas of 2.5 MHz centre frequency. An impulse transmitter was used with an average output power of 35 W and a pulse repetition frequency of 1 kHz. The 16 bit receiver had a capacity of collecting ∼ 1000 traces per second (Table 1). The trace acquisition was triggered by the direct wave between transmitter and receiver. The radar systems Earth Syst. Sci. Data, 6, 331-338, 2014 www.earth-syst-sci-data.net/6/331/2014/ were towed behind snowmobiles at a speed of 5-20 km h −1 , along tracks separated by 2 km. By stacking 3000 traces, a mean trace spacing of 15 m was achieved. Traces were positioned using data from a dual-frequency Global Positioning System (GPS) receiver mounted on the radar receiver sled, 90 m from the common midpoint along the travelled trajectory. The GPS data were processed kinematically using the Canadian Spatial Reference System (CSRS-PPP; Natural Resources Canada, 2013) precise point positioning service (Natural Resources Canada, 2013), which has an estimated theoretical uncertainty of ±0.02 m in the horizontal and ±0.03 m in the vertical. In practice, however, an error of ±1 m in the surface (horizontal and vertical) is expected due to the placement of the GPS antenna relative to the common midpoint of the radar. Several corrections and filters were applied to the radar data: (1) a butterworth bandpass filter, with cut-off frequencies of 0.75 and 7 MHz, was used to remove the unwanted frequency components in the data; (2) normal move-out correction was applied to correct for antenna separation; (3) rubber-band correction was used to interpolate the data to uniform trace spacing; and (4) two-dimensional (2-D) frequency wave-number migration (Stolt, 1978) was used to collapse hyperbolic reflectors back to their original positions in the profile direction. The bed returns were digitized semi-automatically with a cross-correlation picker (Irving et al., 2007). We calculated the ice thickness from the picked travel times of the bed return using a constant wave speed of 168 m µs −1 (discussed further in Sect. 2.3.1). Figure 2 shows an example image of a processed radar profile.

Airborne radar surveys
In addition to the ground-based data, we used two other data sets of subglacial topography provided by the Technical University of Denmark and the NASA IceBridge project, hereafter referred to as the DTU and IceBridge data sets. There are also data in the area collected by the Center for Remote Sensing of Ice Sheets (CReSIS) between 1993 and 2009 ), but we decided not to include these data as they had a lower resolution and provided no significant extended coverage.
The DTU data set in our study area consists of ∼ 3000 km of airborne radar profiles collected using a 60 MHz pulse radar during 2003. Data acquisition took place along flight tracks separated by ∼ 2.5 km, with a sampling rate of 3.125 Hz giving an along-track sample spacing of ∼ 25 m after processing (Christensen et al., 2000). Laser altimetry measurements of ice surface elevation allowed ice thickness to be determined (Forsberg et al., 2001). A southern subset of the DTU data set was published by Ahlstrøm et al. (2005).
The IceBridge data set in our study area consists of in the densely surveyed northern region (Fig. 1). Using a combination of echograms, the ice surface and bed layers were identified and the ice thickness was calculated by subtracting the bed layer from the surface layer. A constant wave speed of 169 m µs −1 was used. The layer tracking of the reflectors was done manually with basic tools for partial automation (e.g. automatic peak detector).

Radar system errors and uncertainty
The data sets were collected using different radar systems with varying specifications and, as a result, have varying vertical and horizontal resolutions (Table 1).

Vertical resolution
The range resolution is the accuracy of the measurement of distance between the antenna and the bed and can be determined from the characteristics of the source pulse (i.e. bandwidth) and the digitization frequency. For the UU data set the range resolution was estimated at 18.8 m. The DTU data set has a range resolution of ∼ 21 m, but the resolution can be significantly better if the signal-to-noise ratio is large (Christensen et al., 2000). The IceBridge data set has an alongtrack range resolution of 4.5 m (Leuschen and Allen, 2010). The use of a constant wave speed (168 and 169 m µs −1 in our case) for the ground-based and airborne surveys is a common method within glaciology, since glacial ice can be assumed to be a homogenous medium (e.g. Lythe et al., 2001). The wave speed, however, can vary spatially depending predominantly on density, for example with the presence of a firn layer and partly on the presence of liquid water in the ice. The profiles were collected in the ablation zone with thin seasonal snow cover; thus, we assume small density variations due to impurities in the ice (e.g. air bubbles). A typical variation of 2 % of glacier ice density (Navarro and Eisen, 2009) gives an uncertainty of ±20 m on the depth conversion (with an ice depth of 1000 m). Variations in the velocity of the radar signal can occur due to varying ice temperature and the presence of inhomogeneities and liquid water in the ice (Drewry, 1975). The effect of inhomogeneities and ice temperature is expected to have a small impact on the average velocity for the whole ice column while a significant water content in the ice can influence the velocity in a substantial way. In most parts of the study area, however, only the basal layer (< 10 % of the total ice thickness) is temperate and therefore contains liquid water, which would give an underestimation of ∼ 0.5 % of the average velocity at 2 % water content (Looyenga, 1965). Errors can also arise when processing the radar signal, since the bed was identified semi-automatically as a reflector on the radar image (Fig. 2).
To estimate picking and positioning errors within each data set and to test the consistency between the data sets we compared the differences (misfits) in the ice thickness and surface elevation estimates between different profiles and data sets at crossover points. The UU data set had a mean crossover misfit in ice thickness of 16.0 m with a standard deviation (σ ) of 20.3 m (based on 159 crossing points). The mean crossover misfit in ice surface elevation was 0.6 (σ = 0.9 m). The DTU data set had a mean crossover misfit in ice thickness of 12.0 m with a standard deviation of 13.0 m (based on 97 crossing points). The mean crossover misfit in ice surface elevation was 8.7 m (σ = 11.3 m). The IceBridge data set had a mean crossover misfit in ice thickness of 18.9 m with a standard deviation of 26.8 m (based on 1909 crossing points). The mean crossover misfit in surface elevation was 11.3 m (σ = 15.9 m). As the crossover analysis within the same data set does not capture systematic errors between the different data sets, and since 10 years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) separate the data acquisition, a comparison between the data sets is essential. When we ran a crossover analysis between all three data sets there was a mean misfit in ice thickness of 19.7 m (σ = 24.6 m). The mean crossover misfit in surface elevation of all the data sets was 8.1 m (σ = 7.5 m). With the inherent accuracy of the radar signal (system specifications and 2 % difference in wave speed), we estimated the total root-mean-squared uncertainty of the ice thickness to be 18.3 m for the UU data set, 18.1 m for the DTU data set, and 16.1 m for the IceBridge data set.

Horizontal resolution
Since the ice is thick, echoes from a large area at the bed will be integrated into the signal both along-track and acrosstrack. Theoretically, the reflected energy from a single reflector does not arrive at the receiver from a single point, but from an ellipsoidal zone of a horizontal plane, called the first Fresnel zone (Robin et al., 1969). The horizontal resolution is also dependent on the vertical variation (roughness) of the bed within the footprint and the errors can be large since the topography and acquisition geometry is not ideal. Given an ice thickness of 1000 m, the first Fresnel zone is 183 m for the UU data set. This is only a theoretical resolution and is improved by the applied 2-D migration along the profiles which collapses hyperbolic reflectors back to their original position. The post-migration resolution of a single reflector is λ/2, where λ is the wavelength of the radar signal at the centre frequency (Welch et al., 1998). Since λ = v/f, where v is the wave speed and f is the frequency, the theoretical best resolution is 34 m for the UU data set (f = 2.5 MHz). This is, however, only valid in the travel direction and does not account for diffraction orthogonal to the profiles. Therefore, large errors in the thickness measurements can still be expected in areas with steep slopes perpendicular to the direction of the profiles. The ideal processing technique in such cases is 3-D migration; however, this could not be applied since this requires a shorter distance between the profiles to avoid aliasing. The DTU data set has not been migrated and has an estimated along-track and across-track resolution of 81 m (Christensen et al., 2000). The IceBridge data set has an along-track resolution of 25 m and across-track resolution of 323 m for smooth terrain and 651 m for rough terrain, where the ice thickness is 2000 m and height above the air-ice interface is 500 m (Leuschen and Allen, 2010).

Assimilation of the data sets
We combined the ground-based and airborne data sets to produce DEMs of ice thickness and bed topography. The measuring interval for the data sets are dense along the profiles, with a data point spacing of 15-25 m, compared with 500-2500 m spacing between individual profiles. As this nonuniform spacing is not optimal for gridding algorithms, we subgridded the ground-based and airborne data sets into a 100 m pseudo-grid to reduce the data density along individual profiles. The subgrid was produced by calculating the median values for the points that fell within the distance of half the grid cell. Before we interpolated the ice thickness data, we added zero ice thickness along the edge of the ice sheet, which we derived from SPOT-5 satellite images acquired in August 2008. We interpolated the subgridded profiles to 250 m resolution in the northern part of the study area (above 66.7 • N), where there was high spatial density of profiles, and to 500 m resolution in the southern part (below 66.7 • N), where the spacing between profiles was the largest. We used a universal kriging interpolation algorithm (e.g. Isaaks and Srivastava, 1989) with a linear drift applied to remove large-scale trends. The interpolation model was based on an anisotropic variogram model, with linear and exponential models representing the spatial variability of the data set. We applied a smoothing nugget effect of 20 m to account for the accuracy of the data points.
We calculated the bed elevation by subtracting the ice thickness from the surface elevation in every grid point. For surface elevation we used the Greenland Ice Sheet Mapping Earth Syst. Sci. Data, 6, 331-338, 2014 www.earth-syst-sci-data.net/6/331/2014/ Project (GIMP) surface elevation model (WGS-1984 ellipsoid;Howat et al., 2014). The GIMP surface elevation model is constructed from a combination of ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) and SPOT-5 DEMs for the peripheral areas of the ice sheet, with a horizontal resolution of 30 m. The root-mean-squared validation error, relative to ICESat (Ice, Cloud, and land Elevation Satellite), is ±10 m for the GIMP data set. We constructed the bed topography data set in this way, instead of using the surface elevation data collected during the radar surveys, to give the compiled data set the same surface reference. To assure that this is a valid step, since the data sets were collected in different years (between 2003 and 2012), we subtracted the GIMP surface elevation at every data-set surface elevation point. The mean difference between the surfaces was calculated to be −5 m with a standard deviation of 10 m.
To assess the error in interpolation we cross-validated the gridded data; this is a common validation technique to see how well a model fits the observed data. By removing one observation from the data set, the remaining data were used to interpolate a value for the removed observation. This process was continued for 1000 random observations in the data; the error is the residual between the observed and the interpolated value (Isaaks and Srivastava, 1989). The standard deviation of the residuals was estimated to be 17 m, and increases with distance from the profiles (Fig. 3a). The data are missing in some sections of the tracks, primarily close to the ice margin (∼ 0 to 20 km; Fig. 3a), where it was difficult to receive a bed signal through fractured ice.
To summarize, the total accuracy of the estimated ice thickness and bed topography depends on: (1) the technical and theoretical capability of the radar systems; (2) uncertainty in the depth conversion; and (3) picking, positioning, and interpolation errors. By assessing all these potential sources of error, we estimate the maximum vertical rootmean-squared uncertainty in the final interpolated DEMs to be approximately ±20 m.

Results
We present DEMs of ice thickness and bed topography of a land-terminating section of the Greenland Ice Sheet at a 250-500 m spatial resolution (Fig. 3c, d). Due to a smooth ice surface and an undulating bed the ice thickness is highly variable (Fig. 3c), with a maximum gridded ice depth of 1470 m and a mean value of 830 m. Consistent with previous studies (e.g. Budd and Carter, 1971;Gudmundsson, 2003), the surface has a secondary component that is a subtle expression of the basal topography. The ice flow direction in the area generally runs from east to west, with a mean surface velocity of ∼ 150 m yr −1 (Fig. 3b; Joughin et al., 2010).
The highly variable subglacial topography (Fig. 3d) resembles the landscape in front of the ice sheet. The bed topography becomes smoother away from the ice margin, consistent with the larger, but lower resolution, bed map of Greenland (Bamber et al., 2013a). The deepest trough lies under the Isunnguata Sermia Glacier and has a minimum gridded elevation of −510 m below the World Geodetic System (WGS)-1984 ellipsoid (−540 m in the raw data) with a maximum difference of ∼ 1000 m between the valley bottom and adjacent subglacial hills. The southern part of the data set (south of 66.7 • N) consists of an area with higher bed elevations and includes the highest subglacial peak of 1060 m above the WGS-1984 ellipsoid (1100 m in the raw data). In comparison with the Bamber et al. (2013a) thickness DEM, our data show a standard deviation difference in thickness of 100 m. The relatively large difference between the two DEMs is possibly caused by our higher data density and higher gridded spatial resolution which creates greater detail in the DEMs. In addition, our interpolation model has been optimized for the region compared to the whole Greenland Ice Sheet. This highlights the need for regionally optimized DEMs when conducting detailed studies of the Greenland Ice Sheet.

Summary and outlook
We have compiled ground-based and airborne radar surveys from various sources to produce ice thickness and bed topography DEMs with high spatial resolution (250-500 m) of a land-terminating section of the western Greenland Ice Sheet. The DEMs cover the whole ablation area (12 000 km 2 ) up to the long-term ELA at ∼ 1600 m a.s.l., ∼ 90 km inland. The bed topography shows highly variable subglacial trough systems, resembling the landscape in the proglacial area. The ice surface is smooth and only reflects the bedrock topography in a subtle way, resulting in a highly variable ice thickness. The southern part of our study area consists of higher bed elevations compared to the northern part.
Further improvement of our maps could include filling in the gaps in the ice thickness and bed topography maps by surveying parts of the ice margin that are highly crevassed. This would demand a low-frequency airborne radar system designed for warm and fractured ice. Surveying the area with 3-D radar tomography (e.g. Paden et al., 2010) would also increase the spatial resolution substantially. Our maps, nevertheless, contain enough detail for a wide range of studies and can contribute to improvements in future ice sheet modelling efforts and studies of ice sheet dynamics and hydrology. The compiled data sets of ground-based and airborne radar surveys are freely available at doi:10.1594/pangaea.830314. The combined data set will be updated when the quality of the data is improved or if new data sets become available.