Skip to main content


  • Access
  • Cited by 21
  • Cited by
    This article has been cited by the following publications. This list is generated based on data provided by CrossRef.

    Hoffman, Matthew J. Perego, Mauro Andrews, Lauren C. Price, Stephen F. Neumann, Thomas A. Johnson, Jesse V. Catania, Ginny and Lüthi, Martin P. 2018. Widespread Moulin Formation During Supraglacial Lake Drainages in Greenland. Geophysical Research Letters, Vol. 45, Issue. 2, p. 778.

    Morlighem, M. Williams, C. N. Rignot, E. An, L. Arndt, J. E. Bamber, J. L. Catania, G. Chauché, N. Dowdeswell, J. A. Dorschel, B. Fenty, I. Hogan, K. Howat, I. Hubbard, A. Jakobsson, M. Jordan, T. M. Kjeldsen, K. K. Millan, R. Mayer, L. Mouginot, J. Noël, B. P. Y. O'Cofaigh, C. Palmer, S. Rysgaard, S. Seroussi, H. Siegert, M. J. Slabon, P. Straneo, F. van den Broeke, M. R. Weinrebe, W. Wood, M. and Zinglersen, K. B. 2017. BedMachine v3: Complete Bed Topography and Ocean Bathymetry Mapping of Greenland From Multibeam Echo Sounding Combined With Mass Conservation. Geophysical Research Letters, Vol. 44, Issue. 21, p. 11,051.

    Carrivick, Jonathan L. Yde, Jacob Russell, Andrew J. Quincey, Duncan J. Ingeman-Nielsen, Thomas and Mallalieu, Joseph 2017. Ice-margin and meltwater dynamics during the mid-Holocene in the Kangerlussuaq area of west Greenland. Boreas, Vol. 46, Issue. 3, p. 369.

    Vincent, Joseph A. and Arnold, Emily J. 2017. Beamforming sensitivity of airborne distributed arrays to flight tracking and vehicle dynamics. p. 1.

    Colonia, Daniel Torres, Judith Haeberli, Wilfried Schauwecker, Simone Braendle, Eliane Giraldez, Claudia and Cochachin, Alejo 2017. Compiling an Inventory of Glacier-Bed Overdeepenings and Potential New Lakes in De-Glaciating Areas of the Peruvian Andes: Approach, First Results, and Perspectives for Adaptation to Climate Change. Water, Vol. 9, Issue. 5, p. 336.

    Guan, Yawen Haran, Murali and Pollard, David 2017. Inferring ice thickness from a glacier dynamics model and multiple surface data sets. Environmetrics, p. e2460.

    Brinkerhoff, Douglas J. Aschwanden, Andy and Truffer, Martin 2016. Bayesian Inference of Subglacial Topography Using Mass Conservation. Frontiers in Earth Science, Vol. 4, Issue. ,

    Aschwanden, Andy Fahnestock, Mark A. and Truffer, Martin 2016. Complex Greenland outlet glacier flow captured. Nature Communications, Vol. 7, Issue. , p. 10524.

    Larour, E. and Schlegel, N. 2016. On ISSM and leveraging the Cloud towards faster quantification of the uncertainty in ice-sheet mass balance projections. Computers & Geosciences, Vol. 96, Issue. , p. 193.

    Karlstrom, Leif and Yang, Kang 2016. Fluvial supraglacial landscape evolution on the Greenland Ice Sheet. Geophysical Research Letters, Vol. 43, Issue. 6, p. 2683.

    Vasskog, Kristian Langebroek, Petra M. Andrews, John T. Nilsen, Jan Even Ø. and Nesje, Atle 2015. The Greenland Ice Sheet during the last glacial cycle: Current ice loss and contribution to sea-level rise from a palaeoclimatic perspective. Earth-Science Reviews, Vol. 150, Issue. , p. 45.

    Young, Nicolás E. and Briner, Jason P. 2015. Holocene evolution of the western Greenland Ice Sheet: Assessing geophysical ice-sheet models with geological reconstructions of ice-margin change. Quaternary Science Reviews, Vol. 114, Issue. , p. 1.

    Khan, Shfaqat A Aschwanden, Andy Bjørk, Anders A Wahr, John Kjeldsen, Kristian K and Kjær, Kurt H 2015. Greenland ice sheet mass balance: a review. Reports on Progress in Physics, Vol. 78, Issue. 4, p. 046801.

    Margold, Martin Stokes, Chris R. and Clark, Chris D. 2015. Ice streams in the Laurentide Ice Sheet: Identification, characteristics and comparison to modern ice sheets. Earth-Science Reviews, Vol. 143, Issue. , p. 117.

    Rignot, E. Mouginot, J. Morlighem, M. Seroussi, H. and Scheuchl, B. 2014. Widespread, rapid grounding line retreat of Pine Island, Thwaites, Smith, and Kohler glaciers, West Antarctica, from 1992 to 2011. Geophysical Research Letters, Vol. 41, Issue. 10, p. 3502.

    Herzfeld, Ute C. McDonald, Brian W. Wallin, Bruce F. Chen, Phillip A. Mayer, Helmut Paden, John and Leuschen, Carlton J. 2014. The trough-system algorithm and its application to spatial modeling of Greenland subglacial topography. Annals of Glaciology, Vol. 55, Issue. 67, p. 115.

    Christianson, Knut Peters, Leo E. Alley, Richard B. Anandakrishnan, Sridhar Jacobel, Robert W. Riverman, Kiya L. Muto, Atsuhiro and Keisling, Benjamin A. 2014. Dilatant till facilitates ice-stream flow in northeast Greenland. Earth and Planetary Science Letters, Vol. 401, Issue. , p. 57.

    Tedstone, A. J. Nienow, P. W. Gourmelen, N. and Sole, A. J. 2014. Greenland ice sheet annual motion insensitive to spatial variations in subglacial hydraulic structure. Geophysical Research Letters, Vol. 41, Issue. 24, p. 8910.

    Morlighem, M. Rignot, E. Mouginot, J. Seroussi, H. and Larour, E. 2014. High-resolution ice-thickness mapping in South Greenland. Annals of Glaciology, Vol. 55, Issue. 67, p. 64.

    Morlighem, M. Rignot, E. Mouginot, J. Seroussi, H. and Larour, E. 2014. Deeply incised submarine glacial valleys beneath the Greenland ice sheet. Nature Geoscience, Vol. 7, Issue. 6, p. 418.




      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        High-resolution bed topography mapping of Russell Glacier, Greenland, inferred from Operation IceBridge data
        Available formats
        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        High-resolution bed topography mapping of Russell Glacier, Greenland, inferred from Operation IceBridge data
        Available formats
        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        High-resolution bed topography mapping of Russell Glacier, Greenland, inferred from Operation IceBridge data
        Available formats
Export citation


Detailed maps of bed elevation and ice thickness are essential for understanding and projecting the evolution of the ice sheets. Such maps are traditionally obtained using airborne radar-sounding profiler data interpolated onto regular grids using geostatistical tools such as kriging. Here we compare three mapping techniques applied to a dense radar survey of Russell Glacier, West Greenland, by NASA Operation IceBridge: (1) radar tomography (RT) processing of the radar data to map the bed elevation, (2) interpolation of radar-derived thickness by ordinary kriging (KR) and (3) reconstruction of ice thickness based on the principles of mass conservation (MC) combining radar-sounding profiler and ice motion data. RT eliminates ambiguities caused by off-nadir reflections, but is spatially limited. KR yields a standard error in bed elevation of 35 m, but large errors (>300 m a−1) in flux divergence when combined with ice motion data. MC yields a comparable performance in bed elevation mapping, and errors smaller than 1 m a−1 in flux divergence. When the number of radar-sounding tracks is reduced, the performance of KR decreases more rapidly than for MC. Our study site shows that MC is capable of maintaining precision levels of 60 m at 400 m posting with flight tracks separated by 5 km.


Bed topography and ice thickness are essential characteristics of glaciers and ice sheets for many glaciological applications, but they remain difficult to measure remotely. This is especially true for ice-sheet coastal sectors and mountain glaciers, where the ice surface is significantly broken up, the ice substrate is near or at the melting point, and scattering and signal absorption are high. At present, bed elevation and ice thickness are mapped using airborne radar-sounding profilers, i.e. sounders that detect radar echoes from the bed directly beneath the aircraft carrier and convert these echoes into bed elevation, surface elevation and ice thickness (Gogineni and others, 2001). The radar-profiler mode (PM) data are then interpolated onto a regular grid using various methods, most commonly geostatistical techniques, such as kriging (e.g. Bamber and others, 2013; Fretwell and others, 2013) or more sophisticated algorithms (e.g. Herzfeld and others, 2011). The absolute precision of these gridded products is difficult to quantify in the absence of independent measurements of ice thickness. This is an issue for ice-sheet numerical models, which employ these data, especially in areas with few data and potentially large interpolation errors.

Russell Glacier is a land-terminating glacier in central-west Greenland, 20 km east of Kangerlussuaq near Søndre Strømfjord (67° N, 50°W), that has not changed significantly in recent years (Pritchard and others, 2009). In 2010 and 2011, NASA’s Operation IceBridge (OIB) conducted a dense radar survey (500 m track spacing) of Russell Glacier (Fig. 1), to evaluate our present-day mapping capability, evaluate the precision of the radar data (e.g. through crossover analysis), compare different techniques of bed mapping and evaluate their errors and potential impact on numerical ice flow models. Here we compare the results obtained using three methods: (1) a radar tomography (RT) processing of the radar data to produce a gridded map of bed elevation, (2) a method of ordinary kriging (KR) applied to radar-derived ice thickness data along the flight tracks to obtain a regular grid of ice thickness, and (3) a reconstruction of ice thickness on a regular grid using the principle of mass conservation (MC) applied to the radar-sounding derived thickness data combined with ice motion data from satellite radar interferometry within an optimization framework. The RT method is a novel, high-resolution processing technique, which requires a dense grid of flight tracks to produce a contiguous map of bed elevation over a large area. The second method is a standard method for producing ice thickness and bed elevation maps from radar profilers. The third method calculates ice thickness values on a triangle mesh and solves for the mass conservation using the finite-element method, while at the same time minimizing the departure of the results from the original radar-derived ice thickness data (Morlighem and others, 2011). These methods are here compared over an exceptionally dense network of radar observations. This setting is ideal for evaluating errors in radar-derived ice thickness, bed elevation and flux divergence and for comparing the performance of various methods. Because this level of mapping is not practical over large areas, we also investigate the impact of reducing the number of tracks used to survey bed elevation and ice thickness to determine the configuration required to achieve a given precision in bed elevation and ice thickness mapping. We discuss the impacts of the results on ice-sheet modeling and conclude with recommendations for bed elevation mapping of glaciers and ice sheets.

Fig. 1. Bed topography of Russell Glacier inferred using the mass continuity (MC) method at 400 m resolution, calculated as surface elevation minus ice thickness. The dashed white line defines the limits of the model domain. NASA OIB flight tracks are indicated as black solid lines. Surface elevation is from the Greenland Ice Mapping Project (personal communication from I. Howat, 2012), including over ice-free zone. Note the agreement between ice-free and ice-covered elevations.

Data and Methods

Ice velocity data

We measure ice velocity using a combination of Advanced Land Observing Satellite (ALOS) Phased Array-type L-band Synthetic Aperture Radar (PALSAR), RADARSAT-1 and Envisat Advanced SAR (ASAR) speckle-tracking data (Rignot and Mouginot, 2012). Data resolution is 350 m, precision in speed is 10 m a−1 and in direction is 1° (Fig. 2). The glacier speed is ∼150 m a−1 at the equilibrium-line altitude, and decreases to a few meters per year at the ice margin. The velocity map reveals that Russell Glacier shares a common source with two other glaciers: Isunguata Sermia to the north, which flows at a comparable speed, and a glacier which we identify as South Russell to the south, which flows at ∼380 m a−1.

Fig. 2. Observed ice surface motion derived from satellite radar interferometry data using speckle tracking (Rignot and Mouginot, 2012) overlaid on a backscatter image from Envisat ASAR 2009.

OIB dataset

NASA OIB collected radar-sounding data at 150 MHz (Multichannel Coherent Radar Depth Sounder (MCoRDS); Leuschen and others, 2011) along flight tracks for the most part perpendicular to the coastline and parallel to the ice flow direction, spanning a region about 80 km east–west by 40 km north–south, with an average track spacing of ∼500 m (Fig. 1). This dense dataset is exceptional by Greenland standards, where most tracks are separated by several km to several tens of km. The radar data were processed in profiling mode (PM) to yield ice thickness and bed elevation with a nominal precision of 10 m for ice thickness (Gogineni and others, 2001). The actual precision of the data varies with location and the quality of the bed echoes. Primary error sources include system electronic noise, multiple reflections (off-nadir) and uncertainties in dielectric properties of ice. Corrections for the firn depth are not always applied or assume a constant value. These errors create spurious bed returns in the trace data that may deviate from the true nadir reflection (Wu and others, 2011). A crossover analysis of all the CReSIS (Center for Remote Sensing of Ice Sheets, University of Kansas) data yields a standard error in thickness of 31 m for Russell Glacier, which is a representative value of the actual precision of the data.

Radar tomography (RT)

We apply RT processing to a portion of the OIB data where the flight pattern is most regular and provides a nearly contiguous coverage of the bed. The processed area, located in the northeastern sector of the survey grid, is 15 km × 30 km in size (Fig. 3e). Unlike conventional radar-sounding profilers, where echoes from the left and right sides are not resolvable, RT relies on multiple antennas and multiple receivers in the across-track direction to differentiate the echoes (Paden and others, 2010; Wu and others, 2011). RT resolves the bottom signal arrival angle and strength, which are then used to derive the bed topography. The measurable swath width depends on the signal-to-noise ratio (SNR), which can vary strongly depending on basal roughness and ice absorption rates. RT produces here fine-resolution (20 m spacing) estimates of bed elevation over a swath width of 1600 m, with a standard noise of 10 m (Wu and others, 2011). We use this map as one of our reference datasets to evaluate the performance of other techniques since it is the only technique capable of dealing with off-nadir errors, which have recently been shown to be significant in the conventional PM radar processing (Wu and others, 2011).

Fig. 3. Bed topography of Russell Glacier with respect to mean sea level, with elevation contours at 50, 200, 350 and 500 m, from (a) Bamber and others (2001), (b) Bamber and others (2013), (c) mass continuity (MC), (d) ordinary kriging (KR) and (e) radar tomography (RT); (f) zoom of (e).

Ordinary kriging (KR)

Kriging combines nearby observations to predict bed elevation (or ice thickness) at other locations (Deutsch and Journel, 1997). This technique has been used extensively for mapping surface elevation, ice thickness (e.g. Griggs and Bamber, 2011) and bed elevation (e.g. Bamber and others, 2013; Fretwell and others, 2013). Ice thickness, H, at a given location (x, y) is assumed to be a linear combination of a subset of observations, Hi :


where λi (x, y) are the kriging weights and is the subset of nearby observations. In our application, includes the closest 200 measurements within a radius of 50 km around (x, y). We employ a particular type of kriging named ordinary kriging, which is a very common type of kriging (Deutsch and Journel, 1997). Ordinary kriging assumes a constant but unknown mean and the kriging weights fulfill the unbiasedness condition:


The weights are calculated to minimize the variance of the estimation error, assuming that ice thickness is a random process for which we estimate the covariance structure. We derive the covariance from a semivariogram assuming an exponential model. We use a nugget effect of 50 m to account for measurement errors, and find a best fit between model and observations for a range of 4.3 km and a sill of 58 km. The data are interpolated on a regular grid at 500 m spacing. Bed elevation is deduced by subtracting H from the digital elevation model (DEM) of the Greenland Ice Mapping Project (GIMP), that combines ASTER (Advanced Space-borne Thermal Emission and Reflection Radiometer), SPOT-5 (Système Pour l’Observation de la Terre) and AVHRR (Advanced Very High Resolution Radiometer) photoclinometry (personal communication from I. Howat, 2012).

Mass conservation (MC)

The mass conservation (MC) method combines ice thickness data with ice velocity to infer a spatially consistent map of ice thickness that fulfills the continuity equation. This method was originally introduced by Rasmussen (1988) and applied to Columbia Glacier, Alaska, USA.

The MC method relies on the depth-integrated continuity equation. Ice thickness, H, must satisfy the following equation:


where is the depth-averaged ice velocity, and is the apparent mass balance, i.e. the sum of surface mass balance, basal melting and thinning rate (Morlighem and others, 2011). Ω is the model domain and Γ the inflow boundary, where measurements of ice thickness, H obs, are required and imposed. The bed elevation is deduced by subtracting H from a surface elevation dataset.

This standard approach yields products that may deviate significantly from the measurements due to errors in the input data (Morlighem and others, 2011). The only ice thickness data that constrain the computation of ice thickness in this method are along the inflow boundary Γ.

Here we employ an improved method, which uses all available thickness measurements within the model domain instead of along the boundary only, and the method takes into account errors in the input data, i.e. ice velocity, apparent mass balance, and ice thickness. We formulate an optimization problem by defining the following objective function:


where T ∊ Ω are the flight tracks where ice thickness, H obs, is measured, and γ is a regularization parameter. The first term of Eqn (4) represents the mismatch between modeled and measured thickness. The second term is a regularizing constraint, which penalizes wiggles in ice thickness. The addition of regularization is essential to infer smooth spatial variations in ice thickness. A large value of γ yields a smoother thickness map that deviates more from observations, whereas a small value of γ produces a good fit with observations but with strong gradients near flight tracks. γ is here kept constant for the entire domain, but it could be spatially variable to better model regions of known high thickness gradients; this aspect ought to be investigated in more detail in the future. We select γ = 105 using an L-curve analysis (Hansen, 2001), which optimizes the trade-off between the two terms that should be controlled.

The depth-averaged ice velocity, , is not calculated from an ice-sheet model but is assumed to be equal to the surface velocity derived from satellite radar interferometry data (Rignot and Mouginot, 2012). The absence of a model is justified for two reasons: (1) Relying on an ice-sheet model would introduce a circularity in the problem since the modeled velocity depends on ice thickness, and MC thickness is constrained by ice velocity (Eqn (3)). (2) Because the MC method relies on a transport equation (Eqn (3)), the solution for H is sensitive to errors in flow direction. Even with the additional support of inversion methods (e.g. MacAyeal, 1993; Morlighem and others, 2010), modeled flow directions are generally less precise than direct measurements, so using a model would actually increase the misfit between observations and calculated thicknesses.

The tolerance interval for depth-averaged velocities, , is derived from measurement errors and from the assumption that surface velocities are uniformly representative of column velocities, which is a reasonable assumption for fast-flowing glaciers but may not be reasonable for slow-moving ice (e.g. around domes or ice divides). We checked this assumption with a higher-order ice-sheet model (Blatter, 1995; Pattyn, 2003). The results show that surface speeds are <3% larger in magnitude than depth-averaged speeds over the entire domain Ω. To allow enough flexibility in the optimization process, we use a tolerance on ice velocity that is larger than the nominal error. We use an asymmetric error for : the upper bound is 15 m a−1 above the observed surface velocity; the lower bound is 5% smaller than the observed surface velocity.

For the apparent mass balance, we assume that the rates of basal melting and surface lowering of Russell Glacier are negligible compared to its average surface mass balance. Surface mass balance from a regional atmospheric climate model RACMO2 (Ettema and others, 2009) ranges from −250 cm a−1 at the ice margin to −50 cm a−1 in the interior of the study domain for the time period 1961–90 (Ettema and others, 2009). The tolerance interval for the apparent mass balance accounts for errors in basal melting rates, surface mass balance and ice-thinning rates. Errors in surface mass balance are at the 10% level (∼20 cm a−1). Errors in melting rates are small (<1 cm a−1). Satellite laser altimetry (Ice, Cloud and land Elevation Satellite (ICESat) and Airborne Topographic Mapper (ATM); personal communication from B. Csatho, 2012) shows that thinning rates reach 80 cm a−1 in this region. We therefore consider a tolerance interval of ±1 m a−1 for apparent mass balance.

Here MC is more sensitive to errors in velocity directions than errors in apparent mass balance. The ice velocity controls almost entirely the mass flux at a given location of the model domain compared to the apparent mass balance, which is significantly smaller than the mass flux from the ice coming upstream.

The spatial resolution of MC is constrained by the spatial resolution of the ice velocity data, here 350 m. We use a finite-element mesh of 25 000 elements ∼400 m in size. We stop the optimization procedure when the objective function, , reaches a value that corresponds to a misfit of <30 m between modeled and measured thicknesses along the flight tracks T, which is consistent with the analysis of crossover errors quoted earlier.


The bed topography derived from the MC method ranges from +700 m for the highest points to 300 m below sea level in troughs (Fig. 1). Bed troughs strongly coincide with areas of fast flow in Figure 2, i.e. ice is channelized along topographic valleys. Conversely, ice flow is slower and thinner on the surrounding topographic ridges and peaks. These features are consistent in location, orientation, size and elevation with valleys and ridges present in the digital topography of ice-free areas, which is available at a higher resolution (90 m) with a higher vertical precision (±10 m). The two independent data sources form an almost seamless reconstruction of the bed across the ice margin, except for a band ∼500 m wide, where the MC method does not apply because the ice velocity drops to values that are too low (<20 m a−1)to be applicable. The remarkable consistency of the two maps provides an independent evaluation of the quality of the MC results.

The standard error in bed elevation from the MC method is 36 m with respect to RT (Table 1) and 34 m with respect to PM. The KR bed topography yields similar results, when compared to both RT and PM. In the bed elevation from Bamber and others (2001), which is commonly used in ice-sheet models, most spatial features present in the KR and MC maps are absent (Fig. 3a). The updated version of this dataset (Bamber and others, 2013, fig. 3b) is more consistent with OIB measurements, but errors in bed topography still exceed 50 m when compared to RT and PM.

Table 1. Comparison of bed elevation obtained using ordinary kriging (KR) and the mass continuity (MC) methods versus the radar tomography (RT) results over a limited domain and versus the original measurements in profiling mode (PM) over the entire domain, for different flight-track spacings: inflow boundary only, 5 km, 2.5 km and 500 m. For reference, the standard error (σH ) between PM and RT is 46 m (right column). The bottom rows give the standard error in flux divergence (σ FD) over the entire model domain, and the maximum flux divergence error (ΔFDmax)

The RT bed elevation map (Fig. 3e and f) is of higher vertical precision and higher resolution than the other maps but it is limited in spatial extent because the ensemble of overlapping RT flight tracks is limited in size and also because the RT results have data gaps in the deeper parts of the troughs. These data gaps are associated with the SNR of the data that exceeds a given threshold selected during RT processing. We choose a conservative threshold here to minimize errors so that the method yields a product relatively free of artifacts, which we use in turn as a reference. With a more aggressive selection of the threshold, the processing would recover additional points in the trough, at the expense of data noise.

Comparing RT with other methods along profile A–B (Fig. 4), we observe significant differences between PM and RT, which we attribute to errors from off-nadir reflections in the PM products. The MC and KR methods yield smoother, lower-resolution results.

Fig. 4. Comparison of bed elevations along profile A–B in Figure 3e, overlaid on a radar echogram (frame 02_20110413_03_024). The bed elevations of mass continuity (MC) and ordinary kriging (KR) are smoother than the profile-model processing (PM) and radar tomography (RT) due to their resolution of 400 m.

The flux divergence is another important quantity to examine during the mapping of ice thickness because ice thickness data are typically combined with ice motion data in ice-sheet models to calculate the ice thickening/thinning rate. The mass conservation equation (Eqn (3)) indicates that the flux divergence should be equal to the apparent mass balance when the glacier is in steady state. The flux divergence should therefore be smaller than ∼1 m a−1. Large values of the initial flux divergence (several m a−1) in a numerical ice-sheet model would force the model to redistribute mass in order to reconcile topography and ice motion, thereby changing the glacier surface topography, and converge to a steady state that may differ significantly from the initial condition (Seroussi and others, 2011). When combined with ice velocity, the KR map yields errors in flux divergence up to 360 m a−1, which is large. Flux divergence errors for Bamber and others’ (2001, 2013) maps exceed 2000 m a−1, probably because of their lower spatial resolution (5 km and 1 km respectively). In contrast, MC yields flux divergence errors smaller than 1 m a−1, which is the selected tolerance interval for the apparent mass balance.

In a second set of experiments, we reduce the number of flight lines used in the mapping to determine how the performance of KR and MC methods evolves when the track spacing increases. We first consider only the measurements that are along the model inflow boundary (Fig. 5c), and run the two methods of KR and MC, using the same numerical settings. The KR method produces a nearly flat bed with no valleys (Fig. 5a), and large errors (Table 1), whereas the MC method yields correct spatial features (Fig. 5b), with a standard error of 110 m. We then use only one-tenth of the flight lines, hence simulating a track spacing of ∼5 km (Fig. 5f), and finally one-fifth of the flight lines, for a track spacing of 2.5 km (Fig. 5i). The results are shown in Figure 5d and g for KR, and Figure 5e and h for MC. Both maps are closer to the observations, and the troughs are more visible as the track spacing decreases. The bottom row of Figure 5 shows the maps obtained when all flight tracks are used, which corresponds to a track spacing of ∼500 m (Fig. 5l).

Fig. 5. Ice thickness of Russell Glacier for different sets of radar observations. The first column (a, d, g, j) corresponds to KR maps of ice thickness; the second column (b, e, h, k) corresponds to MC maps; and the last column (c, f, i, l) shows the corresponding radar tracks (PM) used as observational constraints in each row. In the first row (a–c), we only use measured ice thickness at the model inflow boundary. In the next rows, we use, respectively, one-tenth of, one-fifth of and all radar tracks, corresponding to track spacings of, respectively, 5 km, 2.5 km and 500 m.

When the track spacing increases from 500 m to 5 km, i.e. by a factor 10, the error in MC increases from 35 m to 65 m. The error for KR increases from 35 m to 80 m. The error in the KR method degrades more rapidly when the track spacing is increased than for the MC method, but the error in both methods does not decrease linearly with the track spacing. For flux divergence, with only one track at the inflow boundary, the KR error exceeds 900 m a−1. At 5 km and 2.5 km spacing, the error is still large. At 500 m spacing, i.e. when all tracks are employed, the error drops to 360 m a−1, which is still too large for direct application of ice-sheet models. This problem does not occur with MC, which, by design, provides errors in flux divergence smaller than 1 m a−1 for all experiments.

The estimated maximum potential error in MC bed elevation is shown in Figure 6. The error is low along flight tracks, where the results are constrained to match the input data, and it increases along flow. Maximum errors occur between flight tracks, as expected. The error increases as the track spacing increases. Errors are <50 m over most of the model domain, but exceed 100 m on the sides of the domain where few measurements are available and ice speed is low.

Fig. 6. Estimated maximum potential error in ice thickness derived from the MC method


The bed topography of Russell Glacier is more complex than anticipated given the relatively smooth nature of its surface elevation and the absence of major spatial variations in speed. The results, however, reveal that the flow pattern of the ice and the shape of the bed are strongly related. Ice is preferentially channelized along deeply carved, narrow valleys, which in this case are below sea level. These channels are likely to result from repeated glacial erosion over glacial cycles (Cook and Swift, 2012). When the track spacing is 500 m, KR and MC both yield products with low standard error, but MC maintains much lower errors in flux divergence. At higher track spacing, KR yields smoother products that do not preserve continuous features with preferential orientation due to its isotropic nature, hence the presence of bed channels is less well elucidated than with the MC method and artifacts in elevation are introduced along the bed channels.

Durand and others (2011) suggest that small-scale bed topography features of the order of 1 km in size significantly affect grounding line retreat of tidewater glaciers, and hence recommend digital maps of bed topography at a spatial resolution better than 1 km. Our analysis suggests that a precision of 50 m may be difficult to achieve with the KR method, but possible with the MC method. This result is of practical importance because it is difficult and costly to fly large areas at a 500 m track spacing. In contrast, a 2.5 km track spacing along the coastline has already been achieved along most of central west and northwest Greenland in 2011 and 2012.

Our estimated maximum error in bed elevation for the MC method is conservative because most of the OIB tracks flown on Russell Glacier are oriented in the direction of flow, which is the least favorable configuration for constraining mass transport. A more favorable flight pattern would be to orient flight tracks perpendicular to the flow direction because this constrains input and output mass fluxes better and minimizes the impact of errors in flow direction. With a more optimal configuration of flight tracks, the MC method should provide estimates of bed elevation – and ice thickness – with a precision of 50 m or better at 400 m resolution using a flight-track spacing of only 5 km according to the results of our study.

For higher-resolution mapping, i.e. 10 m spacing or similar to the RT method, higher-resolution ice velocity products are needed, which are achievable, for instance, using TERRASAR-X data (higher spatial resolution). Velocity mapping may also be improved by using the interferometric phase instead of speckle tracking (Joughin and others, 2010). Higher-precision products derived from the interferometric phase would make it practical to extend the MC method beyond the range of fast-flow areas. At present, the MC method is best applied to fast-flow areas surveyed with a track spacing of ∼5 km or better. This is also the portion of Greenland most relevant to high-resolution ice-sheet models projecting the evolution of the ice sheets into the coming decades to centuries.


The detailed study of the bed topography beneath Russell Glacier reveals that radar profilers (PM) yield errors in bed elevation four times larger than commonly thought. Furthermore, the traditional method of KR, with track spacings in the km range, yields products that do not capture important spatial features well if the track spacing exceeds ∼1 km and yield large errors in flux divergence when combined with ice motion data that are untenable for ice-sheet models. The reconstruction of bed topography using a combination of radar data and ice motion data offers significant advantages along coastal sectors. No model assumption is necessary other than neglecting internal shear. Areas where this assumption may be called into question are mostly slow-moving (speed less than ∼50 m a−1). In the case of Russell Glacier, a land-terminating glacier in central-west Greenland, the application of the MC method reveals topographic features that are intimately responsible for the flow structure of the ice toward the ice margin. These features are robust, and consistent with spatial structures found on adjacent ice-free areas. The application of the MC method to the entire ice-sheet periphery therefore promises to be transformative of our knowledge of basal topography beneath glacier ice, which is, in turn, critical for the numerical modeling of the evolution of glaciers and ice sheets.


This work was performed at the Department of Earth System Science, University of California–Irvine, and at the Jet Propulsion Laboratory (JPL), California Institute of Technology, under a contract with the NASA Cryospheric Sciences Program, grant NNX12AB86G. We acknowledge the use of data and/or data products from CReSIS generated with support from US National Science Foundation grant ANT-0424589 and NASA grant NNX10AT68G. We also acknowledge the use of the DEM from GIMP available online. Hélène Seroussi was supported by an appointment to the NASA Postdoctoral Program at the JPL, administered by Oak Ridge Associated Universities through a contract with NASA. The authors thank Neil Ross, Ute Herzfeld and an anonymous reviewer for helpful and insightful comments.


Bamber, JL, Layberry, RL and Gogineni, SP (2001) A new ice thickness and bed data set for the Greenland ice sheet. 1. Measurement, data reduction, and errors. J. Geophys. Res., 106(D24), 33 77333 780 (doi: 10.1029/2001JD900054)
Bamber, JL and 10 others (2013) A new bed elevation dataset for Greenland. Cryosphere, 7(2), 499510 (doi: 10.5194/tc-7-499-2013)
Blatter, H (1995) Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients. J. Glaciol., 41(138), 333344
Cook, SJ and Swift, DA (2012) Subglacial basins: their origin and importance in glacial systems and landscapes. Earth-Sci. Rev., 115(4), 332372
Deutsch, CV and Journel, AG (1998) GSLIB geostatistical software library and user’s guide. Oxford University Press, New York
Durand, G, Gagliardini, O, Favier, L, Zwinger, T and Le Meur, E (2011) Impact of bedrock description on modeling ice sheet dynamics. Geophys. Res. Lett., 38(20), L20501 (doi: 10.1029/2011GL048892)
Ettema, J and 6 others (2009) Higher surface mass balance of the Greenland ice sheet revealed by high-resolution climate modelling. Geophys. Res. Lett., 36(12), L12501 (doi: 10.1029/2009GL038110)
Fretwell, P and 59 others (2013) Bedmap2: improved ice bed, surface and thickness datasets for Antarctica. Cryosphere, 7(1), 375393 (doi: 10.5194/tc-7-375-2013)
Gogineni, S and 9 others (2001) Coherent radar ice thickness measurements over the Greenland ice sheet. J. Geophys. Res., 106(D24), 33 76133 772 (doi: 10.1029/2001JD900183)
Griggs, JA and Bamber, JL (2011) Antarctic ice-shelf thickness from satellite radar altimetry. J. Glaciol., 57(203), 485498 (doi: 10.3189/002214311796905659)
Hansen, PC (2001) The L-curve and its use in the numerical treatment on inverse problems. In Johnston, PR ed. Computational inverse problems in electrocardiography. (Advances in Computational Bioengineering 5) WIT Press, Southampton, 119142
Herzfeld, UC, Wallin, BF, Leuschen, CJ and Plummer, J (2011) An algorithm for generalizing topography to grids while preserving subscale morphologic characteristics – creating a glacier bed DEM for Jakobshavn trough as low-resolution input for dynamic ice-sheet models. Comput. Geosci., 37(11), 17931801 (doi: 10.1016/j.cageo.2011.02.021)
Joughin, I, Smith, BE and Abdalati, W (2010) Glaciological advances made with interferometric synthetic aperture radar. J. Glaciol., 56(200), 10261042 (doi: 10.3189/002214311796406158)
Leuschen, C, Gogineni, P, Rodriguez, F, Paden, J and Allen, C (2011) IceBridge MCoRDS L2 Ice Thickness. National Snow and Ice Data Center, Boulder, CO.
MacAyeal, DR (1993) A tutorial on the use of control methods in ice-sheet modeling. J. Glaciol., 39(131), 9198
Morlighem, M, Rignot, E, Seroussi, H, Larour, E, Ben Dhia, H and Aubry, D (2010) Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica. Geophys. Res. Lett., 37(14), L14502 (doi: 10.1029/2010GL043853)
Morlighem, M, Rignot, E, Seroussi, H, Larour, E, Ben Dhia, H and Aubry, D (2011) A mass conservation approach for mapping glacier ice thickness. Geophys. Res. Lett., 38(19), L19503 (doi: 10.1029/2011GL048659)
Paden, J, Akins, T, Dunson, D, Allen, C and Gogineni, S (2010) Ice-sheet bed 3-D tomography. J. Glaciol., 56(195), 311 (doi: 10.3189/002214310791190811)
Pattyn, F (2003) A new three-dimensional higher-order thermomechanical ice-sheet model: basic sensitivity, ice stream development, and ice flow across subglacial lakes. J. Geophys. Res., 108(B8), 2382 (doi: 10.1029/2002JB002329)
Pritchard, HD, Arthern, RJ, Vaughan, DG and Edwards, LA (2009) Extensive dynamic thinning on the margins of the Greenland and Antarctic ice sheets. Nature, 461(7266), 971975 (doi: 10.1038/nature08471)
Rasmussen, LA (1988) Bed topography and mass-balance distribution of Columbia Glacier, Alaska, U.S.A., determined from sequential aerial photography. J. Glaciol., 34(117), 208216
Rignot, E and Mouginot, J (2012) Ice flow in Greenland for the International Polar Year 2008–2009. Geophys. Res. Lett., 39(11), L11501 (doi: 10.1029/2012GL051634)
Seroussi, H and 6 others (2011) Ice flux divergence anomalies on 79north Glacier, Greenland. Geophys. Res. Lett., 38(9), L09501 (doi: 10.1029/2011GL047338)
Wu, X, Jezek, KC, Rodriguez, E, Gogineni, S, Rodriguez-Morales, F and Freeman, A (2011) Ice sheet bed mapping with airborne SAR tomography. IEEE Trans. Geosci. Remote Sens., 49(10, Pt 1), 37913802 (doi: 10.1109/TGRS.2011.2132802)