A comparison of seismic and radar methods to establish the thickness and density of glacier snow cover

Abstract

We show that geophysical methods offer an effective means of quantifying snow thickness and density. Opportunistic (efficient but non-optimized) seismic refraction and ground-penetrating radar (GPR) surveys were performed on Storglaciären, Sweden, co-located with a snow pit that shows the snowpack to be 1.73 m thick, with density increasing from ∼120 to ∼500 kg m–³ (with a +50 kg m–³ anomaly between 0.73 and 0.83 m depth). Depths estimated for two detectable GPR reflectors, 0.76 ±0.02 and 1.71 ± 0.03 m, correlate extremely well with ground-truth observations. Refraction seismic predicts an interface at 1.90 ± 0.31 m depth, with a refraction velocity (3730 ± 190 ms–¹) indicative of underlying glacier ice. For density estimates, several standard velocity-density relationships are trialled. In the best case, GPR delivers an excellent density estimate for the upper snow layer (observed = 321 ± 74 kg m–³, estimated = 319 ± 10 kgm^–3) but overestimates the density of the lower layer by 20%. Refraction seismic delivers a bulk density of 404 ±22 kgm^–3 compared with a ground-truth average of 356 ± 22 kg m–³. We suggest that geophysical surveys are an effective complement to mass-balance measurements (particularly for controlling estimates of snow thickness between pits) but should always be validated against ground-truth observations.

1. Introduction

The surface mass balance of a glacier is a sensitive indicator of climatic change (e.g. Kaser and others, 2006; Zemp and others, 2009), but campaigns of direct mass-balance measurement can be laborious and expensive (Braithwaite, 1984; Huss and Bauder, 2009). The traditional stakes-and-pits method of measuring mass balance (e.g. Østrem and Brugman; 1991; Kaser and others, 2006 involves monitoring ice melt at ablation stakes during summer field campaigns and digging snow pits during the winter to measure both snow thickness and density. Further methodological shortcomings also arise, for example if ablation stakes sink into the firn during summer melt (Østrem and Haakensen, 1999) or if the network of sample points is too sparse to be interpolated accurately across the whole glacier surface (Paterson, 1994; Barrand and others, 2010). The importance of mass-balance measurements has therefore prompted many investigations into performing them from remote platforms (e.g. Hubbard and others, 2000; Kohler and others, 2003; Hagg and others, 2004; Bamber and Rivera, 2007; Machguth and others, 2012). Nonetheless, provided that they are accurate, direct mass-balance measurements provide a vital ground-truth reference for calibrating remotely sensed observations (e.g. Box and others, 2006; Schuler and others, 2007).

Here we explore the effectiveness of a simple geophysical survey as a complement to a campaign of winter mass-balance records. Geophysical methods, including seismic and ground-penetrating radar (GPR), have broad applications in glaciology for imaging the internal and underlying structure of glaciers and ice masses, but are increasingly used for quantifying glacier physical properties (e.g. density, water content, etc.; Endres and others, 2009; Matsuoka and others, 2009; Booth and others, 2012; Peters and others, 2012). Specialized geophysical methods have been developed for characterizing the properties of snow and firn (e.g. Kinar and Pomeroy, 2007; Hawley and others, 2008; Bradford and others, 2009; Tsoflias and others, 2010; Kruetzmann and others, 2011), but here we consider how effectively a co-located set of efficient GPR and seismic refraction surveys can measure both snow thickness and density. We first investigate the accuracy of snow thickness estimates from the seismic and GPR data compared with a measurement made in a snow pit and then compare the densities obtained from a number of familiar mixing models and empirical relationships with ground-truth observations.

We emphasize that the geophysical data in this paper were not purposefully optimized with the aim of measuring snow density (as described by Harper and Bradford, 2003; King and Jarvis, 2007; Tsoflias and others, 2010; Brown and others, 2012). However, our analysis suggests that even with such non-specialized acquisitions, we can obtain a first-order approximation of snow density and accurate snow thickness and we therefore suggest that geophysical survey methods do provide an efficient complement to a campaign of mass-balance measurements. However, we also acknowledge that surveys designed specifically for measuring snow density could provide significant improvements for reliability.

2. Field Site and Geophysical Acquisitions

Geophysical data were acquired on Storglaciären (Fig. 1), a polythermal mountain glacier in the Kebnekaise massif, northern Sweden (67°54’N, 18°34’W). The glacier has a long archive of mass-balance measurements, continuous since 1945 (e.g. Karlén and Holmlund, 1996; Holmlund and Jansson, 1999; Holmlund and others, 2005; Zemp and others, 2009, 2010). Furthermore, the spatial coverage of those measurements is very dense compared with other glaciers: the coverage of winter and summer mass-balance measurement points is ∼100 and ∼-15km–², respectively (Jansson and Pettersson, 2007). Such detailed study is possible at Storglaciären given the nearby logistic base provided by Tarfala research station.

Fig. 1. Map of Storglaciären and location of geophysical acquisitions (edited from Holmlund and Jansson, 2002). Main grid coordinates are in Swedish grid (RT90), and contours are at 50 m intervals above sea level. Inset shows local grid around location of seismic refraction spread and snow pit. Stars show the locations of the seismic shots at each end of the refraction line (white dotted line, 11.5 m long). The centre of the GPR CMP gather is 2 m east of the western shot, and the grey square shows the position of the snow pit.

The GPR and seismic refraction surveys we present here were acquired in April 2012 in support of a seismic reflection profile located in the glacier’s ablation zone along its centre line (Fig. 1). Seismic refraction data play an important role in the processing of a reflection dataset, since they are used to define static corrections for removing the effect of topographic and near-surface velocity variations (Farrell and Euwema, 1984; Cox, 1999; Yilmaz and Doherty, 2001). In the context of our Storglaciären acquisition, static corrections are required to remove the effect of lateral and vertical variations in the properties of snow cover, which was present across the entire glacier at the time of survey. The purpose of the refraction survey was therefore to measure the seismic P-wave (compressional wave) velocity through the snow cover and uppermost glacier ice.

Refraction surveys could also have been used to investigate snow thickness variations along the whole seismic line, but moving and replanting geophones at such a range of surface positions would have been too time-consuming for the available survey time. Instead, a high-frequency GPR system provided a more efficient (and, at the same time, higher-resolution) means of performing this task (e.g. Dunse and others, 2009). Observations from the two datasets were then tied together using a GPR common midpoint (CMP) survey at the location of the refraction spread.

The geophysical survey coincided with the start of winter mass-balance measurements at Storglaciären. On completing the survey, the opportunity arose to dig a snow pit directly at the site of the geophysical measurements to obtain a ground-truth record of snow thickness and density. This provides a useful reference for calibrating the seismic and GPR, but also allows us to test how reliably a non-optimized (i.e. not designed specifically for measuring snow density) but efficient set of geophysical acquisitions can quantify the properties of the snowpack. We first consider how accurately the thickness of snowpack can be estimated and then investigate the derivation of snow density from GPR and seismic data. Both of these estimates require wave propagation velocity to be measured accurately; a description of our velocity analyses is therefore included in the following sections.

2. 1. GPR data

GPR acquisition was performed with a Sensors and Software pulseEKKO PRO system equipped with antennas of 500MHz centre frequency. In the common offset (CO) configuration (Fig. 2a), the antennas were mounted in a sled (at a separation of 0.23 m) and traces were triggered every 5 cm along the profile using a calibrated odometer wheel. The efficiency of GPR acquisition was such that the 1 km length of the seismic reflection line was acquired in ∼30 min by surveyors walking on foot; however, the GPR system could easily be mounted on a snowmobile to give extremely rapid acquisition. For the accompanying CMP acquisition (Fig. 2b), antennas were positioned either side of the westernmost source location in the refraction line at an initial source-receiver offset of 0.4 m, extending to 4.0 m in increments of 0.1 m. The passage of the survey sled during the CO survey had left a smooth track in the snow; antennas were placed in this track during the CMP acquisition to lessen microscale elevation variations and we estimate that offsets are accurate to ± 1 cm.

Fig. 2. 500 MHz GPR data acquired at Storglaciaren. (a) 15 m section of CO profile. Data processing involves dewow and bandpass filters, time-zero correction, Kirchhoff migration and application of automatic gain control for display. (b) CMP gather. Acquisition is centred on the 190 m CO trace and the trace interval is 0.1 m. (c) Coherence response to (b). Successive responses are to individual half-cycles of the GPR wavelet (Booth and others, 2010). Wavelet period T is measured as the travel time across three successive responses. Two reflections (an internal snowpack layer and the base of the snowpack) and a multiple are identified. Picks of stacking velocity (⊗ symbols) yield the solid trajectories in (b), which are shifted to approximate wavelet first breaks (dashed trajectories in (b)).

Two reflections are apparent in the CO and CMP data, interpreted as an internal horizon in the snowpack (upper event; travel time of 8.4 ns) and the base of the snowpack (lower event; travel time of 16.8 ns). Initially, travel times suggest that the lower event is a multiple of the upper, although application of the extended velocity analysis detailed below shows that this is not the case. Although there are several approaches to conducting CMP velocity analysis (e.g. semblance, t ² – x ² analyses; Huisman and others, 2003), here we use ‘coherence’ since it facilitates the application of an important velocity correction procedure (Booth and others, 2010). We use the term coherence to distinguish from semblance: the two analyses vary only in terms of the duration of their analysis window, with that of semblance extending across 1-11/2 wavelet periods (Sheriff and Geldart, 1999) and that of coherence being one temporal sample only. The procedure we apply is a two-phase approach, first requiring the analysis of (specifically) the coherence response to a dataset and thereafter applying t ² – x ² to time-shifted reflection hyperbolae. The key steps in this method are given here; a complete description is presented by Booth and others (2010).

Most approaches to velocity analysis yield models that are initially only strictly suitable for imaging (e.g. CMP stacking, hence the nomenclature ‘stacking velocity’; Yilmaz and Doherty, 2001) and not for quantifying physical properties (Schneider, 1971). Even in the absence of refraction and anisotropic travel-time effects (Tillard and Dubois, 1995; Becht and others, 2006) velocity errors occur because it is only the first break of a wavelet that expresses true propagation velocity. Semblance and coherence analyses respond only to the ensuing half-cycles of a wavelet (i.e. its non-zero samples), hence the velocity that is expressed is biased systematically slow (Booth and others, 2010). Likewise, time picks located at the peak amplitude of GPR wavelets and analysed thereafter in t ²–x ² analysis are prone to equivalent velocity errors. These errors propagate into estimates of interval velocity (i.e. when velocities are substituted into Dix’s equation; Dix, 1955) and, therefore, into velocity-derived physical properties (e.g. depth, water content). The correction of velocity errors then becomes especially relevant where there is a nonlinear relationship between velocity and a quantity derived from it (e.g. Topp and others, 1980; Fortin and Fortier, 2001).

The responses in a coherence panel correspond to reflection travel times along individual half-cycles of a wavelet (e.g. solid trajectories; Fig. 2c). This only occurs because the analysis window in coherence analysis is a single temporal sample (here 0.2 ns); semblance, with its longer analysis window, would not resolve such individual contributions to the overall response. As such, a coherence pick gives the travel time of the GPR wavelet along an individual half-cycle and furthermore allows the delay between that half-cycle and first break to be estimated. In Figure 2b, both picks correspond to the second response within a package of coherence peaks, suggesting that they both pertain to the second half-cycle of the GPR wavelet. Assuming that wavelet period T is constant with offset (valid given the low GPR attenuation rate in snow; Kruetzmann and others, 2011), the lag Δt between its first break and the peak of its nth half-cycle is

(1)

Consequently, subtracting 3/4T from these coherence-derived reflection trajectories leads to a better approximation of first-break travel times and hence improved velocity accuracy. Wavelet period can also be estimated directly from the coherence panel, since the responses to two half-cycles are separated by T/2 (1.20 and 0.85 ns, respectively, for the upper and lower events in Fig. 2b).

The dashed curves in Figure 2c are the approximations to first-break travel times following the application of -3/4T time shifts to the coherence-derived trajectories. The corresponding travel times are then substituted into t ²–x ² analysis, with output velocities more representative of those expressed by first breaks (Booth and others, 2010). The application of zero-phase deconvolution to our data (e.g. Schmelzbach and others, 2012) would negate the application of time shifts (i.e. peak wavelet amplitude would coincide with first break), but our procedure is essentially a coherence-based approach to wavelet estimation that is nonetheless efficient and effective.

All coherence-derived GPR velocities, and the parameters relevant to this procedure, are listed in Table 1. Following the substitution of these models into Dix’s equation, depths of 0.76 ± 0 . 0 2 and 1.71 ± 0.03 m are estimated for the two horizons; note also that the corrected travel times are no longer suggestive of a primary-multiple relationship between the two events. Uncertainties are based on the resolution of coherence responses evaluated using a Monte Carlo simulation (Booth and others, 2011).

Table 1. Velocity–depth model derived from analysis of GPR CMP gather (Fig. 2). Coherence picks are obtained directly from the coherence panel (Fig. 2c), whereas shifted picks are obtained after the application of time corrections (=–3=4T). Interval velocities are estimated by substituting shifted picks into Dix’s equation (Dix, 1955)

2.2. Seismic refraction data

Our refraction survey comprised a spread of 24 vertical-component 10 Hz geophones pressed into a fresh snow surface and then covered with loose snow to dampen wind noise and shot-generated airwaves. Geophones were installed at 0.5 m intervals (giving a spread length of 11.5 m) and data were recorded on a Geometrics GEODE system. The seismic source was a 14 kg sledgehammer impacting a 0.5 m x 0.5 m plastic plate and we triggered the recording system using a piezoelectric sensor mounted on the hammer shaft. We recorded two shot gathers, with the source placed at the first and then at the last geophone location.

Seismic refraction data (Fig. 3a and b) show good signal-to-noise ratio and impulsive first breaks at almost all geophones. The derivation of a velocity model from seismic refraction data requires direct- and refracted-wave travel times to be identified. Travel-time picks for both records show a change in the gradient of first-break picks at ∼5m offset. First breaks are interpreted as direct waves for the first 4.5 m offset and as the first refraction thereafter. The GPR data (Fig. 2) show negligible dip at the base of the snow cover, hence we consider only one-dimensional variation in seismic velocity.

Fig. 3. (a, b) Seismic refraction data and (c) travel-time analysis. Travel-time picks (triangles in (a) and (b), circles in (c)) correspond to direct waves across the first 4.5 m of the spread and to the first refraction thereafter. The geophone installed at 4 m (grey trace) consistently records noisy data, hence it is excluded from this analysis. Direct-wave intercept times are negative in (c), implying a time delay in the system trigger. The inverse slope of each straight line gives the apparent velocity expressed by the first breaks.

Figure 3c shows the best-fit linear trends to direct- and refracted-wave travel times. In each record, the zero-offset direct wave has a small negative intercept time (–0.41 and –0.35 ms). This is attributed to a trigger error, representing the time taken for the switch inside the piezoelectric sensor to close, although non-causality introduced by the recording system is an additional possibility. In any case, the trigger error does not adversely impact our first-break picking and is removed prior to depth estimation.

The velocity of the direct wave v ₁ is simply estimated using the linear relationship

(2)

where Δx and Δt are the ranges of offset and travel time, respectively, over which the direct wave is observed. Prompted by the expected density gradient in the snowpack, we explored fitting power-law and linear velocity-depth functions (Cox, 1999) to direct-wave travel times, but the extra complexity was unjustified since the observed velocity gradient was negligible.

The velocity of the refracted wave v₂ is also established using Eqn (2), and the depth d to the refracting interface is

(3)

where t_i is the intercept time of the refraction trend on the zero-offset trace (e.g. ≈ 3.3ms on the right-hand side of Fig. 3c). Trigger errors are removed by adding the relevant absolute delay to each t_i prior to substitution into Eqn (2).

The two refraction-shot gathers give consistent velocity models (Table 2), with direct-wave velocity 1000 ± 20 m s–¹ and snow thickness ≈ 1.90±0.3m. The mean inverse slowness expressed by slopes of the refracted arrival is 3730 ± 190 m s–¹. The uncertainty in each quantity is derived from the standard error in the least-squares straight-line fitting, with the possible range of velocities used to establish the range of depth uncertainty. The maximum depth observed in the GPR (1.71 ± 0 . 0 3 m) is consistent with the depth of the seismic refracting layer, at least within the broader uncertainty bounds of the latter.

Table 2. Velocity–depth model derived from analysis of two seismic refraction datasets (Fig. 3)

4. Comparison with Observed Snow Depth

Encouragingly, seismic refraction and GPR methods indicate an interface at similar depth and we now compare this with the snow depth recorded in the snow pit. The pit was dug into undisturbed snow ∼10m south of the midpoint of the GPR CMP gather (67°54.150’ N, 18°34.360’E; Fig. 1). We do not expect that snow thickness will vary greatly between the two measurement points, since no discontinuities on this spatial scale were seen anywhere in the GPR CO profile.

Figure 4 shows the variation of density with depth recorded in the snow pit together with representations of the depth models derived from the geophysical surveys (in Fig. 4 the uncertainty in each velocity and depth is represented by the line thickness). Density measurements were made in the field from snow cores extracted at 5 cm intervals from the vertical wall of the pit. There is the expected increase of density with depth, although the maximum value (500 kgm–³) is observed 0.15 m from the base of the pit (density comparisons are revisited in Section 5). The snow is 1.73 m thick and an ice lens (0.05-0.10 m thick) is observed at its base; the pit terminates 1.83 m below the surface on depth hoar.

Fig. 4. (a) Variation of density with depth in the snowpack and (b, c) velocity–depth models derived from (b) GPR and (c) seismic velocity–depth models for comparison. In (b) and (c), the thickness of the black and grey bars represents the uncertainty in velocity and depth, respectively (transparency increases towards the edge of the precision range). The GPR accurately represents density transitions, although strong fluctuations towards the base of the snowpack are not recovered. For the seismic dataset, the estimate of the snow thickness is accurate within the precision range. The seismic direct wave is assumed to sample to 1.4 m depth, hence the shallower seismic velocity trend is dashed between 1.4 and 1.7 m.

The GPR velocity model gives a very accurate estimate of the thickness of the snowpack, and depths agree within the radar precision range. The internal snowpack reflector correlates very well with a zone of increased density (∼50kgm–³) observed between 0.73 and 0.83 m. No reflection from the base of the ice lens is observed, unsurprising since its 0.1 m thickness is below the resolution of a 500 MHz wavelet travelling at any plausible velocity for ice.

The depth of the seismic refraction is consistent with that of the depth hoar (overestimated by 0.07 m), at least within the precision of the former. The velocity that the refracted wavelets express, ∼3730ms–¹, is typical for slightly porous glacier ice (as observed for the uppermost ice at Storglaciaren by Gusmeroli and others, 2010), hence the refraction could also originate from glacier ice beneath the depth hoar. Note that the snow is shown as a dashed line beyond 1.4 m depth given considerations to the Fresnel volume of the seismic wavelet, as explained in Section 5.

5. Comparison with Observed Snow Density

As with depth, density is derived from wavelet velocity. We consider the experimental error in our velocity analysis to be small since our depth estimates are consistent, although allocating (at most) two velocities to the snowpack is probably an over-simplification of its fine-scale structure (e.g. Bradford and others, 2009). Therefore, we test how well our surveys characterize the first-order properties of the snow-pack, rather than replicating its detailed density structure.

Velocity is usually converted to density via mixing models (e.g. Riznichenko, 1949; Looyenga, 1965; Topp and others, 1980; Endres and others, 2009) or empirically defined relationships (e.g. Wyllie and others, 1956; Kohnen, 1972; Gardner and others, 1974; Tiuri and others, 1984; Sihvola and others, 1985; Fortin and Fortier, 2001). A further source of error is introduced by such an approach, since an assumption made in formulating a mixing model or in assigning an experimental variable may be invalid at Storglaciaren. We therefore investigate a range of widely used relationships between velocity and density. It should be noted that we do not advocate one method over another, but highlight the importance of acknowledging the site-specific variability in snow properties. The density estimate derived from each model is recorded in Table 3.

Table 3. Snow density (kg m^–3) as estimated from GPR and seismic velocities using different mixing models (CRIM and Wyllie time-average) and empirical relationships. Upper-layer values relate to the depth interval 0-0.76 m, in which the reference density is 321 ±74kgm^–3; lower-layer values relate to the depth interval 0.76-1.71 m, in which the reference density is 414 ±44 kg m^–3 Refraction seismic values relate to the depth interval 0-1.4 m, in which the refernce density is 356 ±44 kg m^–3. Values in square brackets show the percentage difference between measured and reference densities

5.1. GPR density derivation

The velocity of a GPR wavelet through porous material is implicitly related to its density. Velocity is explicitly sensitive to the dielectric permittivity of a material, although changes to the porosity in that material influence both the bulk dielectric permittivity and density (e.g. Sihvola and others, 1985). For snow, air-filled porosity causes both dielectric permittivity and density to decrease, whereas water-filled porosity causes them both to increase. Throughout this analysis, we assume that pore spaces in the snowpack are completely filled with air and, as such, snow can be described as a two-phase ice/air mixture.

We also assume that by allocating a single GPR velocity to each layer we observe, the density we estimate is the average of the densities recorded therein. The average ground-truth density recorded within the upper and lower GPR intervals is therefore 321 ±74 and 414±44kgm–³, respectively (the uncertainties are the standard deviation in each depth range).

We first estimate density using the complex refractive index method (CRIM), which is a basic but widely applied mixing model (e.g. Roth and others, 1990; Huisman and others, 2003). Harper and Bradford (2003) formulate a CRIM relationship to link the observed snow velocity v _snow to snow density p_snow:

(4)

where p_ice is ice density and v _air and v _ice are GPR velocities through air (=0.300 m ns–¹) and ice, respectively. At Storglaciaren, Gusmeroli and others (2010) use p_ice = 917 kgm–³ and measure v _ice ≈ 0.17mns^–1 for the uppermost glacier ice. Table 3 shows p_snow values estimated by Eqn (4). The density of the upper layer is clearly well characterized (difference of -0.6%), whereas that of the lower layer is less accurate (difference of +25.4%).

Fortin and Fortier (2001) present a series of empirical relationships (Ambach and Denoth, 1972; Tiuri and others, 1984) that express density in terms of relative dielectric permittivity of snow _{r snow} of the form

(5)

where K is a constant of proportionality, experimentally determined to be 2.0 (Tiuri and others, 1984; Fortin and Fortier, 2001) or 2.2 (Ambach and Denoth, 1972), with p_snow expressed in gcm^-3. Treating snow as a two-phase mixture between air and ice, the only control on bulk dielectric permittivity is the amount of snow within a given volume, since _r for air is 1 and its density is negligible. Assuming low-loss GPR propagation (i.e. that electrical conductivity is negligible), with c=0.300mns-1, Eqn (5) is rearranged as

(6)

with the factor of 1000 introduced to express density in kgm–³. Table 3 shows p_snow for K=2.0 and 2.2 for each layer. In the best cases, the lower K provides a good density match in the upper layer (-5% difference), whereas the higher K improves the match in the lower layer (+14.7% difference). Fortin and Fortier (2001) note that their experimental data have a high degree of scatter, hence K may be poorly constrained, and different values could be appropriate in the two layers of the snowpack, potentially warranting a site-specific study.

Finally, we consider the empirical relationship described in Kovacs and others (1995), which states that

(7)

with the same underlying assumption as the previous relationship. Again, we express ε _{r snow} in terms of v _snow and rearrange for p_snow, and show densities in Table 3. The density of the upper layer is better characterized than that of the lower layer (differences of -2.5% and +18.0%, respectively, with respect to reference densities).

5.2. Seismic-derived densities

The velocity of a seismic wavelet is explicitly related to the density of a material, since density appears in the equation for seismic velocity (Sheriff and Geldart, 1999). Consequently, seismic velocity may be explicitly linked to density variations, although mixing models and empirical relationships remain commonplace due to complicating factors including pore fraction, connectivity and fluid material (e.g. Gardner and others, 1974; Mavko and others, 2009).

As with the GPR interpretation, the depth range over which our seismic velocity estimate averages density should be appreciated. Spetzler and Snieder (2004) define the ‘Fresnel volume’ as the finite volume around a geometric ray path that influences the propagation of a band-limited wavelet. The radius r of the Fresnel volume is:

(8)

where λ is the seismic wavelength and x is the offset at which the wavelet is recorded. For a direct wave, the Fresnel radius is therefore the deepest position in the subsurface that can influence the recorded wavelet (Gusmeroli and others, 2010), although Sheriff and Geldart (1999) state that significant contributions probably only arise within a radius . We measure direct-wave velocities to a maximum offset of 4.5 m, and the wavelength of our seismic pulse is ∼4.4m (frequency 225 Hz, propagating at 1000 ms^-1); the radius of the Fresnel volume therefore increases from 0.45 m when observed at 0.5 m offset, to 1.4 m at maximum offset. We therefore expect that the direct wave will average the physical properties of the top 1.4 m of the snowpack (hence the dashed line for depth will exceed 1.4 m in Fig. 4c), and the average density observed in this range is 356 ± 6 6 kgm^-3.

We first compare a density estimate made using the Wyllie time-average equation (analogous to the CRIM relationship in Eqn (4)), which states

(9)

where is fractional porosity (Wyllie and others, 1956). Values of v _air = 330ms^-1, v _snow = 1 0 0 0 ± 2 0 ms^-1 and v _ice = 3 7 3 0 ± 190ms^-1 are substituted into Eqn (9), and a fractional porosity of 0.265 (i.e. 26.5%) is obtained. Since the density of air is negligible, snow density is then simply approximated as p_snow = (1-Φ)p_ice=674kgm^-3 (assuming again that p_ice = 917 kg m^-3; Gusmeroli and others, 2010). This value is almost double that measured in the snow pit (Table 3). We suggest that this mismatch arises because such CRIM- and Wyllie-type mixing models assume that there is a continuum between highly porous snow and non-porous glacier ice such that v _snow=v _ice when all the air-filled pore spaces are removed. While this is appropriate for GPR velocity, it is inappropriate for the seismic case since snow undergoes further densification before becoming ice. As such, our measured v _ice is inappropriate in this context. The same assumptions are made in the mixing model of Riznichenko (1949), hence we do not explore its application here.

Compared with the GPR case, there are few empirical relationships defined for evaluating snow density from seismic velocity. One such relationship between P-wave velocity and density is proposed by Kohnen (1972), and has been used in several firn and ice settings (e.g. King and Jarvis, 2007; Rege and Godio, 2011). The relationship states

(10)

where p _z and v_z are the density and seismic velocity at depth z, and is more appropriately defined for seismic velocities in ice. We substitute v _snow= v_z= 1 0 0 0±20 ms^-1 and obtain a better estimate of the reference snow density (+13% difference; see Table 3).

Other relationships between P-wave velocity and density exist (e.g. Gardner and others, 1974; Castagna and others, 1985) but are typically established for water-saturated rocks rather than unconsolidated sediment that may serve as an analogue for snow. For example, the system of empirical relationships determined by Carroll (1969) was derived for rocks with densities between 1600 and 2700 kgm^-3 (Wadhwa and others, 2010), much higher than the values observed in our snow pit. When density relationships are derived for unconsolidated material (e.g. Prasad, 2002; Zimmer and others, 2007), they are more generally established in terms of a combination of P-wave and S-wave (shear wave) information (King and Jarvis, 2007; Bradford and others, 2008; Tsoflias and others, 2010). S-waves are insensitive to pore fluid (Mavko and others, 2009), hence the properties they perceive pertain only to the grains (or, in our case, the snow crystals). In contrast, P-waves perceive the bulk density of the whole medium and are therefore sensitively affected by changes in pore fluid (i.e. residual liquid water in the snow could be erroneously interpreted as a high-density anomaly in the snow cover). The most representative characterizations of snow density therefore include both P- and S-wave velocity information (King and Jarvis, 2007; Tsoflias and others, 2010).

For our data, the requisite S-wave information could be inferred from multi-channel analysis of surface waves (MASW) (Tsoflias and others, 2010), which would facilitate analytic (rather than empirical) derivation of physical properties. However, given that our refraction experiment was designed for recording P-waves, surface wave amplitudes are clipped and it is impossible to determine the amplitude of individual frequency components, a fundamental step in implementing MASW. However, clipping is not widespread in the larger seismic acquisition and this approach could be explored in future research.

6. Discussion

In these analyses, we attempted to recover accurate estimates of snow thickness and snow density independently of ground-truth values observed in a co-located snow pit. The more successful of these analyses was the estimation of thickness, specifically from the GPR data: the estimated and observed snow thicknesses agreed to within centimetres. The thickness derived from the seismic data was also consistent with the ground-truth observation although rather less precise than the GPR estimate (expected given the wavelengths involved: ∼4.4 and ∼0.4m in the seismic and GPR cases, respectively).

For the upper layer of the GPR model, all but one of the derivation methods we apply delivers snow density accurate to within ±5%. However, the accuracy for the lower layer is much worse and density is consistently overestimated. For characterizing the bulk properties of the snowpack, our preferred method is therefore refraction seismic: there is overlap between the precision range of the observed and seismic-derived estimate of density, even if it is biased to slightly higher values. However, the discrepancies between the observed and estimated snow properties point towards some significant sources of error in these analyses.

First, we greatly simplify wavelet propagation through the snowpack. In the GPR case, we assume implicitly that the variation of travel time with offset is hyperbolic. While the effects of this are clearly small enough to have no significant impact on our depth conversion, the potential for error is greater in the nonlinear relationship between velocity and density. There are strong density fluctuations (amplitude ∼10 kg m–³) in the deepest 0.2 m of the snow pit (Fig. 4a) and this is likely to introduce strong non-hyperbolic moveout into travel times (e.g. Bradford, 2002), specifically causing GPR rays to refract closer to the vertical given the progressive decrease in velocity (Becht and others, 2006). Conventional hyperbolic velocity analysis, even with the corrections applied here, is therefore unlikely to describe the reflection travel time adequately and a higher-order travel-time description is required (e.g. Grechka and Tsvankin, 1998). Failing this, the offset range of the CMP gather could be restricted such that it does not exceed the depth of the target reflector (i.e. honouring the short-spread approximation; Taner and Koehler, 1969), but non-hyperbolic moveout terms not be completely eliminated and the resolution of the coherence analysis would be adversely affected (Booth and others, 2011). These issues could be overcome using the ray-based modelling approach of Brown and others (2012), where a best-fit velocity model is obtained for numerous CMP-derived velocity functions, although the accuracy of associated density estimates still depends on the suitability of a mixing model (see next paragraph). Seismic travel-time analysis should strictly consider curving ray paths, although we determined that this was not significant for the data interpreted here.

Second, our analysis shows the need to appreciate the suitability of an empirical relationship before applying it to density derivation. We make no suggestion that our observations invalidate the established models, but instead highlight the requirement either for the derivation of site-specific empirical relationships or for ground-truth control in geophysical analysis. As a minimum, a range of empirical relationships should be explored such that the experimental variability between them can be appreciated. We note that empirical velocity–density relationships have also been applied in the opposite sense to our approach (i.e. where a ground-truth density is used to calibrate GPR velocity; Dunse and others, 2009) and we would recommend equivalent caution in such cases.

Finally, we neglect in these models the possibility of a multiphase pore fluid (i.e. both air and liquid water within pore spaces) (Bradford and others, 2009). While it is trivial to add a water-phase term into either the CRIM or Wyllie time-average equations, each of the empirical relationships used here assumes that the snowpack is described by a snow–air blend. In the GPR case, the effect of liquid water would be to reduce the propagation velocity, leading to an overestimate of density. In the seismic case, there is a highly nonlinear relationship between velocity and water saturation, hence the net effect on a density estimate is difficult to predict (Bradford, 2010); nonetheless, the most accurate density estimates would also consider liquid water in the snowpack.

We therefore consider that geophysical methods can complement mass-balance measurements, but should not be applied without ground-truth measurement. The efficiency of GPR acquisitions suggests that they could be effectively applied to tie together thickness estimates made in a sparse array of snow pits (e.g. Machguth and others, 2006), although the benefit of acquiring CMP gathers for velocity control is also clear when considering the resulting accuracy of depth conversions. The additional depth control would therefore serve to improve the glacier-wide interpolation of mass-balance measurements or to extend the observations into areas of particularly sparse control. Seismic acquisitions are much less efficient than GPR, but may provide better density estimates given that GPR velocity is further removed from density. The use of MASW may represent an analytic means of measuring snow properties and we speculate that the ‘optimal’ survey design could be sparse MASW acquisitions that are tied to GPR profiles for extrapolating properties. However, for the datasets interpreted here, reliable density estimates would require validated mixing models and/or empirical relationships.

7. Conclusions

Geophysical acquisition can provide a good complement to a campaign of winter mass-balance measurements, but not without the control of ground-truth data. However, the accuracy of thickness estimates from GPR is encouraging and the method could provide valuable control for spatial interpolations between snow-pit sites or in areas of sparse control, particularly since GPR acquisition is efficient. The measurement of density from either seismic or radar wavelet velocity should be considered a first-order estimate if survey design is not optimized for density measurement and/or the method of converting velocity to density is not calibrated. The most useful complement could therefore be a dense network of GPR profiles tied into seismic and/or GPR surveys that are specifically designed for density measurement.