1. Introduction
Extensive satellite and airborne measurements of the surface velocity, surface elevation and mass of the Greenland and Antarctic ice sheets show that both ice sheets are rapidly and unexpectedly losing mass (Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009; Reference ShepherdShepherd and others, 2012). Further analysis reveals that much of this ice loss is confined to the ice-sheet margins via fast-flowing outlet glaciers and ice streams (Reference Joughin, Smith, Shean and FloricioiuJoughin and others, 2014; Reference Rignot, Mouginot, Morlighem, Seroussi and ScheuchlRignot and others, 2014). These documented changes have potentially catastrophic consequences. The Intergovernmental Panel on Climate Change (IPCC) estimates that sea level could increase by 26–98 cm by the end of this century (IPCC, 2013). The large range in projected sea-level rise can be partly attributed to an incomplete understanding of the forces causing the rapid changes in Greenland and Antarctica. To address this deficiency, more thorough information on the bed topography and basal conditions of fast-flowing regions of the ice sheets is critically needed.
Fast flow is initiated by lubrication of the bed, which causes ice to slide. A lubricated glacial bed increases ice velocities by orders of magnitude compared to ice over a frozen bed, and the ice beds of key outlet glaciers and ice streams in Greenland and Antarctica are below sea level. The precise geometry of the subglacial topography of these glaciers and ice streams determines their sensitivity to changes at their seaward margins (Reference ThomasThomas, 1979). As such, information on subglacial topography and basal conditions, in addition to surface speed and surface elevation, is essential to understanding the processes causing outlet glaciers to speed up and ice shelves to break apart, as well as to improving next-generation ice-sheet models for generating more accurate estimates of the ice sheets’ contribution to sea-level rise.
Satellite remote sensing has greatly advanced our monitoring capabilities, as evidenced by the success of routinely deriving surface velocity from interferometric synthetic aperture radar (InSAR) (e.g. Reference JoughinJoughin, 2002; Reference Joughin, Abdalati and FahnestockJoughin and others, 2004, Reference Joughin, Smith, Shean and Floricioiu2014; Reference Simons and RosenSimons and Rosen, 2007; Reference RignotRignot and others, 2008), measuring surface elevation with radar and laser altimeters to infer mass loss (Reference Davis, Kluever and HainesDavis and others, 1998; Reference ZwallyZwally and others, 2005; Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009) and utilizing gravity to estimate total mass loss (Reference LuthckeLuthcke and others, 2006; Reference Velicogna and WahrVelicogna and Wahr, 2006). Although several concepts for spaceborne measurements of ice thickness are currently being explored, most ice thickness measurements are made with surface-based and airborne radars. The only effective way to measure the ice thickness of fast-flowing glaciers and ice streams with extremely rough surfaces is with airborne radars. Radars for airborne sounding of the ice sheets typically operate at the HF and VHF ranges of the electromagnetic spectrum and normally use wide-beam antennas. At these frequencies, much of the surface of the inland ice sheet appears smooth, and most of the off-vertical signals impinging on the surface, referred to hereafter as surface clutter, are reflected away from the radar so bed reflections and internal layers are clearly visible. Faster-flowing glaciers and ice-sheet margins are characterized by very rough surfaces due to extensive crevassing and contain lossy ice near the bed. Between the extremely rough surfaces and the presence of temperate ice near the bed, which results in large attenuation of signals propagating through this ice, airborne radar sounding of these regions is very challenging. High-sensitivity radars using synthetic aperture radar (SAR) processing and cross-track arrays are required to sound the most dynamic and important parts of ice sheets. SAR processing enables the synthesis of a narrow beam to reduce surface clutter in the along-track direction, and the cross-track array enables advanced array processing to reduce surface clutter in the cross-track direction.
The University of Kansas has demonstrated success in the design of advanced radars for sounding fast-flowing glaciers and ice-sheet margins in Greenland and Antarctica. We successfully sounded three key fast-flowing glaciers in Greenland – Jakobshavn, Helheim and Kangerdlugssuaq – with a very sensitive radar for the first time in 2005 and performed additional measurements during the 2006, 2008 and 2009 field seasons to generate bed maps for these glaciers. We performed these measurements with progressively improved radars capable of SAR and array processing in the along- and cross-track directions, respectively. We processed and analyzed data from these measurements to estimate ice thickness for the most challenging glaciers and combined the ice thickness data with surface elevation measurements reported by others (Reference KrabillKrabill and others, 2000; Reference Marzio, Brenner, Schutz, Shuman and ZwallyDiMarzio and others, 2007) to generate the bed topography. In this paper, we provide examples of radar sounding of two key glaciers: Jakobshavn Isbræ, Greenland, and Byrd Glacier, Antarctica. The results of radar sounding of other glaciers are available through the Center for Remote Sensing of Ice Sheets (CReSIS) website (www.cresis.ku.edu; Reference GogineniGogineni, 2012). We are also reprocessing radar data for other glaciers to generate improved bed maps.
Located on the west coast of Greenland, Jakobshavn Isbræ is the fastest-flowing glacier on Earth. Its surface velocity is reported to have recently reached ~17 km a−1 (Reference Joughin, Smith, Shean and FloricioiuJoughin and others, 2014), and it drains ~7.5% of the Greenland ice sheet. The measurement of this glacier’s ice thickness has been a major challenge due to its extremely rough surface and the presence of temperate ice at its bed. Nevertheless, we have successfully sounded many parts of the glacier, starting from ~66 km upstream and proceeding all the way to the calving front in both the along-flow and across-flow directions. We generated a bed map based on these radar soundings by combining available surface elevation data.
On the other side of the world, Byrd Glacier discharges East Antarctic ice into the Ross Ice Shelf through the Transantarctic Mountains and is the fastest Antarctic ice stream/outlet glacier entering a buttressing ice shelf. Byrd Glacier is buttressed by the world’s largest ice shelf, which, combined with its rock side-walls narrowing to ~21 km, makes it one of the most confined Antarctic ice streams. The center-channel flow speed of Byrd Glacier varies from 650 to 850 m a−1 (Reference Stearns and HamiltonStearns and Hamilton, 2005) and actually increased to ~15% above its average during 2005–07. This has been attributed to the drainage of two upstream subglacial lakes (Reference Stearns, Smith and HamiltonStearns and others, 2008). However, detailed investigations of the dynamics-driven changes to this glacier have been hampered by a lack of detailed bed topography. To address the need for information on detailed bed topography and basal conditions, we dedicated the CReSIS 2011/12 Antarctic field season to conducting measurements. In the next section, we provide a brief summary of the major characteristics of Jakobshavn and Byrd glaciers, as well as the radars used to sound and image the ice sheets. We then present an overview of the radars used for measurements in the instrumentation section. Next, we describe signal-processing techniques and show sample results of the application of SAR and array-processing algorithms for the deepest part of Jakobshavn Isbræ. In the results section, we present radar echograms for both glaciers. We conclude the paper with a summary of accomplishments and future plans to improve the bed maps of Jakobshavn and Byrd glaciers.
2. Background
2.1. Jakobshavn and Byrd glaciers
Both Jakobshavn Isbræ, West Greenland, and Byrd Glacier, East Antarctica, have been extensively studied for decades. However, until recently, detailed bed topography has remained incomplete and, in places, incorrect. The glaciers have several characteristics in common, such as the discharge of large amounts of land-based ice to floating ice, channels with overdeepened beds and negative bed slopes, grounding zones in contact with the ocean, and rough, crevassed surfaces. Overdeepened beds are often associated with transient outlet glacier behavior and can lead to rapid changes in flow speed and grounding zone location. Pronounced changes in flow speed have been observed for both glaciers, with some speed changes for Jakobshavn found to be associated with the thinning and break-up of its floating ice tongue (Reference Thomas, Abdalati, Frederick, Krabill, Manizade and SteffenThomas and others, 2003) and some speed changes for Byrd Glacier associated with a release of subglacial water in the upper catchment region (Reference Stearns, Smith and HamiltonStearns and others, 2008). Overdeepened outlet glacier channels are particularly deep targets for radar sounding of the bed and present the additional challenge of temperate ice, which may exist in the deepest parts of the channel. Along with a rough-ice surface, these characteristics have hindered a complete mapping of the subglacial terrain until recently. During the 1978/79 Antarctic field season, a radar survey flight line along Byrd Glacier was conducted from a low-flying (100 m above the surface) LC-130 aircraft and provided a nearly continuous bed echo (Reference Reusch and HughesReusch and Hughes, 2003). However, it is now clear that ice thickness errors of ~1000 m along portions of the profile were caused by a misinterpretation of the bed echo in these data. Using the ice thickness in Reference Reusch and HughesReusch and Hughes (2003), several studies attempted to model the dynamics of Byrd Glacier, models that will now need to be revised due to the very different bed topography described in this paper. For the Jakobshavn Isbræ channel, similar difficulties were encountered in earlier radar surveys attempting to map the bed topography near the calving front, and more complete and accurate channel bed topography will allow a reexamination of the dynamics of this outlet glacier.
2.2. Radar sounding and imaging
Radars operating over the frequency range 1–1000 MHz are frequently used for sounding ice sheets (Reference AllenAllen, 2008). The use of radars for ice sounding began with the pioneering work of Amory Waite. He first used a radar altimeter operating at 440 MHz to conduct ice thickness experiments on the Ross Ice Shelf in the late 1950s (Reference Waite and SchmidtWaite and Schmidt, 1961; Bogordorsky and others, 1985) and is also credited with conducting the first airborne ice measurements in the early 1960s. Following Waite’s pioneering work, the glaciological application of radars was expanded through significant contributions by Evans and his colleagues at the Scott Polar Research Institute, UK (Reference EvansEvans, 1967, Reference Evans and Gudmandsen1970; Reference EvansEvans and Smith, 1969), and Gudmandsen and his colleagues at the Technical University of Denmark (Reference GudmandsenGudmandsen, 1969). Extensive radar soundings were made by Gudmandsen and his colleagues on the Greenland ice sheet during the 1970s with a radar operating at a center frequency of 60 MHz. These soundings resulted in the very first bed map for Greenland. Most of the radars used for sounding ice during the 1960s and 1970s were incoherent, and analog feed networks were used to form transmit-and-receive antenna beams. The incoherent nature of the radar and its inability to form very low side-lobe antenna beams resulted in poor results over fast-flowing glaciers like Jakobshavn and ice-sheet margins with rough surfaces. To obtain side lobes of 30 dB below the main lobe or lower with a five- to six-element array, it is necessary that the amplitude errors between elements be <0.2 dB and phase errors be <3°. Coherent radars for ice sounding have been developed and demonstrated by several groups over the past 20 years (Reference Legarsky, Gogineni and AtkinsLegarsky and others, 2001; Reference Corr and VaughnHélière and others, 2007; Reference Peters, Blankenship, Carter, Kempf, Young and HoltPeters and others, 2007).
Since 1990, several attempts have been made to sound Jakobshavn and other fast-flowing glaciers, with very limited success. Fully coherent radars operating over the frequency range 1–450 MHz are used for these measurements (Reference GogineniGogineni and others, 2001, Reference Gogineni2012; Reference DallDall and others, 2012; Reference Morlighem, Rignot, Mouginot, Seroussi and LarourMorlighem and others, 2014). We first succeeded in sounding Jakobshavn Isbræ using a radar operating at 150 MHz with a cross-track array, digitizing signals from each element of the array, and processing the digitized signals to synthesize low side-lobe beams in both the along-track and across-track directions in 2005. We continued our efforts to improve our radars’ performance and performed more extensive sounding during the 2006, 2008 and 2009 field seasons.
3. Instrumentation
The early version of the multichannel radar depth sounder was developed for operation during the 2005 field season in Greenland on board a Twin Otter aircraft. The radar was operated with a four-element transmit (Tx) and four-element receive (Rx) array installed on each wing of the aircraft. The transmit array was fed through an analog beam-forming network while the signals from each receive antenna element were multiplexed into a single receiver using a single-pole four-throw (SP4T) switch (Namburi, 2003). The received signals from each element were combined with appropriate weights using a delay-sum beam formation technique for areas with clear bed echoes. These signals were summed and saved for post-processing with a SAR algorithm. However, data from all four receive elements were digitized, summed and saved for additional processing with SAR and cross-track array algorithms to reduce clutter in areas with significant surface roughness. The multichannel radar deployed in 2005 was incrementally upgraded for both sounding ice and imaging the ice/bed interface (Reference LohoefenerLohoefener, 2006; Reference MoralesRodriguez-Morales and others, 2014). The upgrades involved modifications to support operation on long-range and short-range aircraft, as well as for low- and high-altitude measurements. Performance improvements included increasing the peak radar transmit power from 200 W to 800 W, expanding the number of transmit-and-receive antennas and channels, using multichannel waveform generators and digitizers to eliminate the need for analog beamforming and receiver multiplexing, and increasing bandwidth (Reference MoralesRodriguez-Morales and others, 2014).
As discussed below (Section 4), the radars have been extensively used for ice thickness retrieval in Greenland and Antarctica. The radar used to collect data during the 2006 Greenland field season consisted of a transmitter, a five-element transmit-antenna array, a five-element receive-antenna array and five data-acquisition subsystems. The transmitter generates a chirped signal from 140 to 160 MHz with peak power of 800 W. This chirp signal is then applied to the transmit array through a weighted feed network to obtain low antenna side lobes. The radars were typically operated with alternate pulses of short (1–3 μs) and long (10–25 μs) duration at a pulse repetition frequency (PRF) in the 10–12.5 kHz range. The short pulse is used to map internal layers in the top 500–1000 m of ice and sound shallow ice up to 1 km thick. The longer pulse is utilized to map deeper layers and sound ice up to 5000 m thick. Other differences between the short- and long-pulse waveforms include the record window, number of hardware pre-sums and receiver gain settings. The center operating frequency and bandwidth values in Table 1 were selected to match the bandwidth of the antenna elements used in each field season. We operated the radar at a center frequency of 150 MHz with a bandwidth of 20 MHz during most flights conducted throughout the 2005, 2006, 2008 and 2009 field seasons. During a reduced number of flights conducted in 2008, the radar was operated with narrower bandwidth (e.g. 1.5 and 6 MHz) to overcome the effects of radio-frequency interference from the aircraft systems. In 2011, the radar was operated at 195 MHz with a bandwidth of 30 MHz. Detailed descriptions of the early versions of the radar are available in Namburi (2003) and Reference LohoefenerLohoefener (2006). A summary of relevant parameters for each of the above radars is presented in Table 1. Here we provide a brief description of the radars used in Greenland in 2008 and 2009 and in Antarctica in 2011.
Summary of relevant (typical) system parameters for the multichannel radar depth sounder operated on board a Twin Otter aircraft from 2005 to 2011. The transmit power listed for years 2006–11 includes the feed network weights for side-lobe reduction
3.1. System overview
The radar depth sounder deployed to Greenland on the Twin Otter aircraft in 2008 and 2009 consists of three main subsystems: digital, analog transmitter and analog receiver. The digital subsystem consists of a waveform generator and a data acquisition system controlled by a computer. The transmit waveform was synthesized using an arbitrary waveform generator (AWG) in combination with an RF mixer driven by a 120 MHz local oscillator (LO) signal. The AWG produces a chirp signal in the 20–40 MHz range that is up-converted by the mixer. The resulting 140–160 MHz signal at the output of the mixer is filtered, amplified and fed to the different transmit antenna elements in the array (right-wing array) using a passive feed network. Phase matching of the transmit antenna array was accomplished by inserting an appropriate-length transmission line into each array element to compensate for time delays measured with a vector network analyzer (VNA). The 120 MHz LO signal for the up-converter comes from the same phase-locked loop (PLL) source as the sampling clock used for the digitizer. The PLL is referenced to a highly stable 10 MHz signal. The signals from the receive array (left-wing array) are passed through separate analog receiver channels. The receiver subsystem consists of a limiter, blanking switch, low-noise amplifier, variable-gain amplifier and bandpass filter. The limiter protects the receiver from high-power transmitter leakage, the blanking switch provides additional receiver protection during the high-power pulse transmission, and the low-noise amplifier and the second-stage amplifier are used to amplify received signals to the level required for digitization. We typically operated the radar with a total analog receiver gain of ~20–30 dB for capturing signals from the ice surface and shallow ice and 50–60 dB for capturing signals from deeper layers and thicker ice.
The output signals from the receivers are fed to the multichannel data acquisition system, which consists of multiple analog-to-digital converters (ADCs) with 12-bit resolution operated in under-sampling mode. The data from the ADCs are stored on high-capacity SCSI drives. For the 2009 Antarctic season and subsequent campaigns (including Antarctica 2011), the transmit-waveform generator was upgraded to a single module that included a set of eight synchronized direct digital synthesizers (DDS) operated with a 1 GHz clock source to directly generate the RF transmit signal and eliminate the RF mixer. The DDS modules provide the capability of modulating the phase, amplitude, frequency limits and pulse duration of the transmit pulse on a waveform-to-waveform basis. The module also generates a sampling clock for digitization of the received signal, the radar PRF, the transmit waveform pulse width, and the receiver blanking intervals for timing and control. The waveform generated by each DDS is bandpass-filtered, amplified and passed to the transmitter subsection. The amplitude and phase of the transmit waveform supplied to each transmit-array element can be adjusted to synthesize a low-side-lobe antenna beam. The analog receivers, analog transmitters and digitizers were also upgraded. The improved transmitter subsystem consists of high-power amplifiers, transmit/receive (T/R) switches, and a bank of bandpass filters. The output from each high-power amplifier is supplied to each element of the transmit array (right-wing array) through the T/R switch and bandpass filter, which is then used to further reduce out-of-band spurious signals.
On the receive side, the signals from each antenna element are passed to the receiver through the T/R switch and a limiter/low-noise amplifier. A variable-gain receiver stage conditions the signal for the ADCs, and the received signals are digitized with 14-bit ADCs at an operator-selectable sampling frequency between 111.11 and 250 MHz. For the 2011 Antarctic season, the analog receivers and data acquisition subsystems were upgraded to support up to 16 channels and used with a 12-element receive array. Figure 1 shows a simplified block diagram of the system deployed to Antarctica in 2011.
Simplified block diagram of the radar system used in 2011.
3.2. Antenna arrays
In 2005, we operated the radar with eight folded-dipole antenna elements, four under each wing of the aircraft. In 2006 we operated with ten folded dipole antenna elements, five under each wing. In 2008 and during subsequent seasons, we operated with a total of 12 elements, six under each wing. With the exception of 2011, during each year the folded dipoles were oriented along the axis of the aircraft. We used the array elements mounted under one of the Twin Otter aircraft wings for transmit-antenna pattern synthesis and the array elements under the other wing for receive. In 2011, the antennas were upgraded to operate at a center frequency of 195 MHz. The smaller antenna size allowed us to accommodate six folded dipoles oriented perpendicular to the aircraft axis under each wing. This was done to help suppress clutter signals at large incident angles and to mitigate the effect of potential in-cabin radio frequency inference (RFI) sources. The 12 antenna elements were used for receive, with the right-wing array shared for transmit/receive operations. Figure 2 shows the antenna arrays used for measurements in Greenland (2008 and 2009) and in Antarctica (2011).
Photograph of the Twin Otter aircraft equipped with wing-mounted antenna arrays. The inset shows the antenna orientation used in 2008, 2009 and 2011.
4. Field Programs
The multichannel radar depth sounder has been operated extensively in several areas of Greenland and Antarctica on board the Twin Otter aircraft. Table 2 presents a summary of the field campaigns conducted with this type of aircraft between 2005 and 2009 in Greenland and in 2011 in Antarctica. Two survey grids stand out due to their extensiveness: one over Jakobshavn Isbræ and one over Byrd Glacier. Most of the flights over Jakobshavn were completed in 2008 operating mostly from Ilulissat Airport. A reduced number of flights were conducted on the east coast of Greenland that year, operating from Kulusuk. In 2009, the surveys targeting Jakobshavn were complemented by surveys over Nuuk Glacier on the west coast of Greenland, as well as Helheim and Kangerdlugssuaq glaciers on the east coast. The flights over Byrd Glacier were completed during the 2011/12 austral summer season in Antarctica, operating from McMurdo Station. Figure 3 shows maps of the survey lines flown over Jakobshavn Isbræ (2006/2008/2009) and Byrd Glacier (2011).
Flight lines surveyed with the multichannel radar depth sounder over (a) Jakobshavn Isbræ, Greenland, and (b) Byrd Glacier, Antarctica.
Flight-line summary for the 2005–09 Greenland field seasons and 2011 Antarctic season
Flight lines are chosen to meet the science requirements of the mission and are compiled prior to deployment. There have been situations where, after post-flight analysis of collected data, flight lines are modified and possibly re-flown to improve the signal-to-noise ratio (SNR) of the bed echoes. Grid spacings are chosen to adequately sample the expected bedrock topography based on modeling needs. This typically results in lines that are more closely spaced near the glacier channel and that are oriented in both the along-flow and cross-flow directions. Flight lines and grid spacings are also dependent on the aircraft configuration and possibly the requirements of other sensors on board, such as repeating laser altimetry lines for estimating height changes. Flight envelopes are driven by safety requirements, flight endurance and instrument requirements. Although lower altitudes reduce surface clutter, surveys are usually conducted at 500 m above the ice surface. Twin Otter aircraft speeds are ~70 m s−1 to maximize coverage, and bank angles are limited to ±20° to maintain a GPS lock. These are typical operating parameters that have been occasionally modified for specific applications.
5. Signal Processing and Radar Echograms
Radar data are processed with SAR and array-processing algorithms to extract weak bed echoes buried in clutter. SAR uses the forward motion of the aircraft to synthesize the long aperture needed to obtain fine resolution in the along-track direction. The principle of SAR operation and associated processing algorithms are well described in many textbooks and papers (Reference Curlander and DonoughCurlander and McDonough, 1991; Reference Raney, Runge, Bamler, Cumming and WongRaney and others, 1994; Soumekh, 1999). The SAR algorithm we used is based on the frequency–wavenumber (F-K) migration technique to improve the SNR and obtain fine resolution for reducing clutter in the along-track direction. We also processed a few selected files with a time-domain processor to obtain a small improvement in SNR for detecting extremely weak signals (Reference SoumekhSoumekh, 1999). Since the use of SAR for processing radar depth sounder data is extensively described by Reference Legarsky, Gogineni and AtkinsLegarsky and others (2001), Helière and others (2007) and Reference Peters, Blankenship, Carter, Kempf, Young and HoltPeters and others (2007), we provide, in Section 5.2, only a very brief summary of the algorithm we used for processing much of the data over the two glaciers. We used traditional sum-and-delay and minimum variance distortionless beamformers (MVDR) to reduce clutter in the cross-track direction (Reference CaponCapon, 1969; Van Trees, 2001).
SAR and array processing, coupled with radar improvements, allowed us to overcome >110 dB of attenuation, scattering and return losses and sound the deepest parts of the narrow channel under Jakobshavn with a radar for the first time.
5.1. Data conditioning
The first step after loading the raw data is to convert the digital quantization level to a receiver input voltage level. This is done to remove variations between datasets caused by different receiver gain settings and digital acquisition systems. If the received signal from the surface is saturated (generally the case when the receiver is in high-gain mode), the saturated portion of the record is set to zero. This is done by first processing the low-gain data to determine the surface location or using another data source (e.g. a microwave radar or laser altimeter) to determine the surface location. Generally, an amplitude-tapered, linear frequency-modulated chirp is used for the radar’s transmitted pulse.
A matched filter with frequency-domain windowing is used to correlate the transmitted pulse with the received data and ‘compress’ the pulse. The matched filter is normalized so that the choice of windowing and pulse duration does not change the output magnitude; this allows different datasets to be compared. The frequency-domain window, usually a Hanning window, is applied to reduce pulse-compression range side lobes. The pulse-compression filter also compensates for the radar system time delay for each receiver channel individually so that zero time for each channel is the moment when the transmit waveform first begins to radiate from the antenna. After pulse compression, the data are typically converted to complex baseband and decimated to the signal bandwidth to reduce the data volume.
To determine the location of each radar pulse, the data are synchronized to the GPS and inertial navigation system (INS) data using GPS time stamps that are stored with the radar data. The lever arm between the GPS trajectory reference position and each radar antenna is used in combination with the aircraft attitude measured by the INS to convert the GPS trajectory into a trajectory for each antenna.
Using the trajectory data, we create a Cartesian flight coordinate system for each output range line after Reference Wahl and GhigliaWahl and others (1996). A position-dependent coordinate system is necessary because the flight path contains turns (sometimes >360°) that need to be processed. The x –axis is located parallel to the flight direction. The z –axis is a projection of the local elevation vector to the plane orthogonal to the x- axis (this definition works as long as the plane is not flying straight up or down). The y –axis is formed as the cross product of the z – and x –axes to complete a right-handed coordinate system. The squint vector is the vector that points from the radar antenna towards the center of the target scene and is taken to be in the negative z –direction.
5.2. F-K focusing algorithm
After the data are conditioned as described in Section 5.1, we process data with SAR algorithms to improve the along-track resolution by synthesizing a long aperture. The primary algorithm used is a frequency–wavenumber (F-K) focusing algorithm that exploits the fast Fourier transform for computational efficiency. However, we require straight and uniformly sampled data before the fast Fourier transform can be applied. Since the aircraft’s speed is not constant and the trajectory is not straight, the raw data do not meet these criteria without additional processing. We approximate a uniformly sampled dataset by spatially re-sampling in along-track using a sinc kernel with a bandwidth equal to the Doppler bandwidth required for SAR processing. The vertical trajectory deviation from the horizontal straight line required by the F-K migration is compensated in the frequency domain with a phase shift that corresponds to the two-way propagation time over the height difference.
The radar reflectivity function, s (x, z, t), is related to its Fourier transform as
where
We can rewrite the above equation as
where kx is the wavenumber in the x –direction (along-track), kz is the wavenumber in the z –direction, ω is the radar frequency (rad) and v p is the velocity of propagation in air or ice.
The measured field, at z = 0, on the aircraft is given by
Thus the processing involves determining the field at the target location in ice by back-propagating the measured field at z = 0 using
The term e jk z z can be considered a filter with a transfer function H (z) that consists of two terms as
The first filter, H 1(z), propagates the signal from the aircraft at z = 0 back to the ice surface at z = h, and H 2(z) propagates the signal from the ice surface at z = h to the target location in ice at depth d. The inverse Fourier transform of the back-propagated signal yields the desired radar reflectivity function s (x, z).
After focusing, the time-delay shifts that were applied to approximate a straight trajectory are removed to preserve the original geometry for array processing, as the array processing can account for non-uniform sampling.
The F-K focusing algorithm produces a constant along-track resolution; synthetic aperture length L SAR increases with target range, and is given by
where λ c is the wavelength at the center frequency, h a is the height above the ice surface, T is the ice thickness,r,ice = 3.15 is the relative permittivity of ice, and σx is the along-track resolution. For the images in this paper, a resolution of 2.5 m is used. At a flight altitude of 500 m and an ice thickness of 2000 m, this is a 650 m synthetic aperture.
The Nyquist criterion to avoid along-track clutter aliasing in the Doppler spectrum for an isotropic antenna is quarter-wavelength sampling, or 0.5 m at 150 MHz. For the Twin Otter’s nominal speed of 70 m s−1 this translates to a minimum PRF requirement of 140 Hz. The nominal PRF we used is at least 10 kHz. Accounting for the use of two waveforms to capture a low-gain and high-gain sample of the data, the effective PRF is 5 kHz. This far exceeds the Nyquist criterion. The data are pre-summed in hardware before being recorded with a boxcar pre-summing filter. Suppression of unwanted Doppler frequencies will be limited by this boxcar filter. In 2008 and 2009, the pre-summing was set so that a record was registered every 156.25 Hz, which meets the Nyquist criterion. However, the system flown in 2006 could not store data quickly enough and a record was registered every 78.125 Hz. In 2006, clutter at ±90° aliases to nadir, our direction of interest. There are two mitigating factors. First, there is a null in the antenna pattern at ±90° both for transmit and receive arrays. Second, the side-lobe structure from hardware pre-summing has a null close to ±90°. Therefore, for the angles close to nadir, clutter aliasing is unlikely to be a severe problem, even in data from 2006.
5.3. Array processing
As described earlier, the complex data collected from each receive antenna channel are first pulse-compressed then SAR-processed to generate a radar echogram before array processing takes place. The pulse compression and SAR processing operates in the fast-time and slow-time domains, respectively, to increase SNR. SAR also improves along-track resolution. Array processing, or receive beamforming, is performed in the cross-track domain, the third dimension, to further improve the SNR and the cross-track angular resolution to reduce surface clutter. Traditional antenna array theory states that an SNR improvement of 20 log(N) and a beam resolution of λ /L apt (where N is the number of elements, λ is the operating wavelength and L apt is the array aperture length) can be achieved by coherently combining the radar echogram from each receive channel with uniform weights and zero delay for the nadir direction. This is commonly known as the sum-and-delay beamformer (Van Trees, 2001). It is usually used to combine the multichannel data with uniform or Chebyshev weights in the interior with little surface clutter. The Chebyshev weights are amplitude coefficients derived from the Chebyshev polynomials. For uniformly spaced linear arrays, Chebyshev weights offer the minimum side-lobe level for a given beamwidth (Van Trees, 2001). For areas with significant surface clutter (e.g. the downstream areas of outlet glaciers), a parametric beamformer with higher angular selectivity based on the MVDR algorithm is used.
The MVDR beamformer uses data to determine the optimum weights for passing the desired signal without distortion and minimizing the interference and noise power (Van Trees, 2001). The weights for combing data from each element of the antenna array in the cross-track direction are computed using
where w is the weight vector, S (θ0) Þ is the steering vector used to receive signals in the nadir direction, and R c +n is the clutter-plus-noise correlation matrix estimated from data. The number of data samples used to estimate the correlation matrix depends on the decorrelation properties of the signal. In general, we used 50 neighboring samples for the correlation matrix estimation.
To demonstrate the clutter scenario that we are dealing with, we performed a cross-track signal analysis on the along-glacier radar data we collected on 30 May 2006 in Jakobshavn, Greenland (frame 20060530_08). Figure 4 shows the corresponding map and Figure 5 shows the radar echograms obtained with the delay-and-sum beamformer for Chebyshev weights and the MVDR beamformer. The radar was operated with only high-gain waveforms and high-gain receivers for this line. We used the blanking switch to prevent receiver saturation from the surface and near-surface internal layers. The ringing of the switch resulted in moiré patterns.
A map showing the flight line for segment 20060530_08 in red and frame 20090401_05_006 in green.
(a, b) Radar echograms obtained by (a) the sum-and-delay beamformer and (b) the MVDR beamformer (the region under the dashed blue line is where the data are used to calculate the correlation matrix). (c) The corresponding angle-of-arrival estimation result obtained by the MUSIC algorithm (red indicates strong signal return and blue represents weak or no signal; the color bar shows the relative power in dB). The inset in (b) shows the comparison between the antenna array radiation pattern obtained by the MVDR and sum-and-delay beamformers at locations indicated by the white arrows.
The radar echogram is divided into three regions for analysis (Table 3). Region I is the first 30 km in the upstream area where the ice surface is flat and the ice loss is low. Region II is the next 20 km of the flight line toward the calving front, shown between dashed red vertical lines. Over this portion of the line we start observing surface clutter, because the surface is rough and the bed topography begins to vary. Region III is the final 10 km to the calving front, where the ice is warm and the surface is heavily crevassed and extremely rough.
Region for array processing along segment 20060530_08
Comparing the radar echograms, it can be seen that both beamformers demonstrate similar performance in region I, as the bed return is not clutter-limited and both beamformers can provide adequate array gain to retrieve the bed signal. In region II, the MVDR beamformer outperforms the sum-and-delay beamformer, and the bed signal is recovered because the undesired clutter signal is minimized by the MVDR beamformer. Finally, in region III, only the MVDR beamformer provides additional clutter reduction, which enables identification of weak bed echoes after displaying the image on a high-resolution monitor and utilizing image processing to reduce clutter and enhance the weaker echoes.
We also performed an angle-of-arrival analysis using the Multiple Signal Classification (MUSIC) algorithm on the received radar signal (signals below 1 km thick ice from the green dashed line in Fig. 5a), and the result is given in Figure 5c. It can be seen that the signal in region I mostly comes from nadir, while the clutter signal (at large angles) becomes stronger in regions II and III. These results not only coincide with the results obtained from the beamformers, they also provide strong support that the radar performance is limited by the surface clutter near the calving front rather than the volume scattering of the ice.
To illustrate the complexity of sounding Jakobshavn Isbræ, we generated estimates of the combined loss experienced by a wave propagating through both the ice and the reflection at the ice/bed interface, referred to as return loss, for three regions. Figure 6 shows the loss for three regions and ice in the interior. The loss is ~38 dB km−1 for region I, 48 dB km−1 for region II, and 70 dB km−1 for region III; it is ~20 dB km−1 for interior ice. Radars that attempt to sound glaciers like Jakobshavn must overcome >100 dB of combined propagation and return loss and must have combined transmit-and-receive antenna side lobes of 60 dB or lower than the main-lobe peak to obtain discernible bed echoes.
Approximate propagation loss for ice (ice loss + return loss) in three regions of Jakobshavn Isbræ and a comparison to loss for interior ice.
Another example of bed recovery with array processing is shown in Figure 7. We applied three different beamformers – sum-and-delay, MVDR and MUSIC – to a cross-channel line near region II; the flight line is given in Figure 4. It can be seen that the glacier channel can be recovered nicely with both MVDR and MUSIC beamformers.
(a–c) Radar echograms obtained with (a) delay-and-sum, (b) MVDR and (c) MUSIC beamformers. (d) Bed echoes recovered with the MUSIC algorithm.
5.4. Image processing
The bed picks for certain areas of Jakobshavn and Byrd glaciers, especially the region close to the calving front, are very challenging due to the low image contrast and low SNR. To further improve the quality of the images of these challenging areas, an image-processing approach composed of contrast enhancement, adaptive median filters and customized filters was developed to enhance selected areas of the radar echograms to generate discernible bed echoes for estimating ice thickness.
The first step of image processing is to equalize the histogram of the selected region to enhance the contrast. The original histogram around the bed region shows only a small range of gray values. After histogram equalization, the image contrast improves considerably, as the equalization processing maps a small grayscale range to a full contrast range. After optimizing the contrast, the adaptive median filter is applied to reduce the speckle noise and enhance the SNR. The adaptive median filter belongs to the class of nonlinear edge-preserving smoothing filters, and it effectively smooths the data, retaining small and sharp details. In the final image-processing step, we applied a Bas Relief-based edge-detection filter (Reference Weyrich, Deng, Barnes, Rusinkiewicz and FinkelsteinWeyrich and others, 2007) to further enhance the edges and bring out weak bed echoes, as shown in the bottom image of Figure 8.
Radar echograms before and after image processing. The top image shows the original image with a weak bed return. The middle image shows the improved image after applying histogram equalization and an adaptive median filter to the region of interest. The bottom shows the final image after applying a customized filter to the enhanced bed return.
Additional image processing of SAR and array-processed data enabled us to bring out weak bed echoes for the most challenging area near the calving front. Examples of echograms after image processing are given in Figures 8 and 9.
6. Results
In this section, we show several full-bandwidth sample radar echograms to discuss the results for both glaciers. We have posted additional echograms on our website (www.cresis. ku.edu/Jakob_Byrd_echograms) to demonstrate how we generated the bed maps for these glaciers.
6.1. Jakobshavn Isbræ
Figure 10 shows the flight lines used to collect the Jakobshavn data given in this section, and Figures 11–13 show the resulting radar echograms in both the along-flow and across-flow directions, starting below the calving front and extending all the way to ~60 km upstream. Figure 11 shows a radar echogram generated from data collected with an aircraft flying along the glacier as close to a flowline as possible on 30 May 2006. The line extends a few km downstream of the calving front to ~10 km upstream of the calving front. The ice thickness is ~850 m at the calving front, increases to ~1500 m, then decreases to ~1200 m between 5 and 7 km. Beyond this, the ice thickness varies between 1400 and 1200 m. We also observe two distinct echoes separated by 100–200 m in the along-flow line echogram caused either by side-wall echoes and bed echoes or multiple interfaces at the bed. The ice thickness data obtained from across-flowline echograms are shown as blue, red and green triangles and as red circles. In 2006, most of the flight lines were flown at a height of ~500 m above the surface. At this height, it was not possible to obtain clearly discernible echoes, particularly in the first 10 km from the calving front, due to surface clutter heavily masking the bed echoes. We did obtain discernible bed echoes for the same along-glacier line flown on 30 May 2006 at a height of ~150 m above the surface, which was allowed by NASA safety requirements at that time. The results from this line are those shown in the left echogram of Figures 10–12. The reduction of clutter for the low-altitude flight also confirms that the off-vertical surface scatter is the dominant source of clutter, as indicated by our direction of arrival (DoA) analysis. Although changes in ice thickness measured in 2006, 2008 and 2009 might be related to ice thinning, they cannot be confidently attributed to this because of returns from multiple interfaces near the bed. These multiple returns are clearly visible in the along-glacier lines, but are not visible in across-glacier lines. The sloping bed limits the length of a SAR aperture for across-flow lines and reduces the SNR improvement required to sound both interfaces.
Narrow-band radar data echogram from Figure 5 after image processing. Top image shows the echogram obtained after applying histogram equalization and a customized filter (Bas Relief filter). The bed picks are shown in the bottom image, which reveals a complex bed topography with multiple interfaces.
A map showing the flight lines for sample data from Jakobshavn in Figures 11–13.
(a) Along-flow echogram for the first 10 km from the calving front with estimated ice thicknesses from across-flow echograms shown as circles and triangles. (b, c) Sample results for across-flow direction lines. The dotted red and blue lines in (a) and the marked blue lines in (b, c) are the bed picks. The vertical red dotted lines in (b, c) represent crossover lines at the along-flow locations indicated in (a).
Figure 12 shows an along-flowline echogram for the next 25 km of the flight line in conjunction with two across-flow echograms. The ice thickness for this 25 km line increases from ~1200–1500 m at the start to ~2000–2300 m at the end. We again observe returns from two layers near the bed. The returns from the top layers are potentially from the cold/temperate ice interface, and the bottom return is from the temperate ice/bed interface. These signals are only a few (2–5) dB above the thermal noise of the radar. It is not possible to observe returns from both interfaces in the echo-grams of the cross-flow lines, because the slope of the glacier bed limits the SAR aperture length, which in turn limits the SNR improvement that can be obtained by increasing the aperture size, as can be done for along-flowline data.
(a) Along-flow echogram for the next 25 km from the calving front with estimated ice thicknesses from across-flow echograms shown as circles and triangles. (b, c) Sample results for across-flow direction lines.
Figure 13 shows the final 30 km of data collected along the flowline in conjunction with two across-flow echo-grams. The ice thickness varies between 2400 and 2700 m over the 30 km length of the glacier. We observe returns from multiple interfaces for this portion of the line, which is clearly documented in Figure 9. The DoA analysis indicates that the returns largely come from the bed at nadir, and a few echoes come from the side-walls of the narrow channel. We observed weak returns from the bottom interface whenever strong echoes appeared from the top interfaces This indicates that some of these are layered echoes instead of off-vertical returns from the walls of the channel. Many of the across-flow ice thickness measurements are in agreement with returns from the bottom interface.
Same as Figure 12 but for the next 30 km of the flight line.
6.2. Byrd Glacier
We obtained excellent data for most of the Byrd Glacier flight lines shown in Figure 3b, except for a few lines passing through the deepest part of the glacier in the middle. We processed these data for lines close to center with a time-domain SAR processor, and used additional image-processing techniques described in Section 5.4 to bring out weak bed echoes. We also interpolated data from nearby lines with good bed echoes to generate estimates of ice thickness for the missing parts of the line and compared estimated bed picks to verify that the weak returns were correctly identified as ice-bed echoes.
Figure 14 shows the radar echogram ice thickness of three tributaries feeding Byrd Glacier; the corresponding flight line is shown in pink on the inset map. This map also includes blue flight lines in the along-ice-flow direction, over which data were collected. The channel of the middle tributary at this location is 2.5 km deep.
An example of Byrd Glacier RDS measurements in the catchment area. The red segment on the inset map shows the location of the echogram. The main trunk of the glacier is at the bottom left of the map.
We sounded the ice bed along the channel for all flight lines offset by 4 km from the center line; the corresponding radar echograms are shown in Figures 15 and 16a. As shown by these two echograms, the maximum depths of the trench on both sides of the center line exceed 2800 m. To increase radar sensitivity to sound the deeper parts along the center line, we increased the radar PRF from 12 kHz to 17 kHz. As shown by the vertical dashed red line in the radar echogram of Figure 16b, we clearly mapped the ice bed at 3300 m depth, where the steep slope suggests a greater depth downstream. Because we used a higher PRF, the data-recording window was not long enough to sound the deepest part with this setting. We processed the data collected with a lower PRF on earlier flights with a time-domain processor to obtain additional improvement in the SNR; we also processed the resulting echograms with image-processing techniques to detect weak bed echoes. Figure 17 shows the enhanced radar echogram generated from TDP output and reveals that the maximum ice thickness for the deepest part of the channel is ~3623 m. To verify that the weak returns are in fact ice-bed echoes, we estimated the missed part of the center-line ice bed by interpolating the results from cross-channel flight lines, as shown in Figures 18 and 19.
A radar echogram of flight lines along the Byrd Glacier trunk to the south side of the center line.
Radar echograms of flight lines along the Byrd Glacier trunk: (a) to the north of the center line and (b) along the center of the trunk.
Deepest part of the Byrd trunk with TDP and image processing.
(a) Flight segments used for interpolation include FLS 1–4. The red circle indicates the location (80.7108° S, 155.5579° E) of the interpolated maximum depth (3623 m). The blue circle indicates the lowest-elevation location (80.7102° S, 155.5634° E; bottom elevation is 3045 m below sea level, ice thickness is 3623 m). The two locations are very close, and the two circles almost overlap. The location of maximum depth is 120 m upstream of the location of the lowest elevation. (b) Spline interpolated ice bottom elevation profile along FLS 1. The ice bottom at crossover 2 is not visible in the echogram of FLS 1 and is derived from the visible ice bottom echo in the echogram of FLS 3. The ice bottom at crossovers 1 and 3 is visible in the echograms of FLS 1 and 2 and FLS 1 and 4; the depth differences at crossovers 1 and 3 are ~7 and ~41 m, respectively. The correlation coefficient of surface and bed topography is maximized (0.94) with a horizontal offset of ~1710 m, and the red circles (delay markers) show the inflection point in surface topography and the expected inflection point at the bed based on this horizontal offset. The flight path of the radar echogram of Figure 15 is ~100 km and parallel to the center line in (a) with an offset of 4 km, showing the similarity of ice bottom topography compared to the interpolated profile.
7. Bed Maps and Error Analysis
The ice-bed elevation can be derived by subtracting ice thickness from ice surface elevation. Figures 20 and 21 show ice bed elevation (with respect to WGS84) maps for Jakobshavn and Byrd glaciers, respectively, generated by uploading the ice thickness profiles derived using radio-echo sounding (RES) data, ice surface digital elevation map (DEM) and IceFree Mask to software package ESRI ArcGIS 10.0. The maps are gridded to a resolution of 500 m × 500 m. The RES data used to derive the ice thickness profile of Jakobshavn Isbræ include the Multichannel Coherent Radar Depth Sounder/Imager (MCoRDS/I) data from three field seasons (2006, 2008 and 2009) in Greenland. The ice thickness profile of Byrd Glacier is derived from the RES data collected by MCoRDS/I during the 2011/12 Antarctic field season. For ice surface DEM, the data used include NASA Airborne Topographic Mapper (ATM) lidar data and Ice, Cloud and land Elevation Satellite (ICESat) DEM with grid resolutions of 80 m (width) 40 m (along-track) and 1 km, respectively, (http://nsidc.org/data/docs/daac/icebridge/ilatm2/index.html; Reference Marzio, Brenner, Schutz, Shuman and ZwallyDiMarzio and others, 2007). The data used for IceFree Mask have a resolution of 15 m (Reference Howat, Negrete and SmithHowat and others, 2014).
Radar echograms of flight lines across Byrd Glacier trunk: (a) upstream; (b) in the middle; and (c) downstream.
Ice-bed elevation maps of Jakobshavn Isbræ: (a) two-dimensional (2-D) including the catchment area and (b) three-dimensional (3-D) of the main channel.
Ice-bed elevation maps of Byrd Glacier: (a) 2-D and (b) 3-D.
The ice-bed elevation uncertainty consists of ice surface DEM errors and ice thickness estimate errors. The NASA ATM lidar data have an accuracy of 10 cm (Reference KrabillKrabill and others, 2000) and the ICESat DEM data have an accuracy of ±14 cm (Reference ShumanShuman and others, 2006). The difference in propagation time between the ice surface and ice bottom echoes is converted into ice thickness using a constant ice refraction index of 1.77. The ice is assumed to be uniform and no firn correction is applied. It will result in <10 m of error in the ice thickness estimate using the constant ice refraction index and ignoring the effects of firn. We generated our ice thickness estimation error by comparing estimates with existing ice-core data at sites such as NEEM, GRIP, NGRIP and Jakobshavn channel (Reference Dahl-JensenDahl-Jensen and others, 1997; Reference GogineniGogineni and others, 2001; Reference Funk, Iken, Gogineni and TrufferLüthi and others, 2002). We found that our estimates were within ±10 m of the core data.
There are several other potential sources of error in ice thickness estimates, including: (1) the uncertainty in airborne platform geolocation; (2) the uncertainty in range measurement and the manner in which the ice surface and bottom were picked; (3) misinterpretation of off-nadir clutter as the ice bottom in scenarios of strong surface clutter and low SNR of the ice bottom; and (4) interpolation for undetected ice bottoms. The RES data were geo-registered with GPS data with <10 m and 100 m of error in the vertical and horizontal directions, respectively. The vertical location error of the platform will not result in ice thickness estimate error, because the offset exists for both the ice surface and the ice bottom, and systematic error is cancelled. However, the horizontal location error may result in error in ice thickness estimates for areas with rapid changes in bed topography.
We traced the ice surface in the radar echogram automatically, based on the peak power of the strong ice surface returns. We traced the ice bottom by manually entering picks along the visible ice bottom interface, then performing interpolation. The manual ice-bed picks may have errors of one to three range bins based on statistical studies. The overall picking error depends on range resolution, range accuracy and the depth of each range bin. Table 4 lists these parameters for the three radar bandwidths of 10, 20 and 30 MHz that were used during data collection. The range accuracy is derived using an SNR of 3 dB (assuming this is the minimum SNR for discernible bottom echoes). Based on Table 4, the ice thickness error from picking is ~25.63 m for the worst case (range accuracy and range bin errors are added for both the ice surface and the ice bottom).
Radar data parameters related to ice thickness error
Differences in ice thickness at flight-line crossovers are a comprehensive metric for ice thickness error from the above four error sources. By performing crossover analysis, we further carefully investigated locations with unusually large crossover errors and corrected these errors, confining the ice thickness estimate errors and reducing the possibility of wrongly interpreting ice surface echoes as ice bed echoes. Figure 22a and b show the crossover errors of the ice thickness profiles of Jakobshavn and Byrd glaciers, respectively. The mean value of thickness error is ~35 and ~30.6 m, and the standard variation of thickness error is ~40.7 and ~40 m, for Jakobshavn and Byrd data, respectively. Figure 22c shows the percentage or probability of ice thickness error greater than a certain value for Jakobshavn and Byrd data (9040 and 1334 crossovers, respectively). For example, the probability of ice thickness error greater than 100 m is ~7% and ~8%, respectively, for Jakobshavn and Byrd data. We observed crossover errors as large as 300 m mainly in the first 10 km of Jakobshavn near the calving front. These large errors are a result of two factors. One is a lack of clearly discernible multiple interfaces in the across-flowline measurements, and the other is the complex bed topography. SAR processing reduces the along-track footprint, but the cross-track footprint is determined by the radar pulse width and effective antenna beamwidth. The typical effective antenna beamwidth is ~20°, corresponding to a cross-track pulse-limited footprint of 315 m and beam-limited footprint of ~580 m at 2000 m depth. The location of the dominant signal for along-track and cross-track returns is displaced from nadir, resulting in large crossover errors.
Ice thickness crossover analysis: (a) crossover errors of Jakobshavn Isbræ; (b) crossover errors of Byrd Glacier; and (c) percentage as a function of thickness error at crossover.
8. Discussion
In this paper, we have presented bed topography for Jakobshavn Isbræ based on measurements made from 2006 to 2011 with progressively improved radars. The bed topography presented extends from the calving front to >200 km upstream in the eastward direction. The main channel of the glacier extends ~150 km into the interior. Three tributary channels can be identified from the bed topography; two of these tributaries go eastward by >150 km, and one has a northeast to southwest direction and is connected to the main channel at ~75 km from the calving front.
The ice thickness is ~850 m at the calving front. For the first 10 km from the calving front, the crossover errors are large, mainly because of the difficulty identifying multiple interfaces in across-flowline echograms and weak echoes barely visible above the surface clutter. The radar echograms in the along-flow direction show two distinct types of echo. One of these is from the ice/bed interface, and the other, particularly upstream, is from the cold/temperate ice interface as verified by DoA analysis. Temperature measurements in the channel and next to the channel indicate 200–500 m of temperate ice in the upstream area (Reference Iken, Echelmeyer, Harrison and FunkIken and others, 1993; Reference Funk, Echelmeyer and IkenFunk and others, 1994; Reference Funk, Iken, Gogineni and TrufferLüthi and others, 2002). Reference Funk, Iken, Gogineni and TrufferLüthi and others (2002) reported the presence of temperate ice near the calving front and estimated an ice thickness of ~850 m from upturned ice. The ice thickness determined from our radar soundings at the calving front is very close to this estimate. In addition, our measured ice thickness of ~2.7 km is within 100 m of that determined with seismic soundings by Reference Clarke and EchelmeyerClarke and Echelmeyer (1996).
Analysis of radar data indicates that ice-bed echoes are masked by off-vertical surface scatter from an extremely rough surface in the first 10 km from the calving front. Radars with combined cross-track transmit-and-receive antenna side lobes of 60 dB or lower with respect to the main lobe peak are required to reduce the surface clutter to obtain discernible ice-bed echoes. Synthesis of antenna patterns with 30 dB or lower side lobes, particularly with a small number of elements (six to eight) in the presence of amplitude and phase errors caused by aircraft wing flexure and other effects, is a major challenge. As a demonstration, we simulated the effect of amplitude and phase mismatch between array elements on antenna array beam formation. We synthesized weights to obtain side lobes of 30 dB below the main lobe for a perfect six-element dipole antenna array. We also introduced random amplitude and phase errors between elements. Figure 23 shows the simulation results, which demonstrate that amplitude errors must be <0.2 dB (20 mV for an average of 1000 mV) and phase errors must be <3° to obtain antenna side lobes of 30 dB or lower. Although we made significant progress in reducing surface clutter, further reduction of 10 dB or more is required to obtain clearer bed echoes near the calving front than those given in this paper.
Simulated antenna array side-lobe level (SLL) for random amplitude (σw ) and phase (σz ) variations. Amplitude and phase errors have to be in the dark-blue region in order to achieve near 30 dB SLL.
The bed topography of Byrd Glacier contains a deep subglacial trench, which is ~3.62 km below the surface in the main channel. We obtained excellent bed returns for most of the flight lines flown over Byrd Glacier, with the exception of those flown over the deepest part of the glacier. We obtained clear bed echoes for ice with a thickness of ~3.4 km. We used advanced signal- and image-processing techniques to obtain identifiable bed echoes and verified that these echoes are consistent with the ice thickness estimated by interpolating data from lines with good bed echoes. The bed topography shows four main tributaries feeding into the main channel, as well as smaller channels feeding into these tributaries.
We generated bed maps for Jakobhavn and Byrd glaciers from data collected with radars capable of along-track SAR and cross-track array processing. The bed map for Byrd Glacier can be further improved with additional data over the deepest part. The radars must have an additional 10– 15 dB of sensitivity to obtain clear bed echoes over this region. The additional 15 dB sensitivity can be obtained by increasing the transmit power to ~4000 W from the 500 W used during the 2011 field season, as well as by obtaining additional antenna and signal-processing gain of ~6–7 dB. The signal-processing gain can be obtained by processing data with a time-domain processor (TDP) that accounts for changes in the dielectric properties near the bed and the bed slope. The antenna gain can be increased by using additional elements, as is being done with an ultra-wideband radar under development by CReSIS for use on a BT-67 aircraft.
The bed map for Jakobshavn Isbræ can be improved with additional flight lines in the first 10 km of the glacier near the calving front. We will reprocess all existing datasets with a time-domain processor that includes corrections for refraction and surface topography effects and bed slopes and changes in the dielectric properties of ice near the bed. The use of incorrect velocity of propagation can reduce SAR processing gain, and errors as small as ±10% result in a loss of gain (Yilmaz, 1987). Correction for slope and velocity effects requires iterative processing of data. Processing of pulse-compression radar data to synthesize a long aperture in the along-track direction is accomplished with two separate matched filters: (1) one to compress the long-chirped frequency-modulated pulse to a short pulse to improve range resolution; and (2) the other to convolve the received signal in the along-track direction with the complex conjugate reference function of a point target to obtain fine along-track resolution. The matched filter operations are implemented in the frequency domain to reduce computational load with a slight loss in processing gain. Although the use of a time-domain processor can result in a small amount of much-needed improvement in processing gain, the TDP has a much higher computational load.
Our results also show that weak bed echoes near the calving front are masked primarily by surface clutter, i.e. off-vertical scatter from the rough ice surface. The surface clutter can only be reduced by synthesizing a narrow beam in the cross-track direction. Thus, the new data near the calving front must be collected with radars operating on small Unmanned Aerial Systems (UAS) so they can be flown on closely spaced gridlines, as shown in Figure 24. We have already tested a UAS for this application (Reference LeuschenLeuschen and others, 2014) and are planning to deploy it to collect data to demonstrate the concept and generate improved bed maps during the next field season in Greenland.
Illustration of how a small UAV can be flown along closely spaced parallel tracks to collect sounding data over a fine grid. This technique, along with the use of array processing in 2-D, helps to reduce surface clutter.
9. Conclusions
In this paper we have demonstrated the need for radar sounding of fast-flowing glaciers and ice-sheet margins. We described an overview of the radars we used to sound Jakobshavn and Byrd glaciers, as well as other glaciers in Greenland and Antarctica. We presented radar echograms to discuss the challenges associated with sounding and imaging fast-moving glaciers and showed results obtained with improved radars coupled with advanced signal- and image-processing techniques. Finally, we presented the bed topography of Jakobshavn and Byrd glaciers, including a detailed error analysis.
The development of advanced radars coupled with the latest signal- and image-processing techniques has allowed us to sound some of the most interesting glaciers in the world and generate bed topography maps. Without detailed bed topography and information on basal conditions, attempts to model the response of the ice sheets to climate change will remain speculative. The results presented in the paper show that we can derive the topography of fast-moving glaciers like Jakobshavn and glaciers with deep channels like Byrd with radars capable of SAR and array processing. Further improvements in quantifying amplitude and phase errors and reducing their effect on cross-track beam formation would lead to the radar soundings needed to generate fine-resolution bed topography and to map the basal conditions of key glaciers in Greenland and Antarctica.
Acknowledgements
This work was supported by the US National Science Foundation under grant ANT-0424589, by the National Aeronautics and Space Administration under grant NNX10AT68G and by the University of Kansas. We acknowledge the contributions of current and former faculty, staff and students at CReSIS in developing the radars used for data collection and processing earlier datasets. Staff members at University Information Technology Services, Indiana University, Bloomington, IN, USA, also supported data collection and archival and data processing in the field. We also thank the pilots and crews of the NASA P-3 and the pilots and personnel from Kenn Borek Air Ltd for their support during many successful field deployments. We acknowledge the support provided by J. Collins in editing and formatting the manuscript and S. McGill in generating graphics. Finally, we thank the two anonymous reviewers for helpful feedback.
Open access
