1. Introduction
There are fundamental glaciological measures necessary to understand the current status and to predict the future behavior of Earth’s ice sheets. These measures are surface velocity field, surface accumulation rate field, surface topography, basal topography (or ice thickness) and internal temperature. Together they form the basis for estimating glacier mass balance and dynamics through the mass continuity and force-balance equations.
Of these, the surface velocity field is well measured with synthetic aperture radar (SAR) interferometry (Reference JezekJezek, 2003, Reference Jezek2008; Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others, 2010; Reference Rignot, Mouginot and ScheuchlRignot and others, 2011). Surface topography is captured in great detail with airborne (Reference Hofton, Blair, Luthcke and RabineHofton and others, 2008) and satellite altimeters (Reference Schutz, Zwally, Shuman, Hancock and DiMarzioSchutz and others, 2005) as well as from stereo-photogrammetry, which again can be acquired from spacecraft or aircraft (Reference Howat, Smith, Joughin and ScambosHowat and others, 2008). Surface accumulation rate is more difficult to measure remotely with high spatial and temporal accuracy (Reference Munk, Jezek, Forster and GogineniMunk and others, 2003; Reference Arthern, Winebrenner and VaughanArthern and others, 2006), but modeling work is closing that gap in our knowledge (Reference EttemaEttema and others, 2009). Temperature data at depth remain an elusive goal using any technique other than direct measurements in boreholes. Remaining in the list is basal topography, or alternatively ice thickness, which is the subject of this paper.
Ice thickness measurements have been carried out in situ on glaciers and ice sheets using gravity (Reference Bull and HardyBull and Hardy, 1956), radar (Reference Waite and SchmidtWaite and Schmidt, 1962; Reference WaiteWaite, 1966; Reference Bogorodsky, Bentley and GudmandsenBogorodsky and others, 1985) and seismic techniques (Reference SorgeSorge, 1933; Reference Peters, Anandakrishnan, Holland, Horgan, Blankenship and VoigtPeters and others, 2008). Both radar and gravity measurements (Reference Peters, Blankenship, Carter, Kempf, Young and HoltPeters and others, 2007; Reference Diehl, Holt, Blankenship, Young, Jordan and FerraccioliDiehl and others, 2008; Reference Humbert and SteinhageHumbert and Steinhage, 2011) have been completed routinely across large portions of ice sheets using airborne instruments and also as part of NASA’s Operation IceBridge (OIB) program (Reference Koenig, Martin, Studinger and SonntagKoenig and others, 2010). In particular, technologically sophisticated radar systems developed by the University of Kansas (Reference Gogineni, Chuah, Allen, Jezek and MooreGogineni and others, 1998, Reference Gogineni2001; Rodriguez-Morales and others, in press) have been mounted in IceBridge aircraft and these instruments have successfully sounded extensive sectors of the Greenland and Antarctic ice sheets (http://www.nasa.gov/mission_pages/icebridge/).
Up until recently, these instruments have been operated optimally as nadir depth sounders, which is a desirable configuration for minimizing off-nadir clutter signals and increasing the sensitivity to weak basal returns. In 2003, a project was begun to investigate the potential for swath mapping of the ice-sheet base. Several approaches were identified, including radar interferometry (Reference JezekJezek and others, 2011a), but analyses showed that radar tomography was best suited to obtaining robust, high-quality elevation data for the surface and base of the ice sheet as well as intensity images of both surfaces (Reference Jezek, Gogineni, Wu, Rodriguez, Rodriguez and FreemanJezek and others, 2009; Reference Paden, Akins, Dunson, Allen and GogineniPaden and others, 2010; Reference Wu, Jezek, Rodriguez, Gogineni, Rodriguez-Morales and FreemanWu and others, 2011). The data were demonstrated to be useful for understanding glacier dynamics (e.g. separating the effects of mass flux imbalance and melting rates on changes in ice-sheet elevation), as well as for investigating the glacial geomorphology of the modern ice-sheet bed (Reference JezekJezek and others, 2011b).
In this paper, we compare and contrast radar data processed using advanced nadir sounding approaches with radar tomography. The data were collected over Isunnguata Sermia glacier in 2011 as part of OIB. Isunnguata Sermia is located in southwestern Greenland. It was selected by the OIB science team as a focus area for studying the glacier dynamics of a land-terminating glacier (www.bprc.osu.edu/ rsl/IST). As shown here, the two techniques provide powerful tools for characterizing the glacier base and for investigating the scale of structures that control local ice dynamics. We go on to integrate the radar topographies with other OIB datasets to investigate the physical properties of Isunnguata Sermia.
2. Data Collection
Isunnguata Sermia is located on the southwestern flank of the Greenland ice sheet (Fig. 1). The relatively small glacier extends into a now ice-free valley forward of the main ice-sheet margin. The glacier terminates on land and so is less susceptible to rapid and large changes now experienced by marine-terminating glaciers such as Jakobshavn Isbræ.
Data used in this paper were collected by OIB on 13 and 16 April 2011 and along flight tracks spaced ∼500 m apart. The aircraft position was controlled automatically using onboard GPS, and deviations from the planned flight path were usually <50 m. The University of Kansas Multichannel Coherent Radar Depth Sounder (MCoRDS) and the Wallops Flight Facility Airborne Topographic Mapper (ATM) were used to obtain ice thickness and surface elevation data. Additionally, a Sander Geophysics airborne gravity meter was used to obtain free-air gravity anomalies.
3. Radar System
The MCoRDS radar system was operated on each flight segment. The primary MCoRDS instrument parameters are listed in Table 1 (Reference LeiLei and others, 2010).
The radar system was equipped with 15 colinear dipole antenna elements. Seven elements were mounted under the fuselage of the P3 aircraft and four elements are mounted under each wing (Fig. 2). The middle seven antenna elements are used for both transmitting and receiving. The eight side elements are used for receiving only. For the seven transmit/receive channels, the receivers can be opened only after the radar waveforms are transmitted using a transmit/ receive switch.
In order to achieve a dynamic range of over 100 dB in the receiver, two radar waveforms are transmitted. The first waveform has a 1 μs chirp duration and is used to probe the ice-sheet near-surface with a small number of hardware integrations. The second waveform has 10 μs chirp duration, uses more integrations and is designed to capture data from deep or warm ice. When the echoes of these two waveforms are received, a lower gain is applied to the surface waveform echo because the surface echo is strong and a higher gain is applied to the base waveform since the bottom echo is usually weak. In this way a dynamic range of 100 dB can be achieved in the receiver with 14 bit analogue-to-digital converter units.
4. Calibration
For a multichannel system like MCoRDS, phase and amplitude calibration are very important if either beam-focusing or tomography are to be achieved successfully. Data calibration includes interchannel radiometric calibration, antenna geometric calibration and interchannel phase calibration.
In 2011, OIB collected calibration data over sea-ice areas with flight elevation ∼5 km above the sea-ice surface. The flight occurred on 9 April when the sea ice was relatively smooth and provided a good reference for phase calibration. The phase difference between any two channels is normally dependent on the environmental temperature of the radar. However, because MCoRDS is operated at low-frequency VHF, the interchannel phase is a slow function of temperature. After calibration, the swath tomographic method (described below) is used to estimate the sea-ice surface elevation and the image magnitude. Prior to tomographic processing, signals from the left and right sides of the aircraft appear without discrimination (Fig. 3a). After processing, the nadir return is centralized in the image and the left–right side echoes are well separated (Fig. 3b).
5. Nadir Sounding Ice Thickness
The two primary challenges in nadir sounding of ice sheets are obtaining sufficient signal through the attenuating ice sheet and reducing surface clutter and side echoes, which can mask weak basal returns. MCoRDS is designed to overcome these challenges by using multiple antennas that increase the power density on the target and mute side echoes and surface clutter by electronic beam formation and steering. As part of OIB, the University of Kansas processes MCoRDS data and provides the data to the US National Snow and Ice Data Center (NSIDC). Processing procedures are described by Reference Leuschen and AllenLeuschen and Allen (2010) and summarized here for better understanding the comparison between nadir sounding and tomographic results.
Calibration data described above are used to apply channel compensations between each of the antenna phase centers (Reference Leuschen and AllenLeuschen and Allen, 2010). This includes time delay, amplitude and phase mismatches. Range resolution is achieved using pulse compression. Motion compensations for attitude and trajectory lever arm corrections are applied and then the data are SAR processed with an along-track spatial frequency window using frequency–wavenumber (F-K) migration. Data from multiple antennas are combined to enhance the basal return. Two methods are used to do this combination: (1) a simple beam-forming method, which sums radar profiles acquired from each antenna and is analogous to stacking traces; and (2) a minimum variance distortionless response (MVDR) method, which is described by Reference LiLi and others (2013). Both algorithms use multi-looking or spatial incoherent power averaging of a neighborhood of pixels surrounding the image pixel being combined followed by along-track decimation. Lastly, echograms from low and high gain channels are combined to form a single image.
The resulting data can be displayed as individual amplitude waveforms or as intensity modulated displays. The former are useful for carefully picking the travel time of an echo, whereas the latter provide context for manual interpretation of the data. The layer tracking of ice surface and ice bottom reflections are manually driven processes with basic tools for partial automation. The tools used are determined by the operator picking the data and include manual picking and interpolation, tracking the strongest return within a window, and tracking the leading edge of the strongest return. Various post-image processing methods (e.g. multi-looking and detrending) are used to better highlight features in the echogram as needed.
The radar was configured for nadir sounding observations during these experiments. Consequently near-nadir returns were emphasized, whereas off-nadir returns were attenuated. This improved the quality of along-track acquisitions but limited application of swath mode techniques. Nevertheless, as illustrated later, both data-processing approaches yielded useful results.
Operation IceBridge nadir measurements of ice thickness were downloaded from the NSIDC. These were interpolated to an ice thickness grid using kriging algorithms implemented in ArcMAP (Fig. 4). Five neighboring points selected from four equal quadrants were chosen to calculate a prediction map using ordinary kriging. Inspection of the selected points from within the data collection area showed that points were distributed both across and along track. The semivariogram was modeled with a spherical function computed using 100 lags of 100 m each. The prediction map was subsequently converted into a 50 m raster grid. The boundary of the calculated raster was masked to remove obvious interpolation distortions.
6. Radar-Tomography Ice Thickness
Nadir ice-sounding provides one-dimensional thickness measurements of the ice sheets along the flight lines of the radar sounder. The along-track platform position is obtained by conventional GPS; radar timing provides the range measurement. However, at least four locations share these same geospatial characteristics: two points (left/right) each on the basal and surface layers will have the same range and cross-track azimuth. Discrimination among these points requires at least four additional independent measurements. MCoRDS provides these measurements by implementing a three-antenna array with 15 total dipole elements. By using phase-coherent tomographic (or angle-of-arrival) processing, one can discriminate among all ambiguities (basal and surface clutter) to obtain topographic maps of both the base and ice surface, together with radar reflectivity maps.
The formal basis for ice-sheet radar tomography and our implementation of the technique is described by Reference Wu, Jezek, Rodriguez, Gogineni, Rodriguez-Morales and FreemanWu and others (2011). Essentially, we formulate the tomographic basal-ice imaging task as a problem of estimating signal arrival angles. The ice mass has two major interfaces: the upper surface interface, between the air and the ice mass, and the basal interface, between the ice mass and bedrock or basal water. Between the interfaces there are internal layers that originate from slight density changes or from changes in chemical impurities. Because the reflection of internal layers is highly specular, only their close nadir returns are significant but arrive at the receiver earlier than the bottom returns. Therefore the echoes from these internal layers are ignored in our analysis. When the radar signal reaches the air–ice interface, some of the energy is scattered back and some refracts through the interface and continues traveling through the ice mass. Some of the transmitted signal will scatter back towards the receiver from the basal interface. Simultaneous with the nadir echo arriving from the bedrock, the receiver may also detect the surface reflected signal from both the left and right sides of the sensor. Off nadir, there may also be simultaneous echoes from the base originating from the opposite sides of the airplane. In addition, there will be thermal and other noise sources. Radar tomography allows for the separation of these various returns given a sufficient number of independent receiver channels.
In earlier work (Reference Wu, Jezek, Rodriguez, Gogineni, Rodriguez-Morales and FreemanWu and others, 2011), we found that tomographic estimates of ice-sheet surface elevation compared to within 5 m rms of Ice, Cloud and land Elevation Satellite (ICESat) elevation data. We also compared tomo-graphic data with independent nadir-radar measurements acquired orthogonal to the tomography flight directions. The estimates compared favorably to within 14–18 m rms. We noted then that the tomographic data contained artifacts near nadir. This was because the earlier algorithm sought to discriminate between left- and right-side arrivals even given, as is the case at nadir, that there is only one signal. For the data presented here, the processing algorithm was modified to seek only one arrival near nadir. This eliminates the artifacts and, as discussed later, improves the comparison between the coincident tomographic and nadir data to a comparable level.
There was a total of 27 processed tracks of data collected over Isunnguata Sermia on 13 and 16 April 2011. The processing steps are range compression, time-domain azimuth compression and tomographic processing for each track. After the ice thickness is generated, a surface digital elevation model (DEM) compiled from OIB ATM data (discussed below) is used to generate the bottom topography. Finally all the tracks of ice thickness and bottom elevation measurements are mosaicked (Fig. 5). Gaps appear in the map for several reasons including low signal level, shadow and layover. Optimization for nadir sounding also means the off-nadir returns will be preferentially weakened. Nadir tracking, as shown below, fills in some of the gaps because manual tracking can often be accomplished with lower signal levels and discontinuous echoes can be interpolated given the overall context of the segments.
7. Subglacial Terrain
To construct models of subglacial relief, we subtracted the ice thicknesses from a DEM of the surface. We constructed the DEM (50 m post spacing) for the Isunnguata Sermia surface and its surrounds. ATM Qfit data (Reference KrabillKrabill, 2010) were downloaded from the NSIDC for flights that occurred on 13 and 16 April 2011. ATM is a scanning laser altimeter. To reduce the total data volume, only data acquired at the aft, left and right azimuth angles of the scan were selected. The data were converted to polar stereographic coordinates and then interpolated using the same method and parameters as described for the interpolated ice thickness data. The 50 m grid was masked to isolate data over the region of interest and remove obvious boundary-region artifacts (Fig. 6).
The tomographic ice thickness data were subtracted from the surface DEM to construct the basal topography. Similarly, the interpolated nadir-data-derived ice thicknesses were subtracted from the surface topography to construct a second estimate of the base (Fig. 7). The nadir data yield a more extensive estimate of the basal topography than the tomographic result, first because the nadir data are interpolated to a continuous grid, second because the nadir signal is stronger than the off-nadir returns and third because manual interpolation between weak echoes along the flight track can successfully bridge gaps. Both approaches reveal a terrain characterized by undulating highlands to either side of a deep channel that runs northeast–southwest into the main trunk of Isunnguata Sermia. Interpolated nadir data reveal a >600 m deep trough that channels flow northwest towards the westward terminus of the glacier. A 300 m deep channel, visible in both maps, strikes northeast before intercepting the main trunk. The channel intercepts another northwest–southeast-striking 150–200 m deep channel farther upstream.
Basal topography maps estimated from tomography and interpolation of nadir data were differenced in areas of overlap. The mean difference is 4 m with a standard error of 31 m. Profiles of the along-track nadir data are taken as the reference and selected based on continuity of the tomography map over the length of the profile (Fig. 8). The nadir data are compared with the interpolated nadir data and the tomography data along the three nadir sounding tracks (Fig. 9). For the central line and after interpolating data to a common geographic location, the point-by-point rms difference between the nadir data and the tomographic data is ∼16 m, similar to the measure reported by Reference Wu, Jezek, Rodriguez, Gogineni, Rodriguez-Morales and FreemanWu and others (2011). The tomography estimate disagrees by ∼50 m from the other two datasets near 18 km along the northerly line. Disagreement between the tomographic data and the nadir data is often due to errors in tracking the ice bottom in steep terrain. Specifically, when a strong off-nadir reflection occurs, the nadir sounding method may not perfectly reject this side echo and the analyst tracking the ice bottom may incorrectly interpret this as the nadir echo or it may be the only echo available. Tomography better handles this situation because it generally has finer cross-track resolution and properly resolves the side-looking and nadir echoes. In fact, strong side echoes will be better resolved because they are stronger and have a higher signal-to-noise ratio (SNR).
The three basal topography datasets contain spatial wavelengths at the 5–10 km scale with similar amplitudes (Fig. 9). Shorter wavelength (1–2 km) undulations along the flight direction are evident in the tomography and nadir data. Evidence in support of short wavelength topography comes from inspection of ATM-derived elevation data collected forward of the glacier snout (Fig. 9). Fourier transform of the conditioned (zero-meaned, Hanning-weighted and padded) rock surface profile shows that there are spectral peaks at 1.6, 2.7 and 3.6 km cycle−1. Detection of longer wavelengths is limited by the length of the record. Fourier transform of the conditioned tomographic profile along the central line also shows a strong spectral peak at 1.6 km cycle−1. There is a weaker peak at 2.1 km cycle−1 followed by peaks at 3.8 and 10.0 km cycle−1. There is a strong peak at 10 km cycle−1 in the measured nadir data followed by spectral peaks at 3.9, 2.1 and 1.6 km cycle−1; however, the shorter wavelength components are poorly defined. As expected, the interpolated data are smoothed relative to the tomography and nadir data. There remains a strong spectral peak at 1.5 km cycle−1, but the peak at ∼2 km cycle−1 is muted at best and at 3.8 km is nearly hidden in the spectral slope leading up to the expected peak at 10 km cycle−1.
As noted above, the tomographic and interpolated topography estimates are consistent to better than 10 m on average, though the interpolated surface is smoother than the tomographically derived surface. Agreement is a consequence of the exceptionally fine across-track spacing (∼500 m) achieved in this survey. Unlike the continuous surface created by interpolating the nadir data, the tomo-graphic surface is less extensive and contains small gaps, but the tomographic surface better retains short-wavelength roughness in the along-track direction and, as shown by Reference Wu, Jezek, Rodriguez, Gogineni, Rodriguez-Morales and FreemanWu and others (2011), in the cross-track direction as well. In summary and not unexpectedly, different methods of processing the data have advantages and disadvantages. So we rely on the strengths of each method to investigate glacier processes in Section 8.
8. Discussion
Isunnguata Sermia flows through a deep bedrock trough. The glacier terminates in a sediment-filled valley that carries away summer meltwater. A perspective view of Isunnguata Sermia is shown in Figure 10, which illustrates the different environs of the glacier. An August 2010 Landsat-7 image obtained from the US Geological Survey (USGS) is overlaid on the surface elevation model. Flow directions upstream on the glacier are noticeable as weak longitudinal flow stripes embedded within the center of the generally undulating surface. Surface lakes on the ice sheet appear as small dark spots, which correspond well with the positions of surface depressions. Lateral curvilinear bands in the center of the image are exposed layers within the ice in this region of strong surface ablation.
Details of the surface topography are shown in Figure 11. TerraSAR-X SAR interferometry data collected by the German Aerospace Center in support of OIB (personal communication from D. Floricioiu, 2011) in 2009 provide surface velocities (courtesy of I. Joughin). Velocities range from ∼20 m a−1 near the terminus to ∼120 m a−1 in the channel feeding Isunnguata Sermia. In the eastern part of the image, velocities bear to the west more or less orthogonal to the surface contours. In the south central part of the scene, the effect of the deep subglacial channel is seen to reorient the velocities to the northwest and then down towards the terminus of the glacier. Even slow-moving ice along the more stagnant margins flows in a direction orthogonal to the contours, adding confidence to the quality of the velocity dataset.
The velocity field is little influenced by basal topography other than the deepest trough feeding the terminus. This is illustrated in Figure 12 where velocity data are superimposed on a hill-shaded relief map of the surface elevation and on the basal topography, respectively. Many of the northwest- and northeast-striking lineaments that intersect in a checkerboard pattern noticeable in both the subglacial topography and the surface topography have little influence on the orientation of the surface velocity field.
Surface elevation and ice thickness data were resampled to 1 km grids to calculate surface slopes and driving stresses on the order of one ice thickness. Elevation data used here are relative to the ellipsoid. The 2 m geoid undulation, which trends east–west in this area, was ignored in the calculation. Driving stresses are defined in the usual way as
where ρ is the ice density, g is the acceleration due to gravity, H is the ice thickness, h is the surface elevation and s is the distance downslope along the surface (Reference Van der VeenVan der Veen, 1999). Driving stresses range from 10 to 240 kPa (Fig. 13). Driving stress is on average ∼104 kPa. Highest driving stresses (∼220 kPa) are generally associated with the thicker ice. Lowest driving stresses (10–20 kPa) are correlated with small surface lakes, presumably because slopes are low in the depressions where lakes are likely to form.
The relatively uniform orientation of the velocity field, at least in the eastern part of the study area, and the likelihood that lateral forces are small for this relatively slow-moving glacier suggest that a simple laminar flow model is appropriate for describing the gross behavior of the upstream area. The theoretical surface velocity (U s) for the case where basal drag dominates the resistive stresses is
where U b is the basal sliding velocity (Reference Van der VeenVan der Veen, 1999). Using the measured ice thickness, the driving stress as computed above, a warm value for A (2 × 10−17 kPa−3 a−1) and n = 3, the average of the first term is ∼60 m a−1 less than the measured surface velocity. The difference can be explained by attributing ∼50% or less of the surface velocity to basal sliding. Our estimated contribution of basal sliding is higher than that modeled by Reference Brinkerhoff, Meierbachtol, Johnson and HarperBrinkerhoff and others (2011). The acknowledged simplifications in our analysis and also the important differences between our measured basal topographies and those used by Reference Brinkerhoff, Meierbachtol, Johnson and HarperBrinkerhoff and others (2011) together may account for the discrepancy.
Subglacial terrain adjacent to the troughs consists of hilly terrain reminiscent of the deglaciated coastal regions of Greenland. To further compare the proglacial and subglacial terrain maps, we hill-shaded the surface and subglacial DEMs compiled from ATM-interpolated nadir data and radar tomographic data, respectively (Fig. 14). As quantitatively shown earlier using profile data, the tomographic basal topography map seems little different from that of the exposed rock forward of the ice front. The map data in Figure 14 extend this conclusion to two dimensions (2-D). Continuity of the subglacial geology across the glacier terminus is evidenced by free-air gravity anomalies (Reference Cochran and BellCochran and Bell, 2010) (Fig. 15). There are noticeable north–south anomaly gradients, but east–west anomaly gradients across the terminus are weak. Consistent with the topography data, strongest negative anomalies are associated with the glacier troughs and sediment-filled valleys. Similarity between the proglacial terrain and the tomographically derived sub-glacial topography suggests that at least at the scale of the measurements (50 m pixels), erosional processes acting on the bed of the modern ice sheet are largely unchanged from those that operated on the now deglaciated area to the west.
9. Summary
Nadir sounding and tomography processing techniques are both useful for investigating different geophysical aspects of the radar measurements. Nadir sounding data, improved by surface clutter and side echo reduction, can be interpolated to build spatially extensive grids of detailed subglacial topography including sounding of deep subglacial troughs. A limitation of this approach is that the nadir data must be interpolated to build the grid and interpolation inherently smooths the data. Moreover, multiple closely spaced flight lines are required to assemble a high-fidelity grid and only in rare instances, as here, is that logistical requirement achievable. Radar tomographic techniques operating on the same data can be used to directly sample in 2-D and retain high-frequency signals in the terrain. Flight-line spacing can also be increased, making tomography more practical. However tomography requires good SNRs to obtain reliable off-nadir estimates. From past experience (Reference Wu, Jezek, Rodriguez, Gogineni, Rodriguez-Morales and FreemanWu and others, 2011), an SNR better than 0 dB for each channel is desirable for tomography processing.
When applied in this study of Isunnguata Sermia, interpolated nadir data capture an intersecting network of deep subglacial troughs. Upstream of the main glacier trunk, the troughs have little influence on the direction of ice flow, which is governed in about equal measures by basal sliding and creep. At finer scales, the tomography data capture details of short-wavelength topography that compares favorably with the exposed terrain west of the glacier terminus. Consequently the processes currently shaping the bed of this relatively slow-moving land-terminating glacier are likely to be the same as those which once shaped the rocky terrain now forward of the glacier terminus.
Acknowledgements
This research was supported as part of NASA’s Operation IceBridge project and also by NASA’s Cryosphere Program. The University of Kansas also received support from the US National Science Foundation through the Center for Remote Sensing of Ice Sheets. This research was (partly) carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Many helpful suggestions were received from two anonymous reviewers.