An assessment of sub-snow GPS for quantification of snow water equivalent

Global Navigation Satellite Systems (GNSS) contribute to various Earth observation applications. The present study investigates the potential and limitations of the Global Positioning System (GPS) to estimate in situ water equivalents of the snow cover (snow water equivalent, SWE) by using buried GPS antennas. GPS-derived SWE is estimated over three seasons (2015/16–2017/18) at a high Alpine test site in Switzerland. Results are validated against state-of-the-art reference sensors: snow scale, snow pillow, and manual observations. SWE is estimated with a high correspondence to the reference sensors for all three seasons. Results agree with a median relative bias below 10 % and are highly correlated to the mean of the three reference sensors. The sensitivity of the SWE quantification is assessed for different GPS ambiguity resolution techniques, as the results strongly depend on the GPS processing.

1 Introduction

Knowledge of snow cover characteristics is an important basis for climatology, natural hazard forecasting, early-warning systems, and hydro-energy industries. An extensive amount of water stored in snow cover has a high impact on flood development during snow melting periods. High damage is caused worldwide by floods originating from mountain catchment areas. Early assessment of the snow water equivalent (SWE) (Fierz et al., 2009) in mountain environments enhances early warning and thus prevention of major flood events.

Several point-wise measurement methods already exist to continuously determine SWE. SWE is measured usually in situ with manual or automated observation techniques and is expressed in units of mass per area (kg m⁻²) or in millimeters of water equivalent (mm w.e.). Using SWE tubes, a sample is taken out of the snow profile and weighted afterwards leading to SWE. Furthermore, SWE is calculated indirectly based on snow depth and the bulk snow density, measured manually in a snow pit or along a transect (Sturm et al., 2010; WMO, 2008). Both techniques are state of the art and considered to be the most reliable at the moment. They are, however, labor intense, time-consuming, and destructive and have a low temporal resolution. Automated and continuous SWE measurements are provided by a snow pillow or a snow scale (Beaumont, 1965, 1966; Johnson et al., 2007, 2015) and are, however, prone to errors, especially during snow melt events. Additionally, cosmic ray neutron probes and other passive or acoustic instruments indirectly measure SWE (Harding, 1986; Kinar and Pomeroy, 2007, 2015a; Kodama et al., 1979; Rasmussen et al., 2012). Smith et al. (2017) and Kinar and Pomeroy (2015b) provide a detailed summary as well as a comprehensive description of terrestrial SWE measurement techniques. SWE observations based on satellite remote sensing are limited to large plains due to low spatial resolution and problems with steep orography of mountain chains like the Alps. However, accurate and reliable in situ data are still needed for calibration and validation of remote-sensing data (Derksen et al., 2005; Goodison and Walker, 1995; Takala et al., 1995).

Global Navigation Satellite Systems (GNSS) remote-sensing techniques are capable of providing reliable, accurate, efficient, and continuous observations independent of weather conditions. Sub-snow GNSS techniques have been tested lately for SWE estimation (Henkel et al., 2018; Limpach et al., 2013), which is suggested to determine liquid water content (Koch et al., 2014) or considered for avalanche rescue (Claypool, 1997; Olmedo et al., 2012; Schleppe and Lachapelle, 2008). Most studies concentrate on the use of Global Positioning System (GPS) single-frequency signals (L₁, with a frequency of f₁=1575.42 MHz). GPS antennas buried under snow were first used to evaluate the potential of GPS signal reception and positioning performance in avalanche rescue research (Claypool, 1997). Schleppe and Lachapelle (2008) experimentally analyzed the GPS tracking performance under avalanche-deposited snow at two test sites in Canada. The potential of a GPS-based rescue system for victims buried under avalanches was investigated by Olmedo et al. (2012) in the framework of the SICRA project. Steiner et al. (2018) analyzed the GPS receiver behavior from antennas submerged in water and developed a model to estimate the water depth above the submerged GPS antenna based on the path delay of the GPS L₁ signals. Another study deals with the analysis of GPS L₁ signal-to-noise ratio (SNR) for snow property estimation from a ground GPS antenna compared to a GPS reference antenna above the snow surface (Appel et al., 2014). Koch et al. (2014) continuously estimated the liquid water content (wetness) based on GPS L₁ signal strength observations for a seasonal snowpack in the Swiss Alps. Henkel et al. (2018) show exemplarily the possibility of SWE estimation using GPS L₁ carrier phase residuals for the dry snow period during winter 2015/16 at the Weissfluhjoch test site of the WSL Institute for Snow and Avalanche Research (SLF).

Furthermore, GNSS stations above a snowpack are used for SWE estimation based on the interference of direct and reflected GNSS signals (GNSS reflectometry). Thereby, the SWE is estimated using GPS L₁ C/A or L₂c SNR (Jacobson, 2010, 2012; McCreight and Small, 2014). Different GNSS reflectometry methods are used for snowpack characterization (Boniface et al., 2015; Larson et al., 2009; Najibi and Jin, 2013; Ozeki and Heki, 2012; Tabibi et al., 2017; Vey et al., 2016), mainly for snow depth estimation. Cardellach et al. (2012) investigated the potential to remotely sense sub surface snow structures in dry snow areas using bi-statically reflected GNSS signals.

This paper presents a method to estimate SWE based on phase-based differential GPS by using an antenna buried in the snowpack (sub-snow GPS) and a reference station above the snowpack. A model is developed to estimate SWE based on the refraction and path delay of GPS code and phase signals while propagating through a snowpack. A sensitivity analysis on SWE estimation results is carried out, analyzing different ambiguity resolution techniques within the differential GPS processing. Results are validated against the state-of-the-art reference sensors, the snow pillow, snow scale, and manual observations for three full winter seasons at the high alpine test site Weissfluhjoch of the SLF.

Section 2 describes the study site, the geodetic GPS equipment, and the reference sensors. Section 3 derives the model to estimate the SWE from the GPS signal refraction and propagation delay. The SWE of a seasonal snowpack is estimated using buried GPS antennas and the method of investigation is described in Sect. 4. The developed model is applied to a seasonal snowpack over three winter seasons and results are analyzed in Sect. 5. Section 6 discusses the results and Sect. 7 concludes the paper.

2 Study site and instrumental setup

The GPS snow monitoring system was installed during a snow-free period at the SLF test site Weissfluhjoch. The network includes geodetic and low-cost GPS stations, operating permanently since October 2012. The present study focuses on the geodetic station in order to understand the receiver behavior. A follow-up study will investigate the potential of the low-cost system.

2.1 Weissfluhjoch test site

The Weissfluhjoch test site (Fig. 1) is located in the Swiss Alps above Davos at 2536 m a.s.l. The 40×40 m flat test site is obstructed at elevations below 20 ^∘ by mountain faces in the north (N) and west (W) and open to the east (E) and south (S). The open E and S direction and sky visibility above 20 ^∘ in the N and W enhances the GPS satellite visibility and thus the observation geometry. The main winter wind direction for the site comes from the NW and SE (MeteoSwiss, 2014). The test site is equipped with energy supply, internet connection for automated data transmission, and reference data for validation of the GPS snow monitoring system (Fig. 2).

Figure 1 Horizon of the Weissfluhjoch test site (MeteoSwiss, 2014).

Figure 2 Weissfluhjoch test site, equipped with the GPS snow monitoring system, snow scale, and snow pillow. Manual transects are carried out biweekly at the north side as shown by the dashed line.

2.2 GPS snow monitoring system

The GPS snow monitoring system is set up at the Weissfluhjoch test site and consists of a geodetic GNSS system. However, as the present study only deals with GPS signals (for future comparison to low-cost systems), the term GPS is used further for all equipment. The use of a geodetic system allows a better understanding of all snow effects on the observations by minimizing additional effects. A GPS antenna (Leica AS10) is mounted at 5.3 m in height on a pole and serves as reference station. A second GPS antenna (Leica AS10) is installed on the ground next to the reference station (Fig. 2) and is referred to as the sub-snow GPS station. This setup results in a very short baseline for differential GPS processing. The corresponding GPS receivers (Leica GR10) are sheltered in the hut next to the test site, allowing power supply, internet connectivity for remote access, and anytime access in case of receiver problems. Dedicated low-loss GPS cables (Huber+Suhner SPUMA 400) are used to mitigate signal loss caused by the long path below the snow surface of approximately 40 m. The use of the same equipment for both stations (antenna, receiver, antenna cable type, and length) ensures the best possible consistency for the differential GPS processing. As visible in Figs. 2 and 1, a lot of metal surrounding the GPS antennas is present at the test site, generating multipath effects in the GPS signals.

In this setup, the sub-snow GPS antenna is placed on the ground and is buried by snow after the first snowfall already. A GPS antenna is, however, designed for an usual environment of air. The antenna impedance matches the impedance of air in order to avoid refraction effects at the antenna–air boundary. When changing the environment by placing the antennas in snow, an antenna impedance mismatch could occur, influencing the tracking performance (Rao et al., 2012).

As the GPS system was installed in the snow-free period, these observations have been used to determine best possible reference coordinates for the sub-snow GPS station (Sect. 4). Later, SWE estimation results from the GPS monitoring system are validated against state-of-the-art reference data operated by the SLF.

2.3 Reference data

SWE reference data (Marty, 2017) are provided by the SLF from the test site Weissfluhjoch (Fig. 2). A snow pillow (Sommer Messtechnik SP3) and a snow scale (Sommer Messtechnik SSG1000) serve as reference data for SWE estimation. The snow pillow measures the overlying pressure of the snowpack on a fluid-filled bladder. The snow scale uses a weighing surface and load cells to measure the weight of the overlying snow (Beaumont, 1966; Johnson et al., 2007). The sensors are located on the ground, 20 m next to the sub-snow GPS station (Fig. 2). Both sensors are offset corrected to zero each autumn before the first snowfall and acquire data in a 30 min time interval. Additionally, biweekly manual SWE observations are available from snow profiles. Here the water equivalent of the snow cover is obtained from multiple seamless vertical coring using an aluminum cylinder with a cross-sectional area of 70 cm² and a length of approximately 55 cm. The profiles are dug along transects north of the sub-snow GPS station and thus their relative position changes over time. Furthermore, the sum of the measured water equivalents of snowfall is used for comparison purposes. As the latter observations match the manual observations well over the three processed winter seasons, they are not shown in all figures of the present paper for clarity reasons.

The GPS-derived SWE is compared to the snow pillow, snow scale, manual observations, and a combined reference (daily average of manual, snow pillow, and snow scale measurements). The combined reference is used as the three reference sensors are not consistent with each other over the three seasons, which is discussed in Sect. 6.

3 SWE estimation model

A model is developed according to Steiner et al. (2018) to estimate the SWE above the sub-snow GPS antenna. The model is based on the path delay of the GPS signals while propagating through the snow cover. Refraction at the air–snow interface and deceleration of the signal propagation velocity strongly depend on the amount and wetness of snow the signal has to travel through. Similar to atmospheric refraction, the GPS signals received by an antenna under a snowpack are delayed by these refraction effects and lead to a longer electrical path through the snowpack (excess path length).

The excess path length δL_s depends on the bulk refractive index n_s of the present snow, the incident angle ν_a, and the snow depth d:

$$\begin{array}{}\text{(1)}& \mathit{\delta }{L}_{\mathrm{s}}=d\cdot \left(\sqrt{n_{\mathrm{s}}^{\mathrm{2}}-{\sin }^{\mathrm{2}}{\mathit{\nu }}_{\mathrm{a}}}-\cos {\mathit{\nu }}_{\mathrm{a}}\right).\end{array}$$

The derived model is, of course, a simplification as the snowpack consists not of one but of several layers with differing path delays. The model would allow us to estimate the snow depth d if the bulk refractive index n_s of the snowpack is known. Using the attenuation of the GPS signal strength would provide the snow wetness (Koch et al., 2014), which could be converted to the refractive index n_s. However, the bulk density of snow ρ_s and the density of water ρ_w is still needed in order to derive the SWE:

$$\begin{array}{}\text{(2)}& \mathrm{SWE}=d\cdot {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{w}}}}.\end{array}$$

Because the bulk refractive index n_s and the snow density ρ_s are not known, the present study aims to use only the GPS path delays to estimate the SWE directly. For that reason, an assumption is made: the GPS path delay in the layered snowpack is assumed to be nearly independent of snow density and the distribution of liquid water. Water has the highest impact on GPS observations within a snowpack (Steiner et al., 2018). Neglecting snow density and the water distribution, the model simplifies to a single water layer (Fig. 3) with the refractive index of water n_w (Hofmann-Wellenhof et al., 2001). The depth d of the water layer thereby corresponds directly to the SWE (mm w.e.):

$$\begin{array}{}\text{(3)}& \mathit{\delta }{L}_{\mathrm{s}}=\mathrm{SWE}\cdot \left(\sqrt{n_{\mathrm{w}}^{\mathrm{2}}-{\sin }^{\mathrm{2}}{\mathit{\nu }}_{\mathrm{a}}}-\cos {\mathit{\nu }}_{\mathrm{a}}\right)=\mathrm{SWE}\cdot F\left({\mathit{\nu }}_{\mathrm{a}}\right),\end{array}$$

where F(ν_a) can be interpreted as a mapping function to estimate the SWE of a snowpack (Sect. 5.1).

Figure 3 Geometry of signal paths in air (l_a) if no snow is present and in a snowpack (l_s) of depth d. d corresponds to the SWE in the case of the single water layer model assumption.

4 Method of investigation

The present study focuses on the use of GPS L₁ data to allow a comparison to single-frequency low-cost systems in a future study. The GPS L₁ data of the sub-snow station are processed differentially to the base station data. The same equipment is used on both stations, eliminating the impact of antenna phase center offsets and variations. Atmospheric influences are mitigated by the very short baseline of approximately 5 m.

The present study assumes a flat surface at the air-to-snow interface. This is an idealization as the roughness of the snow surface caused by environmental influences such as snow redistribution due to wind would change the signal path delays. However, the roughness effects at the Weissfluhjoch test site are considered to be of a small order and assumed to average out.

All data are processed over 24 h (daily solutions) using the Bernese GNSS software (Dach et al., 2015) with sampling data at 30 s. Observations from the snow-free periods serve as GPS reference measurements without snow effects. Precise coordinates of the sub-snow GPS station are computed as a combined solution of 7 snow-free reference days. The double difference processing is carried out for the very short baseline between the geodetic base and sub-snow station. The GPS L₁ processing utilized ionospheric, clock, and precise orbit products, provided by the International GNSS Service (Dow et al., 2009).

For the days when snow fell on the antenna and the snow depth increased over the winter, the SWE is estimated as a daily solution using the derived model (Eq. 3). The corresponding time stamp is thereby set to noon as it can be interpreted as a daily average. The comparison to the automatic reference sensors (snow pillow and scale) is made on the daily average. Manual observations are usually carried out around 09:00 UTC+1 (CET) in the morning. The time shift in between the manual and GPS observations is neglected in the present study.

4.1 Snow water equivalent estimation

The error induced by the excess path length δL_s (Eq. 3) when the signal propagates through a snow layer above the antenna is introduced in the zero difference phase observation equation as an additional parameter δL_s:

$$\begin{array}{}\text{(4)}& L=\mathit{\rho }+\mathit{\delta }\mathit{\rho }+\mathit{\delta }{L}_{\mathrm{s}}+\mathit{\lambda }N+\mathit{ϵ}.\end{array}$$

Thereby, L is the observed path length and ρ is the range between the sub-snow antenna and a GPS satellite. All known path delays from the sub-snow antenna to the GPS satellites (e.g., ionospheric and tropospheric delays) are included in δρ. Hofmann-Wellenhof et al. (2001) thoroughly describe all known GPS path delays. The unknown number of ambiguities N is contained in λN, with the GPS wavelength λ. The measurement noise is expressed as ϵ.

The developed SWE estimation model was implemented directly into the Bernese GNSS software. The SWE estimation is expected to be highly correlated with other processing parameters such as ambiguities, station height, clocks, and troposphere parameters. The direct implementation permitted a processing, which takes these correlations into account by simultaneously resolving ambiguities and estimating the SWE in one processing step.

By forming the single difference to the base station, the satellite clock and the atmospheric and relativistic effects cancel out at short baselines, in our case approximately 5 m. As the snow only affects the sub-snow GPS station, the error introduced by the snow remains, together with the receiver clock errors, ambiguities, and multipath effects. The receiver clock errors can be eliminated by double difference processing. Multipath effects are not eliminated by differential processing and are not modeled. Therefore, multipath effects remain in the GPS observation residuals and are strongest for low-elevation signals. Due to the GPS constellation repeatability, multipath effects impact the GPS signals analogously every sidereal day. The SWE parameter can be estimated, if the coordinates of the base and sub-snow GPS stations are precisely known and fixed. The daily solution processing provides one SWE estimate per day.

4.2 Sensitivity on GPS ambiguity resolution

The GPS phase observation equation includes the unknown ambiguity term N (Eq. 4). The ambiguity is an integer number by definition as it refers to the counted number of full wave cycles between the receiver and a satellite. However, the estimated ambiguity term is a float number due to instrument biases of the satellite and the receiver. The ambiguity term might be resolved to its integer value during the GPS processing. This ambiguity resolution leads to an accuracy increase, especially for small observation periods below 1 h. The reduction of parameters in the GPS processing enhances the SWE parameter estimation. Successful ambiguity resolution depends on the satellite geometry (high number of observations, long observation periods), the baseline length (ionospheric and tropospheric refraction as well as orbit biases), and the multipath effects (Hofmann-Wellenhof et al., 2001).

Not modeled error sources, such as the snowpack above the GPS antenna, degrade the ambiguity parameters. The subsequent systematic bias in the float ambiguity parameters leads to false resolved ambiguities and, if fixed for further processing, false estimated SWE. Due to the systematic bias in the float ambiguities, the ambiguities cannot just be rounded to the nearest integer. Sophisticated ambiguity resolution methods are used for the processing. Different GPS ambiguity resolution strategies are therefore investigated in order to assess the effect of ambiguity resolution on SWE estimation.

• L₁ float includes no ambiguity resolution. We estimate SWE with float ambiguities.

• L₁ fixed resolves L₁ ambiguities in a first step. We introduce the resolved ambiguities and estimate SWE in a second step (Henkel et al., 2018; Limpach et al., 2013).

• L₁ SWE fixed resolves L₁ ambiguities and estimates SWE in one step.

• L₅ fixed resolves wide-lane (L₅) ambiguities in a first step. We introduce the resolved ambiguities and estimate SWE in a second step. These results are not plotted for visibility reasons; however, they match the behavior of the L₁ fixed solutions.

• L₅ SWE fixed resolves L₅ ambiguities and estimates SWE in one step.

Chosen names of the different strategies are marked in bold, L₁ corresponds to the single-frequency GPS L₁ data, L₅ to the wide-lane linear combination, and fixed to resolved ambiguities. SWE stands for the simultaneous SWE estimation and ambiguity resolution.

Table 1 Regression coefficients and additional statistical values for the comparison of the sub-snow GPS (L₁ float solution) with each reference sensor. “Combined” indicates the average of the manual, snow pillow, and snow scale observations. b and m are the offset and slope of the regression line, cc is the correlation coefficient, MRB is the median relative bias, and n is the number of samples. The combined solution of the seasons 2015/16–2017/18 (in bold) shows the overall result.

The wide-lane linear combination (L₅) is expected to enhance the ambiguity resolution. It is a linear combination of the L₁ and L₂ observations and is defined as

$$\begin{array}{}\text{(5)}& L_{\mathrm{5}}={\displaystyle \frac{f_{\mathrm{1}}}{f_{\mathrm{1}}-f_{\mathrm{2}}}}{L}_{\mathrm{1}}-{\displaystyle \frac{f_{\mathrm{2}}}{f_{\mathrm{1}}-f_{\mathrm{2}}}}{L}_{\mathrm{2}},\end{array}$$

with the two frequencies f₁=1575.42 MHz and f₂=1227.60 MHz. The wide-lane wavelength λ_w of 86 cm is significantly larger than the L₁ wavelength (λ=19 cm), allowing a better separation of the SWE and the ambiguity parameters.

$$\begin{array}{}\text{(6)}& {\mathit{\lambda }}_{\mathrm{w}}={\displaystyle \frac{c}{f_{\mathrm{1}}-f_{\mathrm{2}}}}=\mathrm{86}\mathrm{cm}\end{array}$$

5 Results

Section 5.1 describes the GPS-derived SWE estimation results. The dependency on the GPS processing, especially ambiguity resolution techniques, is illustrated further on.

From the end of April to the beginning of May 2016, no data are available due to a storage failure caused by internet connectivity breakdown. A loss of lock in the middle of April 2018 also resulted in missing data of the sub-snow GPS station. The receiver was restarted at the end of June 2018, however, after the snow had melted completely.

5.1 Snow water equivalent estimation

The SWE is estimated based on an ambiguity float solution (L₁ float) using the derived model (Sect. 3). The time series of the sub-snow-GPS-derived SWE estimations for the 2015/16–2017/18 seasons is shown in Fig. 4a. Generally, the SWE derived from sub-snow GPS corresponds well to the reference sensors over the three seasons. However, the GPS-derived SWE corresponds inconsistently to each of the reference sensors as illustrated by the difference (ΔSWE) of the GPS-derived SWE from the individual reference sensors (Fig. 4b). The GPS-derived SWE fits very accurately to the snow scale during the dry snow period in 2015/16. Later on, at the beginning of the melting season, the GPS-derived SWE fits in between the snow pillow and snow scale observations and fits very accurately to the manual observations.

The sub-snow GPS system seems to overestimate the SWE from January to May 2017, compared to all three reference sensors, but fits best to the snow pillow observations until April 2017. In the following melting period, the GPS-derived SWE corresponds accurately to the snow scale observations. The same behavior is observed until January 2018. Afterwards, the SWE is again overestimated by the sub-snow GPS system. This could be due to an uneven snow distribution, which is not captured by the model yet. Additionally, deviations between the snow pillow and snow scale in wet snow conditions are present at the Weissfluhjoch test site for each melting period. Figure 4c shows the RMS of the sub-snow-GPS-derived daily SWE estimations from the Bernese GNSS software output over the three seasons. The SWE is estimated with a RMS of approximately 1 mm w.e. on average. An increase in the RMS with SWE is seen during each season. The high RMS values in May 2016 are caused by a significant reduced number of observations in this time period (Sect. 7).

Figure 4 Time series of the GPS-derived SWE estimations and the reference sensors (a), the differences among the three reference sensors (b), and the SWE estimation RMS (c) of the GPS float solution for the 2015/16–2017/18 seasons.

Figure 5 shows the regression analysis for the sub-snow-GPS-derived SWE estimations with the reference sensor's SWE observations for the 2015/16–2017/18 seasons. The regression analysis is calculated with respect to (a) the combined reference, (b) manual, (c) snow pillow, and (d) snow scale observations. The regression coefficients and additional statistics are listed in Table 1. Time periods without snow above the sub-snow GPS antenna are thereby excluded.

The sub-snow-GPS-derived SWE is highly correlated with a cross-correlation coefficient (cc) of 0.99 to all reference sensors and 0.98 to the manual observations. Please note that the manual observations are more sparse and thus have much fewer values to compare with (number of samples n, Table 1). The regression slopes (m) underline close agreement between the GPS-derived SWE and the reference sensor measurements. The sub-snow GPS overestimates the SWE compared to the manual observations with an offset (b) of 21 mm. The GPS-derived SWE values are biased by −28, −27, and −31 mm w.e. with respect to the snow pillow, snow scale, and combined reference measurements, respectively.

Generally, the GPS-derived SWE shows a good agreement over the three full seasons to the reference sensors with an RMSE of 70 mm w.e. to the combined reference and the snow scale observations. The RMSE increases to 77 mm w.e. for the snow pillow and 95 mm w.e. for the manual observations. Best agreement is demonstrated to the snow scale with a median relative bias (MRB) of 3.5 % over 664 samples, followed by the snow pillow with a MRB of 10.1 % over 680 samples. The manual observations fit least to the GPS-derived SWE with a MRB of 19.1 %; however only 51 samples are available for this comparison. Nevertheless, the GPS-derived SWE agrees well with the combined reference with a MRB of 8.5 % over 685 samples for the seasons 2015/16–2017/18. A MRB below 10 % is considered a good agreement as each reference sensor is prone to errors, which is discussed in Sect. 6.

The statistics of the comparison to the combined reference are also listed in Table 1 for each individual season 2015/16, 2016/17, and 2017/18 with 204, 299, and 182 samples per season, respectively. The best results are obtained for the 2015/16 season with a correlation coefficient of 0.99, a MRB of 1.4 %, a regression slope and offset of 1.1 and −24 mm w.e., and a RMSE of 52 mm w.e. The 2016/17 season shows the least agreement in terms of the MRB of 12 %. The regression statistics illustrate good agreement as well with a cc of 0.99, a regression slope and offset of 1.2 and −16 mm w.e., and a RMSE of 61 mm w.e. Results of the 2017/18 season agree well with a MRB of 9.6 %, a cc of 1.0, a regression slope and offset of 1.2 and −84 mm w.e., and a RMSE of 97 mm w.e.

Figure 5 Regression analysis of SWE estimation from GPS float solution to (a) the mean of the three reference sensors, (b) the manual observations, (c) the snow pillow, and (d) the snow scale for the 2015/16–2017/18 seasons. The black dotted line represents the 1 : 1 line.

5.2 Sensitivity on GPS ambiguity processing

The SWE estimations derived from the sub-snow GPS L₁ float solution agree well with the reference sensors (Sect. 5.1). The SWE estimations are, however, assumed to depend strongly on the GPS processing, especially the GPS ambiguity resolution. The SWE estimation sensitivity on the ambiguity resolution techniques is therefore further investigated for the strategies described in Sect. 4.2: L₁ float, L₁ fixed, L₁ SWE fixed, and L₅ SWE fixed.

Figure 6 illustrates the effect of the different ambiguity resolution techniques on the SWE time series derived from the sub-snow GPS system. The time series are compared to the combined reference observations (daily average of snow pillow, snow scale, and manual observations). The L₁ float, the L₁ SWE fixed, and the L₅ SWE fixed solutions agree with each other over the 2015/16–2017/18 seasons. All three strategies overestimate the SWE over the three seasons compared to the combined reference. The L₁ float and L₁ SWE fixed solutions thereby fit best to the combined reference.

Figure 6 GPS-based SWE estimation (a) and its RMS (b) from different ambiguity resolution strategies. The black line shows the mean SWE of the three reference sensors. Panel (c) shows the differences of GPS-based SWE estimation from different ambiguity resolution strategies to the mean SWE of the three reference sensors.

Conversely, the L₁ fixed solution deviates significantly from the combined reference in all three seasons, especially during the beginning of the melting seasons and the ablation periods (Fig. 6a and c). The L₅ fixed solution shows the same behavior, although not plotted in the present paper for visibility reasons. Both methods resolve the GPS ambiguities in a first step and estimate the SWE in a second step with introducing the fixed ambiguities as known. These methods underestimate the SWE in all three seasons on the order of a half to a full wavelength (about 19 cm for GPS L₁). The bias is thus assumed to be caused by a false ambiguity resolution. A constant part of the SWE is absorbed by the ambiguity parameters. This can occur if the excess path length caused by the overlying SWE is the same as the GPS wavelength. This behavior is shown for all three seasons for the L₁ fixed solution at a SWE higher than approximately 200 mm w.e. The SWE estimation is thus biased as the false resolved ambiguities are fixed in all following processing steps. Note that the L₅ fixed solution has a very low MRB. This is due to an overestimation of SWE in the accumulation period and a stronger underestimation in the ablation period resulting in low values in average.

The combined estimation of the SWE and the ambiguity parameters in one step improves the results significantly (L₁ SWE fixed and L₅ SWE fixed). The correlations of the SWE and the ambiguity parameters are calculated exemplarily for 1 day with high SWE (approximately 700 mm w.e. on 24 May 2017). The SWE and ambiguity parameters are thereby not correlated with correlations of 0.2 on average, allowing a good separation of the SWE and ambiguity parameters in the combined processing.

Figure 6c shows the RMS of the different ambiguity resolution techniques. The L₁ float and the L₁ SWE fixed solutions estimate SWE most accurately with an RMS around 1 mm w.e. The float solution is thereby more noisy. The L₅ SWE fixed solution has a higher RMS caused by the approximately 5.7 times higher noise of the wide-lane linear combination (Hofmann-Wellenhof et al., 2001). All three strategies show a seasonal trend, correlated to the SWE. The RMS of the L₁ fixed solution is the highest and most noisy. Problems occurred in May 2016, when all solutions except the float solution show high noise in the time series (SWE and ΔSWE) and the RMS. This is caused by a significant reduced number of observations in this time period (Sect. 5.3).

Table 2 Regression coefficients and additional statistical values for the comparison of different sub-snow GPS processing strategies with the combined reference sensors (average of the manual, snow pillow, and snow scale observations) for the 2015/16–2017/18 seasons. b and m are the offset and slope of the regression line, cc is the correlation coefficient, MRB is the median relative bias, and n is the number of samples. The lines in bold indicate the solutions which thereby fit best to the combined reference (L₁ float and L₁ SWE fixed) using only single-frequency signals.

The regression coefficients and additional statistics for the different ambiguity resolution techniques compared to the combined reference are listed in Table 2 for the 2015/16–2017/18 seasons. Detailed results from the L₁ float solution are summarized in Table 1. The sub-snow-GPS-derived SWE from the L₁ float, L₁ SWE fixed, and the L₅ SWE fixed solutions is highly correlated to the combined reference with a cc of 0.99 and a MRB of 8.5 %, 8.0 %, and 11.4 %, respectively. The L₅ fixed solution has a cc of 0.98 and a lowest MRB of 5.4 %. The RMSEs of the L₁ float (70 mm w.e.), L₁ SWE fixed (66 mm w.e.), the L₅ fixed (72 mm w.e.), and L₅ SWE fixed (75 mm w.e.) solutions are approximately equivalent. The slope of the regression analysis is 1.1 for the L₁ SWE fixed and L₅ fixed solutions and 1.2 for the L₁ float and L₅ SWE fixed solutions. The GPS-derived SWE is overestimated with all ambiguity resolution techniques with an offset of −15 mm w.e. (L₅ fixed), −23 mm w.e. (L₅ SWE fixed), −26 mm w.e. (L₁ SWE fixed), and −31 mm w.e. (L₁ float). The L₁ SWE fixed solution is favorable in terms of statistics and agreement with the combined reference, followed by the L₁ float solution. However, systematic and stochastic effects in the GPS processing residuals (Sect. 5.4) are still present for the L₁ float and L₁ SWE fixed solutions in all three seasons. The L₁ fixed solution, with a cc of 0.92, is least correlated to the combined reference and has the highest MRB and RMSE of 23 % and 154 mm w.e., respectively. The regression analysis of the L₁ fixed solution has a low offset of −9 mm w.e. and a slope of 0.8. Fewer samples (587) are a result of an unsuccessful ambiguity resolution.

Figure 7 Properties of the sub-snow GPS processing: sigma a posteriori (σ_post, a), number of ambiguity parameters (b), and number of observations (c) for different ambiguity resolution techniques.

5.3 GPS processing properties

Figure 7 shows different GPS processing properties: sigma a posteriori (σ_post, panel a) of the least-squares GPS processing and number of ambiguity parameters (panel b) for different ambiguity resolution techniques during the differential GPS processing over the 2015/16–2017/18 seasons. The number of observations for the differential GPS L₁ processing are illustrated in Fig. 7c. Around 20 000 daily observations are possible at the Weissfluhjoch test site. A high amount of liquid water in the snowpack attenuates the GPS signals strongly (Steiner et al., 2018). The reduced signal strength results in a frequent loss of lock to the satellites and fewer observations. The reduced number of observations are visible in all 3 years during the melting periods and depend on the liquid water content and the receiver tracking threshold. After the data storage failure in April 2016, the number of observations drops significantly due to the wet snowpack and increases very slowly with decreasing SWE. The antenna connector of the sub-snow GPS antenna was nearly broken in September 2016, causing a high noise in the observation numbers in this time period. The connector was changed and the number of observations was again stable at around 20 000 at the end of September 2016.

Figure 7a shows the sigma a posteriori (σ_post) of the GPS processing from the Bernese GNSS software for the different ambiguity resolution techniques described in Sect. 4.2. The L₁ fixed solution is least accurate with a maximum σ_post of 100 mm. The L₁ fixed solution has a high noise, especially during the beginning of the melting period 2018 and shows strong seasonal trends. The L₁ float and the L₁ SWE fixed solutions agree with each other, with a maximum σ_post of 35 mm in April 2018. A seasonal trend is visible for all 3 years. The L₅ SWE solution shows the lowest σ_post of maximally 10 mm and no significant seasonal trend. A larger number of parameters (more observations of the L₅ solutions due to L₁ and L₂ observations or included ambiguity parameters in the L₁ float solution) reduces the σ_post significantly due to increased redundancy.

Figure 8 Mean (μ) and standard deviation (σ) of the GPS double difference residuals from the L₁ float solution before (black) and after (red) SWE estimation.

Figure 9 Mean (μ) and standard deviation (σ) of the GPS double difference residuals after SWE estimation for different ambiguity resolution techniques.

The number of ambiguity parameters is illustrated in Fig. 7b before ambiguity resolution (L₁ float solution) and after (L₁ fixed, L₁ SWE fixed, L₅ SWE fixed). Approximately 70 ambiguities are set up on average for each day of the 2015/16–2017/18 seasons. The number of ambiguities increases for days with a strong melt. This is caused by the strong attenuation of the GPS signals when a high amount of liquid water is present in the snowpack. The resulting loss of lock to the satellites forces new ambiguity parameters. The number of ambiguities decreases, of course, in the case of a reduced number of observations (Fig. 7c), e.g., in May–June 2016. A higher number of ambiguities are, however, seen in June 2017. A huge number of fragmented observations result in more ambiguity parameters. All ambiguity parameters should be resolved and fixed correctly after the successful ambiguity resolution step in the GPS processing. The L₅ SWE solution resolves almost all ambiguities. This confirms the assumption of better ambiguity resolution using the wide-lane linear combination due to the better separation of ambiguity and SWE parameters caused by the large wavelength of 86 cm. In the L₁ SWE fixed solution, not all ambiguities could be resolved, especially in the beginning of the melting periods. The most significant behavior is, however, shown for the L₁ fixed solution. Thereby, no significant number of ambiguities (around 50 ambiguities) could be resolved during melting periods in 2016/17 or periods of a high SWE above 800 mm w.e. in the 2017/18 season.

5.4 GPS processing residuals

The GPS processing residuals provide information on the applicability of the developed model to the SWE estimation. The GPS processing residuals should be normally distributed around zero without error influences, such as the overlying snowpack. The snowpack above the GPS antenna deteriorates the GPS observations based on the path delay of the GPS signals (Sect. 3). This effect should induce a significant seasonal effect in the GPS double difference residuals (further called GPS processing residuals) if no model is applied in the processing. Figure 8 shows the daily mean (μ) and standard deviation (σ) of the GPS processing residuals from the L₁ float solution over the 2015/16–2017/18 seasons. The mean of the residuals illustrates systematic effects, whereas the standard deviation points out stochastic effects. The black line indicates no applied model and the red line illustrates the residuals if the SWE estimation model is applied. Significant systematic and stochastic effects show up in the mean and standard deviation if no model is applied. Both effects follow the seasonal SWE development with a maximum of around 18 mm in the mean and 300 mm in the standard deviation.

The systematic trend in the mean of the GPS processing residuals is eliminated by applying our derived model. Noise of the order of 1 mm on average is still present, especially in 2018, when the SWE was above average. The stochastic effects in the standard deviation are significantly reduced to 40 mm on average. The model is thus able to correctly estimate the SWE as the effects in the residuals due to the overlying snowpack are significantly reduced. The remaining noise could be due to multipath effects, the assumption of a flat surface at the air–snow interface, the neglected roughness, or layered snowpack (Sects. 4 and 3).

Figure 9 illustrates the influence of the different ambiguity resolution techniques on the daily mean (μ) and standard deviation (σ) of the GPS processing residuals. The L₁ float solution is already described for Fig. 8. The L₅ SWE solution agrees with the L₁ SWE fixed solution in the daily mean (Fig. 9a) with small deviations. Large deviations are present in May 2017. Both techniques model the systematic effect of the overlying snowpack quite well, with a remaining noise around 3 mm. The standard deviation of the L₁ SWE solution performs almost similar to the L₁ float solution for all seasons. The standard deviation (Fig. 9b) of the L₅ SWE solution is approximately 5.7 times higher than the L₁ float solution due to the increased noise level of the linear combination. The residuals for the L₁ fixed solutions are significantly higher, with maximal values reaching 20 mm in the mean and 200 mm in the standard deviation of the GPS processing residuals. About 70 % of the snowpack effect remains over the three seasons compared to the solution without an applied model (Fig. 8). This suggests a wrong ambiguity fixing in the first step, leading to a strong weakness of the derived model to estimate the SWE correctly when using the L₁ fixed processing strategy.

6 Discussion

6.1 Comparison with other measurements

Generally, the sub-snow-GPS-derived SWE estimations correspond well to the reference sensors with a cross-correlation coefficient ranging from 0.97 (manual) to 0.99 (snow pillow, snow scale, combined reference) over all three seasons. All reference sensors are prone to errors and SWE observations are not consistent with each other, resulting in a lack of real ground truth data to compare with. The snow pillow and snow scale have a larger measurement area than a manual point sample, whereas the sub-snow GPS uses signals arriving at different lines of sight within an area of approximately 5 m. The three reference sensors are located at different locations within the test site (Sect. 2) and the manual observations cannot be taken exactly at the same spot each time due to the destructive method. The manual SWE observations were, however, converted to one observer location in the middle of the test site. The snow pillow and snow scale measure at different locations, approximately 20 m west next to the sub-snow GPS station, resulting in small uncertainties.

The snow depth is observed to vary about 10 cm spatially within the test site. This variation is less than 10 % of the seasonal snow depth except at the beginning and end of the three snow seasons. The maximal snow depth was approximately 1.7 m in 2015/16 and 2016/17 and approximately 3 m in 2017/18. A larger snow depth at the sub-snow GPS site could explain a small overestimation of GPS-derived SWE. Ideally, co-located SWE data should thus be used to evaluate the derived SWE estimation results. Although this would reduce the influence of the unequal snow distribution, the sensor biases still remain. Manual SWE observations have an uncertainty of about 10 % due to the within-site variability in snow density (Jonas et al., 2009). Existing ice layers or percolated liquid water complicates the snow sampling (Smith et al., 2017).

Snow-weighting sensors like the snow scale and snow pillow experience measurement artifacts (e.g., bridging effects), especially at the beginning of the melt season. Snow pillows are usually less affected than snow scales due to the larger surface (Beaumont, 1966). Furthermore, meltwater percolates through the snowpack towards the ground, forming a basal liquid water or ice layer, infiltrating the soil or percolating out of the measurement area. These effects can cause SWE over- or underestimation by the snow pillow, snow scale, or manual observations. A basal meltwater layer still delays the GPS signals and could contribute to the overestimation of SWE in the beginning of the melting periods, increasing the bias between the GPS and reference SWE observations.

6.2 Sensitivity on GPS ambiguity processing

The SWE estimation based on the sub-snow GPS system is highly sensitive to the GPS ambiguity resolution techniques in the phase-based double difference GPS processing (Sect. 5.2). A successful SWE estimation is possible by simultaneously resolving the GPS ambiguities (L₁ SWE fixed solution). Ambiguity resolution in a separate processing step (L₁ fixed solution) before SWE estimation with the introduced fixed ambiguities leads to less resolved ambiguities or false ambiguity resolution and thus biased SWE. This is especially the case for SWE higher than one wavelength (λ=19 cm), as the path delay induced by the SWE is partly compensated for by the ambiguity parameters. The error in the SWE can thereby reach 200 mm w.e. Using a larger wavelength of the L₅ SWE solution facilitates the separate SWE and ambiguity estimation, compared to the L₁ fixed solution. However, the measurement noise is significantly increased when using a linear combination and multiple frequencies are required.

The L₁ float and L₁ SWE fixed solutions performed best over all parameters. The SWE is estimated with the lowest RMS, agrees best with the reference sensors (lowest RMSE, MRB, and highest cc), and the fewest systematic and stochastic effects are present in the GPS residuals. The L₁ SWE fixed solution performs slightly better than the L₁ float solution in terms of the statistical values. The L₁ float solution is, however, chosen as the favorable SWE estimation strategy as no ambiguity resolution is required. Thereby, the processing is simplified significantly and false ambiguity resolution is prevented. Moreover, only single-frequency data are needed for the L₁ float or L₁ SWE fixed processing strategies, allowing cheaper equipment and a faster processing.

6.3 Advantages and limitations

The sub-snow GPS system provides a new and promising method for daily SWE quantification. The method is easy to install, requires little maintenance, is nondestructive, and provides automatic observations with a high temporal resolution. Due to the small size, the system could even be installed on mountain slopes. No access is required during the snow period and this method could therefore be applied, for example, in avalanche-prone terrain. Power supply of the GPS receiver can be a limiting factor in these areas. Lossless cables from the sub-snow GPS antenna to the receiver, however, allow GPS receiver and power supply installation outside the avalanche terrain.

A small SWE (10 mm w.e.) could already be quantified from early snowfalls, e.g., October 2016. Precise SWE estimations were possible until 800 mm w.e. in 2016 and 2017. SWE values above 800 mm w.e. (2018) are overestimated by the sub-snow GPS SWE estimation model, compared to the reference sensors. The overestimation could be explained by an uneven snow distribution within the test site, which is not yet captured by the model. The upper limit is not assessed yet. Melted snow percolated to the ground, forming a liquid water layer above the sub-snow GPS antenna that could prevent GPS signal reception if it were deeper than 35 mm (Steiner et al., 2018).

6.4 Representativeness

The present study is carried out over three snow seasons in the Alps at a high altitude. The study should thus be representative for an Alpine snowpack at similar altitudes, as a large number of data are analyzed, including melting seasons. Using longer baselines between the GPS reference station and the sub-snow GPS antenna with a large difference in elevation could complicate the SWE estimations due to tropospheric refraction. The derived sub-snow GPS model should also be representative for a polar snowpack as it is usually dry. Problems could be caused by the GPS satellite distribution with missing satellites in the north. Adding more satellite systems could increase the observations and improve SWE estimation analogously. The presented method is, however, not representative for forested areas, as satellite visibility is very limited and the GPS phase signals are highly attenuated, delayed, or obstructed by the trees. A good GPS visibility is important at all locations and should be considered, especially in narrow mountain areas.

7 Conclusions

The newly developed model is applied to a seasonal snowpack in order to investigate the potential of using GPS L₁ observations from a geodetic GNSS system for daily SWE estimation. Sub-snow GPS is a promising method for point-wise SWE quantification. The snowpack is not destroyed or disturbed due to the automated, continuous, self-sustainable observation method and the effort for installation is relatively small. Remote (online) access is possible and almost no maintenance is required for the small-sized equipment. The presented model enables the direct estimation of SWE if both the reference and submerged station coordinates are precisely known. The use of a single water layer model for SWE quantification is encouraged as the SWE depends on the bulk relative permittivity in the snowpack. SWE could be estimated with a relative bias below 10 % compared to a combination of three independent reference sensors (snow pillow, snow scale, and manual observations). The proposed method successfully estimated SWE over three full seasons, including ablation periods.

The assumption of GPS ambiguity resolution as a critical parameter for GPS observations within a snowpack can be confirmed. False ambiguity resolution biases the SWE estimation by up to 200 mm. It is shown to be important to resolve the ambiguities in one step together with the SWE estimation, instead of separately. In any case, an ambiguity float solution performs better than using false resolved ambiguities.

The promising results of this study encourage the assessment of a low-cost GPS system and the comparison with the geodetic equipment in a next step. The potential of sub-daily SWE estimation, e.g., higher temporal resolution, will be evaluated further on.