Freely available online through the SEG open-access option.

ABSTRACT

The water balance in alpine watersheds is dominated by snowmelt, which provides infiltration, recharges aquifers, controls peak runoff, and is responsible for most of the annual water flow downstream. Accurate estimation of snow water equivalent (SWE) is necessary for runoff and flood estimation, but acquiring enough measurements is challenging due to the variability of snow accumulation, ablation, and redistribution at a range of scales in mountainous terrain. We have developed a method for imaging snow stratigraphy and estimating SWE over large distances from a ground-penetrating radar (GPR) system mounted on a snowmobile. We mounted commercial GPR systems (500 and 800 MHz) to the front of the snowmobile to provide maximum mobility and ensure that measurements were taken on pristine snow. Images showed detailed snow stratigraphy down to the ground surface over snow depths up to at least 8 m, enabling the elucidation of snow accumulation and redistribution processes. We estimated snow density (and thus SWE, assuming no liquid water) by measuring radar velocity of the snowpack through migration focusing analysis. Results from the Medicine Bow Mountains of southeast Wyoming showed that estimates of snow density from GPR (360±70kg/m3) were in good agreement with those from coincident snow cores (360±60kg/m3). Using this method, snow thickness, snow density, and SWE can be measured over large areas solely from rapidly acquired common-offset GPR profiles, without the need for common-midpoint acquisition or snow cores.

INTRODUCTION

In many alpine systems, snow is the dominant form of precipitation (e.g., Serreze et al., 1999). Snowmelt drives surface hydrology in much of western North America, where the annual hydrograph is dominated by high spring runoff driven by snowmelt. The warming climate in much of the western U.S. is altering snow accumulation and melt timing, producing declining mountain snowpack (Mote et al., 2005), higher winter flows, earlier spring runoff peaks, and decreased summer flows (Rood et al., 2008). Improved understanding of snow accumulation, redistribution, and runoff during melting is thus of growing importance to human societies (e.g., Bales et al., 2006).

The water contained in snow is measured as the snow water equivalent (SWE), which has units of depth and is the product of snow density (relative to water) and snow depth. Hydrologic models of runoff from mountain watersheds depend on accurate understanding of SWE at the watershed scale (Jost et al., 2007), which is challenging to acquire. The accurate estimation of snow covered area can be readily obtained through a variety of remote sensing techniques, but determining water content in the pack has heretofore required intensive sampling. The controls on SWE distribution in mountain watersheds are imperfectly understood, but they include net radiation, slope aspect, slope steepness, elevation, canopy cover, and wind redistribution (e.g., Anderton et al., 2004; Molotch et al., 2005; Jost et al., 2007; Rutter et al., 2009; Jonas et al., 2009; Dixon et al., 2014). Common methods for measuring SWE include (1) snow cores (e.g., Dixon and Boon, 2012), which typically provide average snow density for the entire snow pack at a point, or along a snow course; (2) snow pits, which can provide detailed snow density versus depth at a single location (e.g., Elder et al., 1998); and (3) snow pillows, permanent installations that use a weighing sensor and snow depth measurement to determine SWE, such as at the U.S. network of snowpack telemetry (SNOTEL) stations (e.g., Serreze et al., 1999; Andreadis and Lettenmaier, 2006). Snow pits and snow courses are labor intensive and typically rely on interpolation across hundreds of meters or more to characterize watersheds (e.g., Elder et al., 1998; Balk and Elder, 2000; Molotch et al., 2005). Snow depth is variable at the 10 m scale, so that snow cores can be in error if sampling strategies do not account for this variability (Lopez-Moreno et al., 2011). Snow pillows provide valuable time series of snow accumulation and melt, but they are spatially limited, expensive to install and to maintain, and subject to errors from bridging effects (Johnson and Marks, 2004). Remote sensing from satellites has been used to estimate SWE over large areas (e.g., Pulliainen and Hallikainen, 2001; Pulliainen, 2006; Dietz et al., 2012), but these methods often suffer from calibration issues and uncertainties regarding snow grain size and vegetation cover (Foster et al., 2005). More recently, upward-looking GPR data (Schmid et al., 2014, 2015) and reflected signals from campaign-quality GPS receivers (Larson et al., 2009; Larson and Nievinski, 2013; Rittger et al., 2013) have been used to estimate snow depths and SWE with promising results, but these techniques require independent estimates of snow depth or density.

GPR data provide one means of rapidly measuring snow properties over large areas. GPR instruments send an electromagnetic (EM) pulse from the snow surface down through the subsurface and record reflections produced by dielectric permittivity contrasts (e.g., Jol, 2009). Numerous studies have shown that ground-penetrating radar is effective at imaging snow and glacial ice (Frezzotti et al., 2002; Pivot et al., 2002; Harper and Bradford, 2003; Marchand et al., 2003; Marshall et al., 2005; Machguth et al., 2006; Bradford et al., 2009; Lundberg et al., 2010; Kruetzmann et al., 2011; Sundstrom et al., 2012; Gacitua et al., 2013; Sold et al., 2013; Sylvestre et al., 2013; van Pelt et al., 2014). Snow GPR surveys have been conducted from various platforms, including pulled behind an operator on skis or snowshoes (e.g., Bradford et al., 2009), towed behind snowmobiles (e.g., Gacitua et al., 2013), and, more rarely, flown from helicopters (e.g., Sold et al., 2013). In snowpack, GPR data provide accurate measurements of the two-way traveltime of EM waves to reflective boundaries within, and at the base of, the snow layer. Converting traveltimes to depth requires measurement or estimates of the EM velocity. These estimates typically come from either independent measurements of snow density (which bears a simple relationship to radar velocity for dry snow) or an assumed density profile, or they are measured with common-midpoint (CMP) gathers either collected with single-channel (labor intensive) systems at sparse locations or multichannel systems (not labor intensive) at every point in the survey (Harper and Bradford, 2003; Bradford et al., 2009).

We present an approach to estimating SWE over long mountain transects using a novel snowmobile-mounted GPR unit, which allows rapid acquisition of data at the watershed scale. We use migration velocity analysis (MVA) to determine the average radar wave velocity in the snow layer, which can be related to snow density in a straightforward way, assuming dry snow. Our approach is similar to that of Bradford and Harper (2005), who use the MVA of low-frequency (25 MHz) radar data to determine the density of glacial ice. In the present work, we apply MVA to high-frequency (800 MHz) data to estimate average snowpack density; a compilation of snowpack density values versus snow depth compares favorably with snow density measured on coincident cores. We show that the combination of high-resolution snowpack images with MVA provides an independent method for estimating SWE that does not require the collection of snow cores or CMP GPR data.

SETTING

The study site is a partially wooded region of the Snowy Range in the Medicine Bow Mountains of southeastern Wyoming (Figure 1). The area receives approximately 90 cm in annual precipitation, most of which (74%) falls as snow, with peak SWE accumulation typically occurring in early May (USDA, 2015) (Figure 2). Stream systems are dominated by snowmelt, with most runoff occurring in early summer as the snowpack melts. The median number of days this region has snow on the ground is 236, with a mean annual temperature of 0.6°C. Elevations at the study site range from 3100 to 3300 m. Vegetation is composed of mixed forest ecology with grassland and shrub interspersed among spruce, lodgepole pine, and subalpine fir.

METHODS

GPR acquisition

Here, GPR data were acquired by mounting two GPR antennas in front of a Polaris RMK 600 snowmobile on a custom aluminum frame (Figure 3). Mounting the antennas in front of the snowmobile allows for maximum mobility and ensures that measurements are taken on pristine snow (unlike systems in which a GPR unit is dragged on a sled behind the snowmobile). In addition, having the antenna mounted above the snow surface enables recording of the reflection from the snow surface. Data can be acquired at speeds of 6–8m/s, allowing rapid acquisition of data on long transects. There is no need to acquire data on perfectly straight lines. We acquired approximately 15 km of GPR data on three field days (2, 6, and 27 March) in 2014. Typical snowmobiling speeds during acquisition were 4–8m/s.

Radar data were acquired with two common-offset, MALÅ pulse-radar systems, with nominal antenna frequencies of 500 and 800 MHz. We present data from three lines (Figure 1) and only from the 800 MHz antennas, which have superior resolution of snow layering. We recorded approximately 82 ns records with 392 samples per trace (0.210 ns sample interval) on lines 25 and 26, and approximately 80 ns records with 512 samples per trace (0.157 ns sample interval) on line 34. The transmitting and receiving antennas are located in a single unit at 14 cm separation. Acquisition was done by firing continuously in time at approximately 100 traces per second. Data were converted to constant-distance sampling by interpolation using positions recorded once per second by a 12-channel Trimble R8 GPS antenna fixed to the rear of the snowmobile and integrated into the GPR acquisition system. Real-time kinematic corrections were applied to the GPS data after acquisition, resulting in absolute horizontal positioning errors of approximately ±2m after processing. Interpolated trace spacing was approximately 2.5 cm.

GPR data processing

To create migrated images of snow structure (e.g., Figure 4), GPR data were processed through a sequence consisting of the following steps: (1) reset trace time zero, (2) transform traces to equally spaced traces in distance, (3) median filter (global background removal) to suppress antenna “ringing,” (4) time-variable gain correction, (5) mute traces above the snow surface reflection, (6) frequency band-pass filter from 0.3 to 1.8 GHz, (7) f-k migration over a suite of constant-velocity models, and (8) conversion to depth using the best-fitting velocity (usually 0.23m/ns). In addition to these standard steps, additional processes were applied as needed, including interpolation to fill in null traces, trace amplitude equalization, and static shifting of traces to account for topography. Most processing was carried out in the matGPR processing package (Tzanis, 2010), with additional modules (static shifting, trace interpolation) written in MATLAB.

An additional processing sequence was conducted to enable rapid estimation of total snow depth along transects, via automated tracking of the snow surface and ground surface. Prior to automated tracking, we applied (1) trace equalization to the median trace amplitude, (2) calculation of instantaneous amplitude (the envelope of the Hilbert transform), and (3) application of a 50-trace mean spatial filter across the section. This produced a smooth amplitude map on which we tracked amplitudes that crossed a user-selected minimum amplitude threshold, in time windows corresponding to the snow surface reflection and the ground surface reflection on each profile.

Snow cores

On several transects, we established snow courses and took snow depth and density measurements to ground truth GPR estimates of snow properties. These snow measurements were taken following standard methods, and observations of depth, density, and SWE were made at each site. A standard metric Federal snow sampler was used to extract snow cores. Upon arrival at a site, the necessary number of snow tubes was assembled, ice or snow was cleared from the tube, and a tare weight was taken using a digital scale. The tube was inserted to the base of the pack, rotated slightly to extract a small soil plug, the snow depth was noted, and the snow core removed. If a soil plug was not observed or the length of the snow core was less than 75% of the observed depth, the core was discarded and a new observation was taken. Upon successful extraction of a core, the soil plug was carefully removed and the weight of the core was recorded. The SWE was determined by subtracting the tare weight from the recorded weight of the tube and snow. Potential errors in using this technique include oversampling by improperly inserting the tube, overpenetrating the soil, and human error in transcribing observations. Woo (1997) compared snow tube measurements to snow pits and found that the error is less than 5%.

Estimation of radar velocity from migration focusing analysis

Migration is a process that is widely used in reflection seismology and GPR processing to place observed reflections in their correct subsurface position (e.g., Claerbout, 1985; Sheriff and Geldart, 1995). Migration is necessary because diffracted and reflected energy from a subsurface target is recorded at many surface positions along a transect, not just immediately above the reflection point, so that arrivals in a raw data section extend beyond their true spatial positions. For a point diffractor in the subsurface (i.e., an object whose size approximates the wavelength of the radar wave; approximately 30 cm for an 800 MHz wave in snow), the return comprises a hyperbola whose shape depends on the depth of the diffractor and the velocity of the medium above it. Migration of the hyperbola using the correct subsurface velocity collapses the image to a point at its true subsurface position. Because the shape of the hyperbola is sensitive to the radar velocity, it can be used to infer the velocity of the snow. Here, we use MVA to determine the average velocity above diffractors (e.g., Bradford and Harper, 2005).

In our data, numerous diffractions are visible at or just above the ground surface reflection (Figure 5); these likely come from the rocks, logs, and vegetation buried beneath the snow pack. We use these diffractions to estimate radar velocity with the following workflow: (1) Migrate each section over velocities from 0.15 to 0.29m/ns, at intervals of 0.005m/ns. (2) Scan the migrated panels for clear, isolated diffractions. (3) For each diffraction, inspect the migrated panels and select the optimal focusing velocity and the arrival times of the snow surface reflection and the diffraction. (4) Apply the Dix equation (Dix, 1955) to separate the snow layer velocity from the air layer velocity (0.299m/ns), given the height of the antenna above the snow surface.

Calculation of dielectric permittivity and snow density

Radar velocity is sensitive to snow density (e.g., Tiuri et al., 1984; Kovacs et al., 1995; Bradford et al., 2009). In dry snow (no liquid water), the radar velocity depends only on the relative proportions of air and ice and their dielectric constants, and snow density can be determined solely from the radar velocity (Denoth et al., 1984; Tiuri et al., 1984; Mätzler, 1996). If liquid water is present, additional information is required, such as the relative values of the real and imaginary parts of complex permittivity (Bradford et al., 2009). In our analysis, we assumed dry snow (see the justification below under the section “Time-frequency analysis”).

We use the radar velocities determined by migration focusing analysis to determine the snow dielectric permittivity, using v1/εμ0 (Bradford et al., 2009), where ε is the dielectric permittivity of snow and μ0 is the magnetic permeability of free space. We then used an empirical equation (Tiuri et al., 1984) to calculate density from the dielectric constant:  
κ=1+1.7ρ+0.7ρ2,
(1)
where ρ is the snow density in gm/cm3 and κ is the real component of the dielectric constant (ε/ε0, where ε0 is the dielectric permittivity of free space). We then compiled these estimates from many sections to derive a snow density-depth relation.

Time-frequency analysis

The approach for calculating snow density outlined above assumes dry snow. Air temperature data recorded at the Brooklyn Lake SNOTEL station, located approximately 1 km from our transects, show that average air temperatures in the area were typically well less than 0°C during our data acquisition, supporting the dry snow assumption. We also tested for the presence of liquid water by computing time-frequency plots of the GPR data. As described by Bradford et al. (2009), dry snow has a real dielectric permittivity, but liquid water in the snowpack introduces a nonzero imaginary component of the dielectric permittivity. This creates attenuation of the GPR signal in the snowpack, which causes a frequency shift toward lower frequencies in the snow-ground reflection, when compared to the air-snow reflection. Thus, if significant liquid water is present in the snowpack, we should observe a lower central frequency for the ground surface reflection than for the snow surface reflection.

We visualized the frequency content of our sections by calculating time-frequency plots (Figure 6). The plots show the centroid frequency of the returns (e.g., Irving and Knight, 2003), computed using the S-transform (Stockwell et al., 1996). There is no significant difference between the centroid frequency of the snow surface reflection and any deeper reflections, indicating that the snowpack during our data acquisition was dry.

We note that the time-frequency analysis was conducted on processed GPR data. We also analyzed the unprocessed data and found no change in centroid frequency between the snow surface reflection and ground reflection. The processing steps applied here (including median filtering) did not produce substantial shifts in the centroid frequency of the snow surface or ground reflection.

RESULTS

Images of snowpack

The radargrams show detailed images of snow layers through the entire depth of the snowpack (Figure 4). The data in Figure 4 are plotted relative to the height of the antenna, with no corrections for variable topography along the transect. This plotting mode is most efficient for showing the details of snow layering in a single figure, but dips are not necessarily accurate: The snow surface appears flat even though it may change in elevation along the transect. Prominent reflections come from the top of the snow surface and from the snow/ground interface, which generally produces the strongest reflections on the images. Between the snow and ground surfaces are numerous reflections that are laterally continuous and correspond to layering in the snowpack. Although we do not have detailed measurements of snow density with depth on our transects, it is reasonable to infer that subtle contrasts in the snow density, and thus the dielectric constant, at the centimeter to decimeter scale produce the reflections (Harper and Bradford, 2003) (the wavelength of 800 MHz radar waves in snow with a density of 360kg/m3 is 28 cm, implying a Rayleigh criterion quarter-wavelength resolution of 7cm). These density contrasts likely correspond to individual snowfall events and could be produced or enhanced by slight contrasts in the initial snow density, ablation, or wind erosion, icy layers caused by prolonged exposure to sunlight between snowfall events, compaction during burial, or metamorphism during snowpack evolution (e.g., Harper and Bradford, 2003).

Comparison of snow depths: GPR versus snow cores

As a validation check, we compared the snow depth estimates from GPR line 34 to those measured directly on nine snow cores located with 1 m of the line (Figure 7) (core A1 was located 11 m off line and was not included in the analysis). Snow core depths agree with estimated depths from GPR to within 10% (average mismatch of 13 cm in average snowpack of 133 cm). Mismatches are likely due to (1) errors in snow core measurements, (2) errors in GPS positioning of the GPR line, and (3) the use of a constant lateral velocity to depth convert the GPR image, ignoring possible lateral variability in snow density and velocity. The largest misfits are at snow cores A2 and A3, which lie over the steepest ground surface, where small errors in lateral positioning could lead to significant errors in predicted depth from the GPR data. The overall agreement shown in Figure 7 indicates that the GPR data provide a robust estimate of snow depth; indeed, the GPR image offers a much more detailed view of the shape of the snow/ground interface than is possible even from densely sampled snow cores (the prominent subsurface peak at 70 m distance, for example, is missed by the snow cores).

Migration focusing analysis and snow depth-density relationship

Diffractions in migrated images collapse to a “point” when the correct migration velocity is selected. In our data, a “point” comprises a trough-peak-trough (white-black-white) triplet with a time duration that corresponds to the radar source pulse length and a width described by the migrated Fresnel zone (here, approximately 15 cm wide). An example is shown in Figure 5. Here, a clear diffraction is observed centered at 409 m with its apex between 13 and 14 ns two-way traveltime. When the migration velocity is too low, the image is undermigrated, leaving a residual diffraction. This residual diffraction is easy to identify when the velocity is far too low (e.g., 0.20m/ns), but it becomes more subtle because the correct velocity is approached. At 0.24m/ns, for example, the top of the diffractor is slightly curved downward, indicating a velocity that is slightly too low. The images from 0.245 to 0.255m/ns are very similar, and each could be considered an acceptably migrated image (though a slight upward curvature may be detectable at 0.255m/ns). At 0.26m/ns and higher, the upward curvature is obvious, indicating that these velocities are unacceptably high. Thus, in this example, we selected a migration velocity of 0.25±0.005m/ns for a depth corresponding to 13.2 ns two-way traveltime.

Numerous diffractions are seen in our GPR lines; to generate a statistically meaningful data set, we selected 100 of these and identified the optimum migration velocity for each one. Applying the Dix equation yielded snow layer velocity estimates ranging from 0.20 to 0.26m/ns, with an average velocity of 0.229±0.005m/ns. After removing five outliers with anomalously slow velocities that predicted unrealistically high snow densities of >500kg/m3, we applied the Tiuri equation to convert radar velocity to density at each location. We then determined a single snow density versus depth curve for our study area by compiling all snow densities and depths and fitting a second-order polynomial to the data (Figure 8). This yielded a curve of depth-averaged snow density as a function of total snowpack thickness, which shows increasing average density with depth from 1 to 4 m, which then levels off and decreases slightly to 6.5 m depth (this decrease is likely an artifact of the second-order polynomial and is not interpreted to be real). The average snow density of the entire data set is 360±70kg/m3 (where the error bound is one standard deviation; including the five outliers yields 370±90kg/m3). This estimate is essentially indistinguishable from the average snow density measured in the snow cores (360±60kg/m3), with slightly higher scatter (Figure 8).

Estimates of SWE

We applied the density-depth relationship calculated above between the tracked reflections of the snow surface and ground surface to estimate SWE over each transect. The SWE was calculated by tabulating the snow thickness at each location and multiplying by the average snow density to that depth from the relationship shown in Figure 8. To display these data in a more geographically realistic way, we applied a static shift to each depth-converted GPR trace based on an interpolated elevation profile from the GPS data along the transect (Figure 9). The result shows that SWE varies strongly along a single transect; in the case of line 25, SWE varies locally from approximately 0.2 m to nearly 2 m. The SWE variation is driven by changes in snow depth, which varies over a factor of 10 on this transect.

DISCUSSION

Errors

The scatter in snow density inferred from radar velocities in Figure 8 likely has several sources. First and foremost, the error of ±0.005m/ns in the determination of migration velocity corresponds to uncertainty in snow density of ±50kg/m3; this probably accounts for most of the scatter in the data. In addition, several other factors could come into play. Any out-of-plane diffractors (that is, objects lying within reach of the GPR antennas but not directly in the vertical plane of the radargram) would produce diffractions that would focus at the “correct” velocity (the velocity of the medium between the antennas and the diffractor) but would be plotted in the vertical section at an incorrect two-way traveltime. Applying MVA on sections uncorrected for topography could result in distorted hyperbolas that might focus at an incorrect velocity, though we expect this effect to be small because most diffractions extend at most 5–6 m horizontally, and elevation changes over such short distances are small.

The error bounds on the SWE plot (Figure 9) were calculated by applying the maximum and minimum density profiles implied by the standard deviation of ±70kg/m3 in the density estimates from Figure 8. This is likely a conservative estimate of the error bound because any error in snow density estimates from MVA will be counteracted by a compensating error in snow depth. For example, an underestimate of migration velocity of 0.005m/ns will produce an overestimate of snow density by 50kg/m3, but it will also underestimate the snow depth (because a slower EM velocity predicts a thinner snowpack). We consider a 4 m deep snowpack with an average density of 360kg/m3, and thus a SWE of 1.44 m and a radar velocity of 0.23m/ns. If we underestimate radar velocity to be 0.225m/ns, we would infer a snow density of 390kg/m3, an error of 8.3%. However, the velocity error would also underestimate snow depth as 3.9 m, so the estimate of SWE would be 3.9m×(390kg/m3)/(1000kg/m3) or 1.52 m—an error of only 5.5%. This trade-off between snow density and depth estimates will tend to stabilize estimates of SWE from GPR data despite small errors in inferred radar velocity.

Snow stratigraphy

The images produced by the snowmobile GPR system provide remarkable detail about snow stratigraphy that can illuminate processes of snow deposition and redistribution. Figure 10, for example, shows a zone of bright, disorganized reflections near the left edge of the figure that corresponds to a buried tree approximately 3 m high. Our interpretation that these bright reflections were due to a bank of trees was confirmed by examination of high-resolution surface photographs and by a return visit to the site. The depositional patterns of snow layers to the right of the tree indicate that the tree has acted as a snow fence, concentrating snow accumulation downwind of the tree (the dominant wind direction at this site is from the west, or from left to right in Figure 10). Tracking the thickness of successive packages of snow (labeled 1–4 in Figure 10) shows that the locus of snow deposition moves progressively downwind through the winter — the peak thickness of each successive depositional unit moves approximately 12 m downwind. The images show that the mean age of the snowpack changes systematically in the downwind direction behind the tree, even though the total snow depth does not change significantly; immediately adjacent to the tree, the snow consists almost entirely of the first-deposited units 1 and 2, whereas 30 m downwind, layer 1 is absent and the snowpack consists mostly of younger units 3 and 4.

Spatial distribution of SWE

The rapid acquisition of data using a snowmobile, combined with the ability to independently estimate snow density from migration analysis, opens up new possibilities for understanding the distribution of SWE on mountain landscapes.

Histograms of SWE from our data allow a more detailed look at the distribution of SWE in our study area (Figure 11). The overall pattern (Figure 11A) shows an average SWE of 0.61 m on our transects, with a highly skewed distribution, showing a majority of values at SWE <0.6m, but a long tail that extends out to nearly 2 m. The average SWE on the transects is slightly lower than the SWE measured at the Brooklyn Lake SNOTEL station, which reported SWE values of 0.68 m on 6 March and 0.76 m on 27 March (Figure 2). More importantly, the histograms point to systematic lateral variability at scales of approximately 1 km (the distance between lines) because the SWE distributions for the three lines are rather different. Line 26 shows a skewed distribution similar to the average histogram (because line 26 was the longest line and thus dominates the total data), but lines 25 and 34 show different patterns. Line 25, the highest elevation transect, shows a more Gaussian distribution of snow, with a higher mean SWE of 0.92 m. Line 34, the lowest elevation line, shows a bimodal distribution, with peak SWEs of approximately 0.3 and 0.8 m.

Limitations and future work

Because the migration focusing approach relies on the presence of serendipitously located diffractors (e.g., rocks, boulders, bushes, and logs), it provides snow density estimates at unpredictable spatial sampling. Although such objects are plentiful in our study area, the method cannot be applied where diffractors are absent, such as over very smooth ground. Moreover, because the locations of diffractors are unknown in advance, velocity sampling will be unpredictable and random. Despite these limitations, we suggest that the migration focusing analysis will be applicable to many data sets, especially in mountain watersheds, and our results show that it can provide reliable estimates of snow density and SWE without the need for time-consuming CMP acquisition (e.g., Harper and Bradford, 2003). Supplementing the radar velocity estimates with a few well-chosen snow cores adds confidence in the results as well as the ability to select specific sites for more detailed analysis. For studies in which large spatial coverage is important, the snowmobile-mounted GPR is an effective choice.

Given the present data distribution, it is not possible to draw firm conclusions about the causes underlying the contrasting mean snow depths or patterns of distribution (Figure 11). However, the ability to rapidly cover terrain with our system points to a productive line of research where SWE distribution can be measured and modeled at finer scales than ever before. Especially when combined with the details of topography and vegetation cover provided by modern LiDAR data, we will be able to investigate the controls on snow distribution by such environmental variables as aspect, slope steepness, vegetation, elevation, and microtopography.

At present, the principal inefficiency in the workflow is the manual inspection of numerous panels of migrated images to select the optimum focusing velocities. A priority for future work is to automate this analysis by introducing an objective measure of image focusing. This would have the dual advantages of greatly speeding the analysis and minimizing human error in selecting focusing velocities.

CONCLUSION

We have developed a new mode of acquisition of GPR data over snow by mounting GPR antennas to the front of a snowmobile. This system enables rapid acquisition of data on long transects on any terrain that can be safely navigated by a snowmobile. Mounting the antennas above the snowpack, in front of the snowmobile, enables recording of the snow surface reflection and ensures that measurements are taken over undisturbed snow.

Our results demonstrate the following points: (1) GPR data provide high-resolution images of snow stratigraphy, enabling analysis of snow deposition and redistribution processes. (2) The MVA of diffractions produces estimates of snow density that agree with those from snow core measurements and have similar errors. Velocity errors of ±0.005m/ns in the velocity analysis correspond to uncertainties of ±50kg/m3. (3) Mean SWE inferred from our data was 0.61 m, within error of the value measured at a nearby SNOTEL station. In detail, however, SWE varied spatially over a factor of approximately 10, from 0.2 to 2 m, principally due to a strong lateral variability in snow depth.

This work demonstrates the feasibility of quickly measuring SWE over long transects. The success of the MVA approach suggests that, in any GPR data set containing ample diffractions, snow density can be accurately estimated without the need to collect time-consuming CMP-style GPR data or snow core data. Because the spatial variability of SWE is driven more by variations in snow depth than in snow density, we suggest that long GPR transects should be useful in validating models of snow hydrologic processes in mountain watersheds.

ACKNOWLEDGMENTS

We thank E. Traver for field work support and M. Dogan and K. Hyde for helpful conversations. This publication was made possible by the Wyoming Experimental Program to Stimulate Competitive Research and by the National Science Foundation under award no. EPS-1208909.