Freely available online through the SEG open-access option.

ABSTRACT

Snow stratigraphy and liquid water content are key contributing factors to avalanche formation. Upward-looking ground-penetrating radar (upGPR) systems allow nondestructive monitoring of the snowpack, but deriving density and liquid water content profiles is not yet possible based on the direct analysis of the reflection response. We have investigated the feasibility of deducing these quantities using full-waveform inversion (FWI) techniques applied to upGPR data. For that purpose, we have developed a frequency-domain FWI algorithm in which we additionally took advantage of time-domain features such as the arrival times of reflected waves. Our results indicated that FWI applied to upGPR data is generally feasible. More specifically, we could show that in the case of a dry snowpack, it is possible to derive snow densities and layer thicknesses if sufficient a priori information is available. In case of a wet snowpack, in which it also needs to be inverted for the liquid water content, the algorithm might fail, even if sufficient a priori information is available, particularly in the presence of realistic noise. Finally, we have investigated the capability of FWI to resolve thin layers that play a key role in snow stability evaluation. Our simulations indicate that layers with thicknesses well below the GPR wavelengths can be identified, but in the presence of significant liquid water, the thin-layer properties may be prone to inaccuracies. These results are encouraging and motivate applications to field data, but significant issues remain to be resolved, such as the determination of the generally unknown upGPR source function and identifying the optimal number of layers in the inversion models. Furthermore, a relatively high level of prior knowledge is required to let the algorithm converge. However, we feel these are not insurmountable and the new technology has significant potential to improve field data analysis.

INTRODUCTION

Snow stratigraphy, i.e., the thickness and properties of the layers that form the seasonal snow cover, is one of the principal contributing factors considered in avalanche formation (Schweizer et al., 2003). In particular, layering in the form of cohesive slab layers above a weak layer is a prerequisite for dry-snow slab avalanche release. Until recently, snow stratigraphy was determined with manual observations. Conventional methods are destructive. Possible options include snow profile recordings by digging snow pits (Fierz et al., 2009) and measuring the penetration resistance with a snow micropenetrometer (Schneebeli and Johnson, 1998) to derive essential layer properties and snowpack stability (Proksch et al., 2015; Reuter et al., 2015). Furthermore, the layering on a pit wall can be characterized with near-infrared photography (Matzl and Schneebeli, 2006), and the liquid water content of the snow layers can be derived with instruments that measure the permittivity of the wet snow, such as the Denoth capacity probe (Denoth et al., 1984; Denoth, 1994) or the Finnish snow fork (Sihvola and Tiuri, 1986). All these measurements are not only cumbersome and time consuming but can also be performed only at accessible locations. Gathering information from avalanche starting zones is therefore not always possible.

Time-domain reflectometry allows noninvasive continuous monitoring of density and liquid water content (Schneebeli et al., 1998; Waldner et al., 2001). However, operating these sensors in avalanche starting zones is not possible because the installation requires poles that are prone to be destroyed by snow creep or avalanches. If meteorological input data are available, the snow stratigraphy can be modeled (Bartelt and Lehning, 2002; Lehning et al., 2002a, 2002b). The meteorological data can even be extrapolated for any location (Durand et al., 1999; Lehning et al., 2006), but modeling snow stratigraphy on slopes is not straightforward, and the validation data are largely lacking.

In recent years, upward-looking ground-penetrating radar (upGPR) systems have been developed that allow the snowpack to be characterized in a nondestructive manner (Figure 1). This technique has been applied successfully for monitoring the temporal evolution of snow stratigraphy (Gubler and Hiller, 1984; Gubler and Weilenmann, 1987; Heilig et al., 2009, 2010; Mitterer et al., 2011; Okorn et al., 2014; Schmid et al., 2014). Because the radar antennas are buried in the ground beneath the snow, they are able to withstand avalanches without damage. Schmid et al. (2014) demonstrate that snow height, settlement of old snow surfaces, new snow height, and the dry-to-wet transition could be derived from upGPR data if the density could be reasonably estimated. They also derived bulk snow density, bulk liquid water content, and snow water equivalent using external snow height measurements, which is not even needed if, in addition, the attenuation of concurrently recorded global positioning system signals is considered (Koch et al., 2014; Schmid et al., 2015).

A limiting factor of upGPR data analysis performed so far is the restriction to two-way traveltimes determined from the reflected waves. Exploiting the full information content offered by the upGPR data may provide more insights on the internal snowpack structures. This could be achieved using full-waveform inversion (FWI) techniques. This technology was proposed more than 30 years ago (e.g., Tarantola, 1984; Pratt, 1999), but it has only recently received wider attention and has been applied to a broad range of seismic (data sets) and a few GPR data sets (e.g., Ernst et al., 2007b; Virieux and Operto, 2009). Lambot et al. (2002) use a global multilevel coordinate search optimization algorithm combined sequentially with the classic Nelder-Mead simplex algorithm (Nelder and Mead, 1965) to determine the hydraulic properties of partially saturated soils. A local optimization algorithm is also used by Lavoué et al. (2014) in 2D frequency-domain FWI for the simultaneous reconstruction of dielectric permittivity and electric conductivity. Busch et al. (2014) adapt an FWI scheme for on-ground GPR. In addition to the soil electric conductivity and permittivity, this approach required the estimation of the amplitude and phase of the source wavelet. Furthermore, FWI has been successfully applied to estimate permittivity and conductivity in stratified structures such as concrete (Kalogeropoulos et al., 2011, 2013) or layered soils (Busch et al., 2012) and used in 2D crosshole sections (Ernst et al., 2007a; Klotzsche et al., 2010, 2012). It was even possible to characterize high-permeability zones in gravel aquifers with high resolution and in 3D (Klotzsche et al., 2013). Recently, the first applications to cryosphere problems have started to emerge (Babcock and Bradford, 2014, 2015).

Our aim is to explore possibilities and limitations of FWI techniques applied to upGPR data. After a brief introduction to the theory, we present numerical simulations of realistic upGPR scenarios. In particular, we investigate effects of liquid water content in the snowpack, the amount of a priori information required for successful FWI, nonlinear effects caused by strong contrasts of the physical snowpack properties, and the ability to detect thin layers. In this study, we consider only synthetic data, but we have chosen realistic experimental setups as they would be present during field experiments.

THEORY

Model parameterization

The snowpack located directly above an upGPR device often approximates a 1D layered medium, and this geometry allows a 1D formulation in our approach. Our model parameterization is shown in Figure 2. It includes a stack of M layers. The first layer represents the air space between the GPR antennas and a wooden board located on the ground that protects the antenna, and layer 2 represents the wooden board itself. This experimental setup corresponds to that used by Schmid et al. (2014). Layers 3 to M1 comprise the snowpack, and the uppermost layer M is an infinite half-space representing the air above the snowpack. The thicknesses and physical properties of the lowermost two layers and the overlying half-space are assumed to be known (Table 1).

Forward modeling

Electromagnetic wave propagation is governed by Maxwell’s equations, which can be written in source-free regions in the frequency domain in the form of the Helmholtz equations as follows (Ward and Hohmann, 1987):  
2E+k2E=02H+k2H=0,
(1)
where E and H are the frequency-domain representations of the electric and magnetic fields, respectively, and k is the electromagnetic wavenumber, defined as follows:  
k=μ0ε0εrω2iμ0σω¯,
(2)
where μ0 is the magnetic permeability of free space (4π×107VsA1m1), ε0 is the permittivity of free space (8.85×1012A2s4kg1m3), εr is the relative dielectric permittivity, σ is the electric conductivity (106Sm1 [Evans, 1965]), ω is the angular frequency (2πf, f=frequency), and i is the imaginary unit. The real part of k characterizes the electric displacement currents, and the imaginary part represents the conduction currents. For typical upGPR frequencies f of approximately 1 GHz and average snow properties εr2 and σ106Sm1 (Evans, 1965), the ratio of displacement to conduction currents is very large (17,700), i.e., the displacement currents dominate.

Wave propagation in layered media can be expressed with recursive relationships, as described, e.g., by Cory and Zach (2004). For the sake of simplicity, we make a plane-wave assumption. As discussed by Babcock and Bradford (2015), such an approximation is justified at recording distances that are larger than a few wavelengths.

The electromagnetic properties of each layer L (L=1,,M) are contained in wavenumber kL (equation 2). Therefore, the reflection response from a layered medium can be fully described by kL and the layer thicknesses DL. Following Cory and Zach (2004), the overall response can be written by means of a recursive relationship, which is based on the reflection coefficients defined as follows:  
rL1,L=kL1kLkL1+kL.
(3)
Because no downward traveling waves can occur within the uppermost air layer, the algorithm is initialized by setting the reflection response ΓL at the snow surface equal to the corresponding reflection coefficient; i.e.,  
ΓM=rM1,M.
(4)
Then, the reflection response at all layer boundaries is computed using  
ΓL=rL1,L+ΓL+1ei2kLDL1+rL1,LΓL+1ei2kLDL,(L=M1,,2).
(5)
Finally, the reflection response at the receiver antenna is obtained via  
Γ1=Γ2ei2k1D1,
(6)
where Γ1 is the Green’s function for a layered medium as shown in Figure 2. It needs to be multiplied with a realistic source function A(ω) to obtain the frequency-domain radargram W(ω) as follows:  
W(ω)=A(ω)×Γ1(ω).
(7)
A common choice for the source function A(ω) is the Ricker wavelet, defined in the frequency domain as follows:  
A(ω)=ω2ω0eω2ω02,
(8)
where ω0 is the center angular frequency (2πf0). According to the experimental setup used in Schmid et al. (2014), we considered a central frequency f0 of 1.6 GHz. The wavelet was normalized such that the amplitude at f0 was two. The corresponding frequency spectrum is shown in Figure 3. Throughout this paper, we assume A(ω) to be known.
The actually observed time-domain radargram w(t) can be computed by application of an inverse Fourier transformation:  
w(t)=F1{W(ω)}.
(9)
Because we are ultimately not interested in the electromagnetic properties of snow, but in its density ρ and liquid water content θw, we use the constitutive relationship introduced by Looyenga (1965) to obtain these quantities:  
εL=(θw,Lεwβ+θi,Lεiβ+(1θw,Lθi,L)εaβ)1β,
(10)
where εw, εi, and εa are the dielectric permittivities of water, ice (3.18), and air (1.0), respectively, and β=0.5 (Roth et al., 1990); θw,L is the Lth layer’s liquid water content, and θi,L is the ice content for a dry snowpack and is calculated from the dry-snow density ρL and the density of ice (ρi=917kgm3) as follows:  
θi,L=ρLρi.
(11)
The εw is a complex quantity. According to England et al. (1992), Behari (2006), and Hasted (1972), it can be obtained via  
εw=n2+εsε(1+iωτ1)1α+εn21+iωτ2n2=1.8εs=295.681.2283T+2.094×103T21.41×106T3τ1=5.62×1015e3.01×1020/(kbT)sτ2=4.2×1014s,
(12)
where n is the refractive index of water, εs is the static dielectric constant, ε=4.2 is the infinite-frequency dielectric constant, τ1 and τ2 are the relaxation times (the times required for the water molecule to align itself with an applied field), α=0.012 is the spread of relaxation or Cole–Cole spread factor, T is the absolute temperature in K (we assume a constant snow temperature of 273.15 K, i.e., at the melting point for ambient atmospheric pressures), and kb is the Boltzmann constant (1.3806×1023J/K).

In the presence of liquid water (i.e., θw,L>0), εL becomes complex, which introduces a dispersive (frequency dependent) GPR wave behavior (Ward and Hohmann, 1987).

Inversion

As discussed earlier, the snowpack model can be parameterized either in terms of the real and imaginary parts of εL and the layer thicknesses DL or by using densities ρL, liquid water contents θw,L and DL. We tested both options and found that the results were essentially the same. In the remainder of the paper, we will consider a parameterization in terms of ρL, θw,L, and DL. Formally, we merge these quantities into a model vector m of dimension 3M×1 as follows:  
m=(ρLθw,LDL)(L=1,,M).
(13)
The strong nonlinearity of the FWI problems often results in cycle skipping, which causes the inversion algorithm to be trapped in a local minimum. A possible option to mitigate this problem is to initially consider only low-frequency data. During the inversion, higher frequencies are successively added. This strategy favors FWI to be performed in the frequency domain (Pratt, 1999). Therefore, we consider the frequency-domain response W(ω) to form our data space d:  
dK=W(ωK)(K=1,,N),
(14)
where N is the number of data points used.
FWI aims at minimizing the discrepancy r between the observed data dobs and the predicted data dpred, which are computed using a set of estimated model parameters mest as follows:  
r=K=1N(dKobsdKpred)2.
(15)

Compared with typical geophysical FWI problems, which include a large number of model parameters, our 1D problem can be typically casted with only a few to a few tens of parameters. For the solution of the inverse problem, we used the Nelder-Mead algorithm as described in Lagarias et al. (1998) and implemented in the MATLAB function fminsearch. A big advantage of this type of algorithm is the straightforward consideration of additional constraints. For example, we have restricted the liquid water content θw,L to lie within a physically reasonable range (i.e., θw,L<10%), or we can impose that no liquid water is present (i.e., θw,L=0). Furthermore, the upper limit for densities ρL was set to 917kgm3. Finally, the model parameters of the air layer above the snowpack were fixed to the true values (Table 1). The inversion algorithm stops when the change in r and the norm of the model m between two iterations are less than 1×104.

A crucial step in any FWI scheme includes the determination of the initial model parameters with which the inversions are started. The closer the initial model parameters minit are to the true parameters mtrue, the higher the likelihood that the algorithm will converge to the true solution. If the snow properties are monitored continuously with upGPR, then it often can be assumed that the actual properties are not deviating too much from the previously determined parameters, and it makes sense to use mest from the previous experiment as minit for the actual experiment. Therefore, some prior knowledge can be assumed. In fact, we exploit this assumption in our inversion algorithm by supplying damping constraints so that deviations from the prior model are only allowed when positively dictated by the data.

We quantify our state of prior knowledge with the parameter p that ranges from one (no prior knowledge) to four. For all simulations, the initial model parameters for density and liquid water content were randomly selected from a uniform distribution (p=1) or a Gaussian distribution around the true model with standard deviations of σρ,p=100, 50, and 20kgm3 for the density (for p=2, 3, and 4, respectively) and σθw,p=2%, 1%, and 0.5% for the liquid water content (Table 2) for p>1:  
ρLinit=N(ρLtrue,σρ,p2)andθw,Linit=N(θw,Ltrue,σθw,p2),
(16)
where N(μ,σ2) is the normal Gaussian distribution with mean μ and variance σ2. Again, we restricted ρLinit to lie within a physically reasonable range of 50 and 917kgm3 and θw,Linit to lie within 0% and 10%.
For obtaining initial estimates for the layer thicknesses DL, we benefited from the two-way traveltimes of the reflected phases that could be picked on the observed time-domain radargrams wtrue(t). In this study, we used perfect values for the reflection traveltimes TL, calculated for each layer with the forward model. Suppose that we have a set of reflection traveltimes TL, an initial thickness estimate can be obtained using  
DLinit=c0(TLTL1)εL,
(17)
where c0 is the speed of light in a vacuum and εL is calculated using equation 10, where εw is assumed to be a constant value of 87.9 (Mitterer et al., 2011).

We used a combined global search and local minimization algorithm (Kalogeropoulos et al., 2011; Busch et al., 2012). For each simulation, 10 initial models mqinit were generated, and each initial model was refined by minimizing the discrepancy between the observed data dobs and the predicted data dpred (equation 15). If the different solutions were not identical, the solution with the smallest discrepancy r was chosen, thereby assuming that it represents the global minimum.

SIMULATION EXPERIMENTS

Experimental setup

To investigate possibilities and limitations of FWI applied to upGPR data, we designed four suites of simulations to address the following issues:

  1. A priori knowledge requirements: How much a priori knowledge is required to derive internal snowpack structures with a realistic complexity?

  2. Noise and imperfect picks: What is the influence of noisy data and imperfect picks on the simulations results?

  3. Influence of heterogeneity: How do the magnitudes of the contrasts in density and water content between two adjacent layers influence the success or failure of the inversions?

  4. Minimal-layer thickness: What is the minimal resolvable thickness of a layer?

All our simulations followed a workflow as described below:

  1. Calculate Wtrue(ω) and wtrue(t) from a prescribed set of model parameters mtrue using equations 3–11.

  2. Calculate TL from wtrue.

  3. Set fmax to 200 MHz.

  4. Calculate the initial density and liquid water content profiles using equation 16.

  5. Calculate the thicknesses using equation 17.

  6. Perform an inversion using W(ω) for ω<2πfmax as dobs (i.e., all frequencies up to fmax are used in the inversion).

  7. Increase fmax by 200 MHz.

  8. Repeat steps five to seven until fmax=4GHz.

Experiment 1: A priori knowledge requirements

We assumed 30 true synthetic snowpack models mtrue each consisting of six snow layers. The thickness of each snow layer was 0.3 m, and the density of each layer was randomly chosen between 50 and 900kgm3. Furthermore, we considered dry snowpacks, in which the liquid water content of the true models was zero and was kept fixed during the inversions, and wet snowpacks, in which the liquid water content was allowed to vary between 0% and 10%. The algorithm inverted for 16 parameters for dry snowpacks (density and height for six snow layers, the wooden board, and the air above the radar) and for 24 parameters for wet snowpacks (density, liquid water content, and height for six snow layers, the wooden board, and the air above the radar). Initial models were computed by means of equations 16 and 17, varying the amount of prior knowledge expressed by the parameter p between one and four.

Experiment 2: Noise and imperfect picks

Experiment 2 was the same as experiment 1, but white Gaussian noise with a signal-to-noise ratio of 18 dB was added to the observed dobs. In addition, instead of using perfect picks, we perturbed TL using a normal Gaussian distribution with mean TL and standard deviation 0.1 ns, which corresponds to one-sixth of a wavelength. This ensures that 99.7% of the picks are within half a wavelength away from the perfect pick. Such picking accuracy is judged to be realistic for the observed data.

Experiment 3: Influence of heterogeneity

Here, we assumed a true synthetic snowpack mtrue consisting of three snow layers. The thickness of each snow layer was 0.5 m, and the density and the liquid water content of the uppermost and the lowermost snow layer were set to 500kgm3 and 5%, respectively. For the intermediate layer, either the value of density ρ4true was varied between 100 and 900 kg m−3 with steps of 100kgm3 and the liquid water content was set to 5% or θw,4true was varied between 0% and 10% with steps of 1% and ρ4true was set to 500kgm3. The algorithm inverted for 15 parameters (density, liquid water content, and height for three snow layers, the wooden board, and the air above the radar). The initial models minit were calculated again on the basis of equations 16 and 17, and the state of prior knowledge p was varied between one and four.

Experiment 4: Minimal layer thickness

Here, we considered dry and wet snowpacks including three snow layers. The parameters of the uppermost and lowermost layers were chosen to be identical to those in experiment 3. For the middle layer, the true density was 200 or 800kgm3 and the true liquid water was 2% or 8%. We tested all four combinations of positive and negative contrasts in density and liquid water content (only two combinations in case of a dry snowpack). For all models, we systematically varied the thickness of the middle layer between 1 mm and 50 cm and we assumed prior knowledge p=4. Based on our field experience, we think that this level of prior knowledge corresponds to typical uncertainties of field measurements in snow pits. The algorithm inverted for 10 and 15 parameters for dry and wet snowpacks, respectively.

RESULTS

Before the specific results of experiments 1 to 4 are discussed, we illustrate the general performance of our FWI algorithm using one of the six-layer models generated for experiment 1. Figure 4a and 4b shows the true density and liquid water content models (green lines) and the corresponding initial models (black lines). Here, we assumed prior knowledge at a level of p=2. In Figure 4c, the frequency-domain data Wtrue(ω) (green line) and the predicted response Winit(ω) (using minit, black line) are displayed. The corresponding time-domain data wtrue(t) and winit(t) are shown in Figure 4d.

After having included frequencies up to f1=0.2GHz, the model fit improved significantly (the blue lines in Figure 4a and 4b), but the discrepancies between the true and predicted data in Figure 4c and 4d are still substantial, particularly in Figure 4c at frequencies beyond f1. After including frequencies up to f2=0.4GHz, the data fit has greatly improved (Figure 4c), but there are still considerable deviations between the estimated and true model parameters (the cyan lines in Figure 4a and 4b). When frequencies up to f3=2GHz are included, there is an excellent fit of true and predicted data and true and estimated model parameters (the red lines in Figure 4). No significant improvements can be achieved by further adding frequencies up to 4.0 GHz.

Experiment 1: A priori knowledge requirements

It is evident that the performance of any FWI algorithm improves with the closeness of mtrue and minit. In contrast, if the initial model is too far away from the true model, cycle skipping may occur and the inversion algorithm will not converge. With the simulations for experiment 1, we investigated how much prior knowledge is required to let the inversions likely converge. For that purpose, we introduced a measure for the model misfit:  
Φ=1Pi=1P(miestmitruem¯)2,
(18)
where P is the number of model parameters, mest is the inversion result, and m¯ is a mean value for a parameter type (500kgm3 for density, 5% for liquid water content, and 0.5 m for thickness). If the final value of the model misfit Φ was smaller than 0.02, an inversion was considered to be successful. This value was determined empirically. When the Φ value was below this threshold, it was guaranteed that the inversion has converged. Whereas all model parameters are considered for calculating the model misfit, we also evaluated the misfit of the density, liquid water content, and layer thickness parameters separately.

The results are shown with the blue bars in Figure 5c and 5d. In the case of a dry snowpack (Figure 5c) and p=1, only one of the 30 simulation experiments was successful. With p values of three or larger, all simulation experiments were successful. Evaluating the success rates of the parameters separately (not shown) showed that the layer thickness could be delineated more successfully than the density. However, it must be stated that an accurate determination of the layer thicknesses is essential for the overall convergence of the algorithm. If the layer thicknesses could not be resolved accurately, the other parameters could not be resolved either.

If the snowpack was wet (Figure 5d), no simulation experiments were successful without a priori information. With p=2, five of 30 experiments and with p=3, 26 experiments were successful. Even with p=4, one experiment did not converge successfully. As for the dry-snow case, thicknesses could be determined more successfully than the densities (not shown). Results for liquid water content are comparable with those for the densities.

To address the robustness of our findings during experiment 1, we repeated the simulations with models including 3, 9, and 12 layers. The main conclusions were essentially the same for all simulations. Therefore, we judge the results of the 6-layer model, presented in this study, to be representative for a wide range of models.

Experiment 2: Noise and imperfect picks

The results are shown with the green bars in Figure 5c and 5d. They were comparable with the results of experiment 1. With no prior knowledge (p=1), no simulation experiments were successful, neither for dry nor for wet snowpacks. With p values of three or larger, all simulation experiments were successful for dry snowpacks. A prominent difference to experiment 1 was found for p=3 in wet snowpacks (Figure 5d), in which only 13 of 30 experiments were successful compared with 26 in experiment 1.

Experiment 3: Influence of heterogeneity

Figure 6a and 6b indicates the range of variation of density and liquid water content of the middle layer. Success was again defined for Φ<0.02, when considering all model parameters. In all simulations with 10 initial models each, at least one of the initial models converged to the true model. Hence, if we used the combined global search and local minimization algorithm, all simulations were successful. Therefore, we show for each simulation the mean Φ for all 10 initial models (Figure 6c and 6d). For p2 (i.e., the bottom two rows in Figure 6c and 6d), the overall success rate was 43%. In other words, 228 out of the 400 initial models did not converge (Φ0.02). For p>2, only 10 initial models did not converge (success rate of 98%). It is interesting to note that the likelihood of failure seems to depend neither on the magnitude nor on the sign of the contrasts.

Experiment 4: Minimal layer thickness

We investigated the minimal layer thickness that can be resolved, with models exhibiting all combinations of contrasts (Figure 7a and 7b). Here, we were primarily interested in the thin layer’s properties, i.e., its thickness, density, and liquid water content. In contrast to the previous experiments, we display the individual misfits for these three quantities as a function of the true layer thicknesses (Figure 7c–7g).

For the dry scenario (Figure 7c and 7e), the inversions for all model types were successful when the layer thickness was greater than approximately 0.05 m. For the wet scenario (Figure 7d, 7f, and 7g), the results are less conclusive. As for the dry scenario, the layer thickness could always be resolved accurately, but the thin-layer density and liquid water content measurements were prone to inaccuracies. However, the remaining parameters of the models shown in Figure 7a and 7b could always be resolved accurately, regardless of whether the thin-layer properties could be determined or not (not shown).

DISCUSSION

Our simulations demonstrated that FWI with upGPR data is generally feasible, but there are a number of prerequisites and constraints that have to be met. Performing the inversions in the frequency domain using progressively higher frequencies proved to be critical. In fact, additional tests with time- and frequency-domain inversions involving a broad range of frequencies already at an initial stage yielded unsatisfactory results.

Results, shown in Figures 5 and 6, highlighted the importance of the initial model and thus the amount of a priori information available. Based on meteorological data, it is often possible to predict if a dry snowpack can be expected or if the amount of liquid water may be significant. Furthermore, a priori information can be obtained either from results of previous measurements and/or from snow pits. Considering the importance of the initial model, we recommend occasional calibrations with snow pits (e.g., in periods when the avalanche danger is low). It is expected that the amount of information obtained from snow pits can be associated with prior knowledge p4. If an upGPR device is operated continuously and acquires data in relatively short time intervals (a few hours or even less), it can be expected that the snowpack properties do not vary substantially between two measurements. Using results from previous measurements as an initial model will likely be equivalent to prior knowledge p>4.

In this study, we used 10 initial models and chose the solution with the smallest discrepancy r. The inversion could further be improved using more than 10 initial models, which increases the probability that one of the initial models is already close to the true model. A trade-off between the computational time and number of initial models would have to be found.

The misfit in layer thickness was always smaller than the misfit in density or liquid water content. The reason for this finding may be that we did benefit not only from the advantages of frequency-domain FWI, but we also exploited information from the time-domain radargrams by calculating the thicknesses of each layer from picked two-way traveltimes (equation 17). It is expected that picking the reflection signatures from the observed radar traces should be possible, at least for prominent layers.

Thin layers play a key role in snow stability evaluation. Therefore, it is critical to identify these small-scale features. Our simulation results indicated that in a dry snowpack, thin layers of the order of a few centimeters can be identified with upGPR data. Considering typical upGPR wavelengths of approximately 10–20 cm, this is a surprisingly good result. However, in the presence of liquid water, inaccuracies in the density and liquid water content of the thin layer must be expected.

We consider the results presented in this study to be an important first step toward FWI applied to field data. To finally achieve this goal, it is necessary to address a number of additional research tasks:

  1. Here, we assumed the source function A(ω) to be known, but in field applications, this quantity needs to be estimated as well. A possible option is to estimate the source function up-front from the observed field data (Babcock and Bradford, 2014). Alternatively, the source function could be treated as an additional unknown during the inversion (Maurer et al., 2012). Because upGPR measurements include only a single radar trace, such an approach may fail because the problem is underdetermined. Possibly, the problem can be solved by simultaneously inverting for a suite of upGPR measurements and assuming that the source function remains constant during all experiments.

  2. Field data are inherently noisy. Our results from experiment 2 (Figure 5) indicate that satisfactory results can still be achieved in the presence of realistic noise (approximately 18 dB) and picking errors. However, our simulations also highlighted the importance of a priori information in the case of noise-contaminated data. Furthermore, ambient noise is expected to be significant primarily in the higher frequency range. As shown in Figure 4, convergence of FWI can be achieved with low to intermediate frequencies, and high frequencies could even be ignored.

  3. Realistic GPR waves are influenced by the radiation patterns of the transmitter and receiver antennas, and the waves exhibit spherical divergence (i.e., the amplitudes decrease with traveling distance). Both aspects were not considered with our plane-wave approach. Because upGPR considers only vertical incident wavepaths, the antenna radiation patterns do not influence the waveforms. Babcock and Bradford (2015) show that the plane-wave approximation produced similar data as obtained with spreading-corrected simulations using a 2D finite-difference time-domain forward algorithm. The plane-wave approach violates the far-field condition within one to three wavelengths. Because the upGPR setup, as shown in Figure 2, includes a 0.25 m air layer and a 0.05 m wooden board, only the lowermost 0.2–0.3 m of the snowpack may be slightly affected by the plane-wave approximation, which is not judged to be critical.

  4. During our numerical simulations, we always included as many layers into inversion models as we were choosing for the true models, but for the inversion of field data, the number of layers is generally unknown. We performed several tests in which we tried to adaptively add or remove layers during the inversions, but we were unable to establish a robust strategy. Visual inspection of the time-domain radargrams, from which the reflection times TL are picked, offers a good estimate of the number of layers, but more research is required for an optimal parameterization of the models. This is particularly important for the detection of thin layers.

CONCLUSIONS

The upGPR systems potentially offer powerful means for continuously monitoring snowpack characteristics in a nondestructive way. We have investigated whether synthetic upGPR data contain enough information for inferring the density, liquid water content, and thicknesses of each layer of a snowpack by FWI. For that purpose, we have implemented a frequency-domain inversion algorithm that is based on a Nelder-Mead simplex solver, and it allows a priori constraints to be included in a straightforward manner.

Our results indicate that FWI applied to upGPR data is generally capable of delineating a detailed stratigraphy of the snowpack if a good initial model is available. Very good results can be obtained for a dry snowpack, in which one needs to invert for density and layer thickness only. In case of a wet snowpack, the liquid water content needs to be considered as an additional material parameter. Such inversions are feasible, but they require an enhanced a priori knowledge, and convergence of the algorithm is not always guaranteed, particularly in the presence of noise. FWI is also capable of identifying thin layers with thicknesses well below typical upGPR wavelengths, but the physical properties of such a thin layer may not be well resolved in a wet snowpack.

From our simulation results, we conclude that applications to field data are feasible, but additional issues, such as determination of the source function and optimal model parameterizations, need to be addressed.

ACKNOWLEDGMENTS

We thank the associate editor and three anonymous reviewers for their very valuable input, which helped us to improve the performance of our inversion algorithm and finally led to a much improved manuscript. L. Schmid was supported by the Swiss National Science Foundation (SNF 200020_144390).