Interactive comment on “ Uncertainty budget in snow thickness and snow water equivalent estimation using GPR and TDR techniques

The manuscript "Uncertainty budget in snow thickness and snow water equivalent using GPR and TDR techniques" by Di Paolo et al. presents an approach to quantify uncertainties from radar determinations in snow depth and SWE. The authors compare measured two-way travel times (TWT) with density and snow depth measured conventionally and TDR point observations in snow pits. From inclusion of device specific errors, they assess GPR specific uncertainties. Such an approach is not novel but interesting. Several studies before assessed uncertainties in conversion of TWT to derive SWE and snow depth (Lundberg with several studies and Sundstroem et al. 2012 for instance). Due to several severe misinterpretations and mistakes conducted mainly during fieldwork, I consider this manuscript (MS) not being sufficient for publication in


Introduction
Snow water equivalent (SWE) is a quantity computed as the product of snow depth and snow density that expresses the thickness of water that would theoretically result if the entire snowpack was instantaneously melted.In a conventional approach, snow thickness is measured using and hand probe (snow rod) and density is usually measured along the wall of a snow pit sampling and weighing the snow (Kinar and Pomeroy, 2015 and references therein).
Several alternatives to such methods have been proposed and tested by various authors, as comprehensively described by a recent literature review article (Kinar and Pomeroy, 2015).Actually, as two different physical parameters are needed (snow thickness and density) it is difficult to find a single instrument capable to properly measure both quantities, thus a strategy combining two different techniques is usually required.The applied techniques should be both accurate and precise to minimize any possible bias affecting the overall evaluation of the snow water equivalent and to provide values that are affected by the minimum possible uncertainty.Furthermore such techniques should be suitable to cover large areas in a reasonable time and, possibly, should be non-destructive.
One of the most promising geophysical tool to quickly and extensively measure snow thickness is Ground Penetrating Radar (GPR) (e.g., Godio, 2009;Forte et al., 2015), due to the high transparency of snow to radio waves (Annan et al., 1994).GPR provides an electromagnetic image of the subsurface as the result of the interaction between the emitted pulses and the electromagnetic properties of the materials in the subsurface (Jol, 2008).The transmitting antenna (Tx) sends short pulses into the snow and the receiving antenna (Rx) gathers such pulses after they have been propagated downward and reflected back by any dielectric interface (e.g., snow/ice or snow/rocks) present at depth (Annan et al., 1994).GPR, as any other radar, measures the time it takes for the pulses to make a round trip Tx-reflecting interface-Rx, named two-way travel time.To convert this parameter into snow thickness the electromagnetic velocity of the radar pulses should be known.In principle, GPR is capable to measure the pulse velocity if used in specific configurations or if appropriate targets are present in the snow (Eisen et al., 2003;Bradford and Harper, 2005;Bradford and Harper, 2006;Bradford et al., 2009;Godio, 2009;The Cryosphere Discuss., doi:10.5194/tc-2016-267, 2016 Manuscript under review for journal The Cryosphere Published: 1 December 2016 c Author(s) 2016.CC-BY 3.0 License.Gustafsson et al., 2012;Forte et al., 2014;Holbrook et al., 2016).In practice, however, there are many cases where these measurements are unfeasible or unreliable like, for example, when a fixed offset configuration is used, thus the evaluation of the wave velocity requires a support from an independent method.Time Domain Reflectometry (TDR) probably represents the best choice as it is conceptually similar to GPR and can be easily applied to snow studies (Stein and Kane, 1983;Lundberg, 1997;Stein et al., 1997;Schneebeli et al., 1998;Stacheder et al., 2005).
Numerous papers have addressed the issue of snow thickness estimation, density evaluation and snow water equivalent calculation using GPR alone, together with some other classical methods (Annan et al., 1994;Sand and Bruland, 1998;Lundberg et al., 2010;Marchand et al., 2001;Bradford and Harper, 2006;Lundberg et al., 2006;Bradford et al., 2009;Godio, 2009;Gustaffson et al., 2012;Sundström et al., 2012;Forte et al., 2013;Godio and Rege, 2016;Holbrook et al., 2016) or combining GPR and TDR techniques (Harper and Bradford, 2003;Previati et al., 2011;Di Paolo et al., 2015).However, our prospective is quite different as we are interested in understanding how accurate and precise the estimation of the snow water equivalent performed using GPR (supported by another method for wave velocity estimation) can be.Indeed, to our knowledge, very few papers have approached the uncertainty issue in GPR data (e.g., Barret et al., 2007;Lapazaran et al., 2016), even if uncertainty estimation is fundamental when a quantitative analysis is required as, for example, in hydrological and climate change studies.In this work we computed step by step the value and the associated uncertainties of each parameter present in the equations, from the quantities effectively measured by each instrument to the final snow water equivalent values.In particular, because the uncertainties propagate through the calculations, we evaluated the "weight" of each uncertainty on the overall results.Moreover, to avoid any misunderstanding in the definition of measurement errors and uncertainties and to give a solid basis to uncertainty computation, in this work, we followed the NIST (National Institute of Standard and Technology) guide lines on the expression of uncertainty in measurements (JCGM 100, 2008).Our approach allowed us to develop a procedure to accurately calibrate GPR in field conditions, that is, not only calibrate the radar system (which can be properly and accurately done in laboratory), but rather evaluate the overall performance accounting for systematic and random phenomena which can affect the field measurements.

Rationale of the measurement procedure
The goal of a measurement is to establish the numerical value of a quantity; however, as no measurement is completely free of uncertainty, an accurate evaluation of such uncertainty is part of the measurement process.The uncertainty should always be appropriately estimated and declared together with the results of the measurement, otherwise the reliability of the experiment could be questionable (Taylor, 1997;Kirkup and Frenkel, 2006).Furthermore, the consistency (or discrepancy) of independent measurements can only be proven, as will be discussed later, if the uncertainties associated to such quantities are known.
Uncertainty estimation procedure has been debated for a long time among physicist and experimental scientists.When the data acquired during a laboratory and/or a field experiment is represented by a set of numbers, assuming that these data are the outcomes of repeated measurements, the uncertainty evaluation is straightforward, as it can be computed using statistical analysis (Taylor, 1997).However, there is a wide range of cases where the uncertainties cannot be estimated in such a way and different procedures should be adopted.Quantitative geophysical studies, which are aimed at retrieving subsurface physical parameters (e.g., hydrogeophysics) are particularly affected by this type of problem; in addition, in field work it is not always easy to correctly evaluate the various contributions to the overall uncertainty of the retrieved physical parameter.
The main goal of this study is to develop a strategy applicable to field data, to accurately estimate snow thickness (ST) and snow water equivalent using GPR.In theory, GPR is a self-consistent technique in snow thickness estimation as it is capable to measure the signal two-way travel time (twt) and the wave velocity, if specific techniques are employed, like: Common-Mid Point (CMP) (Eisen et al., 2003;Bradford and Harper, 2006;Godio, 2009;Gustafsson et al., 2012), hyperbola fitting (Bradford et al., 2009), migration velocity analysis (Bradford and Harper, 2005;Holbrook et al., 2016) and reflected amplitude analysis (Forte et al., 2014).However, both two-way travel time and velocity parameters are subjected to uncertainties that are difficult to be evaluated, especially in field measurements, therefore they are usually neglected.In particular, the uncertainty on two-way travel time measurement is strongly affected by the unknown "real" position of time zero on the GPR time scale and the unknown "real" shape and time length of the GPR signal after propagation and reflection (Economou and Kritikakis, 2016).It follows that it is not easy to properly pick the first break on transmitted and reflected signals and to compute the relevant uncertainties.On the other hand, the uncertainties on wave velocity measurements have been rarely evaluated (Barret et al., 2007).
In the present work, we collected four independent sets of data (TDR, density, snow thickness and GPR) in three different sites, according to the measurement chain sketched in Fig. 1.We acquired TDR and density data along the wall (vertical direction) of a snow pit and we estimated and compared the wave velocity computed with these two independent methods ("Pit measurements" square in Fig. 1).The consistency among the estimated velocity values, allowed us to apply a calibration procedure to compute the uncertainty on the GPR two-way travel time as follows.We selected the data acquired with the hand probe (HP) for specific snow thicknesses, than using TDR velocity data and assuming a "theoretical GPR" generating and receiving only Dirac delta function signals we converted such data in two-way travel times ("predicted twt estimation" square in Fig. 1).We extracted the corresponding GPR data (i.e., collected in the same positions) and we estimated the uncertainties associated to the two-way travel time measurements under the assumption that, within the uncertainty interval, the wave velocity is constant in the three sites.Consequently, the variations in the two-way travel time can only be ascribed to changes in snow thickness ("HP-GPR comparison" square in Fig. 1).Finally we estimated the depth of the snow and the associated uncertainties combining TDR and GPR data and, from that, the snow water equivalent.These values were compared with those computed using density data collected in the pit ("SWE estimation and comparison" square in Fig. 1).
All uncertainties associated with the snow physical parameters retrieved using this procedure, were calculated applying the most up to date approach to measurement uncertainty estimation, i.e., the "Guide to the expression of uncertainty in measurements" (JCGM 100, 2008).
The test site area is located in the Eastern Italian Alps, between Trentino-Alto Adige and Lombardia regions, in the Ortles-Cevedale group, the largest glacierized group of the Italian Alps (Carturan et al., 2013).In April 2014, we performed a field campaign in the nearby of the Solda glacier (coordinates 46°29′11″ N, 10°35′48″ E), at the base of the highest peaks of the Ortles-Cevedale group: Mt.Ortles (3905 m), Mt.Gran Zebrù (3851 m) and Mt.Cevedale (3769 m), which are aligned in a NW-SE direction; at the time of the survey, the snow cover had reached its maximum thickness, before starting to melt in the following month.The elevations of our test area were comprised between 2600 and 2800 m AMSL.We selected three different sites for our survey, which comprise two almost flat snow covered areas and the tongue of the Solda glacier, which was completely covered by snow during the period of the measurement campaign.

Uncertainty estimation
For many years, in the scientific literature, the terms "measurement error" and "measurement uncertainty" have been used as synonyms.The "Guide to the expression of uncertainty in measurements" (GUM for short) (JCGM 100, 2008) clearly distinguishes these two terms, as described in Fig. 2 for the specific problem of snow thickness measurement.ST true is the true, but unknown, value of the snow thickness to be measured by any available technique; such quantity is named in GUM as the measurand.The value ST meas is the result of the measurement.The error in the evaluation of the snow thickness is unknown as well, as we do know neither the value of ST true nor all possible sources of errors (random or systematic) in the measurement.The best we can do is to estimate the value of the measurand, i.e., ST meas , and compute the associated uncertainty u, taking into consideration all the possible contributions to such value.In this framework, ST meas represents the best estimate of the measurand, i.e., the snow thickness, and the uncertainty interval represents a range of possible values that are consistent with all observations and data and "…that with varying degrees of credibility can be attributed to the measurand" (JCGM 100, 2008).
The evaluation of a physical quantity and its associated uncertainty is more reliable, for obvious reasons, in a laboratory experiment; nevertheless the same rigorous approach can be applied in field measurements, if the uncertainty associated to each quantity accounts for all random and systematic effects that can be identified in the specific measurement procedure.
For each measured parameter two possible approaches can be used (JCGM 100, 2008).i) If the quantity is the result of a series of repeated measurement, the best estimated value of such quantity (ST meas in our case) is the arithmetic mean of the values, and the uncertainty is computed as the standard deviation of mean, assuming that the data set follows a normal distribution.This procedure is named by GUM Type A evaluation of standard uncertainty.ii) If the quantity is not the result of such repeated measurements (e.g., the thickness of the snow can only be measured once) the best estimate of the measurand is the single value ST meas , and the associated uncertainty is computed on the basis of the available information assuming an a priori probability distribution from which the standard deviation can be calculated.This procedure is named by GUM Type B evaluation of standard uncertainty (for details see JCGM 100, 2008).
In this work, we mainly deal with Type B standard uncertainties as only in the procedure to evaluate the GPR two-way travel time uncertainty we performed a series of repeated measurements (Type A standard uncertainty).Indeed, for the uncertainty estimation of the quantities directly or indirectly measured in the snow, we have usually assumed an a priori uniform (i.e., rectangular) distribution.Such a choice is dictated by the lack of knowledge about the parameters controlling the measured quantities in the field.In fact, in the uniform distribution the probability density is zero everywhere except in a particular region where the probability is constant, thus it is assumed that only the lower and upper bounds of the values are known and nothing can be said about the distribution inside this interval (Kirkup and Frenkel, 2006).The use of this type of distribution usually introduces in the calculation relatively large uncertainties but, on the other hand, prevents the underestimation of important effects on the uncertainty computation.For this distribution the standard uncertainty is  =  √3 ⁄ , where a is the half-width of the distribution.Note that in one case we have chosen a different distribution (normal distribution), for which the appropriate standard uncertainty has been computed (see Sect. 4.3).
The standard uncertainties have been used to estimate the combined standard uncertainties   of the various quantities Q, under the assumption that the uncertainties   ,   ,   , … are independent and uncorrelated, applying the following equation: Appendix A reports the values of all standard uncertainties associated to the measured quantities.
Finally, where appropriate, the estimated values of the uncertainties were compared considering two confidence levels, 68 % and 95 % that is, assuming infinite degrees of freedom, computing the expanded uncertainty  =  using a coverage factor  = 1 or  = 2 (JCGM 100, 2008).

TDR measurements
TDR measurements were acquired using a Tektronix 1502C Cable Tester connected to a three-prongs probe with conductors having a diameter of (4.00 ± 0.05) mm and a length  = (30.0± 0.3) cm.The distance between the central and the external conductor is (3.20 ± 0.03) cm and the theoretical characteristic impedance in air (165 ± 2) Ω, calculated according to Ball (2002).The measurements were performed inserting the probe horizontally in the pit wall at various depths.The depth of the was estimated applying the derivative method (see e.g., Mattei et al., 2006) and the wave velocity was computed as: The combined standard uncertainty on the velocity was computed using Eq. ( 1) and the Type B uncertainties on probe length and travel time (  ,   1 and   2 ) estimated as described in Sect.3.1 (uncertainty values reported in Table A1).

Density measurements
Snow density was measured sampling vertical cores nearby the edge of the snow pit, to obtain a continuous density profile vs. depth.Each sample was collected using a cylindrical corer, with a diameter  = (4.700± 0.029) cm and a height  = (30.00± 0.29) cm (i.e., a volume of about half liter).The mass m of the samples was measured using a dynamometer.
The combined standard uncertainty on the density was estimated using Eq. ( 1) and the standard uncertainties   ,   and   associated to the mass of the snow and the height of the corer, respectively.

GPR measurements
GPR profiles were collected along 7 transects for a total of 13 profiles, using a bistatic EkkoPro GPR system (Sensors and Software, Inc) equipped with 250, 500 and 1000 MHz antennas.GPR acquisition was performed in single fold mode with common offset, automatically collecting the data with an odometer at specific step size (5, 10 or 15 cm).A 4 traces stacking was set for all profiles and different time windows (from 70 to 300 ns) were chosen depending on the location of the profile.
As reported in Table 1, all sections (aside from profile AB which has been investigated only with the 250 MHz antenna) have been acquired using 1000 MHz antennas.Profiles CD, OP and QR have also been investigated with the 250 MHz antennas and profiles GH, MN and IL with the 500 MHz antennas (see Table 1).Note that, among all profiles, only CD, GH and OP pass close to the snow pits.
The quality of the radar cross sections was very good for the entire data set.The interface snow/bedrock or snow/ice, as well as some internal layering in the snow pack, are well detectable in each section even without any gain applied, due to the very low signal attenuation of the snow.Figure 3

Snow thickness measurements
Snow thickness was measured along each radar profile using an avalanche hand probe; the measurements were performed after radar data acquisition to measure the snow depth with the same degree of compaction (created by the weight of the antennas).The step size between consecutive measurements varied from 4 to 5.5 m depending on the profile.Snow thickness was ranging along the profiles from 1 to 4 m, depending on topography and snow cover, with an average thickness of about 2.50 m.As the hand probe has marks every 50 cm, we estimated a Type B uncertainty associated to each thickness measurement of 14 cm, that is, the standard deviation of a uniform distribution having 50 cm width (see Sect. 3.1).
Furthermore, we estimated an uncertainty of 15 cm in the spatial position of the hand probe along the profile.

Wave velocity in snow
As discussed in Sect.3.2.1,Eq. ( 2) allows to compute the value of the snow wave velocity using TDR measurements (probe length and travel time) as input parameters.A similar but more elaborate procedure can be applied to estimate the wave velocity starting from snow density values.Indeed, to reach this goal we have first to compute the relevant permittivities and then convert those into velocities.Such conversion is straightforward for the snow as it is non-magnetic and non-conductive, so the wave velocity depends mainly on the real part of permittivity.To compute the permittivity we applied two different relationships commonly used in dry snow studies (Kovacs et al., 1995;Godio, 2009;Previati et al., 2011;Di Paolo et al., 2015).The first one is the Looyenga mixing model between ice and air, which assumes a host material with spherical inclusions, and is given by (Looyenga, 1965): where  is the snow density.We assumed   = (920 ± 10) kg m -3 (as reported by Kwok and Cunningham, 2008) and   = 3.18 ± 0.01 (as measured by Bohleber et al. (2012) in a frequency range 10 MHz -1.5 GHz).
The second one is Robin's empirical equation (Robin, 1975), given by: where  is the snow density expressed in kg m -3 .The permittivities extracted from Eqs. ( 3) and ( 4) have been converted into velocity using the well-known equation: where c is the wave velocity in a vacuum.Again, the standard propagation formula (Eq.( 1)) has been applied to estimate the relevant combined standard uncertainties on the permittivity values estimated using Looyenga and Robin models (   and    ), as well as those on the velocities (   and    ).The velocity values shown in Fig. 4 have been used to compute a "weighted" average velocity of the electromagnetic waves in the entire snow pack.Figure 5 illustrates the procedure used to compute such quantity with TDR (a) and density (b) data.
In particular, for the TDR data, which were acquired at depths   inserting the probe horizontally in the pit wall, we assumed a model of constant interval velocities, that is, a model where the velocities   are constant in layers of thickness , with the sign minus (plus) valid for even (odd) values of i.Thus we calculated the average velocity as: where  is the snow pit depth (local snow thickness).
The combined uncertainty associated to the TDR average velocity was computed applying Eq. ( 1) with   ,    * and    estimated as Type B uncertainties.Note that due to the lack of data below 1.50 m, at site 1 we applied Eq. ( 6) assuming no velocity variation below the last TDR measurement point.
A similar equation has been applied to estimate the average velocity from density data: where   are the heights of the density samples (30 cm) and    are the velocities computed with Eq. ( 5) for both Robin and Looyenga models.The combined uncertainty associated to these average velocities were computed as described in Sect.3.1, with   ,    and     estimated as Type B uncertainties (see Table A1).The results for the three sites are illustrated in Fig. 6 assuming two different confidence levels, that is,  = 1 and  = 2 as coverage factors.TDR values are in red, Robin in green and Looyenga in blue.In agreement with what already shown in Fig. 4, the measurements performed with TDR are the most accurate (about 1 % uncertainty) whereas for the velocity values computed from density data the uncertainties are slightly higher (about 2-3 %).Furthermore Fig. 6 shows that, taking into account the relevant uncertainties, no significant difference in wave velocity is observed among the three sites; this fact confirms that, in terms of electrical properties, the snow pack was rather homogeneous.This condition allowed us to develop a procedure to estimate the overall uncertainty in the GPR two-way travel time and consequently, to accurately evaluate the uncertainty associated to snow thickness and snow water equivalent computed using radar data.

Two-way travel time uncertainty estimation
The overall uncertainty associated to the GPR two-way travel time is the sum of different uncertainties due to system characteristics, measurement procedure and data analysis.These effects are difficult to be separately evaluated, especially in field measurements, thus the best way to assess the uncertainty is to use a calibration procedure.We assumed that the hand probe, which makes a direct measurement of the snow thickness, is the "best estimator" of such parameter and TDR the "best estimator" of wave velocity in the snow.In this way it is possible to compute the overall uncertainty accounting for different contributions like time base stability, antenna coupling, reflector properties and time picking procedure.To apply this method, we used the whole set of snow thickness data (i.e., collected on different profiles) to generate three subsets of measurements, one for each frequency, having the same number of samples.In particular, we have chosen, for each frequency, the thickest value of the snow from which the same number of samples could be extracted.We found 7 hand probe readings at 300 cm depth (250 MHz), 230 cm depth (500 MHz) and 280 cm depth (1000 MHz) (see Table 2); all these snow thickness values were affected, as described in Sect.3.3.2,by the same reading uncertainty of 14 cm.Then, assuming a "theoretical GPR" generating and receiving only Dirac delta functions (i.e., zero width signals) and using the average wave velocity estimated with TDR, we converted the hand probe thickness values in two-way travel time applying the following equation (Lauro et al., 2013): Where h is the snow thickness measured with the hand probe along the profiles, s the antenna separation (38, 23 and 15 cm for the 250, 500 and 1000 MHz respectively) v the average wave velocity measured with TDR in the specific site (〈  〉).
We also computed the combined uncertainty on   using Eq. ( 1) and the Type B uncertainties   and  ℎ associated to wave velocity and snow thickness (see Table 2).
Once defined for each frequency the set of seven "correct" two-way travel time values, we analyzed the radar sections to pick the corresponding traces.Considering that the spatial position of the hand probe is known with an uncertainty of 15 cm, we extracted from each radar section all traces located inside this interval around the probe.Thus, depending on the step size of the radar profile, one, three or five traces were picked and averaged to obtain one radar trace for each hand probe position.To compute the two-way travel time (  ) we applied a cross-correlation procedure to the GPR data (Lauro et al., 2013).For each radar section we defined a reference wavelet under the assumption that it represents the signal emitted by the antenna.The wavelet was extracted from a set of reflected signals collected at the snow/bedrock (first and second site) or snow/ice (third site) interface, in an area where the reflector was flat and sharp.Finally, we computed the cross-correlation between the reference wavelet and the extracted radar traces to pick the starting and ending point of the two-way travel time.
The procedure described above allowed us to compute three sets of predicted two-way travel time   and to extract three sets of measured two-way travel time   related to snow layers of 300 cm, 230 cm and 280 cm, respectively.We treated these (  ,  ) sets as the result of a repeatability (calibration) test where the   quantities were assumed as reference values and the   quantities were considered the data to be tested.Hence we estimated the combined uncertainty associated to the GPR two-way travel time as follows: In Eq. ( 9) the component   (associated to the picking procedure) is a Type B uncertainty, estimated assuming a normal distribution model, as the spectrum of the reference wavelet can be approximated to a Gaussian function having a RMS width   .That is: Conversely   (associated to the measurement repeatability) is a Type A uncertainty, which was computed using the set of seven calibration values as: . Note that Eq. ( 11) attributes to GPR the instability of the   values computed using the "hand probe + TDR" data, thus it might contribute to an overestimation of the uncertainty    .Table 2 summarizes the parameters used for the calibration and the estimated uncertainties for the three antenna frequencies together with the uncertainties associated to   .As shown by the values reported in the table, the component   is always larger than   , being about three times larger at 250 MHz and about twice at 500 and 1000 MHz.As a consequence, this component is the main contribution to the uncertainty on the GPR two-way travel time.Indeed the influence of   on the combined uncertainty (Eq.( 9)) is practically negligible at 250 and 1000 MHz and it is rather small at 500 MHz.Furthermore, it is interesting to notice that the uncertainty    on the predicted two-way travel time is quite large, being remarkably smaller than    only at 250 MHz and becoming even larger than    at 1000 MHz.

Predicted vs. Measured two-way travel time
Applying Eq. ( 8) to the entire dataset of hand probe values we converted depth into two-way travel time, we generated the scatter plots   vs.   for the three frequencies and we estimated the parameters of the least squares linear fit for each frequency.Given the distribution of the profiles/antenna frequencies among the three sites (see Table 1), we obtained three scatter plots rather different in terms of spatial location of the data; in fact the 250 MHz plot was generated combining site 1 and site 2 data, the 500 MHz only taking data from site 2 and the 1000 MHz combining data coming from all three sites.Furthermore, the data collected at 500 and 1000 MHz cover a larger two-way travel time range (10-40 ns) whereas those collected at 250 MHz are mainly distributed between 20 and 40 ns.frequencies, with slope  = 0.82 ± 0.13 at 250 MHz,  = 1.038 ± 0.083 at 500, and  = 1.023 ± 0.048 at 1000 MHz.
Note that the uncertainty associated to the slope is the standard uncertainty of such quantity; that is, the slope has 68 % chance of lying within the uncertainty (Kirkup, 1994).Regarding the intercept we found that the absolute values decrease when the frequency increases, being  = (4.8± 3.6) ns at 250 MHz,  = (−1.0± 1.7) ns at 500 MHz and  = (−0.2± 1.2) ns at 1000 MHz, where the uncertainties associated to the intercepts have the same statistical meaning described above for the slope.Furthermore, to quantitatively evaluate the degree of correlation between N pairs of values, the correlation coefficient r has to be associated to the distribution   (, ) which expresses the probability that the observed data could have come from an uncorrelated parent population (Bevington and Robinson, 2003).Small values of   (, ) imply that the observed variables are likely correlated.In our experiment we found the probabilities  −250 ( = 0.77,  = 105),  −500 ( = 0.95,  = 735) and  −1000 ( = 0.96,  = 157) always much lower than 10 -3 , indicating that all three   −   data sets exhibit a fairly good degree of correlation.

Snow water equivalent estimation
As a last step of the procedure, we computed the thickness of the snow and the snow water equivalent from GPR and TDR data as follows.We used the GPR data collected along CD, GH and OP profiles (i.e., the closest to the pits, see Table 1); as the pits were approximately 1 m long, for each pit and for each frequency we extracted and averaged all traces present in the meter of profile located in front of the pit.The snow thickness ℎ  was estimated as: Where v is the average velocity computed with TDR as described in Sect.4.1.The uncertainty on this value was computed, as usual, applying Eq. ( 1), with   and    the uncertainties associated to the average wave velocity (measured by the TDR) and GPR two-way travel time and assuming a negligible uncertainty on s.Finally, the snow water equivalent was calculated using CRIM (Complex refractive Index Model) (Annan et al., 1994) as: where the factor 0.93 accounts for the ice density reduction with respect to the water.The uncertainty    was computed combining the uncertainties   ,    and  ℎ  according to Eq. ( 1).The snow water equivalent values estimated from GPR data (for the three sites) were compared with those computed for each pit, using the following equation: where   and   are the height and the relative density of the i-th sample in the snow pit.The uncertainty    was computed combining the uncertainty    and    , according to Eq. (1).Table (3) summarizes the results obtained using Eq. ( 13) and ( 14) together with some parameters useful to quantify the discrepancies between the methods.In particular, the accounts for the uncertainties, the discrepancy has been evaluated according to the following formula: where k is the coverage factor.Note that in Table (3) the square root term of Eq. ( 15) is indicated as df (discrepancy factor).
As shown in Table 3, the ∆-parameter is of the order of 1-2 % for all profiles, except for the data collected on profile CD using 1000 MHz antennas, for which ∆= 6 %.For this profile, there is also a weak discrepancy between the two snow water equivalent values when  = 1, whereas the other data do not present any significant difference between the snow water equivalent values, even for the smallest coverage factor.

Wave velocity comparison
As one of the main goals of this work is to highlight potentials and limitations of different combinations of measuring techniques for snow thickness estimation, we start the discussion comparing TDR and density methods as wave velocity estimators.As described in Sect.4.1, we computed the wave velocity using TDR and density data from measurements collected along the wall (vertical direction) of three different snow pits excavated in the same area.As output of our analysis we reported, for each site, one set of values computed from TDR data and two sets of values estimated from density data (see Fig. 4).We actually used density data to test two different models (Robin and Looyenga) from which permittivity and, thus, wave velocity was calculated.As a consequence the two sets of values derived from density measurements cannot be considered totally independent results.In general, the velocities computed with the three methods (Eqs.2-5) are in good agreement as, considering the uncertainty intervals and the coverage factor  = 1, only few values do not overlap.However, the discrepancy among the values totally vanishes if a larger expanded uncertainty is assumed (i.e., choosing a coverage factor  = 2).The small point by point difference among the three data sets collected at each site (see Fig. 4) can be explained by the different sampling procedure used, as the TDR probe was inserted horizontally and the density probe vertically.On the other hand, the most significant differences between points (mainly present in the top portion of the snow pit in site 2 and 3) could be due to the presence of some water in the snow, as the TDR estimates the overall velocity regardless the state of the snow, whereas Robin and Looyenga models are accurate only for dry snow.This statement is also supported by the observation that whenever the measurements do not agree, TDR velocity values are generally lower than those estimated with Robin and Looyenga models, as the presence of liquid water in the snow increases the bulk permittivity (decreases the wave velocity).This finding indicates that some caution should be taken when employing empirical models, as they could introduce a bias in the velocity estimation due to their inability to take into account the presence of liquid water in the snow.The consistency among the three different methods applied to estimate the wave velocity is better seen when the average velocity of the entire snowpack is computed, as described in Sect.4.1.Indeed, assuming a larger expanded uncertainty (coverage factor  = 2) the uncertainty bars overlap nicely, as shown in Fig. 6; nevertheless a partial overlap is still preserved when  = 1.In this case, a small discrepancy is present only between the TDR average wave velocity computed for site 3 and those for site 1 and 2, whereas the velocity values estimated through density still agrees.
TDR measures the snowpack properties in a very similar fashion as GPR, thus it seems logical to rely on it for wave velocity estimation.Nevertheless our study indicates that there are other important reasons to consider this technique the "best estimator" of wave velocity: TDR is accurate (as will be better discussed later), very precise (only 1 % of uncertainty) and suitable for both dry and wet conditions.Therefore in our analyses we used TDR values (assuming  = 1 ) together with hand probe measurements to calibrate GPR data and estimate two-way travel time associated uncertainties.However, our results also demonstrate that density methods can be successfully applied to compute wave velocity if the precision requested is not too high (i.e., the uncertainty needed is above few percent) and the snow is dry.

GPR travel time: calibration and uncertainties
The average wave velocity estimation allowed us to check the lateral homogeneity of the snowpack in the investigated area as, within the uncertainty interval, a similar velocity in the three sites was found.This condition is fundamental to correctly apply a calibration procedure because in such procedure the compared quantities are two-way travel times predicted and physically measured for specific snow thickness values that are sparsely distributed in the entire investigated area.We applied a method (as described in Sect.4.2) which is conceptually equivalent to a laboratory procedure used to check the performance of an apparatus by comparison with a calibrated instrument, assuming that the combination hand probe/TDR was the "best estimator" (i.e., the calibrated instrument) of the two-way travel time.
Before discussing the results of the calibration test it is important to remember that in the data analysis we dealt with combined uncertainties given by two or more uncertainty contributions.In estimating the magnitude of an uncertainty a convenient rule of thumb is to neglect all those terms that are less than 10 % of the largest uncertainty contribution (Bevington and Robinson, 2003).However, we decided to present the results keeping all contributions (see Sect. 4) as they can help to evaluate the weight of each uncertainty on the various quantities involved in snow thickness and snow water equivalent estimation.We start analyzing the retrieved values of the uncertainty    associated to the predicted two-way travel time; we found a constant value, 1.3 ns, for all frequencies (see Table 2) dominated by the uncertainty  ℎ = 14 cm on the snow thickness, being the uncertainty on the wave velocity rather small.In principle, such uncertainty can be reduced if a different procedure to measure the snow thickness is followed; for example, if the part of the probe sticking out the snow is measured with a meter or if a hand probe with a higher sensitivity (finer scale divisions) is used (Kinar and Pomeroy, 2015).
However, it should be bear in mind that large part of the uncertainty on snow thickness is due to the measurement of a very irregular surface; therefore the use of high sensitivity instruments would not significantly reduce the uncertainty.The uncertainty    on the GPR two-way travel time is given by the combination of two contributions (Eq.( 9)): the first one, MHz and 1000 MHz are quite small, confirming that at these frequencies there is essentially no bias in the measured twoway travel times.As a consequence, for such data the estimation of   can be considered quite accurate.Conversely, the intercept value at 250 MHz is significantly larger than those at higher frequency, highlighting the presence of a nonnegligible bias.The reason of this discrepancy is not clear and should be investigated in more detail, even if it may be a consequence of the lack of data regularly distributed along the full range of the two-way travel time values.Something should also be said about the uncertainties on the b-values; in general, these are expected to be larger than the uncertainties on the   values as the process to extrapolate the fitting line back to the y-axis can introduce large uncertainties (Taylor, 1997).As a consequence, the values we found for the three regression lines, and especially for the 250 MHz data, which are clustered at larger two-way travel times, can be considered quite reasonable.

Snow water equivalent
At the end of the overall procedure, once proven that   is reliable and accurate (especially at high frequency), we converted the   values in snow thickness, using the TDR average wave velocity, and we estimated the snow water equivalent from the data collected on radar lines segments close to the three snow pits (as described in Sect.4.4).Such data were compared with those retrieved from density measurements acquired in the pits, which we considered our snow water equivalent "best estimator".As described in Sect.4.4 we evaluated the compatibility between measurements using both ∆parameter, which has been used by other authors, and Eq. ( 15) that represents a more rigorous approach as it also accounts for the uncertainties.In general, we found an excellent agreement between GPR+TDR and density methods for all data except in one case (as discussed below) (see Table 3).In fact, in our analysis ∆-parameter is of the order of 1-2 %, that is, similar or even better than those reported in previous works.For example, Sand and Bruland (1998) found 5 % at 500 MHz as maximum value for ∆-parameter, Lundberg et al. (2000) 5 % at 1200 MHz, Bradford and Harper (2006) 0.9 % at 900 MHz, and Bradford et al. (2009) 12 % at 900 MHz.Furthermore, according to Eq. ( 15), the data reported in Table 3 show that   and   values (aside from one case) are always compatible even with  = 1.Regarding the data acquired on profile CD with the 1000 MHz antennas, we found ∆= 6 % and a weak discrepancy between the two snow water equivalent values when the coverage factor is  = 1 (see Table 3).Such discrepancy was in some way unexpected as the other two measurements at 1000 MHz provided a good estimate of the snow water equivalent (see Table 3) and, more notably, the snow water equivalent value retrieved from 250 MHz data acquired on the same profile, is very accurate.
Indeed, the analysis presented in this work suggests that the 1000 MHz measurements are generally more accurate and precise than those performed with the low frequency antenna, thus the reason of the discrepancy should be sought into the measurement conditions.The source of this error is not clear however it may be in some way link to the antennas properties.
In fact, as well known, the signal amplitude acquired by the GPR antenna is inversely proportional to the frequency; thus in

The
Cryosphere Discuss., doi:10.5194/tc-2016-267,2016   Manuscript under review for journal The Cryosphere Published: 1 December 2016 c Author(s) 2016.CC-BY 3.0 License.three pits was 2.94 m (site 1), 3.33 m (site 2) and 3.21 m (site 3) respectively, and matches in each site the overall thickness of the snow.In particular, in the first site only the top half of the pit was investigated (at irregular spacing) as some connection problems raised during the first day of data acquisition.In any case, in the second and third site, TDR measurements were regularly performed about every 30 cm to match the sampling interval of snow coring (see below), along the vertical length of the pit.The travel time ( 2 −  1 ) associated to the one-way signal propagation in the transmission line illustrates, as an example, the section acquired in the first site along profile AB with 250 MHz antennas.The Cryosphere Discuss., doi:10.5194/tc-2016-267,2016 Manuscript under review for journal The Cryosphere Published: 1 December 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 4
Figure 4 summarizes the results for the three sites.The velocity values extracted from TDR data are single points, blue dots with the associated uncertainty bars, and those computed from the density data are represented by solid segments (having the length of the corer, i.e., 30 cm).The values calculated from Robin's model are in red and those from Looyenga's model in black.The uncertainties associated to these values are represented by dash segments of the same color red and black, respectively.Note that TDR data show the smallest uncertainties whereas those associated to Robin and Looyenga models are comparable.At site 1, as described in Sect.3.2.1, the TDR data are only available for the first 1.50 m of snow cover, whereas in sites 2 and 3 the velocity values extracted from TDR measurements extend down almost to the bottom of the pit.All sites show the same trend: slightly higher velocity values near the surface, a decrease of velocity with depth up to 1.50 m and approximately constant values below such depth.

Figure 7
Figure7illustrates the scatter plots at the three frequencies, the fitting lines and the relevant parameters.All fits have been computed without imposing a zero intercept to check for possible biases.The plots indicate a fairly good linear trend for all

The
Cryosphere Discuss., doi:10.5194/tc-2016-267,2016 Manuscript under review for journal The Cryosphere Published: 1 December 2016 c Author(s) 2016.CC-BY 3.0 License.percentage difference between the estimated values of the snow water equivalent (not accounting for the uncertainties) has been computed as ∆= 100 • |  −   |   ⁄ (to allow a comparison with other reference values) whereas to

The
Cryosphere Discuss., doi:10.5194/tc-2016-267,2016   Manuscript under review for journal The Cryosphere Published: 1 December 2016 c Author(s) 2016.CC-BY 3.0 License.be explained by the low statistics at short two-way travel times as the low frequency data were collected (by chance) mostly where the snow was deeper (20-40 ns range).Regarding the intercepts (b-value), we can observe that those retrieved at 500

Figure 1 :
Figure 1: Schematic of the measurement procedure.

Figure 2 :
Figure 2: Error and uncertainty in snow thickness (ST) measurement according to GUM.

Figure 3 :
Figure 3: Radar cross section acquired using the 250 MHz antenna along trisect AB.The reflections from the snow pack/bedrock 10

Figure 4 : 5 The
Figure 4: Velocities estimated from TDR and density measurements on the wall of the three pits.Blue dots with the uncertainty bars are TDR values, red (Robin's model) and black (Looyenga's model) segments the values computed from density data.Dash segments indicate uncertainty intervals.5

Figure 5 :
Figure 5: (a) Sketch of the positions of the TDR probe   and thickness of the snow layers   * used in Eq. (4).(b) Sketch of density sampling along vertical cores with height   .

Figure 6 :
Figure 6: Average wave velocity computed from TDR and density data in the different sites.Red dots are TDR values, green dots Robin and blue dots Looyenga values.Uncertainty bars are computed using a coverage factor  =  (solid lines) and  =  (dashed lines).

Figure 7 :
Figure 7: Least squares linear fit for the three antennas frequencies and relevant fit parameters.Note the different two-way travel time ranges among the scatter plots.

Table 3 .
Snow water equivalent measured in snow pits and estimated from GPR data using CRIM.