1. Introduction
The estimation of the spatial distribution of the water equivalent of the snow cover (SWE) in mountain catchment areas is presently the most important unsolved challenge in snow hydrology [
1]. In situ data is important for calibration and validation of extensive remote sensing observations and snow hydrological modeling [
2,
3,
4]. The works of [
5,
6] give a detailed description of various ground-based techniques for observing the SWE.
Global navigation satellite system (GNSS) refractometry permits accurate and continuous in situ SWE monitoring [
7,
8]. A GNSS reference (base) antenna is mounted above the snowpack, while a second GNSS (rover) antenna is placed on the ground where it is continuously covered by snow accumulation (
Figure 1). Observed GNSS signals are refracted at the air/snow boundary (and between layers within the snowpack) and delayed while propagating through the snowpack, depending on the SWE above the buried GNSS antenna.
In contrast to manual [
9,
10] point-wise in situ SWE observations, GNSS-refractometry-based SWE monitoring is non-destructive, needs little effort and cost for installation and maintenance, and provides high temporal resolution. Compared to automatic snow depth sensors, the GNSS refractometry method directly estimates the SWE without the need for very scarce and laborious density observations for the conversion of snow depth to SWE [
11]. The spatial footprint of GNSS refractometry depends on the permittivity and the snow depth of the snowpack overlying the buried GNSS antenna and the angle of incidence of the GNSS signals [
12]. For the most accurate GNSS processing, available GNSS signals are tracked from all elevations and azimuths. In case of dry snow, the spatial footprint of GNSS refractometry varies from 0.045
to 41
(0.2 m to 7.2 m in diameter) for snow depths above the buried GNSS antenna of 0.1 m and 3 m. For wet snow, the spatial footprint is much smaller with 0.012
to 11
(0.1 m to 3.7 m in diameter) for snow depths of 0.1 and 3 m. The representativeness of GNSS refractometry using a single setup is thus limited to a point-wise scale.
The temporal and spatial resolution is comparable with automatic SWE observations by buried snow pillows, scales, ground radar, cosmic ray, and acoustic snow gauges [
13,
14,
15,
16,
17,
18,
19,
20]. The main advantage of GNSS refractometry compared to these buried automatic SWE observation sensors is the very low cost and size of the system, as well as the global availability and coverage of GNSS signals, enabling the permanent installation of several distributed buried GNSS sensors at a low cost. Additionally, a mobile GNSS refractometry setup enables distributed SWE measurements at relatively high speeds and low cost, increasing the spatial footprint of this method. Due to the high spatial and temporal variability of the SWE at small to regional scales, several distributed in situ SWE observations from intermediate scales of several hundred meters (e.g., mountain catchment areas) are of high value to complement regional-scale snow modelling and remote sensing products [
21].
Following the idea of [
7,
22], a SWE estimation model was developed by [
12] based on the Global Positioning System (GPS) phase signal path delay due to a present single water layer of depth
d, being equivalent to the SWE above a submerged/buried GNSS antenna. Refs. [
23,
24] demonstrated the successful daily and hourly SWE estimation over three complete (dry and wet snow) seasons using refracted GPS signals. Single-frequency GPS (GPS L1) data collected by high-end GNSS sensors were thereby thoroughly analyzed using the scientific Bernese GNSS processing software. Refs. [
25,
26] used a similar, phase-based, approach for SWE estimation with low-cost GPS L1 sensors. The use of low-cost sensors for hourly SWE retrieval based on GNSS refractometry was independently verified by [
27] for a complete dry and wet snow season. In these studies, the SWE was introduced and estimated as an additional parameter in the GNSS phase observation equation. The impact on the coordinate solution in case of a buried GNSS antenna was demonstrated for the first time by [
24]. A strong correlation (99.7%) between the SWE and the coordinate Up component was revealed, illustrating the feasibility of using the Up component, biased due to the overlying snowpack, for SWE estimation.
In contrast to previous studies, the present study investigates the feasibility of high-temporal (near) real-time in situ SWE determination by GNSS refractometry. Multi-frequency and multi-system GNSS observations from newly available (since 2020) multi-band low-cost GNSS sensors are used to retrieve the SWE based on the estimation of the bias in the coordinate Up component. With this method, no additional SWE parameter is needed in the GNSS observation equation, and the SWE can be directly estimated using standard, open-source GNSS processing software, which significantly simplifies the implementation. For the first time, GNSS refractometry can thus be applied in (near) real time directly in the field using real-time kinematic (RTK) GNSS refractometry in a fixed and mobile low-cost snow monitoring setup.
An overview about fixed and mobile GNSS refractometry snow monitoring setups and available reference data is given in
Section 2.
Section 3.1 summarizes the methodology. The data processing is specified in
Section 3.2.
Section 4 delineates results of post-processing in
Section 4.1 and (near) real-time SWE estimations from a fixed setup in
Section 4.2 and a mobile setup in
Section 4.3.
Section 5 provides a discussion and conclusion.
2. GNSS Refractometry Snow Monitoring Setups
GNSS refractometry snow monitoring systems were set up at two different study sites for investigating the feasibility of GNSS refractometry for SWE estimation using the standard open-source RTKLIB software in post-processing (15 min resolution) and (near) real time. For the fixed, low-cost setup and near real-time processing, the GNSS baseline was estimated with a temporal resolution of 1 Hz and subsequently filtered by a moving median over 24 h and resampled to 10 min. For the (near) real-time processing from the mobile, low-cost setup, the GNSS baseline was determined directly in the field over a 15 min observation window.
The fixed setups were installed on stable ground during the snow-free period to monitor changes in the SWE over time at a specific location. Data were stored or processed in real time. Due to the fixed location, the estimated SWE time series can be filtered subsequently which leads to enhanced SWE results. The fixed setups have a high temporal, but a low spatial resolution due to the point-wise measurements. To increase the spatial resolution, the mobile setup was developed. This lightweight and low-cost mobile setup was designed to enable SWE measurements in snow-covered areas at different locations in near real time. Little effort is needed as the small-sized system fits in a hole drilled into the snowpack, making the measurement additionally fast compared to traditional snow tube samplers where digging snow pits is a labor-intensive prerequisite, especially in deep snowpacks. Due to the availability of the SWE result directly in the field, SWE measurements are possible in different snow depths in the hole to indirectly derive layered density information and/or in different spatial locations to estimate the spatial heterogeneity of the SWE.
2.1. Fixed, High-End Setup for Post-Processing
GNSS refractometry data from the Davos Weissfluhjoch test site at 2536 m a.s.l., operated by the WSL Institute for Snow and Avalanche Research SLF (WSL SLF), are available for the season 2016/17 from a previous study [
23]. The high-end GNSS reference station (Leica GR10 receiver, Leica AS10 antenna) was deployed on a pole at a 5.3 m height. A high-end GNSS antenna (Leica AS10) was mounted on the ground below and connected to a high-end GNSS receiver (Leica GR10). Raw GNSS observations (phase, code, signal-to-noise) are available with a 1 Hz sampling rate for GNSS post-processing. Pressure based snow scale measurements (30 min) and biweekly manual snow tube observations from the test site (10–20 m next to the GNSS sensors) serve as ground truth.
Multi-frequency and multi-system GNSS data collected from this fixed, high-end GNSS refractometry snow monitoring setup is publicly available at [
28] and reprocessed to evaluate the feasibility of post-processed SWE estimation using RTKLIB (
Section 3.2.1 and
Section 4.1).
2.2. Fixed, Low-Cost Setup for Near Real-Time Processing
For the first time, a near real-time (1 Hz) GNSS refractometry snow monitoring system was set up and installed at the WSL SLF Davos Laret test site at 1500 m a.s.l. (
Figure 2). The test site is located in a North–East to South–West facing mountain valley and surrounded by high mountains. Ground-truth data is provided by a close-by snow scale with 10 min sampling rate and weekly manual snow tube observations, both within a distance of 5 m to the GNSS refractometry setup.
The fixed, low-cost GNSS refractometry snow monitoring setup is composed of two Emlid Reach M2 receivers and two u-blox ANN-MB-00 antennas (
Figure 2). These user-friendly instruments allow multi-frequency and multi-system GNSS observations at a low cost, low weight, minimal size, and high precision. The GNSS reference antenna is mounted on a pole at 2.8 m height with a very short GNSS baseline to the sub-snow antenna which is fixed on the ground. The system was installed in November 2021 after the first snow fall and operates in RTK GNSS mode. Estimated real-time (1 Hz) GNSS baselines, in East, North, and Up components, were recorded from November 2021–April 2022 and are publicly available at [
29]. Due to a system outage at the base receiver on 13 April, no baselines are available for the last week of the season during the snowmelt. The near real-time feasibility of GNSS-refractometry-based SWE monitoring using the fixed, low-cost setup is assessed for the most recent 2021/22 season at Davos Laret (
Section 3.2.2 and
Section 4.2).
2.3. Mobile, Low-Cost Setup for (Near) Real-Time Processing
A mobile GNSS refractometry snow monitoring setup was designed and consists of two low-cost multi-frequency and multi-system GNSS antennas (u-blox ANN-MB-00) mounted on an avalanche probe (
Figure 3b,c) and connected to the low-cost multi-frequency and multi-system GNSS receivers (Emlid Reach M2). The antennas were mounted on a 5 mm thick aluminum plate with 8.3 cm in diameter. The plate was attached to the probe by a welded half tube of 3.5 cm length and cable ties, allowing us to move the upper base antenna and place it above the snowpack. The lower rover antenna was mounted at the bottom of the probe so that the antenna could be placed on the ground to estimate the bulk SWE, averaged over an area of several square meters around the antenna. The lightweight setup is portable and delivers in situ SWE estimates directly in the field. To place the rover antenna below the snow, a vertical hole was drilled with an 11.2 cm diameter ice drill (Kogha Arctic). Drilling a hole (
Figure 3a) proved to be significantly faster and less destructive compared to digging a snow pit. Furthermore, snow pits can affect the signal propagation in the cone above the rover antenna which, in turn, deteriorates the SWE estimates.
The system is configured in RTK GNSS mode (
Section 3.2.3) where the base receiver transmits the base position and all GNSS observations to the rover via an internal Wi-Fi. To initialize the system, both antennas need to be above the snowpack. After the rover initialized, the RTK solution state needs to be fixed, indicating that the float phase ambiguities are resolved to integer values. The rover was then slowly lowered to the bottom of the drilled hole. It is important to keep the fixed RTK solution which could be lost by lowering the rover antenna too fast in the snowpack. In this case, the system needs to be elevated above the snowpack until the ambiguities are fixed. Reducing the signal-to-noise (SNR) mask in the receiver configuration improves the RTK solution stability within the snowpack as GNSS signal strengths are decreased while propagating through snow [
12]. The Emlid Reach M2 can be connected via the built-in Bluetooth or Wi-Fi which allows us to control and configure the receivers via a smartphone or a tablet. The receivers are accessed through the Emlid ReachView mobile application, allowing us to read out the baseline and the coordinates of both antennas directly in the field and to configure or monitor the receiver settings and state.
The new mobile, low-cost GNSS refractometry snow monitoring setup was tested for its (near) real-time SWE estimation feasibility (
Section 4.3) at Davos Weissfluhjoch on 26 March 2021 in approximately 2.5 m deep snow close to the fixed, high-end GNSS refractometry setup (
Section 2.1). Two holes were drilled 3 m next to each other and measured three times each for 15 min. A snow profile of 31 March 2021 from the same test site (10–20 m next to the mobile experiment) served as reference. The snowpack settled between both dates without new snow accumulation, influencing only the snow density, not the bulk SWE.
3. Method
In the present study, the SWE of the snowpack overlying the buried GNSS rover antenna was estimated using the GNSS refractometry method based on the bias in the coordinate Up (height) component [
24].
3.1. GNSS Refractometry Based on the Biased Up Component
Introducing the GNSS signal path delay
as an additional parameter in the GNSS phase observation equation (Equation (
1); [
30]) enables us to precisely estimate the SWE by knowing and fixing the coordinates of the base and rover station [
23]:
Thereby, is the observed carrier phase measurement in meter, and is the geometric distance between the GNSS antenna A and a satellite j, impacted by the atmospheric, relativistic, and multipath delays, as well as the satellite and receiver clock biases ; represents the ambiguity phase bias and the measurement noise. The differential GNSS processing to the close-by reference station mitigates the atmospheric (tropospheric and ionospheric) and relativistic effects and the satellite clock errors. The impacts of the antenna phase center offsets and variations are minimized by using similar GNSS antenna types with identical orientations for the base and the rover.
Instead of introducing and estimating the additional SWE parameter (
) in the GNSS observation equation, the height bias of the GNSS rover antenna coordinate due to the snow cover was estimated in the present study. The SWE above the buried GNSS antenna was then derived from the bias in the estimated GNSS baseline Up component. The difference in the station height
h or the Up (height) component of the baseline were strongly correlated with the SWE [
24].
Assuming a single water layer model [
12], the path delay
, introduced by the signal propagation in water is collinear to the mapping function of the GNSS station height (
h) and the receiver clock (
) parameters [
8,
24]:
The path delay depends on the SWE above the buried GNSS antenna, the refractive index of the representing single water layer , and the observation zenith angle z. represents the mapping function to scale the SWE observation bias depending on the observed zenith angles.
The
part of
is almost constant due to
and varies only up to 0.06 between
and
. This part can thus be absorbed by the constant station clock parameter (
).
The
dependent part of
is absorbed by the station height coordinate by 99.9% [
24]. The influence of the SWE can therefore be represented as a linear combination of the station clock and the height coordinate:
The station clock error and thus the constant part of cancels out by double-difference GNSS processing. The change in the SWE above the sub-snow antenna can therefore be derived based on the change in the estimated height component (Equation (4b)) of the rover coordinate. Equivalently, the Up component of the estimated baseline between both antennas can directly be used.
3.2. Data Processing with RTKLIB
The GNSS refractometry method based on the bias in the coordinate Up component [
24] was applied for post and (near) real-time SWE estimation using RTKLIB.
3.2.1. Post-Processing
To test the performance of the open-source standard GNSS processing software RTKLIB, the available seasonal dataset of the fixed, high-end GNSS refractometry snow monitoring setup from Davos Weissfluhjoch (
Section 2.1) was post-processed in RTKLIB (version 2.4.3 b34). The multi-band GNSS observations were processed in time intervals of 15, 20, 30, 40, 45, and 60 min. As no significant deviations in the estimated SWE time series were present for the different time intervals and least noise was observed for the shortest interval, results are only presented here for 15 min processing intervals. Observations are used from all elevations as the use of an elevation mask reduces the SWE estimation accuracy [
24].
Appendix B.1 lists the command used for automated RTKLIB post-processing with the configuration file and python script uploaded to [
31]. Outliers were removed in resulting SWE time series based on a statistical threshold (
) and filtered by a moving median over 24 h. The GNSS-derived SWE results were resampled to match the temporal resolution of the manual (daily) and snow scale (30 min) reference observations.
3.2.2. Near Real-Time Processing
To evaluate the feasibility of near real-time SWE monitoring using GNSS refractometry, the fixed, low-cost GNSS snow monitoring system at Davos Laret operated autonomously in RTK GNSS positioning mode. RTKLIB was also used in real-time applications, as it was implemented in the Emlid Reach (M2) RTK receivers whereby the mobile ReachView application serves, among others, as a web interface for RTKLIB [
32]. The GNSS baseline in East, North, Up (ENU) between the base and rover antennas was measured and logged continuously at 1 Hz sampling rate. The change in the estimated baseline Up component (stored in .ENU log files) was equivalent to the SWE. Similarly to
Section 3.2.1, outliers were removed in the resulting SWE time series based on a statistical threshold (
) and filtered by a moving median over 24 h. As the system was installed after the first snow fall and the true physical baseline could not be measured (
Section 2.2), the
SWE was corrected based on the first manual observation of the SWE. The GNSS-derived SWE results were resampled to match the temporal resolution of the manual (daily) and snow scale (10 min) reference observations.
3.2.3. (Near) Real-Time Processing in the Field
The mobile, low-cost GNSS refractometry snow monitoring setup was also configured in RTK GNSS positioning mode to investigate the feasibility for (near) real-time SWE estimation directly in the field. The base station coordinates were determined and fixed over a 10 min observation interval in the field. For the absolute positioning of the measurement location, a virtual reference station (VRS) of the Swiss Positioning Service (swipos) was used to achieve an accuracy in the centimeter range.
Measurements of the mobile setup were triggered manually and the baseline in East, North, Up components was estimated over a 15 min observation interval. The estimated total baseline was directly displayed in the Emlid ReachView mobile application when connected to the rover receiver. Since the sub-snow GNSS antenna was mounted straight below the base antenna, the total baseline corresponded to the Up (height) component when the probe was held vertically.
The difference between the baseline, estimated with RTK GNSS refractometry to the actual baseline, measured physically on the avalanche probe, corresponds to the bulk SWE above the sub-snow GNSS antenna (Equation (
5)).
with
being the physically measured vertical distance between the base and rover antenna reference points and
the estimated GNSS baseline Up component.
5. Discussion
The present study investigated the applicability of the open-source GNSS post-processing software RTKLIB as well as the (near) real-time feasibility of the GNSS refractometry method using low-cost GNSS equipment. As the SWE induced GNSS path delay is highly collinear with the receiver clock and the height (Up) coordinate, the biased baseline Up component was successfully used for SWE estimation.
The seasonal SWE could be accurately determined over 15 min observation periods in post-processing using the GNSS refractometry method based on the biased GNSS baseline up component. In contrast to previous studies [
23,
24], the open-source GNSS post-processing software RTKLIB was applied to determine the SWE from a fixed, high-end GNSS refractometry snow monitoring setup with a very high temporal resolution (15 min). The software is fast, automatable, free and easy to use, and permits adequate setting options.
Moreover, the SWE was successfully estimated in (near) real time directly in the field with the developed mobile, low-cost GNSS snow monitoring setup. Fixed GNSS solutions are crucial so that the SWE can be estimated with a mean relative bias of 2.5% compared to manual observations. This is contradictory to a previous study which concluded that “resolving the GPS float ambiguities is shown to have no significant impact on the SWE estimation” [
23]. There, the GNSS ambiguities were resolved over longer (daily) processing intervals, leading to an increased stability and accuracy of ambiguity float solutions. Accurate SWE time series can thus be monitored in near real time. However, there are limitations such as GNSS antenna phase center variations not matching with the changed incidence angles due to the refraction in snow. Real-time applications in shallow snowpacks could be limited due to the absolute positioning accuracy of RTK GNSS (2–3 cm). Relative positioning with subsequent filtering is, however, more accurate and could most likely allow an accurate estimation of changes in the SWE over time. Alternatively, the GNSS positioning accuracy is significantly increased (2–10 mm) in case of post-processing and small GNSS baselines, making accurate SWE estimations also possible in shallow snowpacks. Rainfall events are assumed to cause inaccuracies on the daily scale. Such deviations prevent or limit the use of daily GNSS SWE changes as a measure for mass gain or loss. It is noticeable that the noise of the unfiltered GNSS-derived SWE estimates is significantly reduced (approximately by a factor of 10) in case of post-processing compared to the near real-time application. Subsequent filtering of the near real-time SWE estimates is thus crucial to understand if daily GNSS-derived SWE variations are caused by accumulation/ablation or by the uncertainty of the SWE estimates from GNSS refractometry.
The present results show the possibility of SWE monitoring based on 1 Hz RTK GNSS refractometry Up component estimation. In this case, subsequent outlier screening and data filtering is significant for achieving accurate and reliable SWE time series. A highly resolved (10 min) SWE monitoring is very advantageous for new casting applications such as the temporal detection of precipitation or run-off events [
33].
The SWE estimated with GNSS refractometry is currently only representative for a small spatial scale due to the point-wise measurement principle, estimating the SWE over an area of several square meters (up to 41
) above the buried GNSS antenna. The low cost and maintenance, as well as the automatic processing of GNSS refractometry observations and the developed mobile setup allows distributed GNSS refractometry observations to increase the spatial resolution from point-wise to intermediate scales. The presented method is promising due to the small equipment size, cost, and installation effort. Low-cost GNSS systems can be used to monitor the SWE by estimating the biased GNSS baseline Up component between the base and rover antenna in RTK GNSS mode. The near real-time SWE retrieval method is automated, continuous, and nondestructive when installed during the snow-free period. The GNSS snow monitoring system, however, needs to be set up following best possible GNSS practices.
Appendix A summarizes the main points which are important to enable successful GNSS observations.
Up to now, experimental GNSS refractometry snow monitoring systems consist of two GNSS antennas with one antenna buried underneath the snowpack and a close-by antenna which serves as a local reference above the snowpack. Among others, this very short GNSS baseline allows the mitigation of atmospheric GNSS path delays, such as the ionospheric and tropospheric delays. Tropospheric delays increase with higher elevation difference to the reference station and longer baselines. Ionospheric delays increase with longer baselines. They can be significantly reduced by using a linear combination of multi-frequencies. A national permanent GNSS network or virtual reference stations (VRS) could be utilized instead of the local GNSS reference for a distributed SWE monitoring in whole mountain catchment areas. Synergies can be used, and only the buried GNSS station needs to be installed in the snow monitoring site, reducing material, costs, and maintenance and allowing the installation in areas prone to avalanche danger. The feasibility of using long GNSS baselines will be investigated in a future study. The application of RTK GNSS refractometry in shallow snowpacks should be analyzed in more detail due to the limitation in the RTK GNSS positioning accuracy.