Measuring snow water equivalent from common-offset GPR records through migration velocity analysis

Many mountainous regions depend on seasonal snowfall for their water resources. Current methods of predicting the availability of water resources rely on long-term relationships between stream discharge and snowpack monitoring at isolated locations, which are less reliable during abnormal snow years. Ground-penetrating radar (GPR) has been shown to be an effective tool for measuring snow water equivalent (SWE) because of the close relationship between snow density and radar velocity. However, the standard methods of measuring radar velocity can be time-consuming. Here we apply a migration focusing method originally developed for extracting velocity information from diffracted energy observed in zero-offset seismic sections to the problem of estimating radar velocities in seasonal snow from common-offset GPR data. Diffractions are isolated by planewave-destruction (PWD) filtering and the optimal migration velocity is chosen based on the varimax norm of the migrated image. We then use the radar velocity to estimate snow density, depth, and SWE. The GPR-derived SWE estimates are within 6 % of manual SWE measurements when the GPR antenna is coupled to the snow surface and 3–21 % of the manual measurements when the antenna is mounted on the front of a snowmobile ∼ 0.5 m above the snow surface.


Introduction
Many regions of the world are critically dependent on seasonal snowfall for their water resources; accurate estimates of how much water is stored in mountain landscapes are necessary to manage this resource.In the United States, a large network of SNOTEL sites provides continuous information about snow depth, density, and snow water equivalent (SWE) that is used to make water availability predictions (Serreze et al., 1999).While these sites provide valuable information at a site, scaling these point measurements up for basin-or grid-scale estimates can be challenging (Molotch and Bales, 2005).Currently, these data are used to develop empirical relationships between SWE and nearby stream discharge.These predictions are most accurate during average years and may not be reliable during abnormal years (Bales et al., 2006); thus, there is a need to develop new and reliable methods for estimating SWE on a basin scale.
Several previous studies have demonstrated that groundpenetrating radar (GPR) can be used to measure SWE (e.g., Bradford et al., 2009;Tiuri et al., 1984;Holbrook et al., 2016).Tiuri et al. (1984) showed that at microwave frequencies the real part of the dielectric constant for dry snow, which governs the velocity, is almost completely determined by the bulk density of snow.However, when liquid water is present, both the real and imaginary parts are needed to determine the volumetric water content of the snow.The complex dielectric constant can be measured by analyzing both the velocity and attenuation characteristics of the snow (Bradford et al., 2009).In the simplest case of dry snow, bulk density can be estimated directly from radar velocity.Snow depth can be measured from the two-way travel time of the radar pulse between the snow surface and the ground surface and the velocity.SWE can then be calculated as the product of snow density and snow height.
Velocity measurements can be made from the surface in several ways.Common-midpoint gathers (CMPs), where Published by Copernicus Publications on behalf of the European Geosciences Union.
the distance between transmitting and receiving antennas is steadily increased about a central location, provide highly accurate measurements; the two-way travel time to subsurface reflectors is a function of offset and velocity.Collecting CMPs requires separable antennas, and it can be timeconsuming to both collect and process these data.Commonoffset antennas, where both the transmitting and receiving antennas are housed in the same unit at a fixed offset, allow large amounts of data to be collected with minimal effort.Measuring the velocity from common-offset data can be achieved through calibration from measured snow depths, modeling diffraction hyperbola travel times, or migration focusing analysis.
In this paper, we apply the migration velocity analysis (MVA) presented by Fomel et al. (2007) to the problem of estimating radar velocities, and thus snow density and SWE, from 500 MHz common-offset GPR images.After testing the method on a synthetic data set, we estimate SWE from six field data sets.The first two data sets were collected by pulling the GPR along the snow surface, and the remaining four were collected with the GPR antenna mounted on the front of a snowmobile.To validate the method, we compare snow depth, density, and SWE estimates to measurements made in pits and probed depth observations along the profiles.Since our primary goal is to develop a method for quick velocity estimations, we assume that the snow we are measuring is dry.The GPR-derived estimates agree with manual SWE measurements within the estimated uncertainties.

Methods
GPR surveys utilize high-frequency, broadband electromagnetic signals.The signal is generated at the transmitting antenna and propagates in three dimensions at the velocity given by v = c/ √ κ , where c is the speed of light in a vacuum and κ is the real part of the dielectric constant.Signal attenuation is frequency dependent and can be approximated as α ≈ µ 0 κ κ 2 ω, where µ 0 is the magnetic permeability of free space and κ is the imaginary component of the dielectric constant (Bradford, 2007).While both κ and κ are frequency dependent, within the typical frequency range utilized for GPR studies, only κ exhibits strong variations with frequency; in dry snow κ ≈ 0 (Bradford et al., 2009).
When a GPR signal encounters a boundary between subsurface materials with contrasting dielectric constants, some of the energy is reflected back and recorded by a receiving antenna.In this paper, we are specifically interested in targets that have lateral dimensions that are less than the width of the Fresnel zone (the Fresnel radius is given by R f = zλ 2 , where z is depth and λ is the dominant wavelength).These objects scatter energy in all directions and appear on the raw GPR image as hyperbolic events, called diffractions (Landa and Keydar, 1998), whose shape depends on the depth of the object and the velocity of the overlying media (i.e., Claerbout, 1985).The velocity information contained in diffractions can be extracted by fitting hyperbolic curves to the data or by migrating the image until the hyperbola is collapsed to a point or "focus".The latter process is called migration velocity analysis.In this paper, we follow an approach described by Fomel et al. (2007) and develop a semiautomated MVA program in Matlab for the purpose of measuring radar velocities in seasonal snow.The processing flow consists of three steps: (1) separate diffractions from reflections through plane-wave destruction, (2) migrate the filtered images at a range of potential velocities, and (3) use the varimax norm as a measure of diffraction focusing to pick velocities.

GPR data
During February and March 2015, we collected GPR, snow density, and snow depth data in the Medicine Bow Mountains, SE Wyoming.The GPR data were acquired with a Malå pulse radar system with center frequencies of 500 MHz.The data were collected in two ways.In one configuration (lines 1 and 2), we mounted the GPR antenna in a plastic sled and pulled it behind a skier.The unit was set to fire continuously at a rate of 20 traces s −1 and the sample interval on each trace was 0.3223 ns.In the other configuration (lines 3, 4, 5, and 6) the antennas were mounted on an aluminum frame attached to the front of a Polaris RMK 600 snowmobile.The unit was set to fire at a rate of 100 traces s −1 and the sample interval was 0.3181 ns.Mounting the GPR antenna in front of the snowmobile allows us to measure undisturbed snow as well as provide a snow surface reflection, which can be used to analyze the attenuation properties of the snow (Bradford et al., 2009).In both cases, we kept track of our position with a Trimble R8 GPS unit that recorded our location at 1 s intervals.

Snow depth and density data
To validate our snow density and velocity estimates from the GPR data, we manually measured snow depth and densities (Table 1).On lines 1, 2, 3, and 4 we dug snow pits and located them with a handheld Trimble GPS unit.To measure snow densities, we used a 0.001 m 3 , wedge-shaped snow sampler and a scale that is accurate within 5-10 g.We made snow density measurements at 10 cm intervals in the sidewall of the snow pits starting from the snow surface and continuing to the ground.Pit locations were chosen based on the presence of diffractions near the snow-ground interface after viewing the GPR images in the field.On lines 4, 5, and 6 we measured snow depth at regular intervals with a probe.
Probed depth measurements are subject to uncertainties due to uneven ground and deviations in probe angle.We estimate our depth measurements to be accurate within ±5 cm.Snow density observations are subject to over and under sampling and we assign an uncertainty of ±5 g cm −3 .We calculate the average density for each pit profile assigning each snow density observation to a 10 (± 1) cm column of snow and performing a weighted sum.Propagating the uncertainties through the averaging process yields uncertainty estimates of 10-14 % of the averaged value, consistent with uncertainty estimates for snow pit density measurements reported by Conger and McClung (2009).

Preprocessing the GPR data
Prior to MVA we use MATGPR R3 (Tzanis, 2010) to apply several basic processing steps to the GPR data including (1) resetting the trace to time zero, (2) trimming the time window, (3) interpolating the traces to equal spacing using the GPS data, (4) bandpass filtering the data from 100 to 1000 MHz, (5) using a median filter to remove antenna ringing, and (6) scaling the amplitudes by t 2 .

Plane-wave destruction
Plane-wave destruction (PWD) is a predictive filtering method designed to suppress events in a seismic or GPR record having a particular dip (Claerbout, 1992;Fomel, 2002).The GPR image is modeled as the local superposition of plane waves described by the following differential equation (Fomel, 2002): where P (xt) is the wave field and σ (xt) is the local dip.Equation (1) provides the means for predicting a trace in the GPR image from its neighbor as a function of local dip.Fomel's (2002) three-point filter is derived from this equation: where σ is the local dip and the filtering is accomplished by convolving Eq. ( 2) with the GPR image.
Our goal is to suppress continuous reflections that have small dips (such as snow layering and the ground surface) compared to the steeply dipping diffraction limbs.To estimate local dips, we make an initial guess σ 0 for the dip and solve the set of equations for σ .Here, C(σ ) denotes the convolution of the filter with the data (d), C (σ ) is the derivative of the filter with respect to σ (C (σ )d is a diagonal matrix), D is the gradient operator, and ε is a weighting parameter that controls the smoothness of the estimated dip field.Imposing smoothness constraints on the dip field estimate ensures stability in the solution and helps target the reflections in the image, since they generally show higher amplitudes and are more laterally continuous than the diffractions we seek to preserve.The estimated dip field is then used to filter the data.

Migration
Migration is the process that moves reflected and diffracted energy in a seismic or GPR record to its true location in the subsurface (i.e., Claerbout, 1985).The quality of the migration process depends on the accuracy of the velocity estimate.When the correct migration velocity is chosen, diffraction hyperbolas will collapse to a compact focus.With too low a velocity, the hyperbola will only be partially collapsed, while a velocity that is too high will cause the hyperbola to be mapped into a "smile".
For the MVA analysis, we migrate the entire image through a suite of velocities (0.19-0.29 m ns −1 in increments of 0.002 m ns −1 ) using MATGPR's implementation of the Stolt algorithm (Stolt, 1955).The Stolt algorithm performs the migration in the frequency wave-number domain and is computationally efficient.To reduce computational time, we modified the code to perform all the migrations in one function call so that the forward Fourier transform is only performed once.

Velocity picking
After PWD filtering and migrating the data through the suite of velocities, the next task is to use a focusing indicator to pick the image that is optimally focused.Following Fomel et www.the-cryosphere.net/11/2997/2017/The Cryosphere, 11, 2997-3009, 2017 al. ( 2007), we use the varimax norm (V ): where s i is the amplitude of the ith sample and N is the number of samples included in the calculation.V is a measure of the "simplicity" of a signal (Wiggins, 1978).Since the simplest possible signal is a spike and the optimal migration velocity will map hyperbolas to the most compact focus, the maximum V value will correspond to the image migrated with the optimal velocity.To assess the possibility for errors in the migration velocity analysis, we applied our workflow to a synthetic data set generated from diffractors of varying size.We created five synthetic diffractions with a migration velocity of 0.24 m ns −1 .The first four (Fig. 1a) correspond to rectangular objects at 1 m depth with horizontal dimensions of 0.1, 0.2, 0.3, and 0.4 m and thickness of 0.03 m.The fifth corresponds to a circular object with a radius of 0.4 m (close to R f for the 500 MHz ricker wavelet used to generate the diffractions).The corresponding V curves for the windows shown in Fig. 1a are plotted in Fig. 1b.The V curves are peaked at 0.24 m ns −1 for all of the rectangular diffractors, with flatter (less well-resolved) peaks as the horizontal dimension of the diffracting object increases, suggesting a larger uncertainty in the velocity estimate.The peak V value for the circular diffractor is at 0.268 m ns −1 , indicating that curved objects with lateral dimensions close to the size of the Fresnel zone may continue to focus at velocities higher than their true ve-locity.Finally, Fig. 1c shows the V curve for the entire image, which peaked at the correct velocity of 0.24 m ns −1 .This analysis suggests that the peak V value will correspond to the correct velocity if the majority of the diffractions correspond to objects with a radius less than R f .
We choose to compute V in sliding windows that span the entire time section and have a user-defined width.Computing V in this way allows us to incorporate many diffraction events and maximize the likelihood that the bulk of the diffractions satisfy the point diffractor assumption.Moreover, sliding windows offer the potential to capture lateral variability in snow density.After computing V for the entire data set, we choose the maximum V value in each window to get an estimate of the migration velocity.Noise in the filtered image, large diffracting objects, or a lack of diffractions may cause the peak of the V curve to correspond to an incorrect velocity.To reduce the influence of erroneous velocity picks, we smooth the picks in the lateral direction with a boxcar averaging filter the same width as the sliding window.
We use the shape of the upper portion of the V curve to estimate uncertainties in the velocity pick.Comparing the V curves for synthetic diffractions as well as those from our data, we find that V values that are greater than 95 % of the peak value correspond to migrated images that are indistinguishable to the human eye (Fig. 2).We therefore obtain upper and lower bounds on our velocity estimate by finding the minimum and maximum velocities with V values equal to 95 % of the maximum.We use the upper and lower bounds on our velocity estimates to compute upper and lower bounds on all subsequent calculations.

Dix equation
The migration velocity is the RMS velocity of all of the material between the GPR antenna and the diffractor.When the GPR antenna is in contact with the snow and the diffractor is located at the base of the snow, we interpret the migration velocity to be the average velocity of the snow across the width of the diffraction hyperbola.When the GPR unit is mounted on the front of the snowmobile, the signal must pass through the air between the antenna and the snow surface so that the migration velocity is higher than that of the snow.To find the snow velocity from these data, we use the Dix equation (Dix, 1955): where velocity subscripts refer to the migration velocity, the velocity in air, and the velocity within the snowpack and time subscripts refer to the two-way travel times of the snow surface and soil surface reflections.
The Dix equation contains two important assumptions.First, the velocity of the snow must be approximately constant over the width of the hyperbola and, second, the halfwidth of the hyperbola should be small compared to the depth of the diffractor (x z).The diffractions in our data sets are approximately 4 to 5 m wide; thus, we assume that any lateral variations in snow density occur on a larger scale than this.If the second assumption is not valid, then the Dix velocity will be higher than the true velocity, resulting in a density estimate that is too low.The snow depths in our data range from ∼ 1-2 m, which is comparable to the half-width of the hyperbolas.
To determine the minimum snow depth that satisfies the x z assumption, we traced rays from point diffractors at depths ranging from 0 to 5 m through a 0.23 m ns −1 snowpack, representing a snow density of 0.358 g cm −3 (see Sect. 2.7), with a 0.5 m thick air layer between the snow surface and the receiver positions (Fig. 3).For each resulting travel-time curve, we obtained nine different estimates of the migration velocity by performing a least-squares fit to www.the-cryosphere.net/11/2997/2017/The Cryosphere, 11, 2997-3009, 2017 the travel-time data and successively reducing the widths of the hyperbolas from 10 to 2 m in 1 m increments.Using the Dix equation, we obtained estimates of the snow velocity as a function of diffractor depth and hyperbola width (Fig. 4).
The velocity estimates made with the Dix equation approach the true velocity as the diffractor depth increases and the hyperbola width decreases.For hyperbolas that are 4-5 m wide (the average width that we observe in our data), the Dix velocity is within 2 % of the true velocity when the diffractors are about 1.5 m deep, 5 % when the diffractors are about 1 m deep, and 10 % or greater when the diffractors are 0.5 m deep.
We conclude that the use of the Dix equation is justified for diffractors buried deeper than 1.5 m beneath the snow surface.
Although the results of this analysis are only valid for travel-time modeling, the x z assumption may be less severe for migration focusing analysis.Diffraction amplitudes decrease with increasing horizontal distance from the diffractor location; thus, the traces closest to the diffractor have the greatest contribution to the final image, suggesting that the Dix equation may give adequate results for diffractors that are less than 1.5 m deep when velocities are estimated from MVA (we test this with our synthetic data set in Sect.3.1).
To propagate our velocity uncertainty estimates through the Dix equation, we assign an uncertainty of 0.2 ns to our travel-time observations and use Eq. ( 5) along with our ve-locity uncertainty estimates to compute upper and lower bounds on the snow velocity.

Estimating SWE
To estimate SWE from the radar data, we need to know the depth of the snow and the snow density (SWE = z snow ρ snow ).The depth can be found by picking the two-way travel time of the ground reflection and, if applicable, the snow surface reflection and then using the velocity estimate to convert time to depth.Using Eq. ( 1), we convert radar velocity to the dielectric constant (v = c/ √ κ ) and estimate the density of dry snow with the following empirical relationship (Tiuri et al., 1984): where κ d is the dielectric constant and ρ is the density of dry snow.
In this paper, we are primarily concerned with measuring radar velocities and we assume that our data measure the properties of dry snow.The dry snow assumption can be tested from the data by analyzing the attenuation properties of the snowpack (Bradford et al., 2009).Because the real part of the dielectric constant for water (∼ 80) is much larger than that of snow (∼ 1.5-2) and the imaginary part, which describes the attenuation of the signal, is non-negligible (Bradford at al., 2009), the radar signal will travel at a slower velocity and attenuate more rapidly when liquid water is present in the snow.The attenuation coefficient for radar waves in water is frequency dependent (i.e., Turner and Siggins, 1994), with the higher frequencies attenuating more rapidly that the lower frequencies because they go through more cycles per distance traveled.
To test the dry snow assumption, we compare the frequency content of the ground reflection to a reference event (the direct arrival for ski-pulled data and the snow surface reflection for the snowmobile-collected data).We calculate the maximum local instantaneous frequency (Fomel, 2007) within a time window surrounding the event of interest and then average this value across all of the traces in the GPR image.The standard deviation provides an estimate of the measurement uncertainty.We note that at 500 MHz a small shift in frequencies results in a non-negligible volumetric water content.

Data and results
Snow depth, density, and SWE estimates for all of our GPR profiles and pits are summarized in Tables 1 and 2.Here we discuss the processing and describe results for a synthetic data set and two representative field data sets.  1 Lines 2, 3, 5, and 6 are described in the Supplement. 2 Line 3 was located 1.5 m off of Pit 2, disagreement between depth and SWE measurements at this site reflect lateral variations in snow depth. 3RMSE percentages are calculated relative to the mean observed depth along each profile.

Synthetic test
As a first test on the reliability of migration focusing analysis for reconstructing radar velocities, we performed the analysis on a synthetic data set generated with REFLEX software.The synthetic data set was generated using a 500 MHz Kuepper wavelet sampled at 0.0332 ns and traces are 0.01 m apart.
The synthetic model is 50 m long and consists of a 0.5 m thick layer of air overlying a 0.24 m ns −1 layer of snow (corresponding to a density of 0.29 kg cm −3 ) with depths that range from 0.5 to 5.7 m.Beneath the snow is a 0.10 m ns −1 layer representative of soil.Along the snow-soil interface there are 16 diffractors buried at depths ranging from 0.5 to 5.7 m.The purpose of this data set (Fig. 5a) was to test the performance of the Dix equation on velocities estimated from the MVA analysis and, since the migration velocity changes as a function of snow depth, to see if we can resolve lateral variations in velocity.
After applying the PWD filter, the ground reflection was adequately suppressed (Fig. 5b).We migrated the filtered image at 0.002 m ns −1 intervals from 0.18 to 0.28 m ns −1 and measured the optimal migration velocity for each diffractor by computing V in an 8 m wide sliding window (Fig. 5c).We use the Dix equation to convert the migration velocities www.the-cryosphere.net/11/2997/2017/The Cryosphere, 11, 2997-3009, 2017 to the velocity of the snow layer (Fig. 5d).The average of all snow velocity measurements is 0.241 m ns −1 with a standard deviation of 0.002 m ns −1 .
There is no systematic relationship between the velocities recovered and the depth of the diffractor (Fig. 5d).The shallowest diffractor was at ∼ 0.5 m depth and the recovered velocity was 0.232 m ns −1 .The greatest differences between recovered and true velocities were for diffractors at depths of 0.5, 1.6, 2.2, and 2.3 m.Here the recovered velocities were 0.232, 0.247, 0.247, and 0.248 m ns −1 .The shallowest observation underestimates the true velocity, which is the opposite of the effect predicted by our travel-time modeling (Sect.2.6, Fig. 4).The observations for diffractors between 1.5 and 2.3 m all overestimate the true velocity by approximately the same amount.We conclude that the Dix equation is appropriate for snow depths of 0.5 m and greater.
Although the snow in this synthetic model has a constant velocity, the migration velocity changes as a function of the snow depth due to the changing proportions of air and snow in the total travel path.Where the snow is shallow, the velocities are highest and, where the snow is deep, the velocities are low.That the method is capable of resolving lateral velocity variations in this synthetic example is evident in Fig. 5c, where the picked velocities are negatively correlated with snow depth.

Ski-pulled GPR data
We collected two GPR profiles in the ski-pulled configuration on 25 February 2015, in below-freezing conditions.A representative line, Line 1 (Fig. 6), is 74 m long and shows an abundance of diffractions along the snow-ground interface, likely a result of small boulders, and a few isolated diffractions within the snowpack, most likely small trees, bushes, or logs.After interpolation to equal spacing, trace spacing was 0.0362 m.Since the antenna was coupled to the snow, we compare the average frequency of the direct wave to that of the soil reflection to determine whether there is any liquid water present in the snowpack.The average frequency of the direct arrival for every trace in the image along Line 1 is 410 MHz with a standard deviation of 10 MHz, and the average frequency of the soil reflection across the whole line is 457 MHz with a standard deviation of 42 MHz.The soil reflection appears to have a higher-frequency content than the reference frequency, perhaps due to thin-layer "tuning" effects.Since we do not observe a decrease in frequency with travel time, we infer that there was no liquid water present in the snow on this day.
After the PWD filtering step we are left with many diffractions along the ground surface and a few isolated events within the snowpack (Fig. 6b).We compute V in 10 m wide sliding windows and pick the velocity that corresponds to the peak value of V (Fig. 5d, blue line).After smoothing these picks (Fig. 6d, red line) we obtain velocities between 0.237 and 0.276 m ns −1 , with an average uncer-

Snowmobile-mounted GPR data
We collected four GPR profiles in the snowmobile-mounted configuration between 25 February and 17 March 2015.
Here we discuss the processing of a representative profile, Line 4 (Fig. 8), which was collected on the morning of 11 March 2015 in a flat meadow just south of Wyoming Highway 130.This line is 98 m long and shows an abundance of diffractions along the snow-ground interface (Fig. 8).After interpolating to equal spacing, the trace spacing was 0.024 m.
Migration velocities on this line range from 0.237 to 0.277 m ns −1 with an average uncertainty of ±0.01 m ns −1 .The corresponding snow velocities are 0.207 and 0.268 m ns −1 .Here, the exceptionally high velocities are confined to a region between x = 65 and x = 85 m where a number of diffractions from obviously large objects are present (Fig. 8).If we exclude velocity picks from this region, we get a maximum migration velocity of 0.266 m ns −1 and a maximum snow velocity of 0.251 m ns −1 .Estimated snow depths range from 0.7 to 2.1 m with an average uncertainty of ±0.1 m.Estimated snow densities range from 228 to 532 kg m −3 with an average uncertainty of ±50 kg m −3 .Estimated SWE ranges from 0.25 to 0.71 m with an average uncertainty of ±0.09 m.
Snow pits 3 and 4 are located at 50 and 97 m along the profile and showed average snow densities of 379 ± 50 and 360 ± 48 kg m −3 , SWE values of 0.54 ± 0.13 and 0.64 ± 0.13 m, and snow depth values of 1.44 ± 0.05 and 1.8 ± 0.05 m, respectively.The GPR-derived depth, density, and SWE estimates at 50 and 97 m were 1.50 ± 0.08 and 1.91 ± 0.12 m, 389 ± 92 and 394 ± 97 kg m −3 , and 0.53 ± 0.09 and 0.70 ± 0.13 m.GPR-derived estimates for the whole profile are shown in Fig. 9.We also measured 21 snow depths at 5 m intervals (Fig. 9b) along this profile.The RMS error between observed and estimated depths is 0.13 m.
During data acquisition on Line 7, the air temperature was 5 • C, raising the possibility of liquid water in the snow.The average frequency of the snow reflection for every trace in the image is 435 MHz with a standard deviation of 27 MHz, and the average frequency of the soil reflection across the whole line is 464 MHz with a standard deviation of 38 MHz.Again, the frequency content of the soil reflection appears to be higher than the reference frequency.Within the uncertainty bounds there is no resolvable frequency change and we conclude that our dry snow assumption is valid.

Discussion
The primary purpose of this study is to develop an efficient processing flow for measuring GPR velocity and thus snow density SWE from common-offset data that requires a minimum amount of human interpretation.Common-offset GPR data are fast and easy to obtain, and velocity estimates can be made when diffractions are present.However, the common methods of visually inspecting migrated images or fitting curves to diffraction hyperbolas are time-consuming and subject to human error.The migration velocity analysis described in this paper provides an efficient means for extracting velocity information from large GPR data sets.Here we discuss the accuracy and efficiency of the method as well as the level of automation.
To validate the method, we compared estimated snow densities, depths, and SWE to observations made in four snow pits and to 86 probed snow depth measurements.The results are summarized in Table 2 and in Fig. 10.If we exclude the two obvious outliers (Fig. 10a), the RMS error for our depth predictions for the remaining 88 depth observations is 12 % of the mean snow depth observation.The RMS errors for snow density and SWE relative to the mean observed values are 15 and 18 %.Averaging the velocities across the entire line (Fig. 10, red crosses) reduces the difference between predicted and observed depth values to an RMS error of 9 %, suggesting that lateral variations in snow velocity are minimal.Averaging the velocities across the entire line reduces the RMS errors for density and SWE to 8 and 10 %, respectively.
The greatest potential for systematic error in this analysis is the presence of diffracting objects whose dimensions exceed the radius of the first Fresnel zone.The field data offer the opportunity to evaluate the influence of diffractor www.the-cryosphere.net/11/2997/2017/The Cryosphere, 11, 2997-3009, 2017 size on velocity estimates.Line 1, for example, shows four prominent diffractions between 50 and 70 m.The varimax norm has a maximum value at 0.256 m ns −1 , which is the velocity that focuses the two leftmost diffractions (Fig. 6c).The diffractions on the right are clearly not focused because they are caused by an object (most likely a log) with a radius greater than the first Fresnel zone.Because the leftmost two have a higher amplitude than the others, they have the largest influence on the varimax value.Thus, although there are clearly events in the field data that have the potential to give erroneous results, our results suggest that reliable velocity estimates can be achieved so long as the majority of the diffracted energy is related to objects that can be considered point diffractors.One of our main goals was to produce a processing flow that allows for the rapid processing of common-offset GPR data with minimal user interaction.The two most time computationally expensive parts of the processes are the migrations and the varimax calculations.As an example, on a 2016 MacBook Pro with a 2 GHz processor, for the ∼ 100 m long Line 4, performing 51 migrations takes approximately 5 min, the varimax calculation takes about half as long, and the PWD filtering takes a few seconds.The most timeconsuming part of the process is picking the arrival times of snow surface and ground surface reflections.
Although the processing flow is relatively efficient, it does require some user interaction.The PWD method of separating continuous reflectors from diffractions treats the GPR image as the superposition of locally planar waves.Estimating the slope of these waves from the image requires the solution of a regularized inverse problem and the smoothness of the slope-field depends on the choice of regularization parameter.This is the most subjective step of the process, as it may require several attempts to find the optimal smoothness constraints to adequately suppress reflections in the GPR image.However, for our data the majority of the diffractions are located along the ground surface and the internal structure of the snowpack shows dips that closely parallel the ground reflection.A good first guess, and often a good final guess, for the dip field can be computed by picking the arrival times of the ground reflection.Because the ground reflection has to be interpreted to measure snow depth, this strategy can significantly reduce the processing time for each data set.
The data presented in this paper contained an abundance of diffractions located near the soil-ground interface allowing an average velocity for the entire snowpack to be obtained.These events are likely due to small-scale variations in surface topography, rocks, and/or vegetation along the ground surface, which may not be present in all environments.However, we note that mountain watersheds free of vegetation, small undulations in surface topography, and surface rocks are probably rare.Thus, the method may be useful in many regions where a seasonal snowpack exists.

Conclusions
We applied the migration focusing analysis presented in Fomel et al. (2007) to the problem of estimating SWE in seasonal snow.The method was most accurate for the case when the GPR was in contact with the snow, providing GPRderived SWE estimates within 6 % of the manual observation.When the GPR was mounted on a snowmobile, the results were within 12-21 % of the manual observations.

Figure 1 .
Figure 1.(a) Synthetic hyperbolas for four rectangular diffractors with lateral dimensions of 0.1, 0.2, 0.3, and 0.4 m (from left to right) and a round diffractor with radius = 0.4 m (far right; TT2W is the two-way travel time).(b) V curve for windows depicted in (a); V curve colors match the windows in (a).V curves for all four rectangular diffractors show peaks at 0.24 m ns −1 , while the round diffractor is peaked at 0.268 m ns −1 .(c) Varimax curve for the entire image showing a peak at the correct migration velocity of v = 0.24 m ns −1 .

Figure 2 .
Figure 2. Justification for uncertainty estimates.A synthetic hyperbola that is obviously under-migrated (a), migrated at indistinguishable velocities (b-d), and obviously over-migrated (e).(f) The corresponding varimax curve for (a)-(e) showing a peak at the true migration velocity (0.24 m ns −1 ), the shaded area under the curve corresponds to velocities in (b)-(d) and represents varimax values that are 95 % of the maximum.Panels (g)-(l) show the same for a section of field data extracted from Line 1.

Figure 3 .
Figure 3. Ray paths and travel times for point diffractors.(a) A value of 0.5 m of air overlying a 230 m ns −1 snowpack with point diffractors buried at 0.5 m intervals.(b) Two-way travel times for each of the diffractors showing the characteristic hyperbolic shape.

Figure 4 .
Figure 4. Dix velocities for point diffractors as a function of depth for different hyperbola widths.The true interval velocity is 0.230 m ns −1 (red line) and the Dix velocities are shown as black lines.The red dashed line is at 0.234 m ns −1 , which is 2 % greater than the true velocity.

Figure 5 .
Figure 5. Synthetic data set and velocity picking.(a) Synthetic data before filtering.(b) After PWD filtering.(c) Varimax norm for sliding window 8 m wide.(d) Velocities from synthetic data set as a function of diffractor depth.The solid blue line shows measured migration velocities; dashed blue lines show uncertainty bounds.The solid red line shows velocities computed with the Dix equation; dashed red lines show uncertainty bounds.The solid black line shows the true velocity (0.24 m ns −1 ).The light gray region indicates where velocities are within 2 % of the true velocity, and the dark gray region shows where velocities are with 5 % of the true velocity.

Figure 6 .
Figure 6.(a) Raw GPR data for Line 1, (b) GPR data after PWD filtering, (c) diffractions migrated at the mean velocity (0.245 m ns −1 ) for the entire line, and (d) normalized varimax curves for a 10 m wide sliding window.The blue curve shows the peak value for every curve; the red line is smoothed with a 10 m wide box car averaging filter.

Figure 7 .
Figure 7. Line 1 results: (a) density, (b) snow depth (black line), and SWE (blue line) estimates from the GPR data; snow pit data are shown in red.Grayed out region corresponds to areas where velocity picks are unreliable.

Figure 9 .
Figure 9. Line 4 results: (a) density; (b) snow depth (black line) and SWE (blue line) estimates from the GPR data; snow pit data are shown in red.Grayed out region corresponds to areas where velocity picks are unreliable.

Figure 10 .
Figure10.Cross plots of predicted data (horizontal axis) vs. GPR estimates (vertical axis) for all data: (a) snow depths, (b) density, and (c) SWE.Black crosses represent estimates using automatically picked velocities and red crosses represent estimates using the mean velocity for each GPR profile.

Table 2 .
Summary of GPR field data and comparison to manual measurements.