1. Introduction
The snow cover, and in particular its stratigraphy and temporal evolution, is important for many applications including land-surface/climate interactions, water resource management, avalanche forecasting and flood prediction (Reference Barry and GanBarry and Gan, 2011). Currently applied methods for in situ measurements are destructive and do not allow the temporal evolution of the snowpack to be followed. These limitations often result in low spatial and temporal resolution of snow stratigraphy data.
This deficiency can be overcome by continuous point measurements, remote sensing, numerical modeling or a combination of these methods. For example, to derive snow height and snow water equivalent (SWE) for large areas from spaceborne microwave measurements, a combination of ground observations and model data is most promising (Reference Tedesco, Derksen and PullianenTedesco and others, 2012). Radar systems buried in the ground allow for temporally continuous point observations (Reference Heilig, Eisen and SchneebeliHeilig and others, 2010), even in steep terrain where previous instrumentation failed or could be destroyed by avalanches. Modeling has the advantage of not being limited to a certain location (Reference Durand, Giraud, Brun, Mérindol and MartinDurand and others, 1999; Reference Lehning, Völksch, Gustafsson, Nguyen, Sta¨hli and ZappaLehning and others, 2006) as long as the meteorological input data can be extrapolated with reasonable accuracy. However, extrapolating meteorological data is known to be challenging, particularly as a result of snow accumulation or erosion by wind.
For avalanche formation, the snow stratigraphy, i.e. the thickness and properties of the layers that form the seasonal snow cover, is one of the key contributing factors (Reference Schweizer, Jamieson and SchneebeliSchweizer and others, 2003). The existence of a slab and weak-layer combination prone to avalanching can so far only be detected by manual observations (e.g. by in situ sampling of the snow stratigraphy). These point observations are cumbersome, of low temporal resolution and not always and everywhere possible owing to safety concerns. In addition, since the observations are destructive, they do not allow the temporal evolution of snow stratigraphy to be captured at exactly the same location. For assessing the avalanche danger due to wet-snow avalanches, the location, amount and dynamics of movement of liquid water within the snowpack are of particular importance (Reference Mitterer, Heilig, Schweizer and EisenMitterer and others, 2011a). Measurement techniques for the volumetric liquid water content (θ w) in snow are well developed and based on the dielectric properties of snow. Most measurement methods require an open snow pit and thus are destructive (Reference DenothDenoth and others, 1984; Reference DenothDenoth, 1994). Time-domain reflectometry (TDR) allows for non-invasive continuous monitoring of θ w within the snowpack (Reference Schneebeli and JohnsonSchneebeli and others, 1998; Reference Waldner, Huebner, Schneebeli, Brandelik, Rau and DowdingWaldner and others, 2001). However, it is not suited for operational monitoring in avalanche starting zones owing to technical constraints, because the installation of TDR sensors requires poles, which makes the system prone to snow creep or destroyable by avalanches.
Continuous, operationally working measurements are in general limited to flat-field test sites and include only a few snow parameters such as snow height, snow temperature, snow surface temperature and liquid water content at the bottom of the snowpack. All these measurements are not feasible in an avalanche starting zone from where the information is most needed for local avalanche forecasting. Snowpack monitoring in avalanche starting zones requires a system that works independently of the prevailing weather and avalanche conditions, that cannot be destroyed by avalanches and that provides near real-time information on snow stratigraphy, at least on snow accumulation, ablation, erosion and settlement. To fulfil these requirements, a system monitoring the snowpack from below seems most feasible.
The non-destructive recording of snowpack properties can be achieved with radar systems such as frequency-modulated continuous-wave (FMCW) and ground-penetrating radar (GPR). Various measurements with FMCW radar systems from above and below the snowpack were conducted in the past (Reference Gubler and HillerGubler and Hiller, 1984; Reference Schmidt, Gubler and HillerSchmidt and others, 1984; Reference Gubler and WeilenmannGubler and Weilenmann, 1986). Reference Gubler and WeilenmannGubler and Weilenmann (1986) monitored the snowpack evolution and percolating water continuously with an upward-looking system. As they employed a radar system with a frequency range of 3.6–18 GHz, they were not able to give a quantitative value for θ w, as the radar could not penetrate the partly wet snowpack. Recently, Reference Mitterer, Heilig, Schweizer and EisenMitterer and others (2011a) showed that quantifying the liquid water content of wetted parts of the snowpack is feasible with combined radar data and external snow height measurements.
Radar technology has also been applied to determine SWE for large areas (Reference Lundberg, Thunehed and BergströmLundberg and others, 2000; Reference Marchand, Bruland and KillingtveitMarchand and others, 2001). In most cases, SWE was estimated via empirical formulae using the two-way travel time (τ). Under wet conditions this approach suffers from substantial errors when θ w is unknown (Reference Lundberg and ThunehedLundberg and Thunehed, 2000). Therefore, much research was focused on deriving the amount of water within the snowpack using the radar signal. Reference Boyne and GeorgeBoyne and George (1987) quantified volumetric θ w using C-band frequencies under laboratory and field conditions (from above the snow surface) and compared it to θ w values obtained with the dilution method (Reference Davis, Dozier, LaChapelle and PerlaDavis and others, 1985). For laboratory samples and homogeneous snow stratigraphy, results were fairly accurate (error of ± 2.5% at the 95% confidence level). Under natural conditions, however, values of θ w measured with radar were 2–10 times larger than those found with the dilution method. Reference Bradford and HarperBradford and Harper (2006) and Reference Bradford, Harper and BrownBradford and others (2009) derived the complex dielectric permittivity from GPR signals and found good agreement with snow fork measurements (Reference Sihvola and TiuriSihvola and Tiuri, 1986) for snow liquid water contents in the pendular regime (θ w less than ˜4% by volume).
Reference Marshall, Schneebeli and KohMarshall and others (2007) compared FMCW radar signals to snow micro-penetrometer (SMP) (Reference Schneebeli and JohnsonSchneebeli and Johnson, 1998) measurements and standard snow-pit observations. They were able to relate strong reflections to a thin, hard crust. In general, reflections were in good agreement with transitions in penetrating force detected by the SMP.
Reference Heilig, Schneebeli and EisenHeilig and others (2009, Reference Heilig, Eisen and Schneebeli2010) recorded the temporal evolution of snow stratigraphy with an upward-looking GPR. They determined snow height and traced different layers within the snowpack. Transitions from new snow to the old snow surface and various melt–freeze crusts within the snowpack produced distinct reflections which could be followed until the first deep penetration of meltwater.
GPR is also a powerful method in providing laterally continuous information over large distances. Therefore, radar technology was often applied to describe the variability of snow stratigraphy along transects (Reference Gubler and WeilenmannGubler and Weilenmann, 1986; Reference Harper and BradfordHarper and Bradford, 2003; Reference Machguth, Eisen, Paul and HoelzleMachguth and others, 2006; Reference Previati, Godio and FerrarisPreviati and others, 2011).
To our knowledge, in none of the previous studies was radar continuously applied to monitor snowpack evolution over the entire season, or the signal directly related to settlement of internal layers (except for ice in Reference Corr, Jenkins, Nicholls and DoakeCorr and others (2002) and Reference Gillet-Chaulet and KingGillet-Chaulet and others (2011)).
In summary, radar is a promising technology for many applications in snow science, but quantitative results on snow stratigraphy based on radar signals, in particular on the temporal evolution at a specific site, are scarce. We therefore installed upward-looking radar systems with the objective of continuously monitoring the temporal evolution of the seasonal alpine snowpack and deriving snow stratigraphy information from the radar signals. We focused on determining snow height, amount of new snow, snow settlement, liquid water content and SWE from radar signal signatures.
2. Methods
2.1 Weissfluhjoch test site and accompanying data
The measurements were performed at the flat-field test site Weissfluhjoch, which is located above Davos, Switzerland, at 2540 m a.s.l. At the same site, an upward-looking FMCW radar system was installed (Reference OkornOkorn and others, in press). The test site is fully equipped with sensors recording meteorological and snow-cover properties (Reference Mitterer, Heilig, Schweizer and EisenMitterer and others, 2011a; Reference Marty and MeisterMarty and Meister, 2012). We made use of three different snow height sensors (a laser and two ultrasonic gauges), wind, air temperature, snow temperature, water content, radiation sensors, a 5 m2 lysimeter for recording water discharge at the bottom of the snowpack, and a conventional snow pillow to determine SWE of the overlying snowpack. In addition, all the input parameters required to run the one-dimensional (1-D) snow-cover model SNOW-PACK (Reference Bartelt and LehningBartelt and Lehning, 2002; Reference Lehning, Bartelt, Brown, Fierz and SatyawaliLehning and others, 2002a, Reference Lehning, Bartelt, Brown and Fierzb) were measured at the study site.
Settlement of old snow surfaces after new-snow burial was measured utilizing the combined settlement and temperature sensor (SES) described by Reference Fierz and LehningFierz and Lehning (2001). The sensor consists of a lightweight wooden frame (30 cm × 30 cm) clipped onto two vertically fixed wires that allow the sensor’s height to be determined potentiometrically over time. Within the wooden frame a tungsten wire is stretched to record snow temperature (Reference Fierz and LehningFierz and Lehning, 2001). As the sensor is covered by new snow, it settles with the underlying snow cover. During winter 2011/12, six sensors were placed on the snow surface shortly before a significant snowfall.
Conventional manual snow profiles according to Reference FierzFierz and others (2009) were conducted on a bi-weekly basis within 2–8 m of the radar site. Snow density was determined by taking samples of 100 cm3 and weighing them on an electronic scale. For each layer recorded in the snow pit, at least two density samples were taken and averaged. The relative dielectric permittivity of each layer was measured using a capacity probe (Reference DenothDenoth, 1994).
New snow is measured by placing a board on the snow surface (Reference FierzFierz and others, 2009). Every morning at about 08:00, the height of the new snow accumulated on the board since the last observation (HN24) is measured using a ruler. The board is then cleaned of snow and placed flush with the snow surface again.
To derive the snow height automatically from the radar signal (fully automated snow surface picking algorithm; see Section 2.3.2), we used snow height, air temperature and snow surface temperature from a second automated weather station (AWS), WAN7. This station is located within the Steinta¨lli basin (flat field, 2440 m a.s.l.), 2.95 km southwest of the Weissfluhjoch test site.
2.2 Radar set-up and raw data processing
Observation of snow stratigraphy together with total snow height is only possible with sensor systems which are able to penetrate the snowpack and adequately resolve internal layers. GPR systems operating in the high-frequency range provide better resolution of reflecting layers because of their larger bandwidth and smaller wavelengths than low-frequency systems. However, signal attenuation in wet snow strongly increases with increasing frequency above ~ 1 GHz. Therefore, a compromise between penetration depth in wet snow and resolution has to be found. In this study, a commercially available dual frequency GPR system from IDS (Ingegneria Dei Sistemi, Italy) operating at 600 MHz and 1.6 GHz was utilized. According to the manufacturer, the angle of beam spread is 45°.
Since the ultimate goal is to provide information on snow properties from an avalanche starting zone, the system must be able to withstand avalanches without damage. Therefore, we installed our radar system in a wooden box, which was buried in the ground (Fig. 1). The antennas look upward (upGPR) and the radar pulses penetrated the snow from below. Measurements were conducted every 3 hours in dry-snow conditions, and more frequently (every 30 min) as soon as we expected the snowpack to become wet. The measurements were triggered autonomously and data were transferred using an existing internet connection.
Fig. 1. Assembling of radar set-up. (a) The antenna is lifted and lowered by a hoisting device during measurements. (b) The arrangement is located in a wooden box buried in the ground.
Every radar system is subject to noise, which is recorded simultaneously with the desired signal response. For stationary impulse radar systems, the signal response arising from the layered snowpack cannot be distinguished from coherent noise. Therefore, for each measurement we moved the antennas twice vertically (0.12 m) up and down with a hoisting device (Fig. 1). During this time, ~ 1800 traces were recorded. Consequently, horizontal snow-layer interfaces will be imaged in the form of undulating reflections, whereas system-related signals will be visible as horizontal features.
We processed the data in a similar way to that described by Reference Heilig, Eisen and SchneebeliHeilig and others (2010). Contrary to their set-up, the antenna system was now installed beneath the soil surface and not within a snow cave above the soil. The received radar data were further processed to increase data quality. Offsets in the zero line of each radar trace (wow) were corrected utilizing a dewow function for which a running mean over 1–5 ns was calculated. A bandpass filter cut low frequencies (<600 MHz) and high frequencies (>3 GHz) of the spectrum. As described above, the antennas are lifted and lowered twice during data collection. All signals originating from the snowpack showed well-defined up– down curvatures in the radar sections, whereas noise signals generated by the system design remained constant in the time domain, and could be eliminated by applying a background removal procedure on the first 7 ns in travel time. To compensate for geometrical spreading, we applied a gain function, which multiplies raw amplitudes by roughly twice the respective recorded two-way travel time for each trace. Assuming the geometrical spreading is equal to the square root of the area of the corresponding wavefront sector (Reference Granlund, Lundberg, Feiccabrino and GustafssonGranlund and others, 2009) and utilizing the given beam spread of 458, this gain function compensates for divergence losses in signal amplitude strength with increasing travel times for a homogeneous medium. In a further processing step, a static correction on the first reflection (which in our case originates from the wooden plate that covers the radar box underneath the snowpack) was applied. Thus, all reflections originating from the snow cover were horizontally planar, whereas multiples of the direct-wave signals and other interfering signals appeared in an oscillation relative to the vertical movement. The resulting radargram was stacked and averaged over the number of traces. All the measurements recorded over the whole season were then merged in one radar section. As the system was buried in the ground, the surrounding soil caused a permanent signal response. These signals were always constant in the time domain, but with the same up–down curvature as snowpack-related reflections. However, since these reflections remained constant over weeks, we could remove them through a moving-window time filter. All signals at the same sample position over 6 weeks were subtracted.
2.3 Picking the snow surface and internal layers
Changes in density and liquid water content cause partial signal reflections. In particular, transitions from snow to air at the snow surface result in a pronounced reflection. If several changes in density or liquid water content occur within a distance shorter than half of the wavelength λ of a short wavelet (for 1600 MHz:λ /2 ≈ 7.2 cm; Reference DanielsDaniels, 2004), constructive and destructive interferences will appear. Consequently, all signal responses originating from layer interfaces <7.2 cm away in a vertical direction cannot be resolved as single reflections. Additional complications arise from multiples generated by strong signal reflections, which interfere with the primary signals originating from the snowpack. In order to assure a consistent signal interpretation, we chose the following rules for layer picking.
-
1. For each picked layer (e.g. internal layers or snow surface), the location of the pick was set to the maximum absolute amplitude of the first half-cycle of the respective reflection signal (Fig. 2a and c).
-
2. Owing to the fact that the source wavelet consists of three full phase cycles, the signal response from the wooden box (Fig. 1) is partly masked from the source wavelet. Therefore, the first half-cycle (Fig. 2a) is not present in recorded upGPR waveforms. To solve this problem, we picked the zero-crossing before the first detectable positive half-cycle (Fig. 2b) and set this location to a two-way travel time τ = 0.4 ns. This shift corresponds to 1.5 half-cycles of the wavelet, which will set the maximum of the first (masked) half-cycle to τ = 0 ns (Fig. 2a). In addition, the reflection of the wooden box cannot be separated from the bottom edge of the snowpack above the box, resulting in an a priori error of 0.015 m (thickness of the box).
Fig. 2. Exemplary characteristics of a synthetic reflection response originating from the transition from (left) air to wooden board to snow and (right) snow to air (not to scale). The location of the pick was defined as the maximum of the first half-cycle of each reflection signal (a and c). The first two half-cycles are not visible in the wavelet (red dashed). Hence, the zero-crossing after the first two half-cycles (b) was picked and set to a two-way travel time of τ = 0.4 ns. (d) is an example of picking a wrong negative half-cycle (for details see Section 4.1).
We applied two different algorithms for picking the snow surface.
2.3.1. Semi-automated snow surface picking algorithm
The first algorithm requires manual interaction according to the following steps. We used the function ‘phase follower’ of the software package REFLEXW (Reference SandmeierSandmeier, 2010), which always follows the peak of the same half-cycle. If two consecutive traces deviated, we checked whether the height of the snow surface changed due to accumulation, settling or melt. If none of these changes appeared in the recorded weather data, and deviations in the phase sequence occurred (e.g. while surface crusts were persistent or surface melt happened), we neglected phase reversals to avoid arbitrariness. During strong accumulation and melt events, manual intervention was required to reset the follower to the correct phase.
2.3.2. Fully automated snow surface picking algorithm
Alternatively, the snow surface can be determined fully automatically. The transition between snow and air causes the largest peak in the reflected wave due to the large change in dielectric permittivity (as long as dry-snow conditions prevail). In addition, the snow–air transition is the reflection that changes most in the snowpack as snow accumulates, settles or melts. We therefore calculated the difference between two subsequent traces, sample by sample. The snow surface was assumed to be at the position where the difference was the largest. The algorithm works well as long as the snowpack is dry, but it becomes error-prone when a transition from dry to wet snow causes high amplitudes within the snowpack. Figure 3a shows the difference in the amplitudes between the measurements at 03:00 and 06:00 on 6 May 2012. During this time it was snowing. The snow height determined with the semiautomated picking algorithm was 2.53 m, but the largest difference in amplitude (highest peak in Fig. 3a) was at 0.37 m. The automated algorithm would hence strongly underestimate snow height. In order to account for this problem, we used measurements of air temperature (TA), snow surface temperature (TSS) and snow height (HS) to determine whether it was snowing (HS increasing and TA ≈ TSS), the snow was melting (TA > 0ºC and TSS > –0.5ºC) or the snowpack was settling (all other cases). The weather data for this algorithm can be taken from any weather station with comparable precipitation rates allocated at similar elevation in the surrounding area. The algorithm was then refined with a probability function for the snow surface height consisting of two half Gaussian distributions (Fig. 3b). The expected value of the probability functions (red line; 2.66 m) was calculated from the position of the previous snow height (green line; 2.61 m) and the change of the measured snow height Δ HS. The standard deviations of the probability functions corresponded to a search radius upwards (σ up) and downwards (σ down);σ up and σ down are mean accumulation and decrease rates of snow height calculated over several years at WAN7. The three weather conditions, snowing, settling and melting, were treated separately (Table 1). Figure 3c shows the product of the difference in the reflected waves and the probability density function. The highest peak of the curve is at 2.62 m. As the snow surface was defined as the minimum of the first negative half-cycle of the reflection signal (Fig. 2c), the new height of the snow surface is assumed to be at the minimum of the amplitude within a half-wavelength away from the highest peak.
Fig. 3. Method to automatically pick the snow surface. Plots show only exemplary curve patterns (absolute values not to scale). (a) Differences of the amplitudes of the reflected waves between the measurement taken on 6 May 2012 at 06:00 and the preceding measurement (6 May at 03:00). (b) Probability density function for the new snow surface. The location of the previous snow height is denoted with a green line. (c) Multiplication of the two curves. Green line shows the previous snow height, magenta line shows the highest peak and black line shows the new snow height.
Table 1 Size of window for searching next snow surface reflection for various weather conditions
2.3.3. Picking internal snow layer interfaces
Internal layers were picked in a similar way to the procedure of the semi-automated snow surface picking algorithm without any manual interactions.
2.4 GPR-derived snow properties
2.4.1. Snow height and settlement
The height (D) above the ground (snow–soil interface) of picked signals originating from internal layers and the snow surface was calculated using the two-way travel time (τ) of the reflection and the bulk velocity (v bulk) below the reflection:
The velocity in snow depends on the density and the wetness of the snow and therefore differs from layer to layer. We used two models for the velocity profile of the snowpack:
-
1. For velocity profile (i), we simply assumed a constant relative dielectric permittivity for the entire snowpack (Ɛ r = 1.7) as suggested by Reference Heilig, Schneebeli and EisenHeilig and others (2009). The corresponding velocity for dry snow is v bulk = 0.23 m ns–1.
-
2. Alternatively (velocity profile (ii)), we modeled snow densities from recorded weather data utilizing the 1-D snow-cover model SNOWPACK and calculated the bulk velocity v bulk for every signal reflection. For every layer i of the model, the velocity (vi ) was calculated:
(2)
where c 0 is the speed of light in vacuum ρi is the modeled density of the i th layer (kg m–3) and c 1 ¼= 1: 92 × 10–3 and c 2 = 4: 4 × 10–7 are constants proposed by Reference DenothDenoth (1994). For every signal reflection (two-way travel time τ), we needed to know the corresponding layer in the modeled profile. The number n of modeled layers below the reflection was determined by searching the largest possible value for n such that
(3)was fulfilled, where zi is the thickness of the i th layer. The bulk velocity v bulk was then
(4)
2.4.2. Snow density
The bulk density (ρ s) of the dry snowpack was calculated using the radar signal and an external snow height measuring device. Again, we used Eqn (1) to calculate the bulk velocity, where τ is the two-way travel time to the snow surface and D the snow height measured by a nearby laser gauge. Rearranging Eqn (2) and replacing vi by v bulk and ρi by ρ s enabled the bulk density to be determined. The bulk density can only be calculated as long as the snowpack is dry, since Eqn (2) is not valid for wet snow.
2.4.3. New snow height
A conventional instrument for measuring snow height, such as an ultrasonic sensor, only measures the distance to the snow surface. During a snowfall, snow height increases. Since the underlying layers settle while being loaded with new snow, the new snow height is always underestimated, i.e. the amount of new snow cannot be measured automatically. The radar, however, still records the reflection of the old snow surface after it was covered by new snow. We subtracted the two-way travel time of the reflection of the old snow surface from the two-way travel time of the new snow surface. Assuming a wave speed v = 0.274 m ns–1, which for dry snow corresponds to a density of 100 kg m–3, we calculated the new snow height using Eqn (1). The daily new snow height was calculated for each day at 09:00 and summed up for the entire snowfall event.
2.4.4. Dry-to-wet transition and liquid water content
We developed an algorithm to track the penetration depth of liquid water. Since the dielectric permittivity of wet snow (Ɛ r > 2) differs from that of dry snow (Ɛ r ≈ 1.7), a transition from dry to wet snow causes a strong reflection. Amplitudes in the radar signal of such a transition are larger than for any other reflection originating from changes of physical properties. Figure 4a shows the maximum of the amplitudes of each waveform (A max) during winter 2010/11. The snowpack was dry before 24 March 2011 and became wet afterwards.
Fig. 4. (a) Time series of maximum signal amplitude of each waveform. The grey background shows when the snow surface was warmer than –18C. The threshold to find the dry-to-wet transition is denoted black. (b) A waveform measured at a time when the snow was wet (9 April 2011 14:30). Green line shows the snow surface. The dry-to-wet transition (magenta line) is assumed to be at the maximum after the event where the waveform (blue line) exceeds the threshold (black lines) for the first time.
The algorithm starts with the first measurement of the winter and works chronologically through each measurement. For each measurement, the following steps are performed:
-
1. Determining whether the snow is wet or not on a given day: From numerical parameterization, we found that the snowpack is considered to be wet if the highest A max during the day is higher than 110% of A max at 06:00. If the snowpack is classified as dry, no dry-to-wet transition is assigned for that day and the algorithm continues with the next measurement. Otherwise, the algorithm continues with the next step.
-
2. Calculating a threshold to find the dry-to-wet transition: The threshold is the mean of A max for all days so far classified as dry (black lines in Fig. 4).
-
3. Picking the dry-to-wet transition: As defined above, for each picked layer the location of the pick is at the maximum amplitude of the first half-cycle of the reflection signal (Fig. 2). The dry-to-wet transition is set to the maximum of the half-cycle after the waveform exceeds the threshold for the first time (magenta line in Fig. 4b).
Having found the dry-to-wet transition with the algorithm and assuming the snowpack below this transition to be completely dry, the bulk volumetric liquid water content (θ w) of the snow above the transition can be determined. Provided that external snow height data are available (Reference Mitterer, Heilig, Schweizer and EisenMitterer and others, 2011a), θ w is derived from
with β = 0.5 (Reference Roth, Schulin, Flühler and AttingerRoth and others, 1990) and values of the dielectric permittivity of air, ice and water of 1, 3.18 and 87.9, respectively (Reference Mitterer, Heilig, Schweizer and EisenMitterer and others, 2011a). In order to continuously solve Eqn (5), we assumed the bulk dry-snow density (ρ d) to be constant during the melting season. To obtain this value, we averaged the radar-derived bulk density (ρ s) over 100 consecutive measurements before the moment when we defined the snowpack as having become wet. This corresponded to 3–4 days depending on the measurement interval. The bulk dry-snow density ρ d was 319 kgm–3 for the 2010/11 season and 406 kgm–3 for the 2011/12 season. The only unknown in Eqn (5) is the bulk relative effective dielectric permittivity (" eff), which can be determined from the speed of light in vacuum (c 0) and the bulk wave speed in the wet snowpack (v bulk) (Reference DanielsDaniels, 2004):
where v bulk is calculated using Eqn (1), in which D is the thickness of the snow above the dry-to-wet transition (difference between the external measured snow height and the height of the transition) and τ is the difference in two-way travel time between the radar-detected snow surface and the transition.
Values of θ w calculated for different times can only be compared with each other if the location of the transition has not changed in the meantime. Therefore, we slightly modified the above-described algorithm: for the calculation of θ w the transition height was not allowed to move upwards, i.e. even if the snowpack refroze, the transition stayed at the prior position until the snow became wet again. This allowed the temporal evolution of θ w to be followed until the water percolated deeper into the snowpack.
2.4.5. Snow water equivalent
For dry snow conditions, the bulk snow densityρ s was determined in combination with a conventional snow height gauge as described above. Multiplying ρ s by the measured snow height HS yields the snow water equivalent:
Calculating SWE with Eqn (7) assuming that ρ s is given by Eqns (1) and (2) is very sensitive to the values of HS measured by the snow height gauge and the two-way travel time (τ) to the snow surface obtained by the picking algorithm. For instance, we will obtain an over- or underestimation of SWE of 30–40% assuming a deviation in snow height of 10%. If the snow surface picking is erroneous because of a deviation in the two-way travel time of 10%, SWE will deviate by 40–50%. For very low snow heights, uncertainties in the surface pick cause large errors in the calculated relative dielectric permittivity. To correct for these errors, we excluded all values of Ɛ r > 2 for periods when the total snow height was <30 cm and the snowpack was dry. For the rest of the season, all Ɛ > 10 were excluded.
As soon as liquid water is present in the snowpack, it is no longer possible to simply determine the bulk density ρ s with Eqns (1) and (2), as Eqn (2) is only valid for dry snow. Snow has to be considered as a three-phase mixing medium. The bulk density consists of the three contributing components ice (i), air (a) and water (w):
where θ i,θ a and θ w are the bulk volumetric ice, air and water contents (θ i+θ a + θw , = 1) and ρ i, ρ a and ρ w the respective densities (Reference Lundberg and ThunehedLundberg and Thunehed, 2000; Reference Heilig, Schneebeli and EisenHeilig and others, 2009). Since ρ d is the dry-snow density (ρ d =θ i ρ i + θ a ρ a), Eqn (8) reduces to
Owing to the assumption of a constant value for ρ d during the melting season (see above), we underestimate ρ d and therefore underestimate ρ s (though we overestimate ρ w in Eqn (5)). To account for the underestimation, we multiplied ρ w by a factor f:
The factor f was determined using the densities (ρ pit) measured in the snow pits. In order to obtain a robust value for f, we used a leave-two-out cross-validation technique. For each step, two measurements were removed from the training set and the factor f i was found for which the root-mean-square error (RMSE) between measured (ρ pit) and modeled (ρ s′) densities was smallest. For reason of plausibility, ρ s′ was limited to a maximal density of 600 kg m–3. f i was then tested against the two removed measurements and the corresponding errors (e 2i, e 2i+1) were calculated. The iteration was repeated until each possible combination of removing two measurements from the training set had occurred once. The average factor f for the two winter seasons was 3.08 with a standard deviation of σf = 0: 10. The rms of all errors was 16 kg m–3. Accordingly, the SWE for wet-snow conditions was then
Furthermore, we averaged all SWE values over a whole day to account for errors in the snow height measurement and the respective radar picking.
3. Results
3.1 Snow height
Figures 5 and 6 show the continuous radar sections obtained with the 1600 MHz antenna during the winter seasons 2010/11 (Fig. 5) and 2011/12 (Fig. 6). The modeled velocity profile from SNOWPACK was employed (velocity profile (ii) in Section 2.4.1). Red lines indicate the snow height recorded with a conventional snow height sensor (laser range gauge), and green lines display the snow height determined with the radar using the semi-automated picking algorithm. Between 16 December 2011 and 11 January 2012, there were many missing recordings of snow height. We replaced all missing values by measurements obtained with an ultrasonic sensor, which was at the same location as the laser. For dry-snow conditions (before 24 March 2011 and 27 April 2012), the RMSE of the snow height calculated with the radar compared with the laser gauge was about 3.1 cm in winter 2010/11 and 3.9 cm in winter 2011/12 (4.2 cm and 8.1 cm if a constant wave speed of 0.23 m ns–1 was assumed for the whole snow-pack). This is comparable to the accuracy of conventional snow height sensors (Reference Egli and JonasEgli and Jonas, 2009). In fact, the RMSEs between the laser range gauge and a nearby ultrasonic gauge were 3 cm and 11.6 cm for the two winter seasons. The snow height, however, was overestimated by the radar as soon as the snow became wet. The dielectric permittivity of snow increases as soon as liquid water is present; hence Eqn (2) is no longer valid and the snow height is overestimated.
Fig. 5. Processed upGPR data and measured snow height for winter 2010/11. Blue represents negative and yellow positive amplitudes of the reflected wave. For the conversion from two-way travel time to snow height, the wave speed of each layer is calculated from modeled snow densities (velocity profile (ii)). The green line shows the snow height determined with the radar, the red line represents the snow height recorded with a laser gauge, and black circles denote the snow height measured with an avalanche probe directly over the radar.
Fig. 6. Same as Figure 5 but for winter 2011/12. Black lines indicate the evolution of selected internal layers.
The semi-automated and the fully automated algorithm for picking the snow surface are compared in Figure 7. For the fully automated picking algorithm, the weather data from the AWS WAN7 were used. The RMSE of the snow height between the fully automated and the semi-automated algorithm was 9.4 cm for dry and 12 cm for wet snowpack conditions. For the entire 2011/12 winter, 93% of the fully automated picked snow surface values were within ± 10% of the semi-automatically derived snow surface heights (Fig. 7b); 49% of all heights did not show any deviation; when only dry snowpack conditions were considered, the agreement was 66%.
Fig. 7. (a) Comparison between fully automated (black dots) and semi-automated (green line) picking of the snow surface. The red line shows the measured snow height. For the fully automated picking algorithm, the measured snow height of WAN7 (magenta line) is used as input. (b) The red dots show the deviation between the fully automated and semi-automated algorithm (blue lines ± 10%, blue background, wet snowpack).
3.2 Settlement
Figure 8 shows the radar section for a period with a snowfall event in 2012. Just before the snowfall on 20 January 2012, a settlement sensor (SES) was placed on the old snow surface and followed the settlement of this snow layer. The corresponding reflection of the snow layer could be followed even after the snow surface was covered with new snow. When a constant wave speed of 0.23 m ns–1 was assumed for the whole snowpack (velocity profile (i)), the settlement measured by the radar was less pronounced than measured with the SES (Fig. 8a). This can be explained by the fact, that the density of the new snow is low and its wave speed is >0.23 m ns–1. The snow height was therefore underestimated. As settlement of the snowpack proceeded, the snow density increased and the wave speed decreased even below 0.23 m ns–1. The height of the layer was then overestimated. Therefore, the settlement measured by the radar appeared to be too slow. In order to solve this problem, we used the velocity profile (ii) based on densities modeled by SNOW-PACK. The surface reflection now agreed very well with the SES measurement (Fig. 8b). The radar section in Figure 6 was also calculated using the modeled densities. The reflections from old snow surfaces could be tracked until 27 April, when the snow cover became wet. The difference between measured heights of the internal layers and the corresponding reflections in the radar section never exceeded 7 cm for all six SES sensors.
Fig. 8. Processed radar data, measured snow height (red lines) and settlement of selected old snow surfaces (black lines). For the left radar section, the wave speed is assumed to be constant (0.23 m ns–1, velocity profile (i)), whereas for the right radar section the velocity for each layer is calculated from the density profile simulated by SNOWPACK (velocity profile (ii)). The blue arrows indicate the evolution of the old snow surface measured with the radar.
3.3 Snow density
Figure 9a shows a comparison between the bulk densities measured manually and the bulk densities modeled with SNOWPACK for winter 2011/12. The densities were in good agreement (± 1.8%) until 26 February 2012, but subsequently were underestimated by SNOWPACK. In addition, the bulk densities were calculated using the two-way travel time τ measured with the radar and the snow height HS measured using a nearby laser gauge. The green lines in Figure 9a and b show ρ bulk and v bulk. As long as the snowpack was dry (before27 April), the mean error compared with manually measured snow densities was 4.3%.
Fig. 9. (a) Bulk densities and (b) bulk velocities of the electromagnetic wave for the whole snowpack in winter 2011/12. The values are determined with the radar and a nearby snow height gauge (green lines), in manual snow profiles (black circles) and using SNOWPACK (blue lines). The red lines denote the values of wave speed (0.23 m ns–1) and density (357 kg m–3) suggested by Reference Heilig, Schneebeli and EisenHeilig and others (2009) for the conversion from two-way travel time to snow height. The blue background shows when the snowpack was wet.
3.4 New snow height
In winter 2011/12, the new snow amount for nine selected snowfall periods (Table 2) was estimated from radar data and compared with conventionally measured new snow heights. The 24 hour new snow accumulation was determined every morning at 09:00 and summed up for the whole snowfall period. Figure 10 shows the summed-up new snow heights for the nine snowfall periods as well as the snow height increase between two consecutive measurements as measured with the laser gauge. With one exception (period 3), all new snow heights determined by the radar were within ± 10 cm compared with the manually measured new snow amounts. As expected (Section 2.4.3), the snow height difference underestimates the new snow height due to the settlement of the underlying layers.
Table 2 Snowfall events with corresponding mean density (ρ mean) analysed in winter 2011/12
Fig. 10. New snow amount for nine snowfall periods in winter 2011/12. Blue bars show the new snow amount determined with the radar, red bars show the difference between the snow height after and before the snowfall determined with a laser range gauge, and green circles show the manually measured new snow amount. Whiskers show the measurement uncertainties for the whole snowfall period assuming daily errors of ±3 cm and ±7 cm for manual measurements and radar, respectively
3.5 Dry-to-wet transition and liquid water content
In winter 2011/12, meltwater started to percolate at the beginning of March so that the uppermost 50 cm of the snowpack became wet (Fig. 6). Strong multiple reflections indicate when liquid water is present in the snowpack. The wetting of the uppermost layers produced an ice layer as snow temperatures were far below 0ºC after this initial wetting. Once refrozen, this layer caused a strong reflection at ~ 1.85 m. Evidence from snow profiles and modeled snow temperatures underline this fact. The ice layer remained until 27 April, when water percolated through the snowpack down to the soil. When the liquid water content (θ w) increases, Ɛ r increases and thus the snow height is overestimated with our approach. Therefore, even diurnal variations in snow wetness are visible in Figure 6 as the overestimation of the snow height changed during the day and the snow surface appeared to oscillate (note the ripples after 27 April). Whereas in winter 2011/12 the whole snowpack became wet within 1 day, the snowpack became progressively wet from top to bottom over several days in winter 2010/11. Figure 11a shows the advance of the wetting front for the period when the snowpack became wet (28 March to 12 April 2011). The snow below the wetting front was assumed to be dry and the bulk liquid water content of the snow above the wetting front was calculated (Fig. 11b). If the dry-to-wet transition changes, the range in the snowpack for which θ w is calculated changes so that the values of θ w can no longer be compared with each other. In fact, the height of the dry-to-wet transition changed five times during the observed 15 days, as indicated by the gray vertical lines in Figure 11a and b.
Fig. 11. (a) Processed upGPR data and measured snow height in winter 2010/11. The green line shows the snow height determined with the radar, the red line shows the snow height recorded with a laser gauge, and the magenta line denotes the dry-to-wet transition. (b) Volumetric liquid water content calculated from the radar signal. The red crosses show the average liquid water content of the snowpack measured with the capacity probe; the red lines show the range of the liquid water content in all snow layers. (c) Water discharge recorded by a nearby 5 m2 lysimeter.
Therefore, the temporal evolution of θ w is only consistent between two of these lines. Nevertheless, Figure 11b shows that the water content increased during each day, reached its maximum between 16:00 and 17:00 and decreased again during the night. The highest values of water content were measured on 2, 3, 7 and 8 April. The water discharge measured by a 5 m2 lysimeter was largest (Fig. 11c) for these days as well. Furthermore, the bulk liquid water content was calculated for the whole snowpack and compared to the liquid water content measured with the capacity probe in Figure 11b.
3.6 Snow water equivalent
The temporal evolution of SWE during winter seasons 2010/11 and 2011/12 is shown in Figures 12 and 13. SWE was calculated according to Eqn (11) at the upGPR location, (1) derived from snow mass as measured by a snow pillow and converted to SWE, (2) derived from manually measured bulk snow density (ρ pit) times manually measured HS at the pit site or times laser-measured HS, and (3) modeled using SNOWPACK. Since observers were experienced, but observation techniques varied, we assumed an uncertainty of ± 5% for manually derived SWE. To account for spatial variability in snow height, we displayed results at the pit location and at the laser-gauge location. As manually determined SWE is calculated after Eqn (7), differences in snow height between the pit location and the laser site are documented by differences between red circles and black squares in Figures 12 and 13. Especially during winter 2010/11, measured snow heights at the pit locations were >5% larger than at the laser site over almost the entire winter.
Fig. 12. Comparison of upGPR-determined SWE (green) with weight-determined (blue), manually measured (red and black), and modeled (magenta) SWE for winter 2010/11. The blue background shows when the snowpack was wet. (a) The upGPR-determined SWE was averaged per day, weight-determined SWE values were measured utilizing a snow pillow, red circles correspond to SWE measured directly at the pit, and black squares represent SWE calculated from bulk density measured at the pit and snow height measured by the laser gauge. A 5% uncertainty in manual SWE measurements is assumed, and indicated by error bars. (b) Positive and negative deviations from the manual measurements (SWEpit,laser). The dashed line marks the 5% offset.
Fig. 13. Same as Figure 12, but for winter 2011/12.
The winter season 2010/11 was characterized by far-below-average snow height and accordingly below-average SWE at the test site. Furthermore, we faced some problems during early season until mid-January. We measured highly scattered snow heights at the laser site, although no new snow was observed. Therefore, we replaced all laser values until mid-January by measurements recorded with an ultrasonic sensor ~ 20 m away from the upGPR location. Next, the upGPR vertical lift had problems until 3 January 2011, so the snow-surface picking became more difficult and was less accurate for this time period. Owing to these measurement deficiencies, SWE determination for the first 30 days was not reliable; even daily averaged values strongly fluctuated. After 1 February 2011, the SWE derived from the radar signal was in good agreement with SWE calculated using laser-measured HS and, with one exception, was during early season until mid-January. We measured highly scattered snow heights at the laser site, although no new snow was observed. Therefore, we replaced all laser values until mid-January by measurements recorded with an ultrasonic sensor ~ 20 m away from the upGPR location. Next, the upGPR vertical lift had problems until 3 January 2011, so the snow-surface picking became more difficult and was less accurate for this time period. Owing to these measurement deficiencies, SWE determination for the first 30 days was not reliable; even daily averaged values strongly fluctuated. After 1 February 2011, the SWE derived from the radar signal was in good agreement with SWE calculated using laser-measured HS and, with one exception, was always within 5% uncertainty range. Before a significant amount of water was present in the snowpack (24 March 2011), both other methods for estimating SWE (pillow and SNOWPACK) agreed well with the SWE derived from the radar signal. The differences between weighed snow mass converetd to SWE, modeled SWE from meteorological parameters and radar-determined SWE never exceeded 5%. After the snowpack became wet, modeled SWE strongly increased and thereby overestimated measured values, even the directly converted pit measurements. Radar and pillow values did not increase as strongly. UpGPR-determined SWE matched manual measurements even when significant amounts of water were considered to be in the snowpack. The simple parameterization to compensate for offsets due to the assumption that density is constant seemed to work well for the rest of the season. If f +σ f or f –σ f is used in Eqn (10) instead of f, the effect on SWE cannot be displayed in Figure 12 since the deviations are smaller than the line thickness. Weighed snow mass overestimated manually measured SWE calculated with HS from the laser gauge, but matched SWE calculated with HS from the pit. Modeled SWE decreased rapidly after the strong peak in April 2011. Compared with manually determined SWE, RMSEs, peak-SWE offsets (which might not be the seasonal peak SWE as no continuous daily manual observation was performed) and coefficients of determination (R 2) for the snow pillow and upGPR data were very similar (Table 3). The RMSE was twice as large for modeled SWE as for upGPR-derived SWE.
Table 3 Root-mean-square error (RMSE) and coefficient of determination (R 2) of SWE derived from upGPR and snow pillow or modeled with SNOWPACK when compared with SWE derived from measurements of snow height and density in the snow pit. N is number of samples. In addition, the difference to the manually measured maximum SWE is shown (peakmanual). Maximum SWE value determined manually was 585 mm for 2010/11 and 1125 mm for 2011/12.
The winter season 2011/12 was characterized by above-average snow height and accordingly above-average SWE at the test site. During December 2011 and January 2012 it snowed repeatedly and the snow height increased (Fig. 6). For that winter, no significant snow height offset between laser and pit location was recorded. Again for very low snow heights (until mid-December), the radar-based approach had problems in deriving reliable values of SWE. The radar surface pick fluctuated considerably for this period (Fig. 6) as the reflection of the radar box caused interference with the low-lying reflection of the snow surface. Additionally, in late May (when the entire snowpack was wet), upGPR-determined SWE decreased significantly faster than manual recordings in density and HS. Other than that, from mid-January until early March, all three SWE determinations were close to or within the ± 5% accuracy range of manual measurements. From mid-March on, modeled SWE significantly underestimated manual measurements. Right after the snowpack was considered to be wet, SNOWPACK-modeled SWE agreed well (± 5%), but after mid-May values were significantly lower than those obtained with the other methods. Deviations between weighed snow mass and radar-determined SWE were small except for conditions when liquid water was present in the snowpack. When compared with manually derived SWE, the RMSE for the upGPR-derived SWE was slightly lower than for SWE derived from the snow pillow, and again SNOWPACK resulted in ~ 20 mm larger offset. Again, the standard deviation of f in Eqn (10) cannot be displayed in Figure 13 as its effect on SWE is smaller than the line thickness. Peak SWE was distinctly underestimated by all three methods (90 mm by modeling and weighing, 76 mm by radar; Table 3).
4. Discussion
4.1 Error analysis
All radar-derived snowpack parameters depend on a correct pick of the reflection at the snow surface or internal layers and/or no difference in snow height between the location of the radar and the snow height sensor. From these requirements, potential sources of errors arise:
-
1. picking the second negative half-cycle instead of the first one
-
2. uncertainties in pick locations within a half-cycle
-
3. resolution limitations due to the wavelength
-
4. spatial variations in snow height within the test site.
-
1 If the second negative half-cycle (d in Fig. 2) is picked instead of the first one (c in Fig. 2), radar-derived snowpack properties will be either over- or underestimated, depending on the property. For example, for a mean density of 360 kg m–3, picking the second half-cycle would result in an overestimation of density of at least 8% for the range of our snow heights. For snow heights lower than 60 cm, the radar-derived density would be unrealistically high (>500 kg m–3). Since our results (Fig. 9) indicate much smaller deviations (4.3%, Section 3.3), we are confident that our pick is not systemically on the wrong half-cycle. For low snow heights, before 1 January 2012, isolated outliers cannot be excluded.
-
2 Both the semi-automated and the fully automated algorithm automatically correct the pick location to the minimum of the negative half-cycle. For our dataset, surface signals were never clipped owing to too strong amplitudes. Consequently, a clear minimum was always detectable. Uncertainties due to sample resolution are in the range of maximum 0.04 ns and thus negligible.
-
3 Two consecutive layer transitions within the resolution limit of the radar (± 7 cm) cannot be separated. Consequently, small amounts of new snow on top of a high-permittivity layer cannot be distinguished from the old snow surface. For specific time periods during the season, the resolution limit produces significant errors. However, since new snow and settlement of such thin layers are temporary, no error bars for the radar-determined snowpack parameters are presented.
-
4 As the snow height sensor is located at a distance of ~ 5 m from the radar location, differences in total snow height might occur. The RMSE between probed snow heights at the location of the radar and laser-measured snow heights was 5.8 cm for both seasons. Such an uncertainty in snow height induces, for example, a median variation in density of ± 10%. However, deviations of radar-determined density and manually measured density in the snow pits (Fig. 9) were much smaller.
This suggests a less pronounced small-scale variability in snow height (± 4.3%, Section 3.3).
4.2 Snow height
4.2.1. Semi-automated picking algorithm
As long as the snowpack was dry and no velocity profile was used to adjust the wave speed, radar-derived snow height was in good agreement with measured snow height (Figs 5 and 6). The deviations were between 3 cm and 8 cm depending on method and winter season. The largest deviation from conventionally measured snow height occurred between 26 February and 27 April 2012. This deviation cannot be explained by spatial variations in snow height, since the manually measured snow heights directly over the radar matched the snow height measured with the laser gauge (black circle in Fig. 6). The deviation was caused by the large snow density, which was larger than the usually assumed bulk density of 357 kg m–3 (Fig. 9a). Thus the wave speed of 0.23 m ns–1 was too high and snow height was consequently overestimated. Providing model-based velocity profiles (ii) did not significantly improve the accuracy of the radar-determined snow height, since SNOWPACK also underestimated bulk density (Fig. 9a).
4.2.2. Fully automated picking algorithm
In order to provide near real-time snow properties from radar signals to local avalanche forecasters, a fully automatic picking algorithm is necessary. Our approach towards such an algorithm performed well (Fig. 7). Large deviations (i.e. > ± 10%) from the results of the semi-automated picking algorithm occurred at the beginning and end of the season. At the beginning of the season low snow heights prevailed and thus small absolute differences caused high relative deviations. Towards the end of the season, differences were due to the fact that the station WAN7 was free of snow. Therefore, the required input values TSS and HS were not representative for Weissfluhjoch. In general, the automated picking had some problems during snowfall periods. For these cases, the reflection originating from the snow-to-air transition is less pronounced, since the difference in dielectric permittivity becomes small. This was particularly the case for the snowfall event in late January 2012, when the density of new snow was low (75 kg m–3), thereby causing almost no reflections. As a consequence, the algorithm followed the reflection of the previous snow surface. The error prevailed until the snow surface again produced a more pronounced reflection signal (8 February).
In summary, we demonstrated that it is feasible to automatically detect the snow surface by combining the radar signal with meteorological information from an AWS. The station does not need to be close but can be several kilometers away. This fact is important and very promising in view of a future installation of our radar system on avalanche slopes. Without any additional local information, the snow height is reproducible almost within the accuracy of conventional snow height sensors (Reference Egli and JonasEgli and Jonas, 2009). Since the radar is buried, it is not prone to avalanches and snow gliding and might be deployed in potential release areas; it could provide reliable information from where it is most needed.
4.3 Settlement
Detecting internal properties of the snowpack, such as the evolution of old snow surfaces, allows quantifying settlement and new snow height. Settlement of various snow layers was tracked and, compared with an alternative, independent measurement method (SES), the height of internal layers could be determined with an accuracy of ± 7 cm. Until now, settlement of distinct snow layers has had to be modeled based on data from an AWS (Reference Lehning, Fierz, Brown and JamiesonLehning and others, 2004). Measuring of settlement was very difficult (Reference Fierz and LehningFierz and Lehning, 2001) and measurement set-ups were only possible in flat terrain. In addition, the set-up presented in Reference Fierz and LehningFierz and Lehning (2001) was restricted to those snow layers on which a sensor had previously been placed. Knowing the settlement of almost every snow layer now allows the new snow height over a certain period to be calculated.
4.4 Snow density
If an external snow height measurement is available in parallel, the bulk snow density of the dry snowpack can be calculated continuously. The advantage of the upGPR is that it is not biased by bridging effects, as can be the case for measuring systems that weigh the mass of the snowpack above a sensor system (e.g. snow pillow). If the snow height and the bulk density are known, these two parameters can be converted easily to SWE (under dry snow conditions).
4.5 New snow height
The amount of new snow during a snowfall period is a key contributing factor for assessing avalanche danger (Reference Schweizer, Mitterer and StoffelSchweizer and others, 2009), and our results suggest that the radar will provide this property with reasonable accuracy. For the results presented here, there is no clear tendency in regard to the accuracy of our radar-derived new snow height when compared to manual observations (Fig. 10). Over- and underestimations of new snow height are independent of the assumed new snow density of 100 kg m–3 (Table 2). Uncertainties depend on the precision with which the old snow surface can be picked. The exact pick is crucial for a correct estimate of new snow height. Since we can pick the old snow surface with an accuracy of ± 7 cm, we obviously may obtain erroneous values for new snow height. However, we are within the range of our possible resolution, so reducing the error is not possible with our present set-up. Still, the agreement with the observed values was fairly good (Fig. 10).
4.6 Dry-to-wet transition and liquid water content
For a wet snowpack, we monitored the advance of the dryto-wet transition and the bulk liquid water content for the wet part of the snow cover. If the volumetric liquid water content (θ w) constitutes up to 3–4% (Fig. 11b), the water is held by capillary forces; as soon as the water content increases, gravitational forces will govern water movement and water will flow downwards. Figure 11c shows, indeed, that whenever the bulk liquid water content above the dryto-wet transition clearly exceeded 4%, outflow was observed in the lysimeter. This observation agrees fairly well with measurements of residual water content (Reference Coléou and LesaffreColéou and Lesaffre, 1998). The results demonstrate that the radar is capable of calculating the transition from dry to wet snow within the snowpack. Knowing this transition is important for wet-snow avalanche forecasting (Reference Mitterer, Hirashima and SchweizerMitterer and others, 2011b). However, the values of the liquid water content above the dry-to-wet transition in Figure 11b seem unrealistically high. One reason for this is that the thickness of the wet layer (in most cases <20 cm) was hardly larger than the resolution limit of the radar (7 cm). Another source of error might be the difference in snow height between the location of the snow height sensor and the radar. An offset in total snow height of 1–2 cm has a large effect for such thin layers. For these reasons, we cannot reliably measure the absolute values of the liquid water content for such thin layers with the radar. The values derived for the whole snowpack were, however, in a similar range to that measured with the capacity probe (Fig. 11b). The discrepancy between the radar and the capacity probe may be explained by sampling effects, for example caused by an opened snow-pit wall (Reference Mitterer, Heilig, Schweizer and EisenMitterer and others, 2011a).
4.7 Snow water equivalent
The two observed winter seasons were not representative in terms of long-term average SWE. SWE in winter 2010/11 was distinctly below average over the entire season, while during winter 2011/12, SWE was remarkably above average at our test site. However, for both years, upGPR-determined SWE was in very good agreement with snow-pillow results and matched manual measurements within about ± 5% over most of the season. Even for wet snow conditions, the simple parameterization on θ w produced reliable values of the bulk density (Figs 12 and 13). With very few exceptions the density approximation seemed to be reliable. Note that the factor f is not independent of our manually measured SWE values. Despite the fact that the two years differed considerably in terms of SWE, for both years we used the same factor f and obtained similarly reasonable results (Figs 12b and 13b). The factor f seems to be robust against variable amounts of SWE. However, whether the factor would be different for another location or whether its value is generally valid cannot be determined. Performing a cross-validation revealed a mean error of 16 kg m–3. This corresponds, for both years, to an average error of 26 mm w.e.
Especially for shallow snow (early in winter 2010/11), large variations in HS existed, even at this flat test site. Accordingly SWE determined with manually measured density varied for the different locations of snow height measurements. Therefore, whether upGPR-determined SWE is always as accurate as, or more accurate than, weighing the snow masses (e.g. by pillow technology) cannot be determined with reasonable certainty. It cannot be excluded that between the laser and pillow locations a similar difference in snow height existed to that between the laser and pit locations. With distinctly more snow in 2011/12, this spatial variability in snow height was less pronounced. So far, modeled SWE applying SNOWPACK deviated significantly more from manual measurements than weighed or radar-determined SWE. For the two years, SNOWPACK either overestimated (2010/11: below-average year) or underestimated (2011/12: above-average year) measured SWE. For both years, SWE modeled by SNOWPACK deviated by >5% compared with manual measurements.
Applying upGPR combined with a conventional snow height gauge cannot be considered as being more accurate than existing methods such as weighing the snow mass, but enables new possibilities for locations on slopes or more complex terrain. Furthermore, measurements are not biased through bridging effects.
4.8 Potential of upGPR
Although we were able to measure several properties of the snowpack with the present set-up, we still needed external information on the snowpack such as an additional snow height measurement or a modeled density profile. We hence have not yet fully reached our goal of presenting a device that provides independent and relevant information from an avalanche path, since these external data are typically not available (with sufficient accuracy) for the location of an avalanche path. In addition, we can neither verify nor possibly improve layer evolution modeled by SNOWPACK as we still used data provided by the model. The main problem is that we have so far only used two-way travel time () to measure height (D) and density of the layers. As the wave speed in snow (v bulk) depends on snow density, we did not have enough known variables to solve Eqn (1) and had either to make an assumption or get the information from an external source. There are several options to overcome these deficiencies and obtain the necessary information. The two most promising ones are the following. (1) Instead of just placing one transmitting and one receiving antenna beneath the snow, multiple receiving antennas can be placed at different distances from the transmitting antenna. Hence, basically an upside-down wide-angle reflection and refraction (WARR) measurement is carried out. By comparing the change in antenna separation with the time taken for the pulse to return from the reflectors, a value for the velocity for each layer can be calculated. (2) Instead of only evaluating the two-way travel time, the amplitudes of the signal also contain information about the amount of change in dielectric permittivity between two adjacent layers. A jump in dielectric permittivity represents a jump in velocity, which in turn represents a density jump. The velocity profile of the snowpack might be calculated by full-waveform inversion (FWI) (Reference TarantolaTarantola, 1984; Reference Pratt, Shin and HicksPratt and others, 1998; Reference PrattPratt, 1999; Reference SymesSymes, 2008). FWI is an automated method for refining a model by iteratively attempting to match modeled data with recorded data.
For periods when the snowpack is wet, we would still lack information since the wave speed depends in addition on the liquid water content of snow. However, the wetness can be assessed using an expected shift of the frequency content (Reference Quan and HarrisQuan and Harris, 1997; Reference Bradford, Harper and BrownBradford and others, 2009).
Combining the measurements presented in this paper with FWI and frequency analysis might provide quantitative information on snow height, new snow height, snow settlement, liquid water content and SWE from radar signal signatures only.
5. Conclusions
We employed an upward-looking, ground-penetrating radar system (upGPR) to continuously observe and quantify properties of the snowpack during winter seasons 2010/11 and 2011/12. Upward-looking radar sensors, as a nondestructive monitoring technique, enable observation of the temporal evolution of various snowpack parameters over the course of a winter season. Conventional destructive methods (e.g. snow pit) offer only a snapshot in time. As long as the snowpack is dry, snow height can be determined within the accuracy of conventional sensors. We developed an algorithm for picking the snow surface such that snow height can be determined without the need for manual intervention during signal analysis. The radar also offers the unique possibility of following the evolution of internal snow layers. We monitored settling rates of single layers and tracked the reflection of old snow surfaces for up to 4 months after their initial burial. To verify the settling signature within the radar signal, we deployed lightweight sensors attached to potentiometers on the snow surface before every significant snowfall event. These sensors settled with the old snow surface. Settlement curves derived from the radar and measured directly were in very good agreement. As the old snow surface was visible in the radar signal, we were able to calculate the amount of new snow during snowfall events by subtracting the height of the reflection of the old snow surface from the height of the new snow surface. Manual measurements of the new snow amount were generally in fairly good agreement with the values calculated with our radar algorithm.
Using external snow height measurements, bulk density and SWE of the entire snowpack can continuously be determined. In addition, for robust SWE calculations, a training dataset containing manual density measurements collected close to the radar is necessary. Since SWE is not derived by weighing the mass of the snowpack, it is not influenced by bridging effects. Furthermore, measuring SWE with the radar in a slope seems feasible. For wet-snow conditions, we observed the occurrence of strong multiple reflections as well as the daily increase in two-way travel time of reflection horizons. This allows the depth of the transition from dry to wet snow to be determined, as well as the absolute amount of liquid water content for the whole snowpack, the timing of the daily peak in volumetric liquid water content and its decrease due to refreezing during the night.
We demonstrated that it is feasible to quantitatively derive snowpack properties relevant to avalanche formation and monitor their evolution in time. However, in order to derive some of these properties, we still needed additional information such as independently measured snow height or modeled snow density. Hence, the system is not yet able to provide information from avalanche starting zones, since the collection of such data cannot be done in avalanche-prone terrain. In the future, we will therefore evaluate additional techniques (WARR) and apply more sophisticated evaluation methods by combining the algorithms presented in this paper with full-waveform inversion and frequency analysis.
Acknowledgements
We greatly appreciate comments from Edward Adams. Careful and constructive anonymous reviews greatly improved the manuscript. The research project was completed within the framework of the D-A-CH consortium (DFG–FWF–SNF) with the German Research Foundation (DFG) as lead agency. C.M. was funded by the Swiss National Science Foundation (SNF); A.H. was supported by DFG grant EI 672/6; and O.E. by the DFG Emmy Noether programme grant EI 672/5.
Open access
