Joint interpretation of explosive and vibroseismic surveys on cold firn for the investigation of ice properties

Abstract

Two seismic surveys were carried out on the high-altitude glacier saddle, Colle Gnifetti, Monte Rosa, Italy/Switzerland. Explosive and vibroseismic sources were tested to explore the best way to generate seismic waves to deduce shallow and intermediate properties (<100 m) of firn and ice. The explosive source (SISSY) excites strong surface and diving waves, degrading data quality for processing; no englacial reflections besides the noisy bed reflector are visible. However, the strong diving waves are analyzed to derive the density distribution of the firn pack, yielding results similar to a nearby ice core. The vibrator source (ElViS), used in both P- and SH-wave modes, produces detectable laterally coherent reflections within the firn and ice column. We compare these with ice-core and radar data. The SH-wave data are particularly useful in providing detailed, high-resolution information on firn and ice stratigraphy. Our analyses demonstrate the potential of seismic methods to determine physical properties of firn and ice, particularly density and potentially also crystal-orientation fabric.

Introduction

The physical and chemical properties of glacial ice constitute a paleoclimate archive on various spatial scales. Moreover, ice properties also represent a memory of past and present ice dynamics. Investigations of ice properties generally involve drilling an ice core (e.g. Cuffey and Paterson, 2010). This is, however, a labor- and cost-intensive task that captures only point information across a glacier or ice sheet. Geophysical methods, such as radar and seismics, allow remote imaging of the englacial properties and the ice/bed interface over larger lateral distances in significantly less time. Ground-penetrating radar (GPR) measurements, as well as radio-echo sounding, have long been used to trace englacial ice stratigraphy, map ice thickness and deduce basal conditions (e.g. Plewes and Hubbard, 2001; Dowdeswell and Evans, 2004; Navarro and Eisen, 2009). More recently, it has been possible to investigate properties of the crystal-orientation fabric (COF) by radar methods (e.g. Fujita and others, 2006; Eisen and others, 2007; Matsuoka and others, 2009). The sensitivity of seismic waves to variations in wave speed and density within and beneath the ice yields complementary information to that of radar methods, and it is possible to perform a joint analysis of these two geophysical methods in mapping englacial and subglacial properties.

Active seismic methods have been utilized for decades to investigate basal topography and lithology, as well as to map the relationships between seismic velocity and density, temperature and COF (e.g. Robin, 1958; Bentley, 1972; Kohnen, 1972, 1974; Peters and others, 2012). The sensitivity of seismic waves to density variations yields the opportunity to map the structure of the firn pack (Kirchner and Bentley, 1990). Diving waves can be used to derive the density distribution by inversion of the travel-time/offset data (Slichter, 1932; Kirchner and Bentley, 1990). Further, body waves can be used to derive the distribution of the Poisson’s ratios of firn (King and Jarvis, 2007). Seismic reflections can be caused by changes in density or elastic properties, especially COF (Bentley, 1972; Blankenship and Bentley, 1987; Horgan and others, 2008, 2011). Hence, a change in seismic velocity due to a change of orientation in the bulk of the seismically anisotropic ice crystals can lead to reflections from within the ice column. Further investigations include the examination of the damping behavior of seismic waves in ice (Gusmeroli and others, 2010; Peters and others, 2012), investigations of the conditions at the bed for normal incidence (Smith, 1997; King and others, 2004; Smith and Murray, 2009) and the employment of amplitude-versus-offset (AVO) analysis (Anandakrishnan, 2003; Peters and others, 2007, 2008; Booth and others, 2012).

To date, seismic surveys on ice have usually been carried out using explosives, either by drilling holes and deploying charges or by using a detonation cord on the surface (Sen and others, 1998; Johansen and others, 2011). In addition, sources like hammer blows or airguns have been used (e.g. King and others, 1993). In the presence of a thick layer of porous firn, much seismic energy is generally lost in the first few meters of the firn, due to inelastic deformation around the shot location, as opposed to the excitation of elastic seismic waves. Moreover, a considerable amount of energy goes into the excitation of surface waves, often dominating seismic shot records on ice, due to their dispersive nature and high amplitudes. Here a considerable amount of seismic energy is also trapped in the firn column, which acts as a waveguide (Albert, 1998). This effect can be reduced efficiently by drilling holes to reach a depth of higher density within the firn pack; however, this is a time- and energy-consuming task.

Based on the success of using a heavy vibrator as a surface source on an Antarctic ice shelf (Eisen and others, 2010), the survey presented here was set up to investigate the effectiveness of so-called micro-vibrators on shallow, firn-covered ice bodies. A vibrator has several advantages compared with explosives. It can be moved quickly to the next shot location, without the labor- and energy-intensive task of drilling shot holes. The vibrator sweep mostly excites elastic energy, where the energy is distributed over a longer time than with explosives, thus creating smaller instantaneous forces and much less spurious noise. Moreover, the source signal’s frequency characteristic is well known, providing better reproducibility of the shots for stacking.

We present seismic results from two different sources tested on Colle Gnifetti, Swiss Alps, to investigate the potential for joint seismic and radar analysis in detecting ice properties. The seismic impulse source system (SISSY) had limited potential for resolving englacial and clear bedrock reflections. However, strong diving waves allowed a detailed density analysis of the firn pack. The electrodynamic-vibrator system (ElViS) yields a much higher signal-to-noise ratio (SNR), clearly showing englacial stratigraphy as well as bedrock reflections in the processed data. Together with the study by Hofstede and others (2013), which used a heavy vibrator on a local ice dome in Antarctica, we show that seismic vibrator sources with peak forces ranging over two orders of magnitude are efficient alternatives to explosive sources for detecting englacial reflections, given sufficiently high resolution. Altogether, our seismic analyses highlight the ability to utilize seismic data, in the absence of ice-core or radar data, as a means of deriving firn and ice structure, which could potentially be related to firn and ice dynamics.

Field Site: Colle Gnifetti

Colle Gnifetti is a glacier saddle in the Monte Rosa massif on the border between Switzerland and Italy. The plateau lies at ∼4450 m a.s.l. between the summits Dufourspitze and Zumsteinspitze. In comparison to other extensively investigated drilling sites in the European Alps, Colle Gnifetti has very low net snow accumulation caused by wind-driven snow erosion. The net surface accumulation rates in its north-facing flank range between ∼15 cm w.e.a^–1 close to an ice cliff on the eastern side and 50 cm w.e.a^-1 near the saddle point (Alean and others, 1983). The englacial temperature distribution typically lies between -14 and -10°C at 20 m depth and - 13 and -12°C at the bedrock (Hoelzle and others, 2011), which ensures preservation of the stratigraphy over the entire depth range. As Colle Gnifetti falls into the recrystallization/infiltration zone classification (Shumskiy 1964), some typically centimeter-thick melt layers are formed occasionally. Details on the glaciological features of Colle Gnifetti are reported by Haeberli and Alean (1985), Haeberli and others (1988), Eisen and others (2003) and Schwerzmann (2006). GPR measurements revealed an ice thickness of up to 140 m (Haeberli and others, 1988; Wagner, 1996; Luthi, 2000; Eisen and others, 2003). Geodetic surveys and photogrammetric comparisons suggest that the glacier was in a near-steady state throughout the 20th century (Wagner, 1996; Luthi, 2000).

A number of ice-core drillings have been performed on Colle Gnifetti (Schotterer and others, 1978; Wagenbach and others, 1988, 1997; Gabrieli and others, 2011), providing ice-core records on a centennial timescale. In August 2005, a core (KCI) was drilled down to 62 m depth in an area of particularly low accumulation, close to the ice cliff (Bohleber, 2011). The bed was probably not reached, as drilling stopped when the first dirt inclusions occurred in the ice. It is suggested that the glacier bed is likely to be up to 1 m deeper (Bohleber, 2011). The borehole is located at 45.929728 N, 7.876928 E (World Geodetic System 1984 (WGS84), measured in August 2008). Our seismic survey area was centered on this drill site (Fig. 1). We thus have the opportunity to directly compare ice-core and seismic data. A density profile was measured along the KCI ice core at sub-centimeter resolution by means of 7-attenuation (Wilhelms, 1996). These data are later used for comparison with the results of the seismic data, as well as to connect seismic, GPR and ice-core data.

Fig. 1. Acquisition geometry for the SISSY survey at Colle Gnifetti, August 2008. Line 1 is directed from south to north with 34 shots (S1_NS–S34_NS). Line 2 is directed from west to east with 33 shots (S1_EW–S33_EW). Geophone and source locations are marked with triangles and crosses, respectively. In both cases, 24 geophones were placed at 3m distance intervals (G1_NS–G24_NS and G1_EW– G24_EW). The center point (C) is located close to the borehole of the KCI ice core that was drilled in 2005. The inset shows the location of Colle Gnifetti, marked as a filled circle.

Data Acquisition and Processing

Explosive survey

A seismic survey was carried out on Colle Gnifetti in August 2008 to investigate ice properties, particularly the firn-density distribution, and the feasibility of detecting changes in the COF. Two seismic profiles were collected (Fig. 1). Line 1 is located in a north-south orientation along the flow of the glacier, and line 2 in an east-west direction, positioned transverse to the glacier flow. The center point (C) of the two seismic lines is 5.9 m north-northwest of the KCI borehole, providing the opportunity to directly compare our seismic results with the ice-core data.

For each line, 24 vertical geophones (SENSOR SM 4, with a natural frequency of 14 Hz) were placed with a 3 m spacing. Twenty-three shots were detonated between the geophones, while another eight shots, four on each side of the geophone array, were placed with the same shot spacing (3 m) as within the geophone array. Three more far-offset shots were detonated along line 1, all at a 25 m spacing from the neighboring shot (Fig. 1), resulting in a maximum offset of 95 m to the south and 70 m to the north from point C. Only two far-offset shots were detonated along line 2, both at a 25 m spacing and to the west of the central shots (Fig. 1). The positioning of these seismic sources in relation to the geophone lines resulted in a 24-fold survey at point C. The shots were detonated with the SISSY (Druivenga and Fibranz, 2006). The SISSY is a steel tube 1.26 m long, with a total weight of 10 kg, making it easy to handle; a lightweight source, suitable for remote areas, such as Colle Gnifetti. For the shot a Dynergit cartridge, containing 28g of explosive on NC-basis (nitrogen carbon), is screwed into the lower part of the tube, which is then lowered into boreholes of ∼0.5 to 1 m depth. Data were recorded with a Strataview R24 seismograph, with a sampling interval of 30 μs to provide high resolution and a record length of 256 ms.

The SISSY dataset was processed to remove surface-generated noise, in order to investigate englacial and bed reflections. With a source depth of only 1 m, strong surface waves were generated in the frequency range 20-240 Hz. They are visible in all shots, both within and outside the geophone line (Fig. 2, feature B). The first arrivals are near-surface diving waves (Fig. 2, feature A), whose ray paths are shaped by the strong density gradient of the firn pack down to 40 m depth (as observed in the KCI ice core). Further noise in the data occurred due to parasitic resonances of the geophones at near offsets. Due to the strong surface waves and the small offsets, the vertical geophones also have a movement in the horizontal direction when the signal is passing (Faber and Maxwell, 1997). This generates a continuous noise of 205-208 Hz along the near-offset traces, which is expressed by a continuous ringing over the complete recording time (Fig. 2, feature C).

Fig. 2. An illustration of general features from raw data of the SISSY survey, line 1 . Shot 23 is a shot within the geophone line; shot 34 is a far-offset shot with an offset of 25 m to the next geophone. The gain scaling of the two shots is different, chosen such that the amplitudes of the different events are about equally visible. Some dominant features of the raw shots are marked: A – strong diving waves; B– surface waves;C – geophone ringing due to spurious frequencies.

Surface waves, diving waves and geophone ringing are overlaying potentially weak englacial reflections and the bed event at ∼0.04s two-way travel time (TWT), in the time (Fig. 3) as well as the frequency domain. We tested different filters before stacking the data. Frequency filters were not a great success, because the surface and diving waves possess frequency spectra comparable to that of the bed reflection. Some improvements could be achieved by using a frequency/wavenumber fk-filter on the far-offset shots, as the surface wave and geophone noise signals possess a linear moveout, in comparison to the hyperbolic moveout of the ice/bed reflection and any englacial reflections. However, the surface waves are spatially aliased for all other shots, inducing a lot of noise when applying the fk-filter.

Fig. 3. Processed and stacked data from the impulse source SISSY, for line 1 (north-south) and line 2 (west-east), with distances from the center point, C, close to KCI ice core. The bed reflector is visible at ∼0.04s (corresponding to ∼60m depth), marked with arrows. Englacial reflections cannot be unambiguously separated from artifacts, noise and non-reflective events.

Several filters were tested to reduce the strong ringing due to the spurious frequencies. A notch filter induced a lot of additional noise. This could be avoided by a narrow bandpass filter. However, cutting out the frequencies around 200 Hz destroyed the quality of the bed reflection. Killing the near-offset traces, which were affected by the spurious geophone ringing, provided the most efficient means of improving the SNR of the ice/bed reflection. Finally, a wide bandpass filter was used, with a lower limit of 10Hz to suppress noise below the natural frequency, and an upper limit of 450 Hz to reduce high-frequency noise.

The data were stacked with a median stack. This enhanced the quality of the result slightly, as the median instead of the arithmetic mean of the amplitudes is used for the stacking. A considerable amount of low-frequency noise was suppressed by filtering out signals <120 Hz. The bed reflector was thus visible for the common midpoints with the highest fold (Fig. 3). An fk-filter was reapplied to further improve the SNR of the ice-bottom reflection. Nevertheless, no englacial reflections are unambiguously discernible in the data.

Vibroseismic survey

Due to the great success of a vibrator survey in Antarctica (Eisen and others, 2010), a survey with a small lightweight micro-vibrator was carried out at Colle Gnifetti in August 2010. The goal of this repeat seismic experiment was to improve the results of the SISSY survey by using the advantages of a vibrator source: better source control and lower instantaneous forces. ElViS is an electrodynamic micro-vibrator, mounted on a wheelbarrow-like frame, and weighs, including batteries, ∼130 kg (Druivenga and others, 2011). While much heavier than SISSY, it was still feasible to manually move ElViS on the slightly tilted snow surface on Colle Gnifetti. The ElViS vibrator can excite P-waves as well as S-waves by vertical and horizontal motion, respectively, achieving vibration frequencies up to 360 Hz.

The seismic profiles surveyed with the SISSY source in 2008 were repeated with the ElViS source (Fig. 1). For the P-wave survey, a geophone spacing of 3 m with 23 vertical geophones was used. Shots were placed between the geophones, as in the 2008 SISSY survey. However, the source spacing of 3 m was continued outside the geophone line. The extension of shots outside the geophone line was between 75 m in the northerly direction and 30 m in the easterly direction, where space is limited due to a cliff. After some tests, the sweeps were set to 10 s duration, with a linear upsweep of 30-240 Hz for line 1 and 20-160 Hz for line 2. The survey geometry resulted in a 23-fold survey in the middle, decreasing to the ends of the survey line. The record length was set to 11 s, with a sampling interval of 1 ms. Thus the listening time of 1 s was far longer than the TWT of the bed reflection.

For the horizontal shear-wave (SH-wave) configuration, 47 horizontal geophones (10 Hz natural frequency), with geophone orientation across line, were placed at a 1.5 m spacing. The shot locations were again between the geophones and carried out with 1.5 m spacing outside the geophone line for ∼50m in each direction. A 10s linear upsweep of 60-360 Hz was chosen for line 1 and a 30-240 Hz sweep for line 2. Two sweeps of opposite polarity were excited at each location. Due to the higher number of geophones the survey was 47-fold in the middle, with the fold decreasing towards the ends of the line. The record length and sampling interval were set to the same values as for the P-wave survey. Once the survey was set up and working, the measurement itself was quite quick, requiring ∼2 min for two sweeps per shot location.

The processing of data from a vibrator source is somewhat different to that from explosives. As the signal is now a sweep, the first important processing step is the correlation of the recorded data with the source sweep. The signals are compressed to Klauder wavelets, resulting in subsurface signals, similar to those observed by impulse sources. As two shots were acquired at each location these compressed signals were stacked. In the case of the SH-sweeps, the polarity of the second signal was changed. This results in constructive interference of S-wave energy, and destructive interference of P-wave energy of both sweeps (see also King and Jarvis, 2007). The ElViS data also contain dominant diving and surface wave arrivals, though not as strong as in the SISSY data. The data were bandpass-filtered in order to see the bed and englacial reflections. Carefully chosen fk-filters were applied before stacking (Polom and others, in press). Velocities were picked independently for the different wave types and lines for the normal moveout correction and optimal stacking. The derived stacking velocities (i.e. root-mean-square (rms) velocities) were used in a smoothed form for the conversion of TWT to depth for the stacked sections (referred to as depth-conversion velocity below).

Once the stacked data were converted from the travel-time to the depth domain (Fig. 4, results of line 1), a difference in depth was observed between the SH-wave stack and the P-wave stack for the ice/bed reflection. The depth of the SH-wave bed reflector fitted rather well to the 62 m ice-core length, for line 1 and also line 2, while the P-wave bed reflection of line 1 was ∼6m too shallow (8m for line 2). Polom and others (in press) investigate several reasons for the difference in P-wave velocity and discuss different possibilities (e.g. trigger errors, modification of the wavelet away from a zero phase wavelet and errors in the velocity field due to mistakes in the moveout hyperbolas). Trigger errors could be ruled out, as the same equipment was used for the P- and SH-wave surveys. The other possible reasons for incorrect velocity analyses could not be excluded.

Fig. 4. Stacked data from the survey using the ElViS micro-vibrator as a source on line 1 (north-south) in P-wave and SH-wave modes. Distances are shown from the center point, C, close to the KCI ice core. The bed reflector can clearly be seen at ∼60m depth, with coherent englacial reflections visible above. The most prominent reflections are marked with arrows.

Figure 5b compares the depth-conversion velocities of the P- and SH-wave data of line 1 with P- and S-wave rms velocities, respectively, derived from ice-core density data (referred to as ice-core velocity). The S-wave rms velocities were calculated with a -ratio from the P-wave rms velocities. It is clearly visible that the P-wave depth-conversion velocities underestimate the calculated ice-core velocities. The calculated P-wave ice-core velocity is ∼300 m s^–1 faster at the bed reflection than the depth-conversion velocity. The comparison of picked and calculated S-wave velocities is less conclusive, as the S-wave ice-core velocities were calculated from the P-wave ice-core velocities with a fixed Poisson ratio of 0.25 (which is equal to , where v _P and vs are the P- and S-wave velocities). This is not necessarily a true representation of the Poisson ratio in firn. King and Jarvis (2007), for example, derived values for the Poisson ratio in polar firn that are between 0.193 and 0.340, a comparatively large range. However, more reflections, especially in the firn column, could be picked to determine the S-wave velocity profile in more detail.

Fig. 5. (a) Results from the analysis of diving waves from the SISSY survey, used to derive the density distribution in firn with the help of the Herglotz–Wiechert inversion. The measured density data from ice core KCI are shown in black, its smoothed density in green and one standard deviation within the 5 m sliding window in gray (visible down to 30m depth). Calculated density/depth functions from the travel-time/offset data (e.g. feature A in Fig. 2) for lines 1 and 2 are shown, using both the Kohnen and the Robin density/depth formulas. (b) Calculated P-wave and S-wave rms velocity profiles. The P-wave rms velocity profile (magenta) is calculated from the smoothed P-wave velocities (green, 5 m sliding-window), derived from the KCI density-data conversion (black), using Eqn (6) (Kohnen, 1972). Corresponding S-wave rms velocities were calculated for a _v p=v _s-ratio of (dark blue), from the smoothed P-wave rms velocities. For comparison, depth-conversion (DC) velocities of the ElViS P-wave (red) and S-wave (light blue) data of Figure 4 are shown, as well as SISSY rms velocities (purple) of the interval velocities derived from the SISSY diving waves of line 1 (red curve in (a)). The points on the ElViS DC curves indicate the depth at which the rms velocity was picked, which is valid for the interval above this point. To convert the individual velocity profiles from the time domain to the depth domain the corresponding rms velocities were used.

Our comparison between depth-conversion velocities determined from the seismic data and rms velocities calculated from ice-core properties only builds on ice-core densities. Further possible effects not considered by our ice-core velocities could arise from lateral density inhomogeneities or anisotropies of the COF in firn. An uncertainty of 10–15% in deriving depth-conversion velocities from stacking velocities is not exceptional (Hatton and others, 1986; Etris and others, 2001). Especially as we were only able to pick two reflections from the ElViS P-wave data, the depth conversion cannot be as precise as in the case of a smooth profile, such as that calculated from the KCI data. Smaller errors can only be achieved if the depth-conversion velocities can be cross-checked and improved by, for example, vertical seismic profiling data or other borehole results.

To obtain a stand-alone result from our seismic data analysis and allow for an independent evaluation, the depth-conversion velocities were not adjusted to the ice-core velocities. Instead, a direct time/depth adaption (Etris and others, 2001) was applied to shift the seismic P-wave data in the depth domain downwards to align with the bed reflection of the SH-wave depth image. This has a comparable effect to adjusting the P-wave stacking velocities to the ice-core velocities. The overall result of processing the vibroseis data is quite satisfying, as the bed reflector as well as englacial events can be clearly seen (Fig. 4, results of line 1). Details and possible physical origins of englacial reflections are discussed in a later section.

Density Profiles from Diving Waves

Seismic data analysis for velocities and densities

The SISSY dataset, while too noisy to show englacial reflections, provides strong diving waves. We utilize these to determine the seismic velocity and density structure of the firn (Fig. 2). For the density distribution at Colle Gnifetti, the first-arrival diving waves reach ∼30m depth for the largest offset of 120 m. Thus they sample the whole firn pack down to the firn/ice transition zone. These diving waves are used to derive a seismic velocity/depth profile of the firn pack using the Herglotz-Wiechert inversion (Aki and Richards, 2002). We then convert the seismic velocity profile to densities by employing the empirical formulas of Robin (1958) and Kohnen (1972), for comparison with density observations from the KCI core, with which we test the overall performance of the method.

An increasing velocity with depth is mandatory for the application of the Herglotz–Wiechert inversion, as low-velocity zones would lead to incorrect results. Due to the densification process of snow to glacier ice in the dry snow zone, the Herglotz–Wiechert inversion is a suitable method to calculate velocities of the firn pack within the accumulation area of a glacier. Although some melting occurs at Colle Gnifetti, forming ice lenses, these are thin compared to the seismic wavelength, such that the requirement of a continuously increasing velocity is assumed to be fulfilled for the diving waves with respect to the smoothed density distribution.

For the Herglotz–Wiechert inversion, use is made of the fact that the ray parameter, p, equals the reciprocal velocity, v, at the deepest point, Z, of the ray,

(1)

The largest horizontal offset, X(p), can thus be calculated by integrating the ray path over depth, z, up to the deepest depth, Z, corresponding to the ray with the largest offset, X:

(2)

The ultimate goal is to derive z(v) from X(p), which is accomplished by changing the integration variables, thus solving Abel’s problem (Aki and Richards, 2002). This results in

(3)

with offset x. The velocity vD = (∂x/∂t)_D is the gradient of the travel time, t, at the greatest source-to-receiver offset, D, whereas v _A is the apparent velocity given by v _A(x) = x/t.

After picking the travel times, t, of the diving wave for the different offsets, x, a curve needs to be fitted to the data to calculate the velocity, vD. We follow Kirchner and Bentley (1990), who use the Herglotz-Wiechert inversion to derive densities on the Ross Ice Shelf, Antarctica. We fit an exponential function to the picked travel-time/offset pairs of the form

(4)

The parameters a, b, c, d and e are constants. The velocity, vD, is calculated by the derivation of t with respect to x. Thus, velocity/depth data can be obtained from the travel times of the diving waves.

To compare the results of the Herglotz-Wiechert inversion to the densities measured at the KCI ice core, we have to convert the obtained velocities to densities. Two well-known equations exist for this conversion. Robin (1958) derived a linear relationship from laboratory tests, as well as measurements in Antarctica and at the Jungfraujoch, Switzerland (units for density, p, and the P-wave velocity, v _P, are kgm^r-3 and ms^–1, respectively:

(5)

Later, an empirical formula was derived by Kohnen (1972) to calculate densities from velocities,

(6)

with the P-wave velocity of ice, v _i, in ms^–

Comparison with ice-core density data

The density/depth distributions calculated from diving waves, hereinafter called seismic densities, of the near-offset shots are plotted together with the density data from the KCI ice core in Figure 5a. The seismic density profiles for lines 1 and 2 were calculated from the near-offset SISSY data (solid curves in Fig. 5a) and extrapolated through the ice column (dashed curves in Fig. 5a), following the approaches of Robin (1958) and Kohnen (1972) (Eqns (5) and (6), respectively). At first, the densities were calculated using the travel times of all shots per line. This resulted in a rather big difference between the seismic densities in respect to measured KCI densities, as the diving waves collect lateral inhomogeneities of the firn along the travel path over the whole survey area. To reduce this effect, only the travel times for diving waves of the near-offset shots (offset of up to 80 m for shots within 45 m of the center point) were taken into account afterwards. This approach is reasonable, as Alean and others (1983) have highlighted considerable lateral variations in accumulation at Colle Gnifetti; neglecting the far-offset seismic data, the subset samples a region of the firn pack that is more representative of the conditions at KCI. Nevertheless, there are differences in resolution. The ice-core data are point measurements in centimeter resolution in the vertical direction, horizontally averaged over a 10 cm diameter. The seismic result, instead, integrates information over an order-of-magnitude larger area and depth range, so it is also subject to lateral inhomogeneities.

The diving waves for the largest offset used (80 m) reach ∼20 m depth. The velocities and resulting seismic densities were calculated down to this corresponding depth (solid curves, Fig. 5a). Since the velocity is calculated by the fit of the exponential function to the picked data, the exponential function can be calculated to greater depth. Hence, the seismic densities for lines 1 and 2 were calculated down to the bed (dashed curves, Fig. 5a). This is possible, as the curvature of the velocity/depth function is predetermined by the fitted exponential function (Eqn (4)).

The KCI density observations, plotted together with the seismic densities calculated for lines 1 and 2, highlight the broader similarities between the two analytical approaches, particularly at depth (Fig. 5a). As we cannot resolve the small-scale variations of the density with the analysis of the diving waves, we compare the results to the smoothed set calculated with a 5 m sliding window and its standard deviation. The seismic densities from line 1 show very good agreement with the smoothed densities up to 20 m depth, especially between 10 and 20 m. They fit within one standard deviation of the smoothed KCI densities down to 25 m depth. Below, the seismic densities are slightly too low, ∼25 kgm^-3, but still follow the trend of the KCI densities. The seismic densities from line 2 are a bit more off, but still lie within one standard deviation of the average KCI densities down to 15 m. Below this depth they differ by ∼15% of the ice-core density. For both lines the best fit with the smoothed KCI densities is achieved with the results from the Kohnen formula.

There are various possible reasons for the differences between the results of lines 1 and 2. There are slight differences in bed and surface elevation, and thus in the cross-sectional geometry of both lines. Lateral inhomogeneities are likely due to the inhomogeneous accumulation pattern at Colle Gnifetti (Alean and others, 1983), as well as slight inhomogeneities in the temperature distribution, both leading to differences in the densification process and thus density/depth profiles (fig. 2 of Konrad and others, 2013). Another contribution could arise if the bulk COF constitutes a vertically non-symmetric distribution, thus also causing lateral anisotropies in seismic velocities. When we additionally include the far-offset shots in the analysis of the seismic densities, the variations from the KCI densities become rather large, already >15% for 25 m depth for line 1 and without the possibility of extrapolating the density distribution down to the glacier bed for line 2. This makes it difficult to investigate the reasons for the deviation between the line 1 and line 2 results. Altogether, the approach works reasonably well for the firn structure found on Colle Gnifetti, as the obtained densities resemble the smoothed density profile of the firn pack within a relative accuracy of 2-15% for the whole depth, with the best results obtained by the line 1 analysis down to a depth that was really sampled by the diving waves.

In addition to the comparison of ice-core and seismic densities, rms velocities were calculated from the seismic velocities derived from the Herglotz-Wiechert inversion of line 1 (Fig. 5b) for comparison with rms ice-core velocities and depth-conversion velocities from ElViS. For the calculation of SISSY rms velocities, only the interval velocities <2.5m were used to exclude the low-velocity part above this depth, caused by the Herglotz-Wiechert inversion. These rms velocities fit much better to the ice-core rms velocities than the ElViS depth-conversion velocities. The resolution of the SISSY rms velocities is, of course, much higher than that of the ElViS depth-conversion velocities. In the case of the ElViS data a reflection is required to pick a velocity, whereas the rms velocities from the diving waves give a continuously increasing velocity profile. Thus, if the excited diving waves are strong enough for a velocity analysis they can not only reveal the average density distribution of the firn, but can also help to improve the conversion from travel time to depth of the stacked data.

Englacial Seismic Layering

Observations

Whereas the SISSY data did not allow us to detect internal reflections from the firn and ice column, with only a weak ice/bed reflection identified (Fig. 3), it was possible to detect these events from the surveys with the ElViS source (Fig. 4). For the different lines and wave types, different frequency ranges for the sweeps were chosen, resulting in differences in resolution between the surveys. The maximum resolution (quarter-wavelength) can be calculated from the center frequency of the sweeps and the P- and S-wave velocities in ice, of 3900 and 2100ms^-1, respectively. Thus, for the P-wave surveys, we obtain a theoretical resolution of 7 m (maximum frequency 240 Hz) for line 1 and 11 m (center frequency 160 Hz) for line 2; and for the SH-wave surveys 2.5 m (maximum frequency 360 Hz) for line 1 and 4 m (maximum frequency 240 Hz) for line 2. Hence, the SH-wave has a resolution more than twice as good as that of the P-wave with the same maximum frequency content (240 Hz), which is due to the slower velocity of the S-wave in comparison to the P-wave. The resolution is, of course, better within the firn column, where velocities are slower and thus the smallest wavelength starts with ∼2 m for the SH-wave and 6 m for the P-wave. The resolution then decreases over depth with increasing seismic velocities due to increasing densities. Since the ElViS sweeps for line 1 contained higher frequencies, resulting in higher-resolution data for both the P-wave and SH-wave surveys for line 1 than for line 2, only the results of line 1 are discussed for investigating the origin of englacial reflections.

P-wave data from line 1 (Fig. 4) show a first strong reflection at 10 m depth, followed by two weaker reflections. A strong englacial reflector at 30 m depth can be observed, followed by a quiet zone. The strongest reflection in the section is the bed reflector at ∼60m depth, with another strong englacial reflection just above the bed at ∼50m depth.

In the SH-wave section of line 1 (Fig. 4) some obliquely incident signals are seen towards the north side of the profile, where the ice drops towards the Monte Rosa east face. A couple of reflections are observed between the surface and ∼10 m depth and a strong reflection is visible at ∼20m depth. Further down, the reflections are less pronounced, with some laterally coherent signals around 30, 40 and 50 m depth, followed by the strong bed reflection at ∼60 m depth.

Comparison with ice-core and GPR data

An improved understanding of the physical ice properties that produce the englacial seismic reflections observed from the ElViS P-wave and SH-wave surveys can be gained by comparing the seismic data with GPR and ice-core data. In Figure 6, a GPR profile and a single GPR trace are plotted together with the ice-core density data, as well as part of the SH-wave and P-wave sections of line 1. The differences in resolution for these measurements range from sub-centimeter scale for the ice-core densities, to ∼16cm resolution for the GPR data in ice, >4m for the P-wave data and >1.5 m for the SH-wave data. Different events are marked A-F. A comparison of the GPR and the ice-core data was performed by Eisen and others (2003) and Bohleber (2011) to clarify the origin of reflection horizons in the GPR data.

Fig. 6. Comparison of (a) GPR section, (b) single GPR trace, (c) density from KCI ice core with (d) SH-wave and (e) P-wave sections of the seismic survey (north–south) with the ElViS source. The GPR data are plotted with a spherical divergence correction over TWT, so the depth is slightly shifted compared to the other plots. Marked events A–F are discussed in the text. The location of the KCI ice core in the seismic sections is at the center point, C. (GPR records and KCI density from Bohleber, 2011.)

Prior to comparison, the P- and SH-wave sections have to be shifted in depth to obtain a consistent lower boundary of the seismic sections with the ice-core data and GPR sections. There is a time period of 5 years between the drilling of the KCI ice core (2005) and the ElViS survey (2010). During these 5 years, ∼1.75m of snow was accumulated at the surface, measured at the borehole casing. This causes an analogous shift in the depth of some physical properties, such as ice layers and impurities, i.e. they are advected downwards. The GPR data were also corrected for the additional accumulation between 2005 (ice-core drilling) and 2008 (GPR measurement) (Bohleber, 2011). To achieve a consistent lower boundary, the seismic data are shifted down by 3 m, so the bed reflections of the seismic sections fit with the 62 m (±1 m accuracy) of the ice-core and the GPR data. This shift is feasible, as the TWT-to-depth conversion is within an accuracy of 10-15%. Thus we are able to compare the seismic signals with the GPR and ice-core data (Fig. 6).

The depth of event A, at ∼10m, denotes the first clear reflection of the seismic SH-wave data. Strong peaks are visible in this region in the density data. The second event, B, shows a prominent peak in the density data and is the onset of a series of englacial-reflection horizons in the GPR data. At this depth the first strong reflection in the P-wave data is also visible. Reflections near event B exist in both the P-wave and SH-wave data, but in the case of the SH-wave data they are not clearly separated from the reflections of event A.

The region between events C and D contains the firn/ice transition zone with the pore close-off. The density data show no strong peaks in this depth range, but a decrease in density variability. At event C the deepest continuous englacial-reflection horizon in the GPR data can be observed. In the seismic data, a strong reflection is present in the SH-wave profile, that also appears to be subdued, though visible, within the P-wave profile. At event D, just below the firn/ice transition, a distinct reflection is visible in the P-wave profile but is unclear in the SH-wave profile. No corresponding signal can be found in the GPR data. In the depth range 40–50 m, around event E, a rather quiet zone can be observed in the seismic data, where some strong signals are visible in the single trace plot of the GPR data. However, a coherent continuous GPR reflection horizon is missing. Below, near event F, a strong reflector in the P-wave data can be seen, though no distinct reflection is observable in the SH-wave data. The GPR data show a rather blurred subglacial bed reflection. However, the bed is clearly visible in both the SH- and P-wave data (Fig. 6).

Interpretation of englacial reflections

The reflections in the P-wave data, as well as in the SH-wave data, in the region of A and B, seem to be caused by peaks in the density distribution from melt layers and ice lenses within the firn pack. For the region between events C and D we do not find strong changes in ice-core densities over short depth scales, in contrast to the reflections due to density inhomogeneities near events A and B. These observed reflections are due to changes in seismic velocities within the ice, which suggest a change in COF or some change in seismic velocities due to pore close-off in the firn/ice transition zone. Polom and others (in press) were able to derive a change in the Poisson ratio from P- and S-wave velocities for the depth between events C and D, where we observe the strong englacial seismic reflection in the P-wave data. Reflections from the firn/ice transition zone have also been observed in ElViS data from Antarctica (Eisen and others, 2012). The reason for low GPR reflectivity below the firn/ice transition zone remains unclear. We can currently only speculate that signal reduction is because clutter is involved (Konrad and others, 2013).

For the split reflector, ∼2m above the bed in the SH-survey, we suspect that the upper signal is due to some dirt intrusion, which was found when drilling was stopped at 62 m. The second reflection then belongs to the actual bed. The reflection 5 m above the bed in the P-wave data (event F) is more difficult to interpret, especially as there is no counterpart in the SH-section. We cannot connect this reflector with any of the GPR or ice-core signals. The density minimum at ∼51 m is probably an artifact, due to an unnoticed crack in the ice core. The small positive inhomogeneity in the density profile at ∼52 m seems to be real and would cause a reflection coefficient of ∼0.009. However, this density jump should also cause a reflection in the SH-wave data as strong as in the P-wave data. As there is no reflection in the SH-data at that depth, it is not likely that the P-wave reflection is due to the density anomaly. Another possibility for englacial reflections at this depth is again changes in seismic velocity due to changes in the orientation of the anisotropic ice crystals, hence changes in COF over depth. Measurements with a borehole radar carried out during the 2010 field season tentatively indicate that anisotropies might exist in the lowest part of the ice column (Bohleber, 2011). However, as continuous high-resolution in situ measurements of COF do not exist, at present we can only speculate that changes in COF might be a reason for the 50 m reflection within the P-wave section.

Conclusions

Two seismic surface sources, SISSY (impulse) and ElViS (vibrator), were tested to determine their suitability for near-surface investigations on firn and ice. The ElViS vibrator source exhibited clear advantages for detecting englacial layers, as the dominant surface waves and spurious geophone noise were less prominent than with the SISSY impulsive source. We were thus able to clearly observe englacial reflections in ElViS data, demonstrating that a vibrator source is clearly the preferable seismic source in a glacial setting like Colle Gnifetti, where a firn-covered ice mass with a cold recrystallization/infiltration zone is present. However, when the aim of a survey is only to investigate the density distribution of the firn pack, SISSY is an excellent seismic source generating strong diving waves. We demonstrate that diving waves are well suited to deriving the density profile of the firn pack on an Alpine glacier, and thus confirm earlier studies from polar regions. While it is not possible to observe small-scale variations in density within the travel times of the seismic waves, the average density profile we derived fits well with the measured densities from the KCI ice core. Moreover, the analysis of diving waves can help to improve the velocity profile for the time-to-depth conversion of seismic datasets.

We compared englacial reflections from the P-wave and SH-wave data to GPR and ice-core data. In the firn column, we attributed the englacial reflections to density changes caused by ice lenses and refrozen melt layers, which also cause internal layering in the much higher-resolution GPR data. Around the firn/ice transition, where GPR data lack any signals, we find strong reflections in the P- and SH-wave data. Since there are no observable density variations, these deeper seismic reflections are likely due to a change in COF or seismic velocity changes associated with pore close-off in the firn/ice transition zone. Below the firn/ice transition down to bedrock, the GPR data show no coherent signals, while the seismic data contain englacial reflections of considerable strength. With a constant density and no discernible SH-wave reflection, we suspect the P-wave reflectivity below the firn/ice transition zone arises from changes in COF.

We show that radar and seismic data supply comparable but also very different information about the firn and ice column. It is remarkable that in regions with low GPR reflectivity, prominent reflections can be observed in the seismic data, and the converse. Here the different sensitivities of seismic and radar reflections to different physical properties complement each other very well, as seismic waves are sensitive to elastic properties and radar waves are sensitive to dielectric properties of the ice. Hence, joint seismic and radar surveys can improve our understanding of ice properties beyond the present knowledge.