Introduction

Cold ice is relatively transparent to radar, which has been used since the 1950s to map the bedrock beneath glaciers and ice sheets to several kilometres depth. In this connection, interest in the radar properties of the firn layer has traditionally been restricted to correcting radar times in the ice-depth calculation via a time offset.

The last decade has seen an increased number of radar studies of the near-surface, for example, to determine accumulation rates by finding datable internal reflecting horizons in ice or firn (Reference Forster, Davis, Rand and MooreForster and others, 1991; Reference Siegert, Hodgkins and DowdeswellWeertman, 1993; Reference Holmlund, Richardson and P.Holmlund and Richardson, 1995; Kohler and others, 1997; Reference Petrenko and WhitworthRichardson and Holmlund, 1999; Reference Moore and FujitaPalli and others, 2002; Reference Eisen, Wilhelms, Nixdorf and MillerEisen and others, 2003). In addition, the dielectric properties of near-surface snow, firn and ice are important for interpreting synthetic aperture radar (SAR) satellite images (Reference RobinRott and others, 1993), particularly with the emergence of SAR for monitoring ice-sheet elevations and flow patterns.

Radar images of snow or ice can be modelled as a time-domain convolution of an input radar pulse with a time series of reflection and transmission coefficients caused by whatever physical properties produce echoes. Internal reflection horizons within ice or snow masses are thought to be caused by changes in ice chemistry, ice density or crystal fabric (Reference PalliParen and Robin, 1975; Reference Ackley and KeliherAckley and Keliher, 1979; Reference Miners, Hildebrand, Gerland, Blindow, Steinhage and WolffMoore 1988).

In principle, it should be possible to determine snow and ice properties from radar soundings, thus providing better spatial coverage than is possible from snow pits and ice cores. Qualitatively, radar reflections have been shown to correlate well with ice-core electrical conductivity peaks due to acidity layers of volcanic origin (Reference Rott, Sturm and MillerSiegert and others, 1998; Reference Hempel, Thyssen, Gundestrup, Clausen and MillerHempel and others, 2000). Forward modelling, using measured ice-core dielectric properties to generate model radargrams (Reference Miners, Hildebrand, Gerland, Blindow, Steinhage and WolffMiners and others, 1997), has been less successful at replicating prominent reflecting horizons.

In this paper, we forward model a radargram by convolving an idealized input ground-penetrating radar (GPR) monopulse with the complex reflectivity inferred from dielectric measurements made on a shallow core taken directly on the GPR profile, using a simple convolution algorithm and neglecting multiple reflections. The maximum depth considered in this experiment is 10 m, so that the primary source of radar reflections is expected to be permittivity contrasts due mainly to density variations (e.g. Reference Eisen, Wilhelms, Nixdorf and MillerEisen and others, 2003), rather than conductivity contrasts (Reference Fujita and MaeFujita and Mae, 1994; Reference ArconeArcone, 1996; Reference Kohler, Moore, Kennett, Engeset and ElvehøyKohler and others, 1997). Density is highly variable in this regime, both vertically and horizontally, and our objective is therefore to test whether density variations correlate consistently with a radar profile.

Data

We obtained a GPR profile in May 1996 along the centre line of Kongsvegen (78°48’N, 12°59’E), a polythermal glacier near the Ny-Ålesund research station in Svalbard, Norway. Here we use a 30 m long portion of the GPR profile near stake 8, in the glacier’s accumulation zone. The elevation of the site is about 670 m a.s.l., with a mean annual temperature of roughly –10°C.

We used a Geophysical Survey Systems Inc. (GSSI) Spaceborne Imaging Radar-2 (SIR-2) control unit and a GSSI 3102 shielded 500 MHz antenna unit towed behind a snowmobile travelling at a constant speed. Traces were taken at regular time intervals, corresponding to roughly 30 cm intervals, in a 120 ns time window containing 512 16-bit samples. No time-variable gain was applied. The data were post-processed by applying a high-pass filter to remove low-frequency roll in the upper part of the radar profile.

At the study site, a 7.5 cm diameter, 8.8 m long firn core was taken from the bottom of a 1.2 m deep snow pit using a Polar Ice Coring Office (PICO) auger. Air temperatures during drilling were around –10°C. The weight, length and average diameter of each recovered core piece were measured to calculate bulk density (Fig. 1a). Density in the snow pit was measured using snow tubes in the sides of the pit. Since individual core pieces vary from 10 to 50 cm in length, the density measurements are too coarse to show individual ice layers.

Fig. 1. (a) Bulk core density; (b) permittivity and capacitance, as measured (grey line) and after being modified (black (c) conductivity and conductance. Ice lenses 41cm thick are shown by grey lines, and core breaks by dotted lines. The both core depth and the two-way travel time from Equation (7) and density profile.

The core was stored at –20°C for 1 week before analysis. Dielectric profiling was done on the core pieces in their plastic bags using a HP4284 LCR bridge operating at 100 kHz and connected to a 4 cm wide electrode with 140° angle of arc (e.g. Reference MooreMoore, 1993). We measure capacitance and conductance every 2.5 mm along the cores, and then average over 1cm intervals (Fig. 1b and c). It was not possible to make dielectric measurements on the uppermost 1.3 m of the snow-pack because the loose snow did not sit properly in the instrument’s cradle. While core samples were on the cradle, the visual stratigraphy was recorded, noting, in particular, the location of ice lenses.

The measured values of capacitance C and conductance G were converted to permittivity " and conductivity σ using a simple geometric air-capacitance method (e.g. Reference Miners, Hildebrand, Gerland, Blindow, Steinhage and WolffMoore, 1988):

where Cair = 0.97 pF is the geometric air capacitance calculated from measurements on solid ice with known permittivity = 3.17 (Reference Glen and ParenGlen and Paren, 1975), and "0 = 8.854pFm–1 is the absolute permittivity of free space.

Density contrasts produce changes in dielectric impedance through their effect on permittivity (Reference Glen and ParenGlen and Paren, 1975), such that ice lenses correspond, for the most part, to local peaks in the capacitance/permittivity signal (Fig. 1b and c). Conductivity varies in response to changes in ice chemical composition (Reference WeertmanWolff and others, 1995); the drop in conductance values (Fig.1c) at around 3 m depth represents the transition from ion-rich winter snow to the ion-depleted previous summer’s snow.

While the dielectric profiling frequency is considerably lower than the radar frequency, it is commonly accepted (e.g. Reference Glen and ParenGlen and Paren, 1975; Reference MooreMoore and Fujita, 1993) that there are no changes in the dielectric properties of pure ice from roughly 100 kHz up to radio-microwave frequencies. The transition to the high-frequency limit is determined by the frequency range of Debye relaxation processes (Reference Paren, de and RobinPetrenko and Whitworth, 1999). Impurities in ice can increase the transition frequency closer to or even above the DEP measurement frequency. One way to address this problem is to measure the core permittivity at a range of frequencies to determine the plateau value by extrapolation. Our experience with ice cores from several polar regions suggests this is not a problem for Svalbard glacier ice. However, the presence of ionic impurities does raise the relaxation frequency of the composite ice/air firn mixture and plastic/air gap blocking layer capacitor sandwich between the electrodes. This resulted in relatively high permittivity values in the upper core. Given that permittivity and dry-snow density correlate positively, as exemplified by relations of the type

(Reference Richardson and HolmlundRobin, 1975), this means that our permittivity data are at odds with the observed density profile, which, as expected, shows density increasing down-core.

If impurities affected permittivity uniformly throughout the core, then we could simply apply a correction factor to the single-frequency measurements. Permittivity and density correspond most poorly in the upper interval, from 1.3 to 3.1 m, as expected from the drop in conductivity in the ion-depleted snow below 3.1m. We adjust the permittivity signal by applying Equation (2) so that the bulk density implied by the permittivity signal for individual core pieces matches that actually measured. Applying a uniform correction factor for each individual core piece would lead to artificial jumps in the dielectric properties at the core boundaries, so we apply instead an exponentially decreasing correction factor fit to the observed piecewise correction factors. This is defensible on the grounds that the reflectivity coefficients we calculate are primarily affected by the high-frequency variations in permittivity.

Modelled Reflection Coefficients

We assume simple plane-wave propagation and calculate reflection coefficients from the impedance data. We neglect multiple reflections since they are insignificant in firn and ice. For example, the reflection coefficient between solid ice and reasonably low-density firn may be as high as 0.1, and a multi-bounce consists of a minimum of two further reflections back and forth, each with, at most, reflection coefficients of 0.1. This results in a multiple at least 100 times smaller than the primary reflection. In practice, because the ice layers in our cores are thin, multi-bounce is even less important.

We calculate the reflection coefficient between two layers of material of different impedance following Reference Ackley and KeliherAckley and Keliher (1979). In using a single reflection coefficient value, we assume that the dominant response of the snow-pack occurs at the centre frequency of the input wavelet.

The characteristic bulk impedance Z of layer indexed by i is defined by

where is the magnetic permeability of free space, γ_i is the propagation constant given by

and ε_i and σ_i are the relative dielectric constant and conductivity of the layer. The reflection amplitude r_i at a boundary within the medium is given by

where Zi–1 is the characteristic bulk impedance of the i–1 layer, is the surface impedance given by

and hi is thickness of the layer i.

Before the reflection coefficients r _i, which are specified in the core-depth domain, can be convolved with an input radar wavelet, they must be transformed into the time domain. The two-way travel-time vs depth relation is calculated from our data using Reference Richardson and HolmlundRobin’s (1975) expression (Equation (2)) for dry snow with a geometrical correction for the effect of receiving– transmitting antenna separation and taking account of GPR timing calibration constants:

where T is the two-way travel time, T₀ is the effective zero-point for the radar, c is the speed of light, p is the average snow density between the surface and a depth D, and l_a = 18 cm is the antenna separation. We neglect wave refraction within the snowpack since density changes have a negligible effect on wave path length with such small antenna separations.

We can test the relation by comparing the two-way travel times of the reflecting horizon of the previous year’s summer layer (obtained in the 12 km GPR profile obtained along the centre line of the glacier) with depths measured by manual probing into the snowpack (e.g. Reference Kohler, Moore, Kennett, Engeset and ElvehøyKohler and others, 1997). There appears to be good agreement between the two-way travel times of the sounded reflecting horizons and the manual probing depths, using T₀ = 1.8 ns (Fig. 2). The agreement depends, in part, on how representative the core density is for the entire snowpack along the sounding profile; density measurements made in other years along Kongsvegen’s centre line show that there can be some spatial variability (J. Kohler, unpublished data).

Fig. 2. Two-way travel-time vs depth relations given by Equation (7) and density profile measured at stake 8 (thick black line), together with relations corresponding to ice, air and two intermediate density values, and depth to previous year’s summer surfaces determined by manual probing along the glacier’s centre line, plotted against the two-way travel time to the reflecting horizon assumed to correspond to the summer surface (+). The sounding at stake 8 is boxed.

Modelled Radargram

GPRs transmit pulses that are close to monopulses, whose centre frequency is dependent on the tuned antenna and the coupling between the antenna and the surface upon

which it rests. For the SIR-2 500 MHz system on snow, the centre frequency is 400–450MHz. We generated shaped pulses following Reference Miners, Hildebrand, Gerland, Blindow, Steinhage and WolffMoore (1988), using:

where Tl is the pulse length = 1/f where f = ω/2π is the centre frequency of the pulse, chosen to match the dominant frequency observed for each system, β is a phase angle, and t' = t – Tl/2.

The contribution _gi from a particular reflecting layer over an arbitrary time interval ξ is given by convolving the real input wavelet generated by Equation (8) with the complex reflection coefficient r _i:

The dielectric data at 1cm intervals are transformed into the time domain using Equation (7) and resampled using linear interpolation at 0.05 ns intervals. Reflectivity coefficients r_i are calculated using Equation (5). Summing in the time domain the amplitude- and phase-altered wavelets from each reflection layer given by Equation (9) yields a modelled radargram for the entire ice core (Fig. 3c).

Fig. 3. (a) 30 m long 500 MHz GPR profile, with the modelled radargram embedded in the centre, at the location of the (b) “Wiggle”plot of the two radar flanking traces. (c) The modelled radargram computed from core dielectric measurements. Box in upper part shows the input waveform used; dashed lines at 16 and 93 ns demarcate the parts of the radagram unaffected lack of dielectric data, at core top and bottom.

Results

Evaluating goodness of fit by simply comparing modelled and observed radar traces on a sample-by-sample basis (i.e. calculating a correlation coefficient) can go awry for even small offsets in timing, such as those that could arise due to minor errors in the depth-to-travel-time conversion, for example. A simple graphical technique for qualitatively evaluating goodness of fit, often used in artificial seismology, is to insert the modelled radargram into the appropriate space along the GPR profile.We do this, first adjusting the power in the modelled radargram to match that of the appropriate part of the observed traces on either side of the core (Fig.3a).

Figure 3a shows that the modelled radargram does a credible job of reproducing many of the significant reflectors; it is nearly indiscernible from flanking traces in the profile below about 40 ns. Above this time, the visual correlation is not as good, although there appears in both the observed and modelled radargrams at least one distinctive reflector signature, consisting of a –+–+– sequence, at 25–28 ns in the observed trace, and 27–30 ns in the modelled trace.

Discussion and Conclusion

Overall, we find a reasonable correlation between the model radargram and the GPR profile.Where correlation breaks down between the observed and modelled radargrams, this could be due to:

1. (i) timing errors in the radar traces,

2. (ii) error in the input waveform,

3. (iii) errors in the core dielectric data,

4. (iv) incorrectly predicted travel-time vs depth relation,

5. (v) spatial variability in ice dielectric properties on the cm scale.

The first explanation would most likely affect the observed radar profile either by an offset or by a multiplicative factor. The reasonable fit between soundings and picked travel times (Fig. 2) argues for this not being a significant factor.

We tried different model input pulses, as did Reference Miners, Hildebrand, Gerland, Blindow, Steinhage and WolffMiners and others (1997) more exhaustively; in both investigations, varying the pulses did not appear to qualitatively improve the correlation.

The dielectric data are a more obvious source of error, especially given the problems we have in matching permittivity and density with the raw data. As we argue above, though, the latter should affect the lower-frequency content of the simulation, and should not lead to gross errors in placement of reflections. Higher-frequency errors in the dielectric data, on the other hand, such as those due to varying cross-sectional core geometry or apparent changes of dielectric properties across core breaks, would emphasize or suppress reflectors. There is, however, no way to assess this from the data we have available.

Explanation (iii) has merits; judicious stretching of the time coordinates of the dielectric log one way or the other could help line up peaks of the modelled radar signal to match those observed. Stretching could be achieved by applying different density profiles, or by assuming small shifts in the assigned core depths of individual core pieces. Indeed, the correct assignation of core depth is a perpetual difficulty in ice-core work.

We do not reject the above possible explanations outright but tend to favour explanation (v), which we believe is the most likely, as well as the most glaciologically parsimonious, explanation. Radar beam angles are relatively large (Reference ArconeArcone, 1995), so that individual traces comprise returns over larger areas than those sampled by the core. It is clear from the radar profile, if not from the practical experience of digging a shallow snow pit, that layers and thereby reflecting horizons are variable over short distances, even though they comprise information from reflectors integrated over larger spatial regions. It is therefore reasonable that the dielectric properties of sections of the core might be different from the mean properties of the entire snowpack in the vicinity of the coring site. An obvious way of testing this hypothesis, and a future direction for forward modelling of radargrams, would be to collect several to many cores at a site, and take the horizontal average of the dielectric properties to remove the effect of local variations.