## 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 $M\u22121$ 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

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.

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

### Inversion

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 $\theta w,L$ to lie within a physically reasonable range (i.e., $\theta w,L<10%$), or we can impose that no liquid water is present (i.e., $\theta w,L=0$). Furthermore, the upper limit for densities $\rho L$ was set to $917\u2009\u2009kg\u2009m\u22123$. 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\xd710\u22124$.

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 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:

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

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

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?

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

All our simulations followed a workflow as described below:

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

Calculate $TL$ from $wtrue$.

Set $fmax$ to 200 MHz.

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

Calculate the thicknesses using equation 17.

Perform an inversion using $W(\omega )$ for $\omega <2\pi fmax$ as $dobs$ (i.e., all frequencies up to $fmax$ are used in the inversion).

Increase $fmax$ by 200 MHz.

Repeat steps five to seven until $fmax=4\u2009\u2009GHz$.

#### 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 $900\u2009\u2009kg\u2009m\u22123$. 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 $500\u2009\u2009kg\u2009m\u22123$ and 5%, respectively. For the intermediate layer, either the value of density $\rho 4true$ was varied between 100 and 900 kg m^{−3} with steps of $100\u2009\u2009kgm\u22123$ and the liquid water content was set to 5% or $\theta w,4true$ was varied between 0% and 10% with steps of 1% and $\rho 4true$ was set to $500\u2009\u2009kg\u2009m\u22123$. 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 $800\u2009\u2009kg\u2009m\u22123$ 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(\omega )$ (green line) and the predicted response $Winit(\omega )$ (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.2\u2009\u2009GHz$, 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.4\u2009\u2009GHz$, 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=2\u2009\u2009GHz$ 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

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 $\Phi <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 $\Phi $ for all 10 initial models (Figure 6c and 6d). For $p\u22642$ (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 ($\Phi \u22650.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 $p\u22654$. 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:

Here, we assumed the source function $A(\omega )$ 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.

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.

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.

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).