The precision of radar-derived subglacial bed topography: a case study from Pine Island Glacier, Antarctica

Abstract

Recent advances in the measurement of bedforms beneath active ice streams have been made using ground-based grid profiling using impulse radar systems operating with centre frequencies in the 3–5 MHz range. Surveys of Rutford Ice Stream and Pine Island Glacier have shown that features such as mega-scale glacial lineations with topographic relief of as little as 3 m can be traced for many kilometres downstream under more than 2 km of fast-moving ice. In the discussion of these data, it is often asked ‘How is it possible to map such fine-scale topography with such a low-frequency radar’. In answering that question, the key point is the distinction between the precision of a radar range measurement to a single, isolated reflective interface and the ability to resolve the presence of two closely-spaced interfaces of similar reflectivity (commonly referred to as the vertical resolution). This paper will discuss and illustrate this distinction and use the case study of data acquired over Pine Island Glacier to examine the limits of precision of the radar range measurement.

Introduction

Impulsive-source ice-penetrating radar systems are used extensively in glaciology for the measurement of ice thickness and the determination of the internal structure of glaciers and ice sheets (Arcone and others, 1995; Welch and Jacobel, 2003, 2005; King and others, 2009, 2018; Hindmarsh and others, 2011; King, 2011; Kingslake and others, 2016; Lindback and others, 2018; Brisbourne and others, 2019). The systems in use are both commercially produced and academic ‘lab-built’ units, which all share the characteristic that they are classed as ultra-wideband (UWB) radars because the transmitters emit a short-duration voltage spike into the output antenna. In contrast, narrow-band radars emit a controlled burst of frequency-modulated signal which covers a particular frequency range. This discussion is restricted to UWB impulse radars.

The question which this paper addresses is how precisely can the dimensions of subglacial landscape features be determined using an impulsive radar? In this context, the distinction between precision and accuracy can be expressed as follows: precision is the detection limit of changes in bed topography while accuracy is a measure of how well the change in topography can be located in 3D space.

Changes in bed topography are detected by changes in the travel time of the radar pulse reflected from the ice/bed interface as the radar system traverses the profile. How precisely this measurement can be made depends on the characteristics of the radar system and the nature of the reflecting interface. If the ice/bed interface is a single, isolated, specular reflector, then the precision of the travel time measurement is related to how well the radar system resolves the reflected wavelet. Most text books and papers describing ground-penetrating radar state that the distance to an object in the sub-surface can be determined to no better than one-quarter of the wavelength of the centre frequency of the radar in use (Reynolds, 1997; Annan, 2009; Bristow, 2009). Here it will be argued that the commonly-accepted ‘quarter wavelength criteria of vertical resolution’ is the wrong criteria to use when considering the precision of mapping basal topography beneath thick ice and that the measurement using an impulse radar is much more precise than is commonly assumed. One of the few examples of avoiding this pitfall was shown by Conway and others (2009) where the uncertainty in ice thickness was attributed to uncertainties in the wave-speed and picking error. In that paper, the uncertainty in ice thickness was estimated as ±7 m when the quarter-wave resolution was 20 m.

The errors affecting ice-thickness measurements have been treated extensively in a suite of three papers by Martin-Espanol and others (2016), Lapazaran and others (2016a, 2016b). These authors use an even more conservative ‘half wavelength’ criterion for vertical resolution (Lapazaran and others, 2016a, Section 3.1.2) and state that any error in the sampling period can be neglected. However, the majority of the error budget detailed by Lapazaran and others affect the accuracy of the overall ice-thickness measurement but not the precision of the travel-time measurement. I will show that the ability of an ice-penetrating radar to detect a difference in the bed elevation from one point to the next is mainly determined by the digitizer sample period and not the dominant frequency of the radar system.

In the first section, the difference between radar range measurement and vertical resolution will be discussed. This will be followed by the examination of example data from a 3.8 MHz impulse radar acquired over an Antarctic ice stream which will show that subtle bedforms can be precisely mapped beneath thick ice.

Radar Range Measurement Versus Vertical Resolution

The term radar stands for radio detection and ranging and determining the range to a reflecting target is the primary function of any radar system. A separate function of a radar system is the ability to distinguish between closely-spaced targets, which is known as resolution. In this section, I will draw out the distinction between the precision of a radar range measurement and the ability to resolve two closely-spaced reflectors.

Radar range

Radar range is the distance between the radar system and a single reflecting target. It is derived by measuring the flight time of the leading edge of an electromagnetic wave from its origin at the transmitting antenna via a reflection at the target to its arrival at the receiving antenna (Fig. 1). For simplicity, it is assumed that the target is a specular reflector.

Fig. 1. (a) Illustrative trace (sometimes referred to as an A scope) for an impulse radar. The direct or air wave shows that the signal output from the transmitter antennas is a one-and-a-half cycle wavelet. The wave reflected from a positive conductivity contrast at depth beneath the radar retains the same polarity as the direct wave while the amplitude is reduced due to spherical spreading, attenuation over the ray path and the reflection coefficient of the interface. The shape of the reflected wave changes little compared to the transmitted wave because of the small frequency dispersion in ice. (b) Ray paths for the direct and reflected waves.

The criteria that should be considered which control the precision of the radar range measurement are as follows:

1. (1) The timing resolution and jitter of the radar receiver: this is largely controlled by the CPU clock cycle of the controlling computer, which in modern systems is in the GigaHertz range, two to three orders of magnitude greater than the bandwidth of the radar

2. (2) The digitization interval of the recording system: this varies from system to system. In the example data presented below, the digitizer recorded a sample at 10 ns intervals (i.e. a 100 MHz digitizer)

3. (3) The electromagnetic wave propagation speed of the medium between the radar and the target: this parameter determines the propagation time for a given ice thickness. The wave-speed profile is rarely known to a high level of accuracy but it can reasonably be assumed to vary slowly across a survey area. Adjacent radar range measurements along a profile are typically between 1 and 10 m apart. On this scale, the wave-speed profile can be considered invariant. Thus, the wave-speed has a small, generally negligible, effect on precision when mapping detailed bed topography but can have a variable effect on the accuracy on a regional scale.

4. (4) The rise time of the reflected signal: this is a measure of how rapidly the voltage rises at the receiving antenna when the reflection arrives. It is controlled by the rise time of the voltage output of the transmitter, modified by the bandpass filter effect of a dipole antenna and frequency dispersion over the travel path. In the example system used here, a bank of capacitors is connected to the output terminals of the transmitter by a fast switching circuit; as a result, the transmitter output voltage rises at 250 V ns⁻¹ and peak voltage is reached 10 ns after the onset time. This means that a detectable voltage increase takes place in a time short compared to the digitizer interval (Fig. 2).

5. (5) The signal-to-noise ratio (SNR): this determines how easily the first arrival of the reflected wave from the bed can be picked. In general, the SNR for ice-sheet impulse radar data is high because many individual radar soundings are stacked together to produce one trace in the profile. In the case study described below, 1000 soundings were stacked per trace, resulting in a >30-fold reduction in random noise. In practice, noise can never be fully eliminated, but if suppressed sufficiently, interpretation of the reflection onset can be made with confidence. Noise effects are considered further and discussed in the vertical resolution section.

Fig. 2. Oscilloscope trace showing the rise time of the positive side of the transmitter output under no-load condition. The voltage rise from 0 to 2500 V takes place in 10 ns.

Vertical resolution

Vertical resolution is the ability of a system to distinguish between two closely-spaced targets which have similar reflection strength. When reflecting discontinuities are close together, the wave reflected from the second target overprints the wave reflected from the first target (Fig. 3a). The two targets must be a minimum distance apart before the presence of the second target can be deduced. When the two targets are one-quarter of the dominant wavelength apart, an inflection in the combined waveform can be detected.

Fig. 3. (a) A radar reflector exists where there is a step change in dielectric properties. The strength of the reflection or reflectivity is proportional to the dielectric contrast across the interface. When a radar wave impinges on a reflecting interface, the amplitude of the reflected wave is proportional to the reflectivity. When two reflectors with similar contrasts in dielectric properties are located close together, the reflected waves become superimposed. Only when the separation between the reflectors is equal to or greater than one-quarter of the dominant wavelength is there an indication in the shape of the wavelet that two reflectors are present. (b) When there is a single reflector, the onset time of the reflected wave is unambiguous. (c) At an ice/bed interface, the contrast in dielectric properties is very much greater than adjacent intra-ice reflectors. Two 1D models are shown, representing an intra-ice reflection with a reflection coefficient equal to 5% of that of the bed. In the first model, the internal reflection lies a quarter wavelength above the bed, in the second model at 1.25 times the wavelength. The reflection from the bed is blue, from the internal ice horizon is red and the combined return is black. Random noise with a peak amplitude of 80% of the internal reflection is added to the combined return.

At depth in an ice sheet, radar reflections are produced by changes in conductivity which arise from chemical contamination of the ice (Moore, 1988; Eisen and others, 2006). Dielectric profiling of ice cores shows that conductivity changes of between two and five times the background level can occur, often associated with atmospheric fall-out from globally-significant volcanoes. When such conductivity changes occur close to each other in depth, they cannot be distinguished as separate events on a radargram unless more than a quarter wavelength apart because the reflection amplitudes are similar and the wavelets arising from the adjacent conductivity changes interfere (Fig. 3a).

In contrast, the conductivity of water-saturated sediment is between 10 and 100 times the ice value (Cassidy, 2009). This means that the reflection from the ice base of a thick ice sheet is typically 20 times stronger than any adjacent englacial reflection; as a result, the ice base can be considered to be a single reflecting target. This is because reflections from internal ice layers close to the bed have too low amplitude to change the shape of the wavelet reflected from the ice base in a way that makes the time of the reflection ambiguous. In Figure 3c, a simple 1D model of the bed reflection is shown with random noise equivalent to an SNR of 1:25 added to the output trace. The start of the reflected wavelet can still be picked to ± 1 sample. For this reason, the quarter wavelength criteria is the wrong measure of the limitation on the resolution of changes in ice thickness in impulsive source ice-penetrating radar surveys.

This argument applies to clean ice/sediment, ice/bedrock and ice/water interfaces. If there is a debris-rich basal ice layer, the interface may have a gradational change in electromagnetic properties rather than a step change. In this circumstance, the true bed of the glacier or ice sheet may not be detectable by radar.

Example Data from Pine Island Glacier

The example radar data were collected over the main trunk of Pine Island Glacier using the British Antarctic Survey's Deep Look Radio Echo Sounder (DELORES) (King and others, 2009). This is a monopulse radar that is based on designs developed by the University of Washington (Gades, 1998) and St Olaf College (Welch and Jacobel, 2003). The transmitter fired a ±2500 V pulse into the transmit antennas with a 1 kHz repetition frequency. The antennas were fully-damped resistively-loaded wires of 20 m half-dipole length which resulted in a centre frequency of ~3.8 MHz (Fig. 4b). The receiver used a chassis computer with a 100 MHz digitizing card. For profiling, a total of 1000 shots were averaged in the oscilloscope buffer to form each recorded trace in order to improve the SNR. A trigger pulse was transmitted between the transmitter and the receiving digitizer via a fibre-optic cable. The radar equipment was towed behind a snowmobile travelling at ~12 km h⁻¹; therefore, each trace is built up of data acquired over a horizontal distance of ~3 m. Traces were horizontally interpolated such that the final radargram has a trace spacing of 7.5 m, which is the same order of magnitude as the post-migration Fresnel zone (Lindsey, 1989).

Fig. 4. (a) Radargram from a survey of Pine Island Glacier, Antarctica. Trace spacing 7.5 m, total profile length 475 m km. The reflection from the bed of the ice stream is a high-amplitude single wavelet. Undulations in the bed topography of ⩾3 m can be measured. (b) Normalized frequency spectrum of the data in (a), the peak frequency is 3.8 MHz. (c) Wavelet shape.

The example profile (Fig. 4a) shows a 475 m section of the bed of the ice stream. The data were processed using bandpass filtering, a gain function to correct for spherical spreading loss and migration to collapse diffractions to their source points. The reflection from the ice/sediment interface is the prominent event across the centre of the image which has a much higher amplitude than the nearby internal ice reflections. This radargram shows that undulations in the bed topography of >3 m can be reliably detected with this radar system. The shape (or phase) of the wavelet is consistent from trace to trace which allows small changes in the wavelet arrival time to be identified (Fig. 4c). This was confirmed using a simple 2D model (Fig. 5). The model comprised a 3.5 m high bump in the ice/bed interface that was 80 m wide. Model traces at 5 m intervals were computed using the finite difference modelling package in the ReflexW™ (Sandmeier Geophysics) processing software. When the model profile was migrated, the height and width of the bed topographic feature were well reproduced. The model was run with both a 3.8 and a 2 MHz wavelet (Figs 5d and e). The same two-way travel times were picked on both profiles, illustrating that the radar range measurement is independent of the centre frequency used.

Fig. 5. Results of 2D modelling demonstrating that first arrival times are independent of frequency. (a) The model is of a bedrock ridge 3.5 m in height and 80 m in width, placed at a depth of 490 m beneath the ice. (b) On an unmigrated radar profile, the height of the ridge can be determined but not the width. (c) When the profile is migrated with a finite difference migration algorithm, both height and width can be picked from the profile. (d) Profile created with a wavelet with 3.8 MHz centre frequency. Red crosses are the picked zero crossings marking the first arrival (also known as first break or reflection onset). (e) Profile created with a wavelet with a 2 MHz centre frequency. The first arrivals (black crosses) are the same as for the 3.8 MHz wavelet profile.

Precision

In order to test the timing stability of the radar system and the consistency of the transmitted pulse, a static test was conducted by leaving the system running for several minutes while stationary over a section of bed with a clear, high amplitude reflection. The data (Fig. 6) illustrate the following:

1. (1) The transmitter pulse is very stable in shape and amplitude.

2. (2) The triggering process is stable (trigger pulse sent from the transmitter to the receiver over a fibre-optic cable).

3. (3) The digitizer timebase is stable.

Fig. 6. Static test data from Pine Island Glacier shows that the transmitted signal was consistent and that pulse-to-pulse jitter (4.4 ns) was significantly less than the digitization interval (10 ns).

Figure 6 is an extract from a longer record of 10 500 traces. Analysis of all traces showed the std dev. of the time of the positive peak of the waveform was 4.4 ns.

Therefore, the largest factor in defining the precision of the measurement of ice thickness in these data is the digitization interval of 10 ns. This equates to the two-way time through 0.84 m of ice (wave-speed = 0.168 m ns⁻¹). This figure represents the theoretical limit of precision for the radar range measurement for this particular impulse system. In practical terms, such precision could only be achieved with a static measurement, whereas the profile record (Fig. 4a) was obtained with a moving system, with waveforms averaged over a horizontal distance of ~3 m. By visual examination of Figure 4a and adjacent profiles, it is clear that changes in basal topography of ⩾3 m are readily identified and can be traced from profile to profile. These subtle features can be tracked for up to 15 km across the survey area (Fig. 7).

Fig. 7. Bed topography beneath Pine Island Glacier. The survey covered an area 25 km × 18 km, total topographic relief is 423 m. Features as small as 3 m can be traced across the survey area.

Accuracy

We have established that the precision of the DELORES radar system is such that topographic differences of >3 m can be mapped across the bed of Pine Island Glacier. However, the accuracy of the location of the basal topography in 3D space is more difficult to assess. The location of the radar system during acquisition was determined using a dual-frequency survey-grade GPS unit. The GPS data were post-processed using a kinematic precise point positioning service provided by the Canadian Geodetic Survey (https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php). This provided a sub-metre horizontal position accuracy for the midpoint between transmit and receive antennas at the start of the process of stacking individual shots into a single trace. However, the radar transmitter and the GPS receiver were not time-coordinated, so the GPS receiver supplied a position at every UTC second while the stacked trace was written by the receiver at ~1.2 s intervals. Therefore in post-processing, the GPS positions were interpolated to the times that the radar traces were acquired. The interpolation introduces a position uncertainty of ~2.4 m (0.8 s at 12 km h⁻¹).

The accuracy of the vertical position of the bed topography was determined by the accuracy of the GPS-derived position of the radar system; by the surface topography under the radar antennas; and by uncertainty in the radar wave-speed profile between the ice-sheet surface and the bed. Of these, the greatest uncertainty is in the wave-speed profile. No measurements of wave-speed through the firn layer were made during the Pine Island Glacier survey. Accepted practise (e.g. Bingham and others, 2017) in these circumstances is to calculate ice thickness using a constant wave-speed of 0.168 m ns⁻¹ and then add a firn correction of 10 m to allow for the faster wave-speed through the firn. If there were spatial variations in firn properties or thickness on scales of tens to hundreds of metres, then the overall wave-speed profile may change sufficiently to create travel-time variations in the bed reflection that could be interpreted as variations in bed topography. However, it is known from airborne radar profiling (Medley and others, 2013) that snow accumulation thickness varies very slowly across the Pine Island Glacier basin, with gradients on layer surfaces of <1 m km⁻¹ at 20 m depth (Medley and others, 2013, Fig. 2). Variations in firn thickness or wave-speed on the scale of 100–200 m of horizontal distance (the scale required to explain the bed topography observed) are therefore unlikely.

Such an assessment of accuracy becomes relevant when comparing the results of two surveys of the same area to determine changes in ice thickness over time (Davies and others, 2018) but is not relevant to geomorphological studies of, for example, bedform shape and size (Spagnolo and others, 2017). This is because the thickness of ice above the bedforms is of secondary importance compared to the peak-to-trough height of individual drumlins or megascale glacial lineations. That is to say that bedform morphology can be determined using relative trace-to-trace measurements even if the absolute location of the bedform is less well known.

Digital Elevation Models of Subglacial Bed Topography

One of the values of detailed mapping of subglacial topography is the ability to compare contemporary landforms with palaeo examples (Spagnolo and others, 2014) in order to better understand the processes of landform development and landform modification during and after deglaciation. To make valid comparisons using digital elevation models (DEMs) of the different terrains, the DEMs must have similar spatial parameters. In this section, the differences between the common methods of mapping palaeo subglacial landscapes and radar profiling of contemporary ice sheets will be examined.

Extensive palaeo-landform datasets have been recorded from continental shelves around Antarctica and northern high latitudes using swath bathymetry (Dowdeswell and others, 2008; Larter and others, 2009; Graham and others, 2010; Klages and others, 2015). Swath bathymetry systems typically produce DEMs of the seafloor with 20–30 m gridcells. Relict subglacial landscapes have also been mapped onshore using aerial photography, Landsat and SPOT satellite imaging at spatial resolutions of between 5 and 30 m (in both x and y dimensions) (Clark, 1997; Clark and others, 2009; Spagnolo and others, 2014). All these methods utilize a measured data point for each gridcell.

Although some progress has been made on swath-mode radar imaging of ice sheets (Paden and others, 2010), the majority of detailed mapping of subglacial bedforms to date has used 2D profiling techniques (King and others, 2007, 2009, 2016; Bingham and others, 2017). This methodology samples the bed at a high spatial sampling rate along the profile (3–10 m data posting). The spatial sampling in the across-track direction can be adjusted to any desired value by positioning the next acquisition track appropriately; however, practical and logistic considerations quickly become a restriction when trying to reproduce the same spatial resolution across-track as along-track. The example data shown here from Pine Island Glacier were acquired over an ice stream where the subglacial landscape was expected to be dominated by mega-scale glacial lineations (Clark, 1993; Stokes and Clark, 2002), features which are tens to hundreds of metres wide but many kilometres long. In this instance, the radar profile lines were positioned approximately perpendicular to the ice flow with a spacing of 500 m, which allowed sufficient data to be acquired in a short field season to sample the bedforms sufficiently well in both horizontal dimensions to characterize their shape and size.

In order to generate a DEM from 2D profile data with spatial resolution similar to those from swath bathymetry or satellite imagery, the data must be interpolated. In this case, a convergent interpolation was used with a rectangular gridcell of dimensions 25 m in the profile direction and 200 m in the cross profile direction (Fig. 8a). The orientation of the grid with respect to the bedforms is critical (Fig. 8b). If the grid orientation is incorrect, spatial aliasing takes place and the algorithm interpolates an en-echelon series of short bedforms instead of correctly lining up long, continuous ridges. This example shows that setting the parameters for the interpolation of data that are sparse in one direction requires interpretation, in a similar way that hand-contouring does. The aim is to derive a DEM of the subglacial topography that is glaciologically/geologically reasonable and which honours the observed data.

Fig. 8. Effects of grid orientation on 3D data interpolation between widely-spaced lines. Bed locations picked from radar profiles 500 m apart are shown in red. Gridcells of size 200 m by 25 m are shown as thin black lines. (a) Gridcells oriented so as to align similar features from one profile line to the next. (b) Grid-cells oriented orthogonal to the radar profile direction, resulting in unrealistic interpolation.

Therefore, until airborne or ground-based swath radar systems become more generally available for subglacial topographic mapping, the DEMs of subglacial topography will require a greater degree of interpretive interpolation than DEMs of palaeo landscapes. In the particular example given here, the subglacial landscape is dominated by 2D forms which reduce interference by off-nadir reflections on the radargrams. If the bed topography of a survey area has complex 3D shapes, then closer-spaced survey lines and 3D migration techniques would be required in order to produce a precise bed topography.

Conclusion

In this paper, I have been asking the question ‘What size of features in the bed topography of an ice sheet can be detected with an impulse radar operating in the HF frequency band?’ The answer to the question depends on the precision of the system, not its accuracy. The issue of accuracy becomes relevant in answering the different question ‘How far beneath the ice-sheet surface is the bed?’

The underlying assumption is that the bed is a step-change in dielectric properties of very much greater amplitude than any adjacent englacial reflecting horizons. In this circumstance, the quarter- or half-wave criteria for resolving between adjacent reflectors becomes irrelevant.

The precision with which a bed feature can be measured is determined by how well the reflected waveform can be reproduced from one measurement point to the next. This, in turn, dictates what is the minimum detectable change in range between the radar and the bed from one trace to the next. Only when the onset of the reflected wave has shifted by a time greater than the digitization interval can we be certain that the range has changed.

In the example dataset from Pine Island Glacier, it was shown that the transmitted waveform was highly repeatable in both shape and timing. As a result, bedforms of 3 m amplitude can be traced consistently from line to line across a 15 km survey grid using a 3.8 MHz radar system. It was found that the alignment of the interpolation grid was critical to avoiding spatial aliasing when adjacent lines of the grid were spaced far apart.