Introduction
Continuous profiles of temperate glacier depth have been obtained only from the ice surface and have not exceeded 300 m depth using frequencies of 5–30 MHz (Reference Jacobel and AndersonJacobel and Andersen, 1987; Reference Welch, Pfeffer, Harper and HumphreyWelch and others, 1998; Reference Arcone, Lawson, Moran, Delaney, Noon, Stickley and LongstaffArcone and others, 2000). The lower end of this range was prescribed by Reference Watts and EnglandWatts and England (1976) to alleviate the masking of desired signals by scattering from internal water bodies, and their work expanded that of Reference Smith and EvansSmith and Evans (1972) and Reference DavisDavis (1973) by considering realistic dimensions. Frequencies >100 MHz would increase the vertical resolution of bottom and internal structure and allow airborne, high-gain antennas (Reference Arcone and YankielunArcone and Yankielun, 2000) or ground-based antenna arrays (Reference Moran, Greenfield, Arcone and DelaneyMoran and others, 2000) to be used. The possible penetration of hundreds of meters of temperate ice at 100 MHz and above is strongly suggested by the 1400 MHz profiles to 60 m, and possibly 100 m, depths on Black Rapids Glacier, central Alaska, U.S.A. (Reference Arcone and YankielunArcone and Yankielun, 2000), and by the 100 MHz profiles to 190 m depth on Muir Glacier, southeast Alaska (Reference Arcone, Lawson, Moran, Delaney, Noon, Stickley and LongstaffArcone and others, 2000). These results, however, were achieved in ablation zones where crevassing appeared minimal and when snow cover was marginal. The purpose of this paper is to explore some of the limitations that the electrical structure of wet, stratified firn in accumulation zones offers to depth profiling in the 100–200 MHz bandwidth.
Radar profiles of glacier depth or stratigraphy within temperate accumulation zones at 100 MHz or higher have only been tried when the snowpack and firn was still cold (Reference Kohler, Moore, Kennett, Engeset and ElvehøyKohler and others, 1997). In late spring or summer, radar performance could be degraded by transmission losses through wet layer interfaces, and scattering losses from drainage structures and crevasses. These losses could become significant because of the strong contrasts between the refractive indices n of ice (1.8), water (9.4), wet snow and firn (1.8–6. 7) or air (1.0), and the likelihood of irregular interfaces. These indices represent the ratio of the speed of radio waves in air to that in the materials. They are not complex quantities near 100 MHz, so the relative dielectric permittivity, ε, is a real quantity and equals n 2. Scattering, manifested by hyperbolically shaped diffractions, is prominent where the sizes of the inhomogeneities (e.g. the width of a crevasse) compare with an in situ radar wavelength and where individual scatterers are spaced far enough apart to preclude multiple scattering. Unfortunately, the in situ wavelengths that might be needed near 100 MHz (1–2 m) to resolve firn strata in the thicker temperate accumulation zones may also compare with the widths of drainage features and crevasses.
The main objective of this paper is to derive refractive indices that occur in firn layers undergoing significant melting and drainage. By doing so, the transmission losses through the firn can be estimated. I derive these indices by analysis of airborne radar diffractions in profiles recorded over the Bagley Ice Field, southern Alaska, in early summer 1994, when melting was well underway, where yearly accumulation rates can exceed 3 m, and where extensive crevassing was caused by a then-present surge (see Appendix). I used helicopter-borne short-pulse radar with a pulse bandwidth centered near 135 MHz in expectation of both resolving firn layers and gaining some measure of penetration that appeared possible from earlier studies near this frequency (Reference Arcone, Lawson and DelaneyArcone and others, 1995). For the analysis I first use a high rate of stacking to suppress the diffractions and reveal the depth of the stratigraphy. The airborne-recorded diffractions themselves in the unstacked data are then used to derive refractive indices within the firn to about the maximum depths seen and then at much greater depth within a small basin. I verify my analytical approach with a medium of known index, that of water in a supraglacial lake; sampling and direct observation of the firn was precluded because it was not possible to land. I then compare profiles of crevassed and uncrevassed areas to examine causes of the diffractions.
My analysis of airborne-recorded diffractions appears new. Classically, dielectric permittivities are derived from diffractions for ground-based radar surveys by standard move-out analysis (Reference Jezek, Bentley and CloughJezek and others, 1979; Reference Clarke and BentleyClarke and Bentley, 1994; Reference Arcone, Delaney, Noon, Stickley and LongstaffArcone and Delaney, 2000; Reference Murray, Stuart, Fry, Gamble and CrabtreeMurray and others, 2000), but this approach has not been tried with airborne data, as I do here. With airborne data, dielectric properties of snow, ice and/or water layers have only been derived from the relative amplitudes of closely spaced reflections (Reference ArconeArcone, 1991).
Radar System
The radar comprised a Geophysical Survey Systems, Inc., SIR Model 10a+ control unit with real-time video display, a model 3107 prototype antenna transducer unit, cables, Exabyte digital tape recorder, and battery power supply. The model 3107 contains the transmitter (42 W peak power) and receiver electronics and separated antennas. The wavelet spectrum was centered near 135 MHz. The wavelet lasted about three cycles, of which the central 1.5 cycles are strongest. The 3 dB beam-widths of the two-way radar directivity patterns in the two principal radiation planes are each approximately 70° in air (Reference Arcone, Delaney and PerhamArcone and others, 1986), which narrows to about 38° after refraction into ice and slightly less for wet snow.
The antenna unit was carried in a cargo net about 6–14 m above the glacier surface (Fig. 1). The antenna elevation is discernible from the data. The antenna polarization never remained in a constant direction because the unit swiveled (yawed) during profiling but was never observed to swing. Clutter and indirect glacier reflections from the aircraft were detectable despite the built-in antenna shielding that is intended to reduce the intensity of radiation transmitted above the antennas. Such events will be identified in the profiles later.
Fig. 1. Depiction of model 3107 antenna transducer placed in a cargo net, and ray diagram of events recorded in a profile. T and R indicate transmitter and receiver antennas. The ray labels are explained in the text. The wedges depict crevasses.
Data Handling
Recording and positioning
Data were recorded at 25 traces s−1, a density of 1024 8-bit samples per trace, and a time range of 550 ns for the trunk profile; and at 12 traces s−1, 2048 16-bit samples per trace, and a time range of 1200 ns for the cross-glacier profiles. The trace-acquisition rate was limited by the sample and hold technology required to convert these high frequencies to the low-kHz audio range for 16-bit digitization. Time-variable gain (TVG) and high- and low-pass frequency filtering were applied before recording. The TVG suppressed the direct coupling between antennas, approximately compensated for the signal strength lost to spherical wavefront spreading, and kept all signals within the voltage level required for efficient digitization.
Electronic markers were recorded at 1 min intervals on the radar data, and waypoints were recorded on a Magellan (Model MX300) global positioning system (GPS) simultaneously. Helicopter (Aerospatiale 350A) speed (10–15 m s−1) was determined from the GPS positions and the trace recording rate, and antenna altitude (generally 6–14 m) from the time delay between the direct coupling between the antennas and the surface reflection. I determined the flight path from the GPS positions plotted on a topographic map (U.S. Geological Survey (USGS) 1 : 63 360-scale, Bering Glacier Quadrangle Series, 1985; based on aerial photos of 1972–73). These distances were then used to normalize the number of profile traces per kilometer to compensate for variations in flight speed. The positional accuracy of the GPS was about 30 m, but distance and average flight speed accuracy are estimated at 95% because of the long distances covered.
Profile processing
The processing alleviated unwanted reflections received directly from the aircraft, and prepared the profile for diffraction analysis and stacking. The first events in a profile (Fig. 1) are an artificial, start-of-trace (sot) signal wired directly from the transmitting antenna, T, to the receiving antenna, R, followed by the direct air coupling (dc) between T and R. The direct reflections (dr) from the helicopter were removed with a highpass horizontal filter because they occur at constant time delays across the profiles. This filtering did not significantly affect glacier surface (gs) and subsurface (gss) reflection and diffraction (gd) events because the changing helicopter altitude caused them to arrive at a non-constant time delay. The diffractions appear nearly hyperbolic on the radar profiles and can only be of glacial origin. The time-variable gs reflection was then normalized (i.e. the time delay was made constant within the record) with a seismic-type, statics removal algorithm (Reference SheriffSheriff, 1980). This process makes the dc horizon trace the aircraft altitude. Indirect glacier reflections (dashed line, gri) may include the ice surface and helicopter in their propagation paths. These events conform with either the surface reflection or the sot horizon after surface normalization, are virtually impossible to eliminate, and appear as coherent but unwanted reflections in the final profile.
Diffraction analysis
I derived effective thicknesses (d eff) and refractive indices (n eff) for layers defined between diffractions in the unstacked profiles by matching hyperbolic approximations of the profile diffractions with theoretical ones generated by layered models. I chose diffractions which had the steepest asymptotic slopes so that they probably originated beneath the antennas, and which also had time delays to their peaks that were similar to those of reflections revealed by the stacked profiles. Almost all diffractions along any horizon had the same slopes, which also indicates a steady flight speed. I approximated progressively deeper diffractions with one-layer model hyperbolas defined by a refractive index, n m, and a model depth, d m, both of which were determined by the hyperbola shape and the time delay to the hyperbola apex. The accuracy of this matching process was better than 2%. For example, a visual distinction can easily be made between a hyperbola generated for n m = 1.69 and one for n m = 1.71. Starting with the approximation of the first subsurface diffraction, I then used trial and error to generate a matching theoretical diffraction for a two-layer model of air over firn. This process then determined an effective refractive index, n eff, and corresponding thickness, d eff, for the firn layer above the diffraction apex. I computed the average of the squares of the differences between the time delays of the model and matching hyperbolas for many pairs of n eff and d eff and used the minimum for my solution. The average minimum difference, or error, between the approximation and the model for 30 (high indices) to 130 (low indices) distance points was always a very small fraction of a per cent for the best match. The pair values thus generated were then fed into a three-layer model, etc.
The matching procedure used the simple ray theory of geometrical optics. For the simplest case of a two-layer model (air above a snow layer of refractive index n 2), the time delay, t s, for a theoretical diffraction recorded in air from the bottom of the layer at a given offset is
where 
, d
2 is the depth of the second layer, h is the antenna altitude, s
1 + s
2 is the sum of the ray path lengths, and x
1 + x
2 is the total offset of the antennas from the scatterer (Fig. 2). The altitude delay is 2h/c, so that the statically corrected value, t
c, of t
s is then
Fig. 2. Ray-optics depiction of strong refraction towards horizontal above a high dielectric permittivity half-space. The rays exiting from the wet snow become horizontal at the critical angle θ c = sin −1 (1/n). For example, if n = 3.6, then θ c = 16°.
Equation (2) generates the diffraction that I matched in form to the hyperbolic approximation of the diffraction in the profile record for various values of n 2 and d 2.
Stacking, deconvolution and display
The purpose of the stacking was to reveal stratigraphy not seen in the recorded profiles and estimate the time range to the deepest firn reflections. The profile stacking is performed on the recorded profile, so it coherently adds events that occur at the same time delay, such as closely spaced diffraction peaks, and suppresses diffraction asymptotes. It also enhances the signal (S) to random noise (N) ratio according to the formula S/N = 10 log M (Reference Wright, Bradley, Grover and RedmanWright and others, 1994), where M is the number of traces stacked. I also applied minimum-phase predictive deconvolution (Reference SheriffSheriff, 1980) using a predictive lag equal to one cycle of the wavelet to improve horizon resolution for the stacked profile of the main trunk. Although ground-penetrating radar (GPR) wavelets are mixed phase, the algorithm effectively reduced horizon resonance by about one half-cycle. Traces showing wavelet improvement by stacking, however, are taken from the recorded data before deconvolution.
The profiles are displayed in a line-amplitude format that uses a non-linear gray scale that emphasizes weaker events and allows subtle stratigraphy to be viewed. White bands within an event indicate positive signal amplitude, and black bands indicate negative. Trace number is translated into profile distance. The echo time delay, t, in nanoseconds (ns) is translated into depth, d (m), with the simple echo-delay formula
where c is the speed of electromagnetic waves in a vacuum (0.3 m ns−1) and the factor of two accounts for the round-trip propagation path. Equation (3) applies to reflections from flat interfaces or to diffractions from point scatterers and edges oriented orthogonally to the profile direction. I use an effective n value for the profile depth scales, but I account for variations in n in my profile interpretations.
Temperate Firn
Structure
Temperate firn is highly stratified in density to a vertical resolution of centimeters and to depths of meters, with total depth usually depending on elevation (Reference SharpSharp, 1951; Reference Ambach and EisnerAmbach and Eisner, 1966). Firn generally has a dry density ρ = 400 kg m−3 near the surface, reaching an average value of about 700 kg m−3 (Reference SharpSharp, 1951; Ambach and Eisner, 1963) by a few meters depth. Within the percolation zone, firn is intermittently layered with wet and saturated snow and with ice. Vertically channeled drainage (Reference BensonBenson, 1996) can spread horizontally along less permeable, fine-grain snow (Reference ColbeckColbeck, 1974) or along existing ice layers to form saturation layers (Reference Oerter and MoserOerter and Moser, 1982). The water can then refreeze to form ice layers, or lenses, of a few meters width or less. Water undergoing normal gravity drainage within fresh snow will not exceed about 9% by volume (water volume at about 14% of the porosity) (Reference ColbeckColbeck, 1973); 3–4% is typical for drainage in dense firn at ρ = 700 kg m−3. The vertical drainage channels will freeze also and form ice glands (Reference BensonBenson, 1996). Drainage will occur around ice lenses and throughout firn when all the firn reaches 0°C (Reference MillerMiller, 1963) and will form a saturated firn aquifer (Reference FountainFountain, 1989; Reference SchneiderSchneider, 1999) above the impermeable glacier ice. Firn is a competent material throughout which crevasses may extend (e.g. Reference SharpSharp, 1988).
There are no reported observations of firn depth on the Bagley Ice Field; other information suggests that it could exceed 10 m at the elevations of these profiles. Annual rainfall equivalents reported at the nearest station, at Yakataga, on the coast, are about 400 cm. The transect elevations in the accumulation zone were about 1050–1400 m. Firn depths of 13 m (Reference SharpSharp, 1951) and 18–21 m (Reference MillerMiller, 1963) were measured at elevations of 1790–2130 m at Seward Glacier (128 km east of the mid-Bagley), and 34 m near 1250 m on Taku Glacier (528 km southeast) by Reference MillerMiller (1963). These larger values are consistent with experimental studies of the firn-to-ice transition by Reference Kawashima and YamadaKawashima and Yamada (1997). Reference PatersonPaterson (1994) remarks that the 13 m depth found by Reference SharpSharp (1951) appears to be exceptionally “rapid”.
Dielectric properties
The contrasts in dielectric properties between air, ice and water make firn a strong source of diffractions and reflections. The dielectric permittivity of dry snow is virtually lossless in the microwave region (Reference CummingCumming, 1952; Reference MätzlerMätzler, 1996) and will not change near 100 MHz. Cold water of low conductivity (<0.001 S m−1) is also a virtually lossless dielectric near 100 MHz because its relaxation frequency at 0°C is near 9 GHz (Reference Stogryn and DesargantStogryn and Desargant, 1985); I assume no lowering of this frequency from capillary tension because snow and firn grain-sizes are too large.
Since ε is real near 100 MHz for water, air and ice, I used the refractive index, 
. I interpret firn composition by relating neff to the volumetric fractions, v, and refractive indices of the firn components with the simple, empirical complex refractive index method (CRIM; Reference Annan, Cosgrove, Sigurdsson and RedmanAnnan and others, 1994) such that
where subscripts a, i and w represent air, ice and water, respectively, n i = 1.78, n a = 1.0, and n w = 9.38, that of free water at 0°C. Equation (4) then gives n eff = 1.60 for dry firn (ice plus air) at p = 700 kg m−3, which is very close to the data of Reference CummingCumming (1952). If only 1% of v a is replaced by water, then n eff = 1.7; for 6%, n eff = 2.1. The free-water value would occur for a water-filled crevasse, or a lake, as was observed in the ablation zone, and either water or air makes hidden crevasses a probable source of strong diffractions.
The expected n eff of near-surface snow or firn at saturation and for a water content of 9% is shown in Figure 3, along with values for other water contents. Reference ColbeckColbeck (1973) points out that saturation by imbibation (capillary tension only) falls about 8% short of volumetric filling. Even for ρ = 800 kg m−3 (v i = 87%), n eff = 2.8 at saturation because of the high value of water. The n eff values important to this study should cover a firn density of 500–800 kg m−3, which is the range covered by old snow to mature firn (Reference SharpSharp, 1951; Reference Ambach and EisnerAmbach and Eisner, 1966; Reference FountainFountain, 1989).
Fig. 3. Empirically derived dependence of the refractive index of firn upon ice density at different volumetric water contents.
The use of a mixing formula requires that the lesser components (e.g. air and water) be distributed within the matrix (ice) in units much smaller than the in situ wavelength of the whole medium. This criterion is satisfied given the grain-sizes that occur in firn. Equation (4) well predicts the permittivity of mineral, water and air mixtures at values to saturation, such as soil (Reference Topp, Davis and AnnanTopp and others, 1980). Empirical calibrations for 10 MHz–1 GHz and water fractions less than about 10% (Reference Tiuri, Sihvola, Nyfors and HallikainenTiuri and others, 1984; Reference DenothDenoth, 1989; Reference Kendra, Ulaby and SarabandiKendra and others, 1994) give values within 5% of the CRIM approach but unreasonable values at water contents approaching saturation for typical firn densities of 700 kg m−3.
Site Description and Transect Location
Figure 4 shows a simplified plan of the Bagley Ice Field, profile transects 1–3, and an elevation profile for transect 1 based on USGS 1 : 63 360 maps of 1985 (compiled from aerial photography of 1972–73). The profiles were recorded on 23 June 1994, during a surge of the Bagley Ice Field and Bering Glacier. Transect 1 traversed the main axis of the ice field, and the north–south transects extended over small basins and exposed bedrock ridges overlooking the west Bagley. Transect 1 traversed one of several supraglacial lakes within the confluence of the Bagley and the upper Bering trunk. Melting snow covered most of the ablation zone along transect 1, and extensive drainage networks fed the lakes. Pre- and postsurge elevations (Reference Molnia and PostMolnia and Post, 1995) suggest that the surge surface gradients should closely approximate the presurge condition; surge drawdown is estimated at 25–100 m (Reference Sauber, Plafker, Molnia and BryantSauber and others, 2000). The firn line was probably transported down-glacier about 0.5 km on the east side and less on the west, judging by the surge speeds measured by satellite interferometry (Reference Fatland and LingleFatland and Lingle, 1998). The snowline is higher on the east side, where the afternoon sun exposure is more favorable.
Fig. 4. Simplified plan for Bagley Ice Field with radar profile transects and elevation profile for the main transect. The dark areas are rock, and the small ellipse indicates one of several lakes. The average east-side surface slope is 0.6°.
Along the main trunk, regularly spaced, discolored creases in the snow marked subsurface crevasses; some crevasses were open, as discussed by Reference Lingle, Post, Herzfeld, Molnia, Krimmel and RoushLingle and others (1993) for the previous year”s surface. During profiling, the softer texture of the surface and the lack of any visible signs of crevassing marked the onset of the snowfields. However, crevasses revealed in the later USGS photography of the accumulation zones must have become exposed when snow bridges collapsed during the summer.
USGS vertical aerial photography of the Bagley Ice Field (ICP photographs 94V4-205 to 94V4-226, GeoData Center, University of Alaska Fairbanks), obtained on 7 September 1994 (Fig. 5), covers most of the area over which the profile was recorded. Unfortunately, the photography did not cover the westernmost 7.5 km of the transect. The flow direction along the east Bagley was nearly parallel with the transect, especially between 0 and 11 km, where I later interpret subsurface diffractions to be from crevasses. The flow direction along the west Bagley was nearly at right angles to the transect at 34 km, in the ablation zone, progressively decreased to about 30° by 42 km in the accumulation zone, and then was probably parallel to the transect by 50 km. The ablation zone contained a complex network of intersecting crevasses. The east and west Bagley 1994 snowlines appear to have intersected transect 1 at about 9 and 38 km, respectively.
Fig. 5. Reduced aerial vertical photography of Bagley Ice Field, with a simplification of the profile transect superimposed (top), and deconvolved 225-fold stack of the east–west profile (bottom). Bering Glacier heads off the southwest corner of the photograph, and Tana Glacier off the northwest. The large arrow in the profile indicates a surface multiple reflection, and the small vertical arrow at the east end of the horizontal scale locates the traces shown later in Figure 9. The profile vertical exaggeration is about 425 ×, and the depth scales are based on effective refractive indices. The interpolation of the transect location on the photographic distortion is not exact (north arrow is locally positioned).
Results
Transect 1: total profile
A 225-fold dynamic stack of the 50 km east–west profile is shown in Figure 5. The large ratio of elevation change to strata depth prevents elevation corrections (shown later) from being presented at this scale. The surface normalization process plus the zero time offset of the dc signal in the unprocessed profile leaves 400 ns of glacial data from the originally recorded trace durations of 550 ns.
The stack forms nearly continuous horizons from diffractions and masked reflections, theoretically provides 24 dB of random-noise suppression, and allows subtle strata and changes in strata dip to be seen. It also destroys the depth profile (shown later) of a supraglacial lake located at about 30.5 km. The radar strata appear to dip up-glacier because of the normalized surface. In reality, they dip down-glacier, as seen next. The depth scales are based on average values of refractive index that are calculated next. The stacking also reveals a horizon (large arrow in Fig. 5) that occurs throughout the profile and parallels the surface reflection. This horizon was received indirectly from a part of the helicopter, most likely the rotor, because it is revealed after surface normalization. Its removal with horizontal filtering would also destroy the stratigraphic horizons. Its occurrence indicates a wet, reflective snow surface.
Transect 1: east accumulation zone
Figure 6 shows a 45-fold stack of the east-side stratigraphy, along with matching photography and an interpretation of horizon depth based on the diffraction analysis discussed next. The depth scale of the radar profile is based on an effective value of refractive index for the first 233 ns (derivation given later). Horizon a is visible from at least 0 to 14 km and is probably the bottom of the seasonal snow. Horizons b and d begin between 12 and 13 km and deepen toward the east end. Horizon c cannot be followed farther west than about 7 km. The deeper d and e horizons are less defined than the others. Horizon e projects to the surface farther down-glacier, and its events extend from about 220 to 250 ns at the east end; it is about 13 m deep from 0 to 2 km distance. The scale of the line interpretation suggests that all horizons merge smoothly with the surface. However, the exaggerated vertical scale of the radar profile reveals that horizon e has a more unconformable contact with the surface, as any contact down-glacier from the equilibrium line should have in the ablation zone. All horizons project to the surface several kilometers down-glacier from the 1994 snowline seen in the photograph.
Fig. 6. A 45-fold stack of the deconvolved radar stratigraphy of the eastern side (middle) matching vertical photography with the profile transect superimposed (top), and stratigraphic interpretation (bottom). Labels are just beneath the horizons they represent. The west end of horizon b, about 4 km west of the snowline, is best seen by viewing the profile from the edge. The bright, erratic horizon at the bottom of the profile is an artifact of the surface normalization.
In contrast with the stacked profile, the unstacked profile (Fig. 7) is dominated by diffractions. The profile shows four of the five labeled horizons of Figure 6 but far less definition for horizons a, b and c. The diffractions occurring near the 4.5 km distance were analyzed because they occurred at the lowest helicopter altitude (10 m) on the east side and would therefore give best accuracy for values of n eff and d eff for the overlying layers. These diffractions peak at 85, 106, 153 and 233 ns delay and occur very near the surface of horizons b, c and d and within e of Figure 6, the latter interpreted to be near the deepest firn at this location. The five-layer (air is included as a layer) model derived from these diffractions to calibrate depths within the east accumulation zone is shown in Figure 8. The helicopter speed was 15.1 m s−1, so 1.97 traces represents 1 m of distance. The modeling yields equivalent homogeneous half-space values of n m = 1.08, 1.17, 1.31 and 1.46 for each deeper layer, with a margin of error for each value of about ±0.02. By matching, the n eff values for each of the layers are 4.1, 2.9, 2.1 and 2.1, and at successive thicknesses of 3.1, 1.1, 3.4 and 5.7 m for these time delays. I interpret these indices to represent a progressive decrease of water content with depth into the firn. The n eff and d eff values (4.1, 3.1) above horizon b also represent the layer above a, for which I could not find a solution. Similarly the values for the last layer (2.1, 5.7) represent firn from below horizon e that extends up to horizon d. An integrated, equivalent one-layer value of n eff = 2.5 results for the firn overburden above 233 ns, and I used it for the depth calibration for the profile in Figure 6. This value was derived from the total time delay and the total of the layer thicknesses.
Fig. 7. An unstacked profile segment from the east side near 4.4 km distance, where the antenna altitude is 10 m. The arrows indicate the analyzed diffractions (hyperbolically shaped events). Labels a–d indicate the same horizons seen in Figure 6.
Fig. 8. Hyperbolic approximations (solid lines) of profile diffractions from the east and west side segments from the main trunk profile, and two-to five-layer model matching diffractions (dashed lines labeled n m ) that generate values of n eff and d eff for each layer.
Figure 9 compares undeconvolved traces from the unstacked and 45-fold stacked profiles at the same location near 0 km. The surface reflection in the unstacked profile shows the form of the transmitted wavelet, which is reproduced after stacking. Within the top 150 ns the unstacked profile shows relatively weak subsurface wavelets, whereas the stacked profile shows clear events above the noise level. The phase progression of their half-cycles (+ − + −), which are numbered for comparison, indicates a decreasing refractive index with depth because it is opposite that of the surface reflection (Reference Arcone, Lawson and DelaneyArcone and others, 1995; Reference ArconeArcone, 1996) and so is consistent with the indices derived. This phase structure results from integration of diffraction peaks, as much as it might result from masked reflections.
Fig. 9. Comparisons between undeconvolved traces from the stacked and the unstacked profiles. The traces are located at the same distance on the east-side accumulation zone (indicated by the vertical arrow on the horizontal scale in Figure 5). Events sot (far left) and dc are distorted by the post-processing The numbers in the lower trace track the same half-cycles within events and reveal subsurface phase reversals relative to the surface reflection at 0 ns.
Transect 1: west accumulation zone
Figure 10 shows a 45-fold stack of the west-side stratigraphy, along with the limited matching photography and an interpretation of horizon depth based on the diffraction analysis. The depth scale of the radar profile is based on an n eff for the first 280 ns. Almost all the horizons are far stronger and better defined than on the east side. The profile shows no stratigraphy where the transect skirted the edge of the snow between about 38 and 40 km. The strata then begin within the more continuous snow coverage to the west at about 41 km. Between 41 and 45 km the more intense horizons appear more conformable with the surface than do those on the east side. Between 44 and 47 km an oblique view of the profile more clearly reveals that the deepest horizons extend to about 23 m depth and could be below the firn and into the ice at this elevation. There is also an apparent break in the stratigraphy at 49 km, where horizon b ends and b′ begins.
Fig. 10. A 45-fold stack of the deconvolved radar stratigraphy of the western side (middle), the available matching vertical photography with the profile transect superimposed (top), and stratigraphic interpretation (bottom). The labels are just beneath horizons they represent. Subtle stratigraphy can be seen to 23 m depth by viewing the profile obliquely.
The last 300 m were recorded at the lowest helicopter altitude, 6 m, along the whole profile. The unstacked version (Fig. 11) clearly shows only the a, b and c horizons. Horizon a, near the surface, mainly contains reflections, but there are also faint and very wide diffractions that give index values <1. This phenomenon is generated by linear features obliquely oriented with respect to the profile line (Reference Clarke and BentleyClarke and Bentley, 1994), so these events are probably hyperbolically shaped reflections from within the seasonal snow. They appear too weak to be from crevasses. The diffractions analyzed occur at 105, 161, 233 and 280 ns delay, with the latter interpreted to be near the deepest extent of the firn at this location. The first two time delays correspond with horizons b and c, but the latter are not associated with any apparent horizon in Figure 11. The aircraft speed here was 10.0 m s−1, and 2.77 traces represents 1 m of distance. Modeling of these diffractions yields n m = 1.44, 1.64, 1.64 and 1.68 (Fig. 8). The n eff values for each of the layers are 4.5, 2.5, 1.7 and 2.0 and at successive thicknesses of 3.5, 3.3, 6.4 and 3.7 m. This puts horizon a at a depth of about 1 m at the west end. As on the east side, I interpret the first three indices to represent a progressive decrease of water content with depth into the firn. I obtained an equivalent one-layer value of 2.3 for the firn above 280 ns and used it for the depth calibration of the profile in Figure 10.
Fig. 11. An unstacked profile segment from the west side near the 50 km distance, where the antenna altitude is 6 m. The arrows indicate the analyzed diffractions. Labels a−c indicate the same horizons seen in Figure 10.
Transect 1: supraglacial lake
A supraglacial lake of about 0.5 km2 (Fig. 4) occurred near 31 km along the main trunk profile. All of the supraglacial lakes drained after the outburst flood of 27 July 1994 (personal communication from B. F. Molnia, 1994]). The white area in the photograph of Figure 12 is interpreted to be remnant frozen snow in part of the drained lake basin. In the lake profile (Fig. 12) I normalized the first faint reflection surface to provide the datum plane; the wavy character of the stronger near-surface reflections shows that surface skims of ice were present. The profile shows two apparent lake-bottom horizons. The deeper and stronger horizon (b) is interpreted to be the ice bottom with a maximum depth (Equation (4)) of about 3.8 m. This horizon is composed of diffractions (Fig. 13), which I interpret to be responses to water-filled crevasses. The weaker horizon a occurs earlier and is interpreted to be the top of the (submerged) snow layer because it is continuous with the snow layer reflections along the east and west approaches to the lake. Its weaker strength is caused by a weak contrast between the permittivity of the saturated snow and that of the water above. At each end of the profile, horizon a merges with the surface reflection and cannot be resolved. There are also apparent sub-bottom diffractions (events c) whose depth can be determined from the left-side depth scale for ice. All other events except the surface reflection are multiple reflections between the helicopter and the highly reflective surface (Figs 1 and 13).
Fig. 12. Profile of the supraglacial lake, and aerial photograph of 7 September with the profile transect superimposed. An n = 9.38 for water at 0°C was used to calibrate lake depth. The labeled horizons and events are interpreted in the text.
I tested my diffraction analysis with an ice-bottom diffraction because the refractive index of water at 0°C (9.3–9.40) is well known. I used a bottom diffraction from an unstacked segment of this profile (Fig. 13) to derive a refractive index for the composite structure of water and saturated snow. The saturated snow is assumed to have ρ = 500 kg m−3 (n = 5.2). The peak of this diffraction occurs after a 207 ns delay in the water plus a 30 ns delay in the saturated snow. The antenna altitude above this test diffraction is 18 m (and ranges to 20 m above the lake), which is unusually high for the entire transect but provides a further test of the analysis. The minimum in the average least-squares error between my hyperbolic approximation to the profile diffraction and my model diffraction occurs at n eff = 9.9, with erratic behavior between 9.3 and 10.3. This result compares fairly well with the n eff of 8.8 that I calculated from the refractive indices given above for water and saturated snow and from the time delays in these materials above the diffraction. The average errors are very high compared with those discussed next for a different case but at the same altitude. Therefore, the discrepancy is probably caused by the matching procedure and is discussed later.
Fig. 13. Detail of the lake-bottom profile (top) and a least-squares error analysis (bottom) to derive the refractive index and depth of the water above the diffraction indicated by the arrow. The water surface reflection is normalized to 0 ns. The antenna altitude above the surface is profiled above the water-depth profile. The pilot increased altitude over the lake for fear of immersing the antennas. The three nearly horizontal events at about 1, 2 and 3 m below the surface reflection are indirect reflections between the water surface, antennas and helicopter.
Transects 2 and 3: peripheral snowfields
Transect 2 started at the west end of transect 1, crossed the west-side accumulation zone and then crossed three snow basins above the main trunk (Fig. 4). Transect 3 crossed two other basins. All the basins appeared unaffected by the surge because of the lack of near-surface diffractions in their profiles and no visible crevasses during profiling. Although the photography did not cover these basins, no crevasses can be identifed in any of the adjacent basins farther east. The transect profiles are shown in Figure 14. The slower recording rate (12 traces s−1) and short lengths of the transects prevent stacking from effectively improving the stratigraphy across the main trunk, which extends to 2.2 km of transect 2. I interpret the deep horizons to be bedrock because they all culminated in a ridge or bedrock exposure. The lack of prominent events below the near-surface layering and above the bedrock indicate nearly homogeneous conditions.
Fig. 14. South–north profiles, which traversed peripheral snow basins above the western Bagley. Events A are the traces of the direct coupling, which now give the antenna altitude after surface normalization. The glacier surface is the faint horizontal line at 0 ns. Events B are diffractions from the bedrock. Events C are near-surface layering
In contrast with the main trunk profile, the major nearsurface horizons within five profiled basins are composed of continuous reflections, while the deepening subsurface horizons are formed mainly by diffractions (Fig. 15). The lack of diffractions within the layering and the similarity of the reflected wavelets to that of the surface reflection (Fig. 16) indicate sharp transitions in n. The horizon nearest the surface in Figure 16 shows a phase reversal relative to that of the surface and the other reflections. This reversal indicates lower-n material below the horizon and higher-n material above.
Fig. 15. Profile section (top) within transect 3 that crossed a snow basin, and elevational perspective (bottom) based on the 1972 USGS topographic contours. The bedrock is defined by the sloping horizon of diffractions. The depth scales are based on an analysis of the diffractions indicated by arrows.
Fig. 16. Detail of near-surface layering of the basin profiled in Figure 15. A phase reversal, relative to that of the surface reflection, in the first subsurface horizon (arrow) at about 5 m depth indicates the top of a lower-density layer. The phase reversal can also be seen in the trace (wavelet marked a), which was recorded at the position marked by the vertical arrow beneath the profile.
I analyzed the two diffractions indicated by arrows along the slope in Figure 15. The peak of the earlier one is centered at trace 443, where the antenna altitude was 17 m, and occurs at a time delay of 166 ns. The asymptotes are visible to 20 m horizontal distance from the peak, and it is the only discernible diffraction along any of the smooth reflections. Modeling yielded a sharp null for the average least-squares fitting error (Fig. 17) at the values d = 8.9 m and n eff = 2.80. The peak of the later diffraction, centered at trace 283 at an antenna altitude of 14 m, occurs at a delay of 838 ns. Its asymptotes are visible to 47 m from the peak. It is one of many diffractions at depth, all of which my analysis gave similar values of n eff. Using the value n eff = 2.8 for the 5 m (93 ns) of stratified overburden at this distance, the solutions for a second layer (Fig. 17) are d eff = 43 m and n eff = 2.6, which give a total depth at this point of 48 m. A one-layer equivalent model at this location gives n eff = 2.7, which is consistent with the slightly higher values determined for the near-surface layers and suggests that the n eff = 2.8 value at trace 443 also applies to the layers at trace 283. Using these refractive indices for calibration, the maximum depth reached is about 58 m. At this depth the signal-to-noise ratio is about 1–2 dB, which implies that very few meters of further penetration would be possible to interpret.
Fig. 17. Average least-squares errors between approximations to the profile diffractions indicated in Figure 16 and layer models of these diffractions for combinations of n eff and d eff . Values of d eff are indicated along the curves. The curve derived for trace 443 is for a one-layer model; that for trace 283 is for the bottom layer of a two-layer model whose top layer has n eff = 2.8.
Discussion
Stratigraphic interpretation
There is no obvious stratigraphy after about 150 ns in the unstacked sections of Figures 7 and 11, but it does exist to >220 ns in the same stacked sections and to >300 ns elsewhere in the stacked profiles. Therefore, the stacking significantly improves the stratigraphic picture and shows an approximation of the classical conception of snow and firn structure on a temperate glacier (Reference MüllerMüller, 1962; Reference PatersonPaterson, 1994; Reference BensonBenson, 1996). In this structure, by mid- to late summer the previous winter’s snowpack melts back to, or just above, the equilibrium line, above which a wedge of snow remains that thickens up-glacier. At an elevation much greater than reached by these profiles, the limit of melting, or the late-summer dry snowline, occurs. By the following year this wedge of snow (now firn) has moved down-glacier, and a new wedge forms near the same snowline. Thus, after several years a series of layers has formed, each of which, in earlier succession, intersects the ablation zone surface farther down-glacier from the snowline.
When Figures 6 and 9 are viewed obliquely, it is apparent that as horizon a is followed down-glacier it moves toward the surface. Therefore, I interpret a to represent the bottom of the 1993/94 snowpack, which still covered the ablation zone during profiling but cannot be resolved from the surface reflection there. On the east side there is no consistent down-glacier progression of horizon intersections with horizon a. Horizon b intersects a at about 12 km, c fades out by 7 km, d merges with b at about 12 km, and e projects to the surface beyond b. On the west side the down-glacier progression better conforms with the classic picture. There, horizon b appears to intersect a at about 46 km, and all succeeding horizons extend progressively farther down-glacier to about 41 km.
I have no reason to interpret the deeper horizons to represent annual ablation surfaces. However, their continuity and near-conformity with each other imply that they represent isochronal processes. Given the large contrasts in n that can develop between saturated and unsaturated layers, they may represent nearly saturated fine-grained layers. Stratigraphic breaks or discontinuities, such as the end of horizon c at 7 km (Fig. 6) and the transition from horizon b to b′ at 48 km (Fig. 9), show that these horizons may not be continuous on the scale of about 10 km. In addition, the apparent resonance seen in the deeper horizons shows that a higher-frequency pulse is needed to improve their resolution.
Firn refractive indices
Refractive index values above that of pure ice, 1.78, indicate that water was present. For the unlikely case of total saturation (Fig. 3), n should vary between 3.6 (ρ = 700 kg m−3) and 6.7 (ρ = 300 kg m−3). The highest value I calculated is 4.5 for the top layer on the west side of transect 1. From Figure 3 a saturated top layer of ρ = 400 kg m−3 would have n = 6; an n = 4.5 at 400 kg m−3 gives a v w = 0.33 using Equation (4). The one-layer n eff values are 2.5 to 13 m depth on the east side of the main trunk and 2.3 to a depth of 17 m on the west. Therefore, it is unlikely that a significant part of the firn was at saturation. For practical comparison, the average porosity of firn of temperate South Cascade Glacier, Washington State, U.S.A., in late summer was measured to be 0.15, of which about 60% was water (Reference FountainFountain, 1989; Reference Fountain and WalderFountain and Walder, 1998). These measured values give ice, air and water volumetric contents of 0.85, 0.06 and 0.09, respectively, and a theoretical refractive index of 2.4 (Fig. 2), which agrees with my one-layer equivalent values. Similarly, water contents of 0.04–0.09 (equivalent volumetric) were determined for firn at a snowpatch site in Japan (Reference Kawashima, Yamada and WakahamaKawashima and others, 1993).
The higher values near the surface of the main trunk are consistent with the extensive supraglacial drainage network that led through the accumulation area to the lakes. They also suggest possible saturation sub-layers. The lack of any especially strong reflection or high value of n at depth suggests that a water table had not yet developed within 13–17 m depth. The value n eff = 2.6 for the deeper material of the basin is surprisingly high. It can be achieved with ρ = 820 kg m−3, v i = 0.89 and v w = 0.11, which implies that the deeper material was dense firn at saturation, whereas nearice index values were reached by only 4–7 m depth on the main trunk. Since firn at this density should be impermeable, the water would result from internal melting and not drainage. Reference NicollinNicollin and others (1992) measured temperate-ice wave speeds in the Western Alps that give n = 2.2, and Reference Murray, Stuart, Fry, Gamble and CrabtreeMurray and others (2000) meased englacial speeds in Iceland that give n = 2.0 using an offset survey. Murray and others’ surface method of measuring speed is certainly more exact than my airborne approach. As discussed below, my approach could cause an 11 % overestimate in n.
I have computed n values for firn at only two locations along the main trunk and one location in a snowfield. Certainly there was variation for other locations, but the computation process is not automated and so is labor-intensive for any one location. Nevertheless, the n values on the east and west sides of transect 1 are similar in value and structure, considerably above that of ice nearer the surface, and approach that of ice at depth. In addition, my choice of diffractions where helicopter elevation was lowest assures locations where accuracy should be best.
Two important sources of error in the diffraction analysis are the accuracy of the aircraft speed and in finding a hyperbola to match the profile diffraction. The GPS measurements appear accurate to 30 m, and the pilot held steady speeds over all sections analyzed so that the conversion of traces to distance must be accurate to within a few per cent. A 5% error in aircraft speed is a 5% error in distance determination and would result in a theoretical 5% error in n eff for a hyperbola at 15 m depth and matched to ±20 m. Matching of successive surface hyperbolas over a section 2–3 km long at near-constant altitude showed standard deviations in n m of only about 2%, which indicates the consistency of the air speed. The matching of profile diffractions appears to be the largest source of error because both the lake-bottom diffraction (Fig. 13) and the shallower basin diffraction (Fig. 15) were recorded at nearly the same altitude (17–18 m), yet the lake-bottom result shows a greater least-squares matching error, and the result for n appears to be in error by about 11%. The cause is clearly the abbreviated extent of the lake-bottom diffraction, which makes it hard to match and is caused by the large contrast in n between water and air. Given that this is the strongest contrast possible, the maximum estimated accuracy is then about 11%.
Also shown in Figure 8 are the horizontal distances (vertically oriented, finely dotted lines) over which I was able to match the profile diffractions. The longer distances occur on the west side, where deeper hyperbolas were analyzed. It is apparent that the layer model diffractions deviate significantly from the profile approximations beyond these distances and always have less time delay for a given distance from the origin. Therefore, a more sophisticated, multi-layer algorithm for generating exact matches to the profile diffractions, instead of using model hyperbolic approximations, would improve accuracy in matching.
Given the contrasts in n between air and firn, and within the firn, and assuming transmission into solid ice beneath the deepest firn layer along the main trunk for which I calculated a refractive index, the corresponding two-way transmission losses are only 4.5 and 5.7 dB on the east and west sides, respectively, as calculated using the plane-wave, Fresnel interface transmission coefficients (Reference WaitWait, 1970). Most of this loss occurs through the surface. In compensation, the added antenna gain caused by refractive focussing across the surface (Reference Bogorodsky, Bentley and GudmandsenBogorodsky and others, 1985) is 4–5 dB. This means that reflections from these layers of wet firn should not have limited penetration severely, and this is borne out by the basin profile. However, a saturated layer at the firn–ice transition would have severely limited penetration (Reference Jacobel and AndersonJacobel and Anderson, 1987).
Origin of diffractions
The high values of refractive index that occurred in the nearsurface stratigraphy of the small basins imply that drainage structures such as vertical conduits and ice lenses should have been present. However, the lack of diffractions within these horizons suggests that such structures were not a strong cause of scattering at 100–200 MHz. Therefore, the diffractions along the main trunk probably originated from buried crevasses and, in particular, from the bottom of their snow bridges because their integration by stacking forms horizons.
The existence of so many interpretable, hyperbolic-type diffractions shows that single scattering, as opposed to multiple scattering, prevailed. Hyperbolic diffractions would not occur for multiple scattering because the scattering-path lengths are then erratic and not necessarily directed back to the receiving antenna. Single scattering makes no requirement on the size of the scatterers, only that they be spaced several wavelengths apart (Reference Van de HulstVan de Hulst, 1957). This is consistent with shallow and wider snow-bridge origins rather than with deeper, convergent crack origins for the diffractions. It is also consistent with crevasse spacing, which appears on the photography to be no closer than about 20 m on the east Bagley. In addition, all the analyzed diffractions must have originated from either point sources beneath the antennas or edges oriented orthogonally to the flight-line for such reasonable values of n to be obtained; off-center diffractions quickly widen with distance from the transect and give unreasonably small values of n. This orientation is reasonable since I chose my diffractions from the far eastern and western ends of the transect, where crevasse strikes were orthogonal, or very nearly so, to the transect.
Conclusions
The main losses near 100 MHz must be caused by scattering because transmission losses through firn interfaces do not appear severe. The loss process appears to be single scattering, manifested by the interpretable diffractions, but I have no way of estimating this loss because the inhomogeneities are oddly shaped and not small compared with a wavelength. The masking of interfaces by so many diffractions validates the hypothesis of Reference Watts and EnglandWatts and England (1976), but for firn and not just englacial ice. The appearance of distinct and deeper horizons after stacking mainly resulted from constructive interference between diffraction peaks and destructive interference between diffraction asymptotes.
The reasonable agreement between actual and calculated refractive indices for the supraglacial lake, and the reasonable values obtained for wet firn, demonstrate that airborne diffraction analysis can be used to estimate the depths to firn horizons and the water content of the layers. The values also imply that about 9% might be a good estimate of the maximum average water content of draining firn. Newer radar systems employing real-time digitization of the received radio-frequency energy will have far greater trace acquisition rates. Given better differential GPS control on flight speed and higher-trace-density profiles, more exact multi-layer algorithms could be used to improve modeling and calculation of refractive indices.
The penetration of about 58 m achieved in the basin, and the low firn transmission losses calculated for the main trunk, lead to the conclusion that far greater penetration could be achieved in uncrevassed areas with higher-gain 100–200 MHz antennas, high rates of stacking, and low-noise, real-time digitization. With the system used here, it appears possible to profile to the depth of the deepest firn in glaciers with large annual accumulations, especially if they are not severely crevassed.
It appears difficult to interpret any of the horizons below the first as representing annual surfaces because melting and refreezing events within any one year can cause reflecting horizons. The lack of continuity seen in a few horizons suggests that annual changes in the elevation of the snowline along any chosen transect would further hinder this possibility. It does seem reasonable, however, that annual accumulation for a previous year could be measured at the end of winter by following the snow–ice horizon along the ablation zone and into the accumulation zone.
Acknowledgements
Support for this research was provided by a subcontract to U.S. National Science Foundation grant OPP9319872, from the State University of New York (SUNY) at Oneonta, and by the U.S. Army Cold Regions Research and Engineering Laboratory (CRREL) Project 4A24161102AT24, work unit AT24-SS-E10. I gratefully acknowledge the field support of P. J. Fleisher of SUNY Oneonta, Colonel P. K. Bailey of Homer, AK, and A.J. Delaney of CRREL, and the comments of H. Rott (Scientific Editor), J. Kohler and an unnamed reviewer, which greatly improved the paper.
Appendix the Bering Glacier Surge of 1993–95
The Bering Glacier–East Bagley Ice Field is the longest glacier system in North America. Bering Glacier has undergone several major surges in recent history (Reference Muller and FleisherMuller and Fleisher, 1995). The most recent began about May 1993, ended around July 1994 and then continued again in 1995 (Reference FleisherFleisher and others, 1995, Reference Fleisher, Cadwell and Muller1998; Reference Muller and FleisherMuller and Fleisher, 1995). I obtained these profiles during the late stages of the first surge, which also affected the Bagley Ice Field (Reference Herzfeld and MayerHerzfeld and Mayer, 1997; Reference Fatland and LingleFatland and Lingle, 1998; Reference Sauber, Plafker, Molnia and BryantSauber and others, 2000). Thus, because of the season, there was a crevassed surface beneath the annual snow layer and, because of the year, there was displacement of the historical accumulation zones down-glacier, especially on the east side of the Bagley Ice Field (Reference FatlandFatland, 1998).
I acquired my radar profiles in an attempt to establish the extent of hidden crevassing and thus to determine the extent of the surge along the Bagley Ice Field (Reference Arcone, Delaney and FleisherArcone and others, 1994). I chose the 100–200 MHz bandwidth as a compromise between higher frequencies such as 500 MHz, which resolve crevasse structure (Reference Arcone, Delaney, Noon, Stickley and LongstaffArcone and Delaney, 2000) but might not penetrate wet snow, and lower frequencies, which increase penetration but might not resolve crevasse structure. Although my main profile is 50 km long, this proved inadequate because interferometric images obtained by satellite-based synthetic aperture radar provided a continuous profile of axial speed over a far greater distance (Reference FatlandFatland, 1998; Reference Fatland and LingleFatland and Lingle, 1998). In addition, the late-summer photography cited in the text revealed many crevasses in the accumulation zone and farther east than my profile extended. Therefore, the crevasse diffractions themselves and their transformation into continuous reflections motivated this study.
