A new bed elevation dataset for Greenland

. We present a new bed elevation dataset for Greenland derived from a combination of multiple airborne ice thickness surveys undertaken between the 1970s and 2012. Around 420 000 line kilometres of airborne data were used, with roughly 70 % of this having been collected since the year 2000, when the last comprehensive compilation was undertaken. The airborne data were combined with satellite-derived elevations for non-glaciated terrain to produce a con-sistent bed digital elevation model (DEM) over the entire island including across the glaciated–ice free boundary. The DEM was extended to the continental margin with the aid of bathymetric data, primarily from a compilation for the Arctic. Ice thickness was determined where an ice shelf exists from a combination of surface elevation and radar soundings. The across-track spacing between ﬂight lines warranted interpolation at 1 km postings for signiﬁcant sectors of the ice sheet. Grids of ice surface elevation, error estimates for the DEM, ice thickness and data sampling density were also produced alongside a mask of land/ocean/grounded ice/ﬂoating ice. Errors in bed elevation range from a minimum of ± 10 m to about ± 300 m, as a function of distance from an observation and local topographic variability. A comparison with the compilation published in 2001 highlights the improvement in resolution afforded by the new datasets, particularly along the ice sheet margin, where ice velocity is highest


Introduction
The bed elevation and ice thickness of the Greenland ice sheet are important boundary conditions for numerical modelling. Surface ice velocity is roughly proportional to the fourth power of ice thickness (Paterson, 1994) and errors in the latter can, therefore, introduce substantial errors in modelled velocities for the present-day or future evolution of the ice sheet. Bed and surface geometry can be used to determine hydraulic potential and, hence, subglacial hydrological pathways (e.g. Wright et al., 2008), while elucidating subglacial topography can also provide insights on the origin and genesis of landforms at the bed (e.g. Young et al., 2011). For these and other reasons, a large number of airborne field campaigns have been deployed over Greenland during the last decade with a key objective of obtaining ice thickness measurements. The last major compilation of these data, for deriving bed elevations, was undertaken more than a decade ago (Bamber et al., 2001b) and did not, therefore, include the more recent and extensive field campaigns. In particular, in recent years there has been a focus on acquiring data over the fast-flowing outlet glaciers that fringe the margins of the ice sheet and which are responsible, in part, for the recent acceleration in mass loss observed on the Greenland Ice Sheet Howat et al., 2011;Howat and Eddy, 2011). It is now well established that some of these marine-terminating outlet glaciers can respond rapidly, and with large amplitude, to changes in the force balance at the bed, lateral margins and/or calving front and, as a consequence, attention has been focussed on modelling their past (Nick et al., 2009) and future behaviour (Nick et al., 2012). For such applications, accurate basal geometry is critical for determining rates of grounding-line migration and potential pinning points (Favier et al., 2012). Here, we take an approach that is aimed at maximising the resolution and utility of the basal topography in these key, marginal sectors of the ice sheet: areas that are often challenging for conventional ice penetrating radar (IPR) systems. As a consequence, we employ a range of techniques for improving the bed representation in these areas and interpolate the data at two different resolutions, which are then merged into a single product. This product is intended to be dynamic such that, as new data become available, they will be incorporated in new releases. In addition, we intend to update the datasets with new methods, such as the use of a mass conservation model, that improve the interpolation (Morlighem et al., 2011). This paper describes the methods used, and presents the first release of the products.

Airborne datasets
In contrast to the previous compilation (Bamber et al., 2001b), bed elevation is the interpolated parameter rather than ice thickness. The latter is derived from the difference between the bed and surface elevation. This means that the bed elevation varies smoothly and realistically across the land-glaciated boundary. To determine the bed elevation, data from a number of airborne IPR missions have been collated, alongside new datasets for the unglaciated bedrock. At present, we have collated ice thickness data from seven sources, each of which is described below and detailed in Table 1. Figure 1 shows the spatial distribution of the different data sources and is provided at higher resolution in the Supplement.

Center for Remote Sensing of Ice Sheets (CReSIS)
The majority of the data included in this compilation are from a series of instruments developed and flown by CReSIS at the University of Kansas (Gogineni et al., 2001). These data were collected between 1993 and 2009, with those from 1993-1999 being identical to the data used previously (Bamber et al., 2001b). Between 1993 and2002, data (Bamber et al., 2001b). CReSIS00 (blue) includes all data derived from CReSIS instruments between 2000 and 2012.
2009, the Multi-Channel Radar Depth Sounder (MCRDS), and since 2010, the Multi-channel Coherent Radar Depth Sounder (MCoRDS) were flown. Since 1999, around 65 % more data were available compared with the previous compilation and from 2006, more effort has been spent on focused campaigns with dense grids over individual outlet glaciers.
The MCoRDS instrument operates over the 180-210 MHz frequency range with a 10-30 MHz adjustable bandwidth and multiple receivers developed for airborne sounding and imaging of ice sheets. Multiple receivers allow cross-track surface clutter to be suppressed so that relatively weak bed echoes can be retrieved. MCoRDS has been flown on the NASA P-3 and DC-8 aeroplanes. Aircraft navigation was  1996, 2004, 2010Nixdorf et al. (1999  from kinematic GPS and a precision laser altimeter was usually mounted coincidentally. Based on a comparison with ice core sites, the vertical accuracy of the thickness measurement was estimated to be ±10 m (Gogineni et al., 2001), but it is worse in areas where the ice/bed interface is ambiguous or complex (e.g. in hilly terrain or when off-nadir scattering obscures the desired basal return) and larger cross-over errors are commonly seen in these regions. Data processed to contain location and ice thickness are available from CReSIS and were used in the main here. Geolocated radar echo strength profile images (often known as radargrams) were used only to verify the removal of some data which appeared to be unphysical.
In a region close to the outlet of Jakobshavn Isbrae, CRe-SIS produced a 125 m posting bed DEM, which includes all data collected in the region by them between 1997 and 2007. This grid was used in place of the original CReSIS data as they reprocessed all the data in the main channel and collocated the data with coincidental lidar surface elevations. Additionally, they included ASTER data in bare rock areas and fjord soundings to complement the airborne data.
Since 2010, the MCoRDS instrument was flown as part of the NASA Operation IceBridge (OIB) programme (Studinger et al., 2010). OIB is designed, primarily, to provide airborne data to fill the gap between the end of the ICESat record in 2009 and the launch of ICESat-2, which is scheduled for 2016. MCoRDS operates on all flights where ice thickness measurements can be retrieved. Several of the existing dense grids over rapidly changing outlet glaciers were re-flown as well as tracks over previously unsurveyed glaciers. These data significantly improve coverage particularly where narrow, fast flowing outlet glaciers were previously unsurveyed around the north-west and south-west coasts of Greenland. The gridded datasets include flights from the 2012 OIB season, which increases the total coverage by some 30 % (Fig. S2 in the Supplement).

Alfred Wegener Institute (AWI)
The AWI airborne instrument is capable of penetrating 4 km of ice with better than 50 m vertical accuracy and 3.25 m along-track sampling . It operated at 150 MHz transmitting bursts of 60 ns and 600 ns duration. In earlier years a combination of GPS and inertial navigation was used and, since 1997, differential GPS was employed. Horizontal navigation errors are up to ±100 m. Data were collected in [1996][1997][1998][1999]2004 and 2010, operating out of the NGRIP camp site in central Greenland (Nixdorf and Goktas, 2001) and from coastal airstrips at Station Nord  and Qaanaaq. These data provide dense coverage in central Northern Greenland (Fig. 1).

Warm Ice Sounding Explorer (WISE)
WISE is an airborne sounder designed to measure ice thickness in areas of warm and fractured ice. It is based on the MARSIS planetary sounder used on Mars. It operates at a 2.5 MHz centre frequency with a monopole antenna with navigation using conventional GPS. It was operated by the NASA Jet Propulsion Laboratory (JPL) on an Air Greenland Twin Otter as an IPY deployment. Data were collected in 2008, 2009 and 2010 in marginal areas. WISE provides useful coverage for some marginal areas of southwest Greenland.

Pathfinder Advanced Radar Ice Sounder (PARIS)
PARIS flew in 2009 and was operated by Johns Hopkins University on an Operation Ice Bridge flight (Raney, 2010). PARIS successfully demonstrated high altitude soundings with a delay Doppler design. The along track sampling of the data was 250 m with a vertical accuracy of 12.5 m.

High Capability RadarSounder (HICARS)
HICARS is operated by the University of Texas, Institute for Geophysics (Peters et al., 2007(Peters et al., , 2005 and was flown 2011 and 2012 as part of a joint UK-US, NERC/NSF funded project called Greenland Outlet Glacier Geophysics (GrOGG). It is a 60 MHz phase coherent pulsed radar with a 15 MHz bandwidth. It has been flown extensively in Antarctica including surveys over the Thwaites Glacier catchment and surveys of large sectors of East Antarctica (Young et al., 2011). New algorithms employed with HICARS allow a horizontal resolution of less than 5 m and vertical resolution of ±10 m.

Technical University of Denmark (TUD)
A small amount of data are included which were collected in the 1970s by the TUD using a 60 MHz echo sounder, analogue recording and inertial navigation. The quality of these data is, in general, substantially poorer than those acquired more recently so they are only used when there are no other data sources within 50 km ( Fig. 1; Bamber et al., 2001b).

Surface elevations
All the airborne data were provided as ice thickness measurements but we require bed elevation. Thus, a surface elevation estimate is required. For most of the CReSIS and OIB campaigns, the flights also operated the NASA Airborne Topographic Mapper (ATM): an airborne laser altimeter. ATM scans the surface and is resampled to a horizontal spacing of 50 m with data averaged into 80 m diameter platelets. ATM was not always available and so our secondary source of surface elevations was from a DEM of the whole of Greenland. The surface DEM used was produced as part of the Greenland Ice Mapping Project (GIMP). It has a horizontal resolution of 90 m and is a multi-sensor DEM derived from data collected between 2000 and 2009. This DEM was created from MODIS, AVHRR, ASTER, SPOT and RADARSAT datasets merged with the ASTER GDEM and the (Bamber et al., 2001a) DEM vertically co-registered using ICESat data. Validation against ICESat data, indicated vertical errors of ±5 m on the ice sheet and ±7 m for the unglaciated margins (Howat, personal communication, 2012).

Bathymetry
Numerical modelling over long timescales, such as glacialinterglacial cycles, or spinning up the thermodynamics in an ice sheet model requires basal geometry that extends out to the continental shelf: i.e. as far as the maximum glacial extent, which can reach to the shelf edge several hundred kilometres from the present-day ice limit (Evans et al., 2009;Dowdeswell et al., 2010;Cofaigh et al., 2013). To achieve this requires inclusion of bathymetric data. Here we used the International Bathymetric Chart of the Arctic Ocean (IBCAO) v3 (Jakobsson et al., 2012). This is an interpolation of various bathymetric data from the entire Arctic Ocean and a DEM for Greenland, Ellesmere Island and Iceland. It was supplemented with additional data from soundings in the Jakobshavn fjord, which were included in the CRe-SIS Jakobshavn grid described above. Other changes to this dataset, and the reasons for these changes, are as described in Sect. 3.2.3.

Methods
All the airborne data were transformed onto a polar stereographic projection with standard parallel at 71 • N and a central meridian of 39 • W. All invalid data, as defined by the instrument teams, all data outside Greenland, and all data with ice thickness less than 0 m were removed. Rather than interpolating ice thickness, which has a discontinuity at the ice sheet margin, we interpolate bed elevation, as this varies smoothly across the ice edge. By combining the high resolution surface topography with the bed data we can create a more realistic ice margin. Where CReSIS thickness data were acquired within 2 days and within 1 km of an ATM surface elevation, the ATM estimate is used to convert to bed elevation. For 18 % of ice thickness estimates from CReSIS when ATM was also flown, no surface elevation estimate was recorded. In those cases, and for all other datasets, the GIMP DEM was used. This introduces a potential error in the derived bed elevation if there has been a change in surface elevation (dh/dt) between the acquisition of the ice surface and thickness data.
For 74 % this is not relevant because simultaneous LIDAR data (such as ATM) are available. For the remainder, the GIMP DEM was used. The exact time stamp of GIMP is unclear as it was derived from a mosaic of images. For Jakobshavn, Helheim, and Kangerdlugssuaq glaciers, however, it was based on imagery from 2007, and for Petermann glacier from 2003 (personal communication; I. Howat, 2012). As the largest dh/dt values are found over these outlet glaciers, a correction was applied in these areas. Dh/dt values were taken from an ICESat based estimate for Greenland, covering the period 2003-2009, interpolated to 1km resolution (Hurkmans et al., 2013). Annual dh/dt values were added for years between the IPR acquisition date and the GIMP time stamp for the area of interest. We assumed GIMP is representative for approximately mid-2003 and 2007 and using the month of the IPR measurement, the appropriate fractions of the dh/dt values for the GIMP year (2003GIMP year ( or 2007 and the IPR year were taken into account. Of the 26 % of data points for which GIMP was needed, 22 % were located over one of the four glaciers mentioned above and have been corrected, therefore, for dh/dt. For most of the remaining 4 % (e.g. nearly all the AWI data), the data are in the interior where dh/dt values are at the few cm a −1 level.
To interpolate the bed, it is necessary to delineate glaciated from ice-free terrain. A land surface mask was created by merging a number of data sources. The coastline of Greenland was smoothed to a 1 km resolution from the Danish Ministry of Environment (formerly KMS) 1 : 2 500 000 scale vector maps of the coast. The Canadian Arctic and Iceland were separated from the Greenland coast by identifying distinct polygons. The ice sheet and periphery ice caps were identified using a binary mask (Howat and Negrete, 2013) produced from a combination of Landsat 7 panchromatic band imagery from July-September 1999-2001 and The Cryosphere, 7, 499-510, 2013 www.the-cryosphere.net/7/499/2013/ RADARSAT-1 SAR amplitude images from autumn 2000. Data were provided at 180 m horizontal spacing and reprojected onto the polar stereographic grid at 1 km resolution. Grid cells with over 50 % ice cover were considered to be glaciated. Ice shelves were categorised separately from the rest of the ice sheet. The existence of an ice shelf was determined by the presence of a grounding line from Interferometric Synthetic Aperture Radar. Grounding lines were provided by Eric Rignot (JPL) and are a sub-set of those previously published (Rignot et al., 1997). The ice shelf front was determined by a number of means. In areas with ATM coverage, the ice front was located by a step in the elevations. This was confirmed by ice front locations from the KMS maps.
In the absence of ATM data, the KMS ice fronts were used alone. On Storstrommen and Ostenfeld glaciers, no front was present in maps or ATM data but InSAR confirmed the presence of an ice shelf. In these cases, a small shelf was added based on Google Earth imagery. No attempt was made to determine the bed elevation of peripheral glaciers (i.e. separated from the ice sheet) where no thickness data exist (see Fig. 1) and, in this case, it is the glacier surface elevation that is obtained, if resolved at all at 1 km, i.e. ice thickness is zero.

Data editing
The data were gridded with 5 km postings and a 3 standard deviation filter was applied twice to remove elevation outliers. This removed 0.3 % of data points. Visual inspection indicated that the filter had removed noisy data but a small number of anomalous measurements remained. These were tracks which were not picking the ice/bed interface but appeared to be tracking an internal layer. A coarse filter was applied whereby a data point was removed if the bed elevation deviated by more than 500 m from the previous estimate (Bamber et al., 2001b). This was only applied in areas where there was previous data coverage and not in areas with large relief or areas of high surface velocity > 100 m yr −1 ). Visual inspection also led to the removal of several other tracks from the CReSIS data after examination of the echograms. In all, 98.6 % of the data from the various campaigns were determined to be over ice and of sufficient quality to be included.

Interpolation
All data which passed the quality checks were locally averaged into a quasi-regular 1 km and 2.5 km resolution grid, which reduced the disparity in along and across track spacing of the data. The resolution of the two grids was chosen based on the data density as indicated in Fig. 1 and Fig. S4 in the Supplement. A 1 km posting DEM results in about 20 % of grid points containing data in areas where the acrosstrack spacing is greater than 20 km. For sectors of the interior, across-track spacing can be more than 50 km (Fig. 1).
In general, however, bed gradients are smaller in these areas and a lower resolution is adequate for capturing the large scale relief. In areas identified by the land/ice mask as being unglaciated, surface elevation data from the GIMP DEM were included in the quasi-regular grid. The 1 and 2.5 km quasi-regular grids were interpolated to regular grids using ordinary kriging. The GSLIB library (Deutsch and Journel, 1997) was used to calculate variograms. Separate variograms for the two resolutions were calculated and an exponential function was fitted to the first 100 km of each variogram using a nonlinear least squares scheme. These variograms were used to interpolate the quasiregular grids using a nugget of 50 m to take account of uncertainty in the airborne data. A maximum of 50 quasi-regular grid points were considered and the maximum search radius was set to 250 km so that a result was obtained everywhere. The 2.5 km grid was bilinearly interpolated to 1 km resolution. This avoids artefacts produced from interpolating sparse data at 1 km, and results in a grid at a single horizontal posting, which provides a simpler data structure compared with a multi-resolution or nested grid approach. The two grids were combined, with the higher resolution one being used for all areas which had sufficient data density (Fig. S4 in the Supplement) and the lower resolution grid used elsewhere. The merging was done using a Hermite basis function of width 20 km across the boundary.
Ice thickness was derived by subtracting the bed from the GIMP DEM for every grid point defined as ice covered in the mask. The minimum ice thickness at the margin was set to 50 m. Wherever the thickness was less than this, the bed was lowered to be 50 m below the GIMP surface elevation. This is necessary because thickness data around much of the margin of the ice sheet does not exist (Fig. 1) and therefore, subglacial bed elevations cannot be determined directly in these areas but must be interpolated close to a discontinuity in thickness and surface elevation.

Ice shelves
Nine ice shelves are present around the Greenland ice sheet according to our land mask. The surface elevation for these was, in general, found by combining GIMP, ATM and IPR data, which had been converted to surface elevation using the assumption of hydrostatic equilibrium. Some ice shelf gridpoints contain no elevation data from any source, in which case the nearest neighbour interpolation from other ice shelf elevations was used. For Nioghalvfjerdsbrae and Zachariae Isstrom, only ATM and airborne elevations were used as GIMP elevations were not in agreement with the other data sources.
Ice shelf thickness was calculated from surface elevation using the assumption of hydrostatic equilibrium (see Griggs and Bamber, 2011 for a full description of the method). A constant firn density correction of 10 m was used throughout. An ocean water density of 1027 kg m −3 and an ice density of 917 kg m −3 were chosen. The same parameters were used to convert the airborne ice thicknesses to surface elevations and www.the-cryosphere.net/7/499/2013/ The Cryosphere, 7, 499-510, 2013 back again, which was done to increase data coverage. In areas where the surface elevation was less than 10 m above sea level or where the assumption of hydrostatic equilibrium is invalid, surface elevations were interpolated from thicker neighbouring data points.

Mass conservation
In many areas, there were no airborne data within a few kilometres of glacier termini. If there are unglaciated regions closer to the terminus than the closest airborne data, they dominate the interpolated elevation at the terminus, artificially raising the bed elevations. This was corrected for in two ways. In four areas, we had a flight line that crossed the fast-flowing region of the glacier. In these cases, we took the bed elevations at the airborne data points along with the known velocity (Joughin et al., 2010) and the 30 yr average modelled surface mass balance ) and calculated the expected bed elevation based on the principle of mass conservation. The direction of flow of the ice from the location of the airborne data to the terminus was determined from the velocity vector. A more sophisticated approach has been developed and demonstrated on 79 North Glacier (Morlighem et al., 2011) and we intend to incorporate results using this method in future releases of the product. Figure 2 shows the result of correcting the bed elevations using this method for one outlet. In Fig. 2a the original elevations are shown and Fig. 2b shows the bed after the application of the technique. A region of 7 km length is altered from the position of an airborne track inland, following the region of faster flowing ice as shown in Fig. 2c. It is clear that after this change, thick ice is now able to flow along the trough where previously a ridge was mapped due to the inclusion of data from locally higher unglaciated terrain close to the terminus.
In a number of other locations, the same situation arises but either there are no velocity measurements, the velocity vectors make the ice flow out of the region of fast flow, unsurveyed tributaries join the main flow seaward of the airborne data, or the airborne measurements do not cross the entire fast flow channel. In these cases, our mass conservation scheme cannot be easily implemented and instead, linear interpolation was used to remove the artificial ridges at the terminus. In all cases where data were changed from the original interpolated values, a mask is provided which notes the change made, the reason, and the original interpolated value.

Bathymetry interpolation
The interpolated bed elevations were merged with the bathymetry of the fjords and oceans around Greenland (Jakobsson et al., 2012). A smoothing distance of 3 km was used on the seaward side of the coastline to merge  data are particularly high and appear to follow the ice surface rather than the fjord bottom (which is at about 500 m below sea level) . In Jakobshavn fjord, we replaced the IB-CAO bathymetry with that from the CReSIS grid described in Sect. 2.1.1. In other areas, the IBCAO interpolation contains a fjord which is much shallower than the bed elevation at the grounding line or glacier terminus. For example at Hagen Brae, the fjord just seaward of the grounding line was 200 m shallower than the bedrock elevation at the grounding line. To ensure realistic ice flow pathways, we lowered the bathymetry in the vicinity of the fjord for regions affected in this way to create a smooth surface without a discontinuity at the land/ocean margin. Where ice shelves are present, a similar approach was undertaken involving interpolation of the bed elevation at the grounding line seaward to the first measured IBCAO value. A minimum depth for the cavity beneath the ice shelf of 10 m was imposed, distal from the grounding line.
We have not included any new bathymetric data as this is beyond the scope of this study. In the affected areas, we interpolated between the grounding line/glacier terminus and beyond the fjord mouth using triangulation. The results were smoothed over 2 km around the triangulated values. End points for the interpolation were chosen to ensure that ice flows out of the fjord. The approach is somewhat subjective and results depend upon the choice of endpoints and can produce steep gradients in the across fjord direction beyond the lateral region of adjustment. As a consequence, we provide a mask with the data which indicates where the bathymetry has been changed in this manner.

Results and discussion
The final bed DEM was referenced to the EIGEN-GL04C geoid (Forste et al., 2008) and is shown in Fig. 3 alongside the ice thickness grid. All the figures plotting bed elevation are referenced to the geoid, which has a range of around 10 m in the north-west to 65 m in the south-east. In addition to the bed DEM, we also produced grids of ice thickness (Fig. 3b), surface elevation, error maps of surface and bed elevation (Fig. 7), the land surface mask, the geoid-ellipsoid separation, a mask showing changes made to the bed elevation postinterpolation, the bed elevation and ice thickness without any intervention, a mask showing the data sources used for the ice shelves and a grid of the number of IPR data per grid cell. The bed dataset includes detail which was not visible in previous compilations and improves the representation of many features previously observed. This improvement is most noticeable in areas of relatively high relief close to the ice margin and, in particular, where dense grids were flown (Fig. 1, and Fig. S1 in the Supplement). In the Jakobshavn catchment ( Figs. 4 and 5), there is a dendritic channel system extending for about 325 km from the current grounding line into the interior almost as far as the ice divide (Hoch et al., 2011). It seems likely that this is a palaeo-fluvial feature that predates ice cover in Greenland and may be important for subglacial water routing. The presence of a deep trough extending into the interior has been previously reported, based on SAR processing of MCoRDS data (Hoch et al., 2011). We make a direct comparison between the new bed DEM and the previous compilation (Bamber et al., 2001b) in two regions to illustrate the improvements made. Figure 4a shows the Jakobshavn region in the new dataset (see black box in Fig. 3 for location) and Fig. 4b shows the same region in the older dataset. First, it is apparent that the true resolution (as opposed to the grid spacing) of the new DEM is significantly improved throughout the region, better characterising the undulating terrain to the north of the trough. Secondly, the deep trough under the main fjord is present. The previous compilation showed almost no evidence of this trough due to (i) a lack of bed returns in this area and (ii) the resolution of the older grid at 5 km. The width of the trough in the new DEM is 3-4 km and the region of fastest flow coincides fairly well with the location of the deepest ice (Fig. 4c). The trough in the new DEM is 1366 m below sea level at its deepest point compared to a maximum depth over the entire region of 556 m below sea level in the older dataset. The main trough of Jakobshavn Isbrae is not continuous in the new dataset, disappearing around 100 km on the x-axis and reappearing at about 140 km. This does not imply that the The Cryosphere, 7, 499-510, 2013 www.the-cryosphere.net/7/499/2013/ trough is discontinuous, but only that there are insufficient data to confirm the trough's presence or otherwise in this region. Figure 4d and e show a similar comparison for a 200 km by 500 km area along the north-west coast of Greenland (red box in Fig. 3). In the older DEM, there is some evidence of bed troughs but they are wide and not well aligned with the areas of fast flow seen in Fig. 4f. Figure 4d shows numerous troughs, all aligned with areas of fast flow. The wide unconstrained minimum seen in Fig. 4e at a depth of 438 m below sea level has now become a deep, narrow trough under a region of fast flow with a maximum depth of 1219 m below sea level. This improvement is due mainly to the increased data coverage and, to a lesser extent, the higher resolution of the DEM. Similar comparisons can be made elsewhere and the full fidelity of the dataset becomes apparent when examining smaller regions of a few hundred kilometres in extent. To further illustrate this we have produced shaded relief plots of the Jakabshavn and Kangerdlugssuaq regions in Fig. 5a and b, respectively. Figure 5a indicates the small-scale structure of the dendritic network of troughs extending inland and the complex bathymetry near the mouth of the glacier. The effect of lowering the bathymetry near the mouth of the fjord is evident by the steep slopes at the lateral margins of the fjord (around −2150 Northing, −500 y-axis). The adjustment made ensures a smooth surface in the ice flow direction but not in the across flow. Without additional bathymetry data, further improvements to these adjustments will be difficult and somewhat arbitrary.
The impact of coarser flight track spacing on the topographic detail can be clearly seen at −2200 to −2260 Northing, −150 to −50 on the y-axis. This region possesses a feature that, likely, would be deeper and narrower if adequately resolved, like the trough adjacent and just north of it. The Kangerdlugssuaq fjord is another region where the IBCAO bathymetry had to be lowered by several hundred metres to avoid an artificial "cliff" at the glacier terminus. The fjord in this region is up to 700 m below sea level in our dataset, which is only partially captured by the original IBCAO dataset. Inland, the trough that the glacier follows is around 1500 m in depth and in places less than 10 km in width. This type of extreme topographic relief requires both dense data sampling (Fig. 1) and appropriate resolution interpolation (1 km in this region) to adequately define it.

Error assessment
For most applications, a reliable estimate of the uncertainty in the DEM is essential. To determine this we considered (i) the random error in the thickness observation and (ii) the impact of interpolation. The former was assessed from a trackto-track difference analysis of the data (a combination of repeat track and quasi-crossover differences). We consider the differences between any two measurements obtained at different times within a 50 m area. This has the advantage of including along track differences for repeated flights. In total, 24.1 million differences were calculated. The data were split by campaign to calculate both inter-and intra-campaign differences. These data, strictly, only provide information on repeatability and not systematic biases but, because we are considering inter-campaign differences, we believe that biases due, for example, to timing, navigation or radar calibration errors will be captured in this analysis. Only biases common across all campaigns (such as a common error in the radio wave velocity in ice used) will not be seen. Table 2 summarises the differences for each set of campaigns grouped by instrument and/or institute. Figure 6a shows the histogram of cross-over and along-track differences, which has a symmetric distribution with a bias close to zero. We estimate the random error as 1/ √ 2 (because a difference contains two observations) of the standard deviation of the inter-campaign difference of the nearest data to a grid point. In the case of TUD data, where there are no differences, we use the standard deviation of all inter-and intracampaign values.
Next we consider the error due to interpolation. This increases with distance from an observation but is also a function of the properties of the underlying surface. We estimate this uncertainty using a bootstrap approach for two classes of bed topography: coastal and interior.
To differentiate these two classes, we calculated the standard deviation of the bed elevation in overlapping 50 km boxes and used this as a measure of basal roughness (Fig. S3 in the Supplement). From visual inspection, a standard deviation threshold of 170 m was used to distinguish the two classes of bed. For these two zones, all CreSIS data from 2000 onward were used in a bootstrapping approach to determine the effect of interpolation as a function of distance from an observation. The difference between the interpolated and observed elevation as a function of distance was used as a measure of the interpolation error. The standard deviation of these differences for the two zones is shown in Fig. 6b. Exponential curves were fit to the differences and are also plotted in Fig. 6b. These curves were used to determine the error due to interpolation as a function of distance.
The two sources of error were combined in quadrature to produce a map of the uncertainty in the bed elevation (Fig. 7). Standard deviation of difference between bed elevations created using all data and using a sub-set of the data in "marginal" areas with a roughness standard deviation greater than 170 m (asterisks) and the "interior" region where the roughness was below this threshold (diamonds) (see Fig. S3 in the Supplement). An exponential fit to the data is shown as solid ("margins") and dashed ("interior") lines.
As expected the error increases significantly with distance from an observation and the largest errors occur in coastal regions where extrapolation, rather than interpolation, was needed. For the peripheral ice caps where, in most cases, there are no airborne data, the bed elevation is poorly constrained and these areas are included only for completeness. Better estimates of the ice thickness in these areas could be estimated using an ice-surface area/volume scaling approach (Bahr et al., 1997).

Conclusions
A large volume of high quality new data have become available since the last comprehensive dataset of ice thickness in Greenland was compiled (Bamber et al., 2001b). We improve on the earlier compilation in several respects. Most importantly, we have included extensive new data sets acquired by several different groups over the last decade. A significant effort has been made, during this period, to sound the bed of fast-flowing outlet glaciers which were either missed or proved challenging targets for the previous generation of IPR systems due to high attenuation and clutter. These new datasets also provide improved coverage of the, previously, sparsely surveyed interior. The greater coverage and dense network of flight lines in many coastal areas, and some in--8•10 5 -6•10 5 -4•10 5 -2•10 5 0 2•10 5 4•10 5 6•10 5 land regions, (Fig. 1 and Fig. S1 in the Supplement) supports a grid spacing of 1 km, five times higher than the previous compilation. As a consequence, many basal features, in particular basal troughs containing outlet glaciers, are now properly resolved in the bed topography. In the ice sheet interior, flight lines are less dense but the topography is generally smoother and, therefore, a coarser resolution is adequate. For convenience, the bed DEM is provided at a single posting of 1 km alongside a grid indicating whether the value is interpolated or based on observation. An error map for the DEM was also calculated and indicates areas where additional data would be particularly useful. The estimated volume of the ice sheet is 2.96 × 10 6 km 3 compared with 2.93 × 10 6 km 3 obtained previously (Bamber et al., 2001b). In all, 22 % of the ice sheet bed is below sea level and, accounting for this and the thickness of the firn layer but excluding any glacioisostatic adjustment, we estimate that the ice sheet has the potential to raise global mean sea level by 7.36 m were it all to melt. The data presented here represent a major advance in our knowledge of the topography of the bedrock of Greenland. However, there are several areas where data are currently lacking such as along parts of the central and north eastern margins of the ice sheet (Fig. S1). While the interior of Greenland is relatively smooth, there are still many areas where the distance to the closest observation is more than 50 km, resulting in uncertainties in bed elevation exceeding ±100 m and missing potential short wavelength relief that is evident, for example, inland from Jakobshavn Isbrae (Fig. 5a). Although this may be less important for numerical modelling, such detail provides valuable insights into the genesis of the subglacial landforms and the geomorphology of the bed. Recent results suggest reduced uncertainty and significantly increased spatial detail can be obtained from radar tomography (Jezek et al., 2011;Paden et al., 2010), while mass conservation approaches have also shown promise for poorly sampled outlet glaciers (Morlighem et al., 2011). Application of these techniques could significantly improve current mapping by decreasing the amount of interpolation needed over deeply incised outlet glaciers. Bathymetry seaward of the glaciers and beneath ice shelves is currently poorly characterised in some areas and non-existent in most. This is a major gap in our current knowledge and requires further effort by the community. New data are being acquired over Greenland continuously and we intend to issue new releases of the products when warranted. The complete set of grids, metadata and documentation are available in netcdf and geotiff format from the lead author (JLB). Users will be notified of new releases as they become available.