Inversion for the density‐depth profile of polar firn using a stepped‐frequency radar

Translating satellite measurements of ice sheet volume change into sea level contribution requires knowledge of the profile of density as a function of depth within the ice sheet and how this profile changes over time. This paper describes an interferometric method of inverting ground‐penetrating radar returns for the profile of firn density as a function of depth. The method is an interferometric implementation of the common‐midpoint approach, performed using a stepped‐frequency, phase‐sensitive ground‐penetrating radar. By recording the phase difference of returns with a range of antenna separations, the different path lengths through the firn allow recovery of a smoothed representation of the density profile. This density model is characterized by three parameters: surface density and two decay lengths for porosity, each operating over a different density range. Our results suggest that the stepped‐frequency radar used here can accurately recover differences in two‐way traveltime and produce useful estimates of the density profile. In a test of the method performed at Summit station in Greenland, the recovered density‐depth profile agreed with independent density measurements from an ice core and a neutron probe to within 6% root‐mean‐square error.


Introduction
[2] Satellite radar and laser altimeters can be used to measure changes in the volume of the Greenland and Antarctic ice sheets [Wingham et al., 2006;Pritchard et al., 2009]. Before these altimetric measurements can be used to assess the contribution to sea level from the ice sheets, volume changes must be translated into mass changes [e.g., Zwally et al., 2005]. This requires knowledge of the density profile in the upper layers of the ice sheet and how it changes over time [Braithwaite et al., 1994;Arthern and Wingham, 1998;Zwally and Li, 2002]. Changes in the rate of compaction have been observed using strain meters installed in boreholes [Arthern et al., 2010], but this has been limited to locations where shallow ice cores have been drilled.
[3] The density profile is also needed for ground-based radar surveys of snow accumulation [e.g., Arcone et al., 2004]. Ground-penetrating radar used over ice sheets record reflections from horizons beneath the surface. It is usually assumed that laterally continuous reflecting horizons mark snow layers deposited at the same time, referred to as isochrones. This assumption has been validated using ice cores and has formed the basis for surveys of the rate of snow accumulation over the ice sheets. Often, shallow ice cores are drilled along the route and used to date the isochrones. Once the layering has been dated at a core site, the isochrones can be followed for many hundreds of kilometers. The shallow cores also provide estimates of the density-depth profile, and this is used to estimate the speed of propagation of radio waves through the snowpack and to compute the mass of snow lying between any two dated isochrones.
[4] On any particular traverse, it may be difficult to drill enough shallow cores to fully sample the variations in firn density. In this paper, we develop a remote sensing capability to provide estimates of the density profile away from shallow core sites. The method is an interferometric version of the common-midpoint approach used in seismic processing [Yilmaz, 1987] and often applied to radar returns from polar ice sheets and glaciers [e.g., Robin et al., 1969;Macheret et al., 1993;Murray et al., 2000;Eisen et al., 2002;Nath and Vaughan, 2003;Bradford et al., 2009;Brown et al., 2012]. The basis of the method is that returns collected with

Rx1 Rx2
Rx1 Rx2 Tx2 X1 Tx1 X2 Figure 1. Inverting for density by collecting returns with different antenna separation at a common midpoint. The difference in path length introduces a phase difference between the two returns. By measuring this phase difference accurately and comparing it with predictions from a theoretical model, the density profile can be estimated as a function of depth. different separation between transmit and receive antennae, as shown in Figure 1, experience different two-way traveltimes, and this provides a means of inferring the subsurface density profile. The method relies on the sensitivity of the two-way traveltime to the velocity of radio waves in the upper layers, which itself depends upon the density.
[5] The general problem of inferring the subsurface velocity and density from traveltime observations is ill posed. Here the problem is regularized by assuming the density profile is approximated well by a model with three parameters: one describes the surface density and two others, each operating over a different density range, describe the rate of increase of density with depth.
[6] Aside from the density model, the main difference in this study from previous glaciological applications of the common-midpoint method is the use of a stepped-frequency radar [Corr et al., 2002] and the use of interferometry to detect differences in traveltime (see Neidell and Taner [1971] for similar approaches). The stepped-frequency radar affords accuracy and temporal stability of frequency and phase, allowing differences in the two-way traveltime as a function of antenna separation to be determined from interferometry of repeated measurements. We show here that this allows the parameters that approximate the shape of the density profile to be recovered to useful accuracy. The interferometric approach described here should be useful for accumulation surveys, for validating models of firn compaction, and for assessing the sea level contribution from ice sheets.

Theory
[7] If different antenna separations X 1 and X 2 are used to record reflections from a particular layer of snow at depth z beneath the surface, as shown in Figure 1, each antenna will record a signal that has propagated along a different path. The two reflections will arrive after different delays 1 and 2 . Assuming both pulses are transmitted with the same phase and experience the same phase change upon reflection, the phase difference, (X 1 , X 2 , z), is given by in which f is the center frequency of the radar and = 2 -1 . In this section, we derive an expression for .
[8] The propagation delays depend upon the speed of propagation, v(z). This varies with depth, z, depending upon the relative dielectric permittivity, (z). From Kovacs et al. [1995], the permittivity depends upon density and is approximated well by the relationship: where k = 8.45 10 -4 kg -1 m 3 is a constant, c is the speed of light in free space, and (z) is the density as a function of depth.
[9] Because of refraction, the angle of propagation with respect to vertical, Â (z), also depends upon the profile of permittivity, (z), through Snell's law: where Â 0 is the incidence angle, i.e., the angle of propagation in free space just above the surface. Defining s Á sin Â 0 , the two-way traveltime can be approximated for small s as follows: [10] Following Castle [1994], we have retained all terms except those proportional to s 6 and higher powers. We have also defined the following for notational convenience: [11] The antenna separation can be approximated in a similar way as [12] Squaring equation (6), neglecting terms proportional to s 6 and higher, and then eliminating s using equation (4) give [13] Solving this quadratic, we obtain the following form [e.g., Castle, 1994]: [14] From equation (1), the relative shift in phase [15] Equation (9) is derived for small angles from vertical. Nevertheless, it includes all terms except s 6 and higher powers and therefore represents a better approximation to wider angles than the normal moveout theory that neglects terms proportional to s 4 and higher powers [e.g., Eisen et al., 2002]. Furthermore, equation (9) also tends to the correct limit at large angles, i.e., z/X ! 0: If s is the permittivity just beneath the surface, z = 0, we have [16] This corresponds to the change of phase with antenna separation for a ray reflected at shallow depth and traveling almost horizontally just beneath the surface.
[17] We adopt equation (9) as our model for how phase difference depends upon antenna separation. To make use of equation (9), we still need to evaluate D 1 (z) and D 2 (z). These quantities depend upon the permittivity profile, through equation (5), and this in turn depends upon the density profile through equation (2).

Modeling the Density Profile
[18] In this paper, we parameterize the density profile using a simple model in which porosity decreases exponentially with depth. Because density typically increases with depth much more rapidly near the surface, we use a different value for the decay length in each of two density regimes, separated at a critical density c = 550 kg m -3 . The first decay length L 1 applies in the low-density regime, Ä c . In this density range, grains are free to rearrange themselves by sliding [Alley, 1987]. The second decay length L 2 applies in the high-density regime, > c , where grains are jammed and compaction occurs more slowly [Herron and Langway, 1980]. Under this assumption, the density profile can be parameterized as follows: where z c is the depth at which the critical transition density c is reached, s is the surface density, and i = 917 kg m -3 is the density of pore-free ice.
[19] Based on this model and equation (2), we can perform the integrals in equation (5). The results are [20] Wherever (z) appears in these expressions, it can be evaluated using equations (11) and (2). The permittivities s , c , and i correspond to densities s , c , and i , respectively. The function G( ) is defined as [21] For any given choice of s , L 1 , and L 2 , the combined system defined by equations (2), (7), (11), (12), and (13) can be solved iteratively using Newton's method to provide depth z for any combination of two-way traveltime and antenna separation. Evaluating D 1 (z) and D 2 (z) allows phase differences to be computed using equation (9).

Experimental Setup
[22] On 25 May 2009, an experiment was performed at the Greenland Summit using a stepped-frequency radar to detect the difference in phase between returns collected at different antenna separations. A stepped-frequency VHF radar with a center frequency of 314 MHz and a bandwidth of 142 MHz was used for this experiment. The radar transmits and receives at a sequence of 7095 frequencies, equally spaced at 20 kHz intervals across the bandwidth. The amplitude and phase of the return signal are then reconstructed in the traveltime domain by taking the inverse Fourier transform.
[23] Returns were collected at 21 different antenna separations ranging from 6 to 46 m, in increments of 2 m. The transmit and receive antennae were each moved outward by 1 m between shots to maintain both at equal distance from a common midpoint. The point of reflection at any particular depth is therefore assumed the same for all 21 recorded returns.

Results
[24] Figure 2a shows the return power plotted as a function of two-way traveltime and antenna separation. The highest power reflections come from near the surface at small antenna separation. The power decreases with depth and with antenna separation. The arrival time of the direct  Figure 2b) and measured (in Figure 2c). Phase differences are plotted at separation X 1 and are relative to the return collected at separation X 2 = X 1 + X, with X = 2 m in all cases. The measured phase differences were unwrapped by adding or subtracting multiples of 2 to give the closest match to the theoretical phase difference. The arrival time of the direct wave is plotted in white.
wave through air is plotted as a white line in Figure 2a. The power before the arrival of the direct wave is much lower than the power reflected from beneath the surface, indicating a high signal-to-noise ratio. Discrete reflecting horizons are apparent, and many appear consistently across a range of antenna separation.
[25] Figure 2b shows the phase difference predicted by the theory as a function of traveltime and antenna separation. The plot shows the phase difference between a return with separation X 1 and the next return with separation X 2 = X 1 + X, with X = 2 m. The phase difference is small for the deeper reflections, as the two paths are then almost parallel. Phase difference is larger in magnitude for shallower returns or for returns with larger antenna separation.
[26] Figure 2c shows the phase difference measured using the radar. There is an ambiguity in the measured phase, because it always lies between 0 and 2 . Before plotting, the measured phase differences were unwrapped by adding or subtracting multiples of 2 to give the closest agreement between the theoretical phase differences and the measurements. Even allowing for this unwrapping procedure, the theory agrees closely with the observed phase differences. The measured phase differences suffer from some noise, especially where the power in one or other of the returns is very low, since this makes it difficult to unambiguously detect the phase of the return signal. Restricting the comparison to cells where power level in both returns is greater than -40 dB gives a root-mean-square discrepancy of 0.31 radians.
[27] Figure 3 shows how the root-mean-square discrepancy between the measured and theoretical phase difference varies when the surface density, s , and the two decay lengths, L 1 and L 2 , are varied. Again, only returns that both had return power greater than -40 dB were used, giving   Figure 4. Density recovered from common-midpoint inversion (solid black line), from ice core data (red points), and from a neutron density probe (blue points). The 90% confidence interval, estimated from Monte Carlo simulations, is shaded. 660 estimates of phase difference for this comparison. The approximate spread of traveltime and antenna separation spanned by these returns can be identified from Figure 2a. The root-mean-square error is minimized by the parameter choices s = 280 kg m -3 , L 1 = 27 m, and L 2 = 42 m, and these parameters were used to generate Figure 2b. Figure 4 shows the density profile that corresponds to these parameters.
[28] To characterize the uncertainty in the density inversion, we performed a Monte Carlo assessment of the errors. Subtracting the measured and theoretical phase differences shown in Figure 2 produces a population of residual-phase differences. An ensemble of synthetic data sets was generated by "bootstrap" resampling of these residual-phase differences [e.g., Press et al., 1992]. Our inversion method was then applied to each synthetic data set, to gauge the variability in the recovered parameters, and in the recovered density profile. Roughness of the surface can introduce a systematic variation in the phase difference by affecting the path length for all signals to or from a given antenna. To account for this systematic component, the error in the measured phase difference was modeled as a bias plus a random component. The bias was assumed different for each pair of antenna separations used (i.e., for each pair of X and X + X). In the Monte Carlo calculations, this component was sampled at random from the population of the residualphase differences averaged over two-way traveltime. The random component was sampled from the population of residual-phase differences after subtracting the average over two-way traveltime. Figure 4 shows the 90% confidence interval derived from these Monte Carlo simulations. These indicate relative uncertainty in retrieved density of approximately 5% within the upper 100 m. From equation (2), this corresponds to relative uncertainty in velocity of about 2%.
[29] As an independent test of the inversion, densities derived from an ice core are also plotted (data digitized from Alley and Koci [1988]). Density data collected using a neu-tron probe [Hawley et al., 2008] are also shown. The close agreement between the densities from the neutron probe and the ice core data suggests that there have not been major temporal changes in density in this region and that the density profile is similar from location to location in the region of the Greenland Summit. The densities from our inversion agree with the ice core and neutron probe data to within 6% root-mean-square error.

Discussion and Conclusions
[30] We have described a method of inverting for firn density as a function of depth, using a phase-sensitive, stepped-frequency VHF ground-penetrating radar. Our results suggest that the radar recovers differences in twoway traveltime to an accuracy of about 0.3 radians in terms of phase. This compares with an accuracy of about 2 radians reported for a similar glaciological study by Eisen et al. [2002]. The assessment of uncertainty given by Eisen et al. [2002] may be pessimistic, and higher accuracy could probably be obtained by processing of time-domain data. For example, Brown et al. [2012] and Bradford et al. [2009] assume that relative velocity errors are about 2% in their inversions. It would be interesting to apply our inversion method to pulsed radar data to see how it compares with the traveltime inversion method used by [Brown et al., 2012]. Our interferometric approach does not require manual picking of coherent reflections, so it may be better suited to automated processing of large data sets. The accuracy and temporal stability of frequency and phase afforded by the stepped-frequency system would not necessarily be obtained with all ground-penetrating radars.
[31] Our approach applies to larger offset to depth ratios than the normal moveout equation used by Eisen et al. [2002]. Nevertheless, there are a number of other assumptions inherent in our modeling of the phase differences. We have assumed horizontal layering, an isotropic medium, and lateral homogeneity of density. To relax some of these assumptions, more sophisticated approaches to the inversion could be applied. Ray-based approaches such as reflection traveltime inversion [e.g., Brown et al., 2012;Harper and Bradford, 2003] can allow for lateral inhomogeneity, anisotropy, larger angles, and dipping reflectors, provided that a high-frequency, geometric optics assumption is valid. A variety of migration algorithms [e.g., Fisher et al., 1992;Hill, 2001] can compensate for off-nadir reflections from lateral inhomogeneities, reducing noise and establishing the physical depth of any particular reflection. Applying the common-midpoint approach to migrated returns (prestack depth migration) can provide information about lateral variations in velocity through iterative refinement of a velocity model [Leparoux et al., 2001;Pipan et al., 2004;Bradford, 2006;Bradford et al., 2009]. We have assumed that the dielectric discontinuities between layers are not strong enough to give significant multiple reflections. To include the effects of multiple reflections would require full waveform inversion [e.g., Virieux and Operto, 2009]. This approach seeks a model of the subsurface properties that best reproduces the observed amplitude and phase of waveforms using adjoint-based inversion methods.
[32] Another assumption inherent in our method is that the density model is capable of reproducing the approximate shape of the density profile. Numerous density profiles have been recovered from dry-snow regions of Greenland and Antarctica [Herron and Langway, 1980], and these typically exhibit shapes similar to that shown in Figure 4. However, for regions subject to melting, percolation, and refreezing, the density profile is likely to be more complicated, and the particular density model and the scattering assumptions used here may not be appropriate.
[33] It would be tempting to try to resolve the density profile at higher depth resolution. However, even with the three-parameter model used here, the problem is not well conditioned. This is clear from the middle panel of Figure 3 where there is a shallow valley in the root-mean-square mismatch between model and data. This shows that a different combination of surface density and the decay length parameter L 1 can produce almost as good an agreement with the observations as the best fit. If data are subject to higher noise levels or if fewer common-midpoint separations are available, it may only be possible to obtain a two-parameter model based on surface density and one representative decay length, or even a one-parameter model based on a single decay length, in which case an independently specified estimate of surface density would be needed.
[34] The interferometric method should allow density to be recovered as a function of depth more rapidly than shallow drilling. The system deployed at our study site was not designed to recover density continuously along a radar traverse, but a modified version of the method could be applied in a mobile or airborne system [e.g., Bradford et al., 2009]. This would provide continuous profiles of density and layering, giving the accumulation rate more accurately. A moving system with fewer antenna separations available may not be able to recover the full three-parameter model at each common-midpoint location, but there would be more opportunity for averaging data, so the potential of such a system should be investigated further. For mobile systems, pulsed or frequency-modulated continuous wave (FMCW) radars have the advantage of collecting data more rapidly, so it would be worth investigating how our inversion method performs when applied to FMCW or pulsed returns. To make more precise observations while traveling, it would be beneficial to survey the location of transmit and receive antennae using kinematic global positioning [e.g., Bradford et al., 2009], especially if sastrugi or other surface undulations are present.
[35] The interferometric nature of the method should make it possible to study changes in snow compaction by comparing the phase of returns from different years, as well as different antenna separations. This would reveal how changes in the compaction of near-surface layers is affecting satellite observations of surface elevation change. In particular, observations of strain rates that are presently only made in shallow boreholes [e.g., Arthern et al., 2010] could be extended over a wider area. Correcting the satellite observations for temporal changes in firn compaction should provide better estimates of how the large ice sheets of Greenland and Antarctica are contributing to sea level rise.