Travel-time modeling of VSP

In a hypothetical VSP experiment the first prominent arrival is the direct wave, which propagates directly from the source to the in-hole receiver. If the source is located at some offset, x, from the borehole (Fig. 3), the path of the propagating direct wave is oblique; otherwise, when the source is at the position of the borehole and the borehole is straight, the path is vertical and in this case the VSP may be called zero-offset VSP (Reference Sheriff and GeldartSheriff and Geldart, 1999). In layered media the propagation of seismic waves obeys Snell’s law: ray paths refract at the interface between two different layers, according to the incidence angle and the magnitude of the speed contrast. Such behavior causes rays to bend at interfaces.

In common geological situations seismic speeds in the subsurface vary over a very large range (<300 to >2000m s–¹; Reference Miller and XiaMiller and Xia, 1998; Reference Sheriff and GeldartSheriff and Geldart, 1999). Since speed contrasts are thus very high, and ray curvature is not negligible, a bent-ray approach is required to model VSP travel times (Reference Moret, Clement, Knoll and BarrashMoret and others, 2004). A bent-ray approach should also be applied in polar ice sheets or in the accumulation area of mountain glaciers, where the presence of a strong density gradient leads to P-wave speeds <1500ms–¹ in snow, rising to ∼3000ms–¹ in firn and reaching ∼3800 ms–¹ in ice (e.g. Reference King and JarvisKing and Jarvis, 2007). Such a significant speed gradient is not present in the ablation area of glaciers, where snow and firn are absent. Consequently for the ablation area of Storglaciären we can model ray paths using the negligible-ray-bending assumption (Reference Schuster, Johnson and TrentmanSchuster and others, 1988) and ray paths can be modeled as straight.

The travel-time model is structured using the physical properties illustrated in Figure 2. An 85 m VSP is modeled with 1 m depth-step sampling of the direct wave (the first-arriving compressional wave that, following a near-straight ray path, propagates between the source and the receiver). In our model the shot position is located at distance x from the borehole top. Receivers are located at depths z_i, where subscript i indicates the depth of the layer. By considering a VSP experiment of a simple two-layer case (Fig. 3) the travel times, t, of the direct wave recorded at positions z ₁ and z ₂ are given by

and

respectively. The path lengths, AE = EC, are easily found by trigonometrical relationships on the triangle ABC (Fig. 3).

The propagation angle, a, changes according to receiver depth. Thus for each z

Similarly, AE, the travel path in each layer, will diminish as a increases:

Thus the travel time of the received P-wave at depth z, t_z , is simply given by the sum of the travel times in each layer of different speed, v(z):

where n is the total number of layers considered in the model.

Field experiment

We collected four orthogonally orientated walkaway VSPs in an 80m deep, water-filled borehole located in the upper ablation area of Storglaciären in summer 2008 (Fig. 1). The borehole was drilled using the Stockholm University hot-water drill (Reference Jansson Pand NäslundJansson and Näslund, 2009). A walkaway is simply obtained by moving the shot position progressively away from the borehole. Seismic energy was generated with a sledgehammer and received by a 12-channel, 1 m spacing, 10 Hz hydrophone string, which was lowered down the borehole to different depths. The borehole was entirely water-filled. The array was lowered down-hole so the travel times of the direct wave could be estimated at a 1 m depth spacing. The shot position was located using a tape measure. The area investigated by the VSP experiment was generally flat or had very little topography. A flat area was primarily chosen because it facilitated drilling operations. Some topography (up to ∼2 m height difference in 60 m distance) was, however, present in the along-flow direction. In other words, the shot located 30 m west of the borehole was ∼2 m higher than the shot located 30 m east of the borehole (Fig. 1). The effect of this topographic gradient in the along-flow direction is discussed below.

Sample data for a VSP with source located 3 m away from the borehole are shown in Figure 4. The collected seismic profiles are of good quality. They show clear P-wave direct arrival first breaks, and are rich in other arrivals (i.e. converted-wave direct arrivals; englacial and subglacial P-wave reflections; tube waves, a ubiquitous noise source in VSPs) that have varying dominant frequencies as revealed by bandpass filtering (Fig. 4). Here we use only the direct arrivals first breaks to compare seismic travel times with our a priori model. The other seismic signals observed in the VSP records and the full velocity inversion will be the subject of future investigations.

Fig. 4. Sample VSP data collected at Storglaciären. The source (hammer and plate) was located 30 m from the borehole. (a) Filtered data (cut off at 150 and 800 Hz, cosine ramps between 150–300 and 400–800 Hz). (b) Sketch of the main events: P and converted P arrivals (dashed curves), tube waves (dotted curves) and internal reflections (solid curves).

Modeled P-wave speed structure in the ablation area

The modeled P-wave speed through the 80 m ice column imaged (Eqns (1) and (2)) varies by ∼100ms–¹ (Fig. 5, asterisks). The lower speed (∼3700ms^–1) is found in the upper part of the ice column, due to the presence of air within the cold ice. The transition from cold ice to temperate ice is marked by a ∼10 m s–¹ decrease in speed (Fig. 5). The seismic CTS at Storglaciären is thus characterized by a 0.3% decrease in speed, ∼10 times lower than the decrease in speed observed at the radar CTS (Reference Gusmeroli, Clark, Murray, Booth, Kulessa and BarrettGusmeroli and others, 2010a). The air content and its decrease with depth dominates the overall modeled seismic speed structure (Figs 5 and 2b). The temperature-induced decrease of speed from the cold to the temperate side of the CTS is very small (5 m s–¹; Fig. 5, curve T). In the bulk part of the ice column, speed gradually increases with depth, reaching the highest value of 3760 m s–¹ at 80 m depth.

Fig. 5. P-wave speed structure of the ablation area of Storglaciären inferred from the physical properties illustrated in Figure 2. The final speed model used to predict VSP arrivals (asterisks) was calculated from Eqn (2) by considering temperature, T, air and water. The speed structure calculated by temperature only (Eqn (1)) is shown with crosses. Combinations of T, water and T, air are also indicated with circles and dots, respectively.

VSP travel times

The arrival times of the propagating P-wave at different receiver depths were obtained from four orthogonally orientated VSP lines (Fig. 1). In this way, potential azimuthal variation in P-wave speed (e.g. ice anisotropy) can be detected, since four different propagation directions were sampled (up-flow, down-flow and across-flow in two opposite directions). We varied the source-offset distance, x, in 5 m increments from 5 to 30 m, and collected six VSPs for each line. Here we only show the differences in modeled and measured travel time for a representative 30 m offset, since it sufficiently summarizes the main findings of our experiment.

Results from VSP arrival times for the four lines (Fig. 6 for a 30 m offset VSP) show that the modeled P-wave structure generally matches the field experiment. This is especially evident for the two transverse profiles (LINE1 and LINE3; Fig. 6) as the misfit between modeled and measured travel time for these two lines is within 0.2 ms (Fig. 6). This misfit is within the travel-time uncertainties for the P-wave direct arrival first picks, which is approximated as S _t = 0:25 ms (twice our sampling rate of 0.125 ms). _t can be approximated to twice the sampling rate, because the sampled first break might actually be arriving one or two samples earlier. The uncertainty in travel time causes speed uncertainties which vary with the travel path. Our resultant P-wave velocity uncertainty, _v, is a function of the uncertainties in our observed travel time, t, and the ray-path length, d, for a given source-to-receiver arrangement,

Fig. 6. Example of a misfit plot (difference between modeled and measured arrival time) for a 30m offset VSP acquired at Storglaciären.

and the error in speed varies inversely with distance travelled, d. v is very high (between ±400 and ±100ms–¹ ) in the upper 10 m of the ice column, but it reduces to ±50 m s–¹ in the bulk of the ice column. We can thus be confident that the measured travel times at depths >10m in the across-flow direction agree with the speed structure indicated by asterisks in Figure 5, within a margin of ±50 ms^-1.

In the along-flow direction, the observed direct arrival travel times show a deviation of up to 1 ms from our speed model, the up-glacier travel times being consistently slower than the model predictions (LINE4; Figs 6 and 7). The misfit between our model and the seismic observations is indeed higher for these two profiles and well above the 0.25 uncertainty in our travel-time picks, with a deviation of up to ±1.0ms below 40 m receiver depth (Fig. 7). These results suggest apparent anisotropy in the flow direction; P-wave speed appears to be lower when propagating along-flow, whereas it is higher when propagating in the opposite direction. Apparent anisotropy is also confirmed by looking at the mean misfit measured in each survey (Fig. 7).

Fig. 7. Mean misfit for the 24 VSP surveys acquired. The error bars represent the standard deviation of the misfit for each survey. The gray area is 0 ± _t where _t is our mean standard deviation (0.18 ms) which is about twice the sampling rate (_δt = 0:25 ms); a reasonable error for our travel-time estimates. LINE1 and LINE3 are stable on the value of 0 ± _t , whereas those waves sampled in LINE2 and LINE4 are slower and faster, respectively. Theout-of-trend data point in LINE1 is believed to be caused by poor data quality for that particular VSP. Each line has six data points, which correspond to the six walkaway positions covered (5, 10, 15, 20, 25 and 30 m).

Is the straight-ray approach safe and appropriate?

The ice of the ablation area is homogeneous, old metamorphic glacier ice, already densified. We are thus correct to assume an a priori absence of a significant density gradient that leads to an equally important velocity gradient. However, our choice of a straight-ray approach is simplistic and violates Snell’s law (Fig. 3). Figure 5 demonstrates that seismic velocities can vary by ∼-100ms–¹ in the upper 20-30 m of the ice column; a potentially important velocity gradient thus exists and whether or not this impacts our analysis needs to be determined.

We thus undertook further analysis; based on the model observed in Figure 5, we compared modeled travel times following the straight-ray assumption and utilizing a framework that bends rays for each layer and does not violate Snell’s law. The results show that there is no practical difference between the two (Fig. 8) and that the straight-ray approach described above is applicable for the ablation area of a glacier where the near-surface speed gradient is minimal.

Fig. 8. Comparison between the straight-ray and the bent-ray travel-time modeling of a VSP. Dots and solid lines indicate straight-ray and bent-ray model, respectively. There is no practical difference between the two for the small velocity gradients observed in our study.

Seismic speeds and physical properties of glacier ice

Seismic speeds in temperate ice are important because they can be used to infer the microcrystalline water content and fabric, which are fundamental parameters in modeling ice viscosity. Reference Navarro, YuYa and BenjumeaNavarro and others (2005) and Reference Endres, Murray, Booth and WestEndres and others (2009) inferred water contents as high as 2% by inverting low seismic speed for water content in two-phase (ice/water) models. Our model and its agreement with the VSP travel times suggest that the P-wave speed in the cold surface layer is relatively low, ∼3700ms–¹. This relatively low value is also supported by the speed analysis of the refraction surveys collected in the same area (Reference Gusmeroli, Murray, Jansson, Pettersson, Aschwanden and BoothGusmeroli and others, 2010b). Reference GusmeroliGusmeroli (2010) provides further agreement, measuring an average P-wave speed of 3617 ±57 m s–¹ in the uppermost ice. Here we suggest that this relatively low seismic speed is primarily caused by the air content of the ice. Such low seismic speeds have traditionally been interpreted as being typical of temperate ice with some volumetric percentages of microcrystalline water (Reference Benjumea, YuYa, Navarro and TeixidóBenjumea and others, 2003; Reference Endres, Murray, Booth and WestEndres and others, 2009). A simple (two-phase) interpretation of speeds as low as 3600 m s–¹ could lead to the wrong conclusion that the upper 20 m of the ice column are composed of temperate ice with a volumetric water content of ∼3%. This value is calculated using Eqn (2) for a water/ice mixture. Other two-phase models, such as the Riznichencko model (Reference Benjumea, YuYa, Navarro and TeixidóBenjumea and others, 2003), do not differ significantly.

Instead seismic speeds measured and modeled in our study area suggest that relatively low P-wave speeds (e.g. 3600-3700 ms–¹) are to be expected for cold ice in the ablation area. These low speeds are not due to the presence of water but, instead, to the presence of small percentages of air. As observed in polar ice masses, v _P increases because of increasing pressure which decreases bubble size and air content. Although dominant, the magnitude of the air-content-induced change is, however, very small (only 100ms–¹ between the lowest and the highest speed in Fig. 5). An even smaller, perhaps undetectable, change (only 10ms–¹) is caused by the low water contents measured at Storglaciären. We believe that the change in vP due to water content might be hard to detect given the uncertainty often observed in many glaciological seismic speed estimates (tens of m s–¹).

Equation (2) indirectly incorporates the densification with depth using air content. We thus consider the bulk density of the polycrystalline glacier ice to be primarily conditioned by the volumetric air content. Reference Kohnen and BentleyKohnen and Bentley (1973) pointed out that a universal density/speed function applicable to different glaciological settings might be difficult to obtain. Most of our knowledge of density/speed relationships comes from observations in polar firn (Reference Kohnen and BentleyKohnen and Bentley, 1973; Reference King and JarvisKing and Jarvis, 2007), where density changes up to 0.5 gcm^{- 3} lead to P-wave speed variations of ∼1500 m s–¹ (Reference King and JarvisKing and Jarvis, 2007). In the ablation area of a mountain glacier (e.g. this study) there is no firn and the density gradient is significantly lower.

The necessity of including air content when inferring water content from radio-wave speed estimates appears to have been established in recent radioglaciological literature (Reference Gusmeroli, Murray, Barrett, Clark and BoothGusmeroli and others, 2008, Reference Gusmeroli, Murray, Jansson, Pettersson, Aschwanden and Booth2010b; Reference Bradford, Nichols, Mikesell and HarperBradford and others, 2009). The reason is that the presence of air increases the bulk radio-wave speed and thus causes underestimation of water content. From a seismic perspective, this behavior is quite the opposite; the presence of air within the ice matrix decreases the bulk speed, meaning that water content can be overestimated if air content is ignored. The glaciological consequences of this are important; by ignoring the presence of air, bulk speeds in the water-free cold ice of the ablation area of Storglaciaren suggest water content up to 3%.

We stress that the P-wave structure modeled in this study is valid only for the ablation area of the glacier. In the accumulation area the presence of snow and firn might lead to speeds considerably lower than those reported in Figure 5. Additionally, we have very little knowledge about the water content of the ice in the accumulation area. Reference Schneider and JanssonSchneider and Jansson (2004) found values as high as 3% and these values are likely to be significant. Acquiring borehole geophysical data in the ablation area of glaciers is easier than working in the accumulation areas, where the occurrence of snow-covered crevasses and the presence of porous firn result in logistical challenges that must be overcome to improve our ability to use ground-based geophysics to detect ice properties.

Azimuthal variation in seismic results

Results from our VSP surveys (Figs 6 and 7) might suggest an up-/down-glacier oriented seismic anisotropy. P-waves propagating from west to east (down-flow) propagate more slowly than those waves propagating in the opposite direction (Fig. 6). This behavior is consistent with P-wave arrival times, which are higher than modeled in the down-glacier direction but are lower than modeled in the up-glacier direction (LINE2 and LINE4 in Figs 6 and 7). In this subsection we briefly evaluate whether or not intrinsic anisotropy is observable in our data. Seismic P-wave speed is strongly anisotropic in single-crystal ice (see Reference Gusmeroli, Pettit, Kennedy and RitzGusmeroli and others, 2012, for a recent review). When the wave propagates parallel to the c-axis it propagates ∼5% faster than when it propagates parallel to the basal plane (Reference BentleyBentley, 1972; Rothlisberger, 1972; Reference Blankenship and BentleyBlankenship and Bentley, 1987; Reference Anandakrishnan, Fitzpatrick, Alley, Gow and MeeseAnandakrishnan and others, 1994; Reference Gusmeroli, Pettit, Kennedy and RitzGusmeroli and others, 2012). In a polycrystalline aggregate of glacier ice this seismic anisotropy is possible when the majority of crystals are aligned along certain preferential directions. A response such as the one observed in our dataset (slower wave from up-glacier) can be explained by a fabric style which is not typically observed in glacier ice. Another factor that could potentially contribute to along-flow differences is variations in the percentage of water. A slower wave from up-glacier implies higher water content up-glacier from the center of the VSP experiment. Reference Pettersson, Jansson and BlatterPettersson and others (2004) found distinct patterns of water content at the CTS, with lower and higher water content on either side of the glacier center line. We thus may expect across-glacier rather than along-glacier variability. Additionally the horizontal and vertical variability documented by both Reference Pettersson, Jansson and BlatterPettersson and others (2004) and Reference Gusmeroli, Murray, Jansson, Pettersson, Aschwanden and BoothGusmeroli and others (2010b) still lies within (/>_w<1%. For these reasons we do not believe that the observed azimuthal variation can be explained with different cj) _w.

Inclinometry data from the borehole show very little deviation from the vertical. The maximum deviation is <1 m and the main deviation is towards the northeast. The upper 30 m of the borehole exhibits a deviation of <0.5 m; a maximum deviation of 1 m is observed at the base of the borehole (80 m depth), resulting in a 0.01 ms change in travel times from a truly vertical borehole. Deviation of ∼1 m has a negligible effect (∼0.01 ms) on the arrival times. We used differential GPS (dGPS) to measure differences in topography of ∼2 m for the greatest offset distance covered in our study (30 m). It thus follows that a receiver located at depth z is in reality located at depth z + 1 m for the wave propagating down-glacier and at depth z – 1 m for the wave propagating up-glacier. For z = 40 m we obtain travel-time differences of ∼0.5 ms between the two directions. Topography might thus explain most of the observed apparent travel-time anisotropy, because singular profiles (Fig. 6) and mean misfits are <1 ms (Fig. 7). It is therefore likely that the observed misfits are due to topography instead of real physical anisotropy.