Skip to main content Accessibility help


  • Access
  • Cited by 1



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Vertical seismic profiling of glaciers: appraising multi-phase mixing models
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Vertical seismic profiling of glaciers: appraising multi-phase mixing models
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Vertical seismic profiling of glaciers: appraising multi-phase mixing models
        Available formats
Export citation


We have investigated the speed of compressional waves in a polythermal glacier by, first, predicting them from a simple three-phase (ice, air, water) model derived from a published ground-penetrating radar study, and then comparing them with field data from four orthogonally orientated walkaway vertical seismic profiles (VSPs) acquired in an 80 m deep borehole drilled in the ablation area of Storglaciären, northern Sweden. The model predicts that the P-wave speed increases gradually with depth from 3700ms–1 at the surface to 3760ms–1 at 80m depth, and this change is almost wholly caused by a reduction in air content from 3% at the surface to <0.5% at depth. Changes in P-wave speed due to water content variations are small (<10 ms–1); the model’s seismic cold–temperate transition surface (CTS) is characterized by a 0.3% decrease downwards in P-wave speed (about ten times smaller than the radar CTS). This lack of sensitivity, and the small contrast at the CTS, makes seismically derived water content estimation very challenging. Nevertheless, for down-going direct-wave first arrivals for zero- and near-offset VSP shots, we find that the model-predicted travel times and field observations agree to within 0.2 ms, i.e. less than the observational uncertainties.


Polycrystalline glacier ice near the melting temperature (∼0°C) contains μm-scale water inclusions at the three-grain junctions (Nye and Frank, 1973; Mader, 1992). This water significantly increases the strain rate of temperate glacier ice (Duval, 1977), which has motivated numerous studies targeted at measuring microcrystalline water content in glaciers (e.g. Vallon and others, 1976; Macheret and others, 1993; Murray and others, 2000; Pettersson and others, 2004; Bradford and Harper, 2005; Navarro and others, 2005; Bradford and others, 2009; Endres and others, 2009; Gusmeroli and others, 2010b).

Polythermal glaciers are composed of both cold and temperate ice, and represent a perfect experimental setting to investigate changes in microcrystalline water content. This is because the boundary between cold and temperate ice, the cold-temperate transition surface (CTS; Hutter and others, 1988; Pettersson and others, 2003), is a prominent englacial feature. Across the CTS a jump in water content (from near-zero to single-digit percentages) and hardness (stiffer cold ice to softer temperate ice) is typically observed (Aschwanden and Blatter, 2005; Gusmeroli and others, 2010a). Such an increase in water content on the temperate side of the CTS lowers the propagation speed of both radar and seismic waves. The speed of these signals can therefore theoretically be used to infer microcrystalline water content.

Temperate ice also contains air, mainly in the form of bubbles (Raymond, 1976), and is more appropriately described as a three-phase medium containing solid polycrystals, liquid microcrystalline water and air in bubbles. Such a three-phase approach has recently been demonstrated to be particularly crucial in terms of water content estimates from radar speeds (Bradford and Harper, 2005; Bradford and others, 2009; Gusmeroli and others, 2010a). When the air content of temperate ice is neglected (i.e. when only a two-phase ice/water model is used) water content estimates are typically overestimated (Bradford and Harper, 2005; Gusmeroli and others, 2008). This sensitivity is, however, still unclear from a seismic perspective. Several studies (Navarro and others, 2005; Endres and others, 2009) have interpreted variation in seismic speeds as an indicator of water content using a two-phase (ice/water) approach; however, these studies have not addressed the effects of air content (responsible for bulk density; Kohnen and Bentley 1973) and ice temperature (Kohnen, 1974). Thus the effect of air content on seismic waves propagating in a three-phase mixture of ice, water and air remains to be investigated.

Here we investigate the vertical distribution of compressional-wave (P-wave) speeds using a three-phase framework, by forward modeling the effects that known air content, temperature and water content values have on seismic waves. We simulate the P-wave structure of a vertical section in the ablation area of the polythermal glacier Storglaciaren, Sweden (Fig. 1), using a high-resolution (1 m depth-step) density and water content profile measured from borehole radar (Gusmeroli and others, 2010a). We then compare the modeled speed structure with the results of a series of four orthogonally orientated walkaway vertical seismic profiles (VSPs) acquired in the same area in the summer of 2008.

Fig. 1. Map of Storglaciären showing the location of our VSP experiment. (Coordinates Swedish Grid RT90 1616415–7536885.) The inset shows the survey geometry, with the borehole in the middle of the cross. Four walkaway VSP lines were obtained by moving the shot position progressively away from the borehole. The tick-marks on the line indicate shot positions. The approximate position of the equilibrium-line altitude is indicated with a dashed curve.

Storgläciaren: Physical Properties of the Study Area

Storglaciaren (Fig. 1) is a small (3 km long), well-studied valley glacier located in the Kebnekaise massif, Sweden, just north of the Arctic Circle (e.g. Jansson, 1996; Holmlund and others, 2005). Storglaciären is entirely temperate in the accumulation area, whereas it has a polythermal structure in the ablation area. The polythermal structure at our study site is composed of a cold, permanently frozen, surface layer overlying a temperate core (Hooke and others, 1983; Holmlund and Eriksson, 1989; Pettersson and others, 2003). The water content in the temperate ice on Storglaciären varies spatially (Pettersson and others, 2004) but is generally small and homogeneous with depth, typically 0.6–1% (Gusmeroli and others, 2010a). Higher water contents (possibly up to 2%) are suspected for deeper, near-bed ice (Aschwanden and Blatter, 2005, 2009; Gusmeroli and others, 2010a); these higher values result from heat of deformation, which melts some ice and increases the water content.

Gusmeroli and others (2010b) studied the physical properties of a vertical section of the upper ablation area of the glacier (Fig. 1), and obtained vertical profiles of temperature, air content and water content for the upper 80 m of the ice column. A thorough description of the experiments conducted to obtain the englacial measurements is given by Gusmeroli and others (2010b); here we present a brief synopsis of the results.

The temperature profile, measured by thermistors installed within the ice, shows that the CTS is at 21 m depth at this location (Fig. 2a). Ice temperature at the surface is below the melting point and increases with depth. The air content (Fig. 2b) was measured using the near-surface (1m depth) radio-wave speed (Gusmeroli and others, 2010a). This shallow estimate gave an initial value which was subsequently modeled using a pressure/porosity model based on the theory of Bradford and others (2009). The modeled air content profile (Fig. 2b) shows a decrease of density with depth due to air-bubble compression. Gusmeroli and others (2010b) used this air profile to correct the radio-wave speeds. These corrected speeds can be inverted for water content (Fig. 2c). The result shows that water content values measured from cross-borehole radar at Storglaciären are relatively homogeneous with depth in the upper part of the temperate ice, with a tendency towards higher values at greater depths. The mean water content measured within the ice column is 0:6 ± 0 :3%. In this study we simulate the slightly higher water contents at depth with a three-layer water content structure consisting of an upper layer (21-60 m depth) with 0.5% water content and a lower layer (61-85 m depth) with 0.6% water content. The three profiles in Figure 2 represent the starting point of our seismic forward modeling. Although the water content model here is a simplified version of that of Gusmeroli and others (2010b), it sufficiently describes its most important characteristic: a homogeneous layer with only a slight increase in water content at greater depths. These physical properties (Fig. 2) will be used to obtain the P-wave speed structure within the glacier.

Fig. 2. Models of the physical properties of the study area. (a) Temperature distribution with depth. (b) Air content distribution with depth. (c) Water content distribution with depth. Models in this figure are derived from the field measurements described by Gusmeroli and others (2010a).

Modeling the P-Wave Speed Structure

Propagation speed of compressional P-waves within glaciers depends on temperature, air content and water content. Cold temperatures, low air content and low water contents result in high seismic speeds. Based on in situ observations Kohnen (1974) derived the temperature dependence of P-wave speed as


where T is ice temperature (°C). v P(T) is often considered a reference value for cold, water-free ice. If glacier ice was not a multi-phase material, then v P(T) would be the speed structure within the subsurface. However, since temperate glacier ice can also be composed of air (Fig. 2b) and water (Fig. 2c), a three-phase medium should be considered. Gusmeroli and others (2012) reviewed suggestions for temperature correction published since Kohnen’s (1974) work, and concluded that they are all relatively similar. For example, Helgerud and others (2009) suggested a gradient of –2.7. The difference between Helgerud and others (2009) and Kohnen (1974) is almost negligible for the small temperature range used in this study and therefore the Kohnen (1974) relationship is considered adequate.

In polar ice sheets, P-wave speed increases with depth, mirroring the increase of density, p (Kohnen and Bentley, 1973). This behavior reflects the reduction of the relative proportion of air in the two-phase firn. A simple, well-established method to estimate the bulk speed in a multiphase medium is Wyllie’s time-average equation (Wyllie and others, 1958). The bulk propagation speed in the medium will be proportional to the relative volumetric concentrations of the components. Thus, by defining vPa and v Pw as the P-wave speed in air (330 ms–1; Sheriff and Geldart, 1999) and water (1450 ms–1; Sheriff and Geldart, 1999), respectively, we can obtain v P(T, a, w), the bulk speed in the three-phase glacier ice,


where Φa and Φ w are volumetric air and water fractions, respectively. Both vPa and vPw may vary with depth in the borehole, due to temperature and pressure change (Mackenzie, 1981; Sheriff and Geldart, 1999); however, these variations are ∼2ms–1, which is insignificant for our study. We thus assume constant v Pa and vPw. Equation (2) includes the strong dependence of seismic speeds on density (Kohnen and Bentley, 1973).

Vertical Seismic Profiling (Vsp)

Borehole VSP is an established technique in exploration seismology which allows quantification of subsurface seismic speed (Sheriff and Geldart, 1999). In typical geological applications heterogeneity in the stratified subsurface results in high P-wave speed contrasts, which allow differentiation of material properties. Speed contrasts in glaciological applications are smaller than in rock applications, so an additional aim of this study is to examine the usefulness of VSP experiments in alpine glaciers. In VSP a series of receivers, located in a borehole, record seismic energy generated at the surface at some variable distance from the borehole (Fig. 3).

Fig. 3. Schematic representation of a VSP experiment using the negligible-ray-bending assumption. The shot position, at offset x = is indicated by a star located at A. Receivers, indicated by circles, are located in the borehole at depths z 1, z 2 and z 3. A, B, C and E are the trigonometric points used in the calculations.

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 (Sheriff 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–1; Miller and Xia, 1998; Sheriff 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 (Moret 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–1 in snow, rising to ∼3000ms–1 in firn and reaching ∼3800 ms–1 in ice (e.g. King 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 (Schuster 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 zi, 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 1 and z 2 are given by




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, tz , 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 (Jansson 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–1 (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–1 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 (Gusmeroli 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–1; 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–1 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–1 ) in the upper 10 m of the ice column, but it reduces to ±50 m s–1 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).


For down-going direct-wave first arrivals, for zero- and near-offset VSP shots, we find that the travel times predicted by the model agree with field observations to within 0.2 ms, i.e. less than the observational uncertainties (Fig. 7). This is particularly true for the VSP surveys acquired perpendicular to the glacier flow, where we found negligible differences between modeled and measured P-wave arrival times. Potential along-flow seismic anisotropy was also observed because waves propagating down-glacier arrived ∼0.8–1.0 ms after those propagating up-glacier. In this section we discuss the two main issues that arose in this study: (1) the relationship between seismic speeds and physical properties in the ablation area of a polythermal glacier and (2) the presence of apparent along-flow seismic anisotropy.

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–1 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. Navarro and others (2005) and Endres 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–1. This relatively low value is also supported by the speed analysis of the refraction surveys collected in the same area (Gusmeroli and others, 2010b). Gusmeroli (2010) provides further agreement, measuring an average P-wave speed of 3617 ±57 m s–1 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 (Benjumea and others, 2003; Endres and others, 2009). A simple (two-phase) interpretation of speeds as low as 3600 m s–1 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 (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–1) 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–1 between the lowest and the highest speed in Fig. 5). An even smaller, perhaps undetectable, change (only 10ms–1) 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–1).

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. Kohnen 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 (Kohnen and Bentley, 1973; King and Jarvis, 2007), where density changes up to 0.5 gcm- 3 lead to P-wave speed variations of ∼1500 m s–1 (King 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 (Gusmeroli and others, 2008, 2010b; Bradford 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. Schneider 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 Gusmeroli 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 (Bentley, 1972; Rothlisberger, 1972; Blankenship and Bentley, 1987; Anandakrishnan and others, 1994; Gusmeroli 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. Pettersson 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 Pettersson and others (2004) and Gusmeroli 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.

Conclusions and Outlook

We have studied the seismic P-wave structure of the upper ablation area of Storglaciaren, and found that the seismic velocity structure of the ice is largely dominated by the presence of air. P-wave speed typically increases with depth from 3700 m s–1 at the surface to 3760 m s–1 at 80 m depth; this change is almost entirely explained by a reduction in air content from 3% at the surface to <0.5% at depth. Changes in P-wave speed due to water content are very small (∼10ms–1). The seismic CTS is characterized by a 0.3% decrease in speed, about ten times lower than the decrease in radar-wave speed observed for the radar CTS (Gusmeroli and others, 2010a). These results highlight the importance of including air content when inferring physical properties of temperate glacier ice; if air is neglected, a two-phase (ice/water) interpretation of relatively low P-wave speed in cold, water-free ice could lead to the conclusion that the ice was temperate with some volumetric percentage of water.

Our modeled arrival times match the observed arrival times to within 0.25 ms in the transverse direction, corresponding to a seismic speed uncertainty of 50ms–1. This uncertainty in constraining the seismic speed is high for detecting water content. Future studies should focus on reducing uncertainties in seismic speed determination as well as applying a joint inversion of radar-seismic datasets, as in Endres and others (2009). We believe that future application of VSP to glaciological problems, perhaps combined with the use of explosives to increase the signal-to-noise ratio and increase sampling frequency, might enhance important advances in glacier geophysics. These advances include a detailed analysis of the shear waves collected here (as they are more sensitive to water than P-waves), and ongoing efforts to seismically infer englacial temperature (Gusmeroli and others, 2010a; Peters and others, 2012), model ice rheology and, ultimately, obtain reliable datasets that can be useful for modelers when predicting the future dynamics of mountain glaciers and ice sheets.


A.G. is supported by the Alaska Climate Science Center, funded by Cooperative Agreement No. G10AC00588 from the United States Geological Survey. The contents of this paper are solely the responsibility of the authors and do not represent the official views of USGS. Fieldwork was undertaken while A.G. was a postgraduate student at Swansea University, and was funded by the Jeremy Willson Charitable Trust, the Consorzio dei Comuni del Bacino Imbrifero Montano dell’Adda, the Percy Sladen Fund, the Mount Everest Foundation, the British Society for Geomorphology, the Quaternary Research Association, the Dudley Stamp Memorial Fund and the Earth and Space Foundation. A.G. thanks E.C. Pettit for support, discussions and comments on an earlier version of the manuscript. Partial financial support for P.J. was provided by Swedish Research Council grant 621-2007-3752B. B. Reinardy and D. Hjelm helped in the field. Logistics from the Tarfala Research Station were fundamental for this study. Queen’s University Belfast let us use its seismic system. R.S. Sturm provided assistance in revising the figures. Discussions with and help from A.D. Booth and B.E. Barrett are gratefully acknowledged. Careful reviews from H. Horgan and an anonymous reviewer improved the manuscript.


Anandakrishnan, S, Fitzpatrick, JJ, Alley, RB, Gow, AJ and Meese, DA (1994) Shear-wave detection of asymmetric c-axis fabrics in the GISP2 ice core, Greenland. J. Glaciol., 40(136), 491496
Aschwanden, A and Blatter, H (2005) Meltwater production due to strain heating in Storglaciären, Sweden. J. Geophys. Res., 110(F4), F04024 (doi: 10.1029/2005JF000328)
Aschwanden, A and Blatter, H (2009) Mathematical modeling and numerical simulation of polythermal glaciers. J. Geophys. Res., 114(F1), F01027 (doi: 10.1029/2008JF001028)
Benjumea, B, YuYa, Macheret, Navarro, FJ and Teixidó, T (2003) Estimation of water content in a temperate glacier from radar and seismic sounding data. Ann. Glaciol., 37, 317324 (doi: 10.3189/172756403781815924)
Bentley, CR (1972) Seismic-wave velocities in anisotropic ice: a comparison of measured and calculated values in and around the deep drill hole at Byrd Station, Antarctica. J. Geophys. Res., 77(23), 44064420 (doi: 10.1029/JB077i023p04406)
Blankenship, DD and Bentley, CR (1987) The crystalline fabric of polar ice sheets inferred from seismic anisotropy. IAHS Publ. 170 (Symposium at Vancouver 1987 – The Physical Basis of Ice Sheet Modelling), 1728
Bradford, JH and Harper, JT (2005) Wave field migration as a tool for estimating spatially continuous radar velocity and water content in glaciers. Geophys. Res. Lett., 32(8), L08502 (doi: 10.1029/2004GL021770)
Bradford, JH, Nichols, J, Mikesell, TD and Harper, JT (2009) Continuous profiles of electromagnetic wave velocity and water content in glaciers: an example from Bench Glacier, Alaska, USA. Ann. Glaciol., 50(51), 19 (doi: 10.3189/172756409789097540)
Duval, P (1977) The role of the water content on the creep rate of polycrystalline ice. IAHS Publ. 118 (Symposium at Grenoble 1975 – Isotopes and Impurities in Snow and Ice), 2933
Endres, AL, Murray, T, Booth, AD and West, LJ (2009) A new framework for estimating englacial water content and pore geometry using combined radar and seismic wave velocities. Geophys. Res. Lett., 36(4), L04501 (doi: 10.1029/2008GL036876)
Gusmeroli, A (2010) Polythermal glacier dynamics at Storglaciären, Arctic Sweden, inferred using in situ geophysical techniques. (PhD thesis, University of Swansea)
Gusmeroli, A, Murray, T, Barrett, B, Clark, R and Booth, A (2008) Correspondence. Estimates of water content in glacier ice using vertical radar profiles: a modified interpretation for the temperate glacier Falljökull, Iceland. J. Glaciol., 54(188), 939941 (doi: 10.3189/002214308787779942)
Gusmeroli, A, Clark, RA, Murray, T, Booth, AD, Kulessa, B and Barrett, BE (2010a) Seismic wave attenuation in the uppermost glacier ice of Storglaciären, Sweden. J. Glaciol., 56(196), 249256 (doi: 10.3189/002214310791968485)
Gusmeroli, A, Murray, T, Jansson, P, Pettersson, R, Aschwanden, A and Booth, AD (2010b) Vertical distribution of water within the polythermal Storglaciären, Sweden. J. Geophys. Res., 115(F4), F04002 (doi: 10.1029/2009JF001539)
Gusmeroli, A, Pettit, EC, Kennedy, JH and Ritz, C (2012) The crystal fabric of glacial ice from full-waveform borehole sonic logging. J. Geophys. Res., 117(F3), F03021, 1–13 (doi: 10.1029/2012JF002343)
Helgerud, MB, Waite, WF, Kirby, SH and Nur, A (2009) Elastic wave speeds and moduli in polycrystalline ice Ih, sI methane hydrate, and sII methane-ethane hydrate. J. Geophys. Res., 114(B2), B02212 (doi: 10.1029/2008JB006132)
Holmlund, P and Eriksson, M (1989) The cold surface layer on Storglaciären. Geogr. Ann. A, 71(3–4), 241244
Holmlund, P, Jansson, P and Pettersson, R (2005) A re-analysis of the 58 year mass-balance record of Storglaciären, Sweden. Ann. Glaciol., 42, 389394 (doi: 10.3189/172756405781812547)
Hooke, RLeB (2005) Principles of glacier mechanics, 2nd edn. Cambridge University Press, Cambridge
Hooke, RLeB, Gould, JE and Brzozowski, J (1983) Near-surface temperatures near and below the equilibrium line on polar and subpolar glaciers. Z. Gletscherkd. Glazialgeol., 19(1), 125
Hutter, K, Blatter, H and Funk, M (1988) A model computation of moisture content in polythermal glaciers. J. Geophys. Res., 93(B10), 1220512 214
Jansson, P (1996) Dynamics and hydrology of a small polythermal valley glacier. Geogr. Ann. A, 78(2–3), 171180
Jansson Pand Näslund, J-O (2009) Spatial and temporal variations in glacier hydrology on Storglaciären, Sweden. (SKB Tech. Rep. TR-09-13) Swedish Nuclear Fuel and Waste Management Company, Stockholm
King, EC and Jarvis, EP (2007) Use of shear waves to measure Poisson’s ratio in polar firn. J. Environ. Eng. Geophys., 12(1), 1521 (doi: 10.2113/JEEG12.1.15)
Kohnen, H (1974) The temperature dependence of seismic waves in ice. J. Glaciol., 13(67), 144147
Kohnen, H and Bentley, CR (1973) Seismic refraction and reflection measurements at ‘Byrd’ station, Antarctica. J. Glaciol., 12(64), 101111
YuYa, Macheret, Moskalevsky, MYu and Vasilenko, EV (1993) Velocity of radio waves in glaciers as an indicator of their hydrothermal state, structure and regime. J. Glaciol., 39(132), 373384
Mackenzie, KV (1981) Discussion of sea water sound-speed determinations. J. Acoust. Soc. Am., 70(3), 801806
Mader, HM (1992) Observations of the water-vein system in polycrystalline ice. J. Glaciol., 38(130), 333347
Miller, RD and Xia, J (1998) Large near-surface velocity gradients on shallow seismic reflection data. Geophysics, 63(4), 13481356 (doi: 10.1190/1.1444436)
Moret, GJM, Clement, WP, Knoll, MD and Barrash, W (2004) VSP traveltime inversion: near-surface issues. Geophysics, 69(2), 345351 (doi: 10.1190/1.1707053)
Murray, T, Stuart, GW, Fry, M, Gamble, NH and Crabtree, MD (2000) Englacial water distribution in a temperate glacier from surface and borehole radar velocity analysis. J. Glaciol., 46(154), 389398 (doi: 10.3189/172756500781833188)
Navarro, FJ, YuYa, Macheret and Benjumea, B (2005) Application of radar and seismic methods for the investigation of temperate glaciers. J. Appl. Geophys., 57(3), 193211 (doi: 10.1016/j.jappgeo.2004.11.002)
Nye, JF and Frank, FC (1973) Hydrology of the intergranular veins in a temperate glacier. IASH Publ. 95 (Symposium at Cambridge 1969 – Hydrology of Glaciers), 157161
Peters, LE, Anandakrishnan, S, Alley, RB and Voigt, DE (2012) Seismic attenuation in glacial ice: a proxy for englacial temperature. J. Geophys. Res., 117(F2), F02008 (doi: 10.1029/2011JF002201)
Pettersson, R, Jansson, P and Holmlund, P (2003) Cold surface layer thinning on Storglaciären, Sweden, observed by repeated ground penetrating radar surveys. J. Geophys. Res., 108(F1), 6004 (doi: 10.1029/2003JF000024)
Pettersson, R, Jansson, P and Blatter, H (2004) Spatial variability in water content at the cold–temperate transition surface of the polythermal Storglaciären, Sweden. J. Geophys. Res., 109(F2), F02009 (doi: 10.1029/2003JF000110)
Raymond, CF (1976) Some thermal effects of bubbles in temperate glacier ice. J. Glaciol., 16(74), 159171
Röthlisberger, H (1972) Seismic exploration in cold regions. CRREL Monogr. II-A2a
Schneider, T and Jansson, P (2004) Internal accumulation in firn and its significance for the mass balance of Storglaciären, Sweden. J. Glaciol., 50(168), 2534 (doi: 10.3189/172756504781830277)
Schuster, GT, Johnson, DP and Trentman, DJ (1988) Numerical verification and extension of an analytic generalized inverse for common-depth-point and vertical-seismic-profile traveltime equations. Geophysics, 53(3), 326333 (doi: 10.1190/1.1442466)
Sheriff, RE and Geldart, LP (1999) Exploration seismology. Cambridge University Press, Cambridge
Vallon, M, Petit, JR and Fabre, B (1976) Study of an ice core to the bedrock in the accumulation zone of an Alpine glacier. J. Glaciol., 17(75), 1328
Wyllie, MRJ, Gregory, AR and Gardner, GHF (1958) An experimental investigation of factors affecting elastic wave velocities in porous media. Geophysics, 23(3), 459493 (doi: 10.1190/1.1438493)