Can GNSS-R Detect Abrupt Water Level Changes?

Global navigation satellite system reflectometry (GNSS-R) uses signals of opportunity in a bi-static configuration of L-band microwave radar to retrieve environmental variables such as water level. The line-of-sight signal and its coherent surface reflection signal are not separate observables in geodetic GNSS-R. The temporally constructive and destructive oscillations in the recorded signal-to-noise ratio (SNR) observations can be used to retrieve water-surface levels at intermediate spatial scales that are proportional to the height of the GNSS antenna above the water surface. In this contribution, SNR observations are used to retrieve water levels at the Vianden Pumped Storage Plant (VPSP) in Luxembourg, where the water-surface level abruptly changes up to 17 m every 4-8 h to generate a peak current when the energy demand increases. The GNSS-R water level retrievals are corrected for the vertical velocity and acceleration of the water surface. The vertical velocity and acceleration corrections are important corrections that mitigate systematic errors in the estimated water level, especially for VPSP with such large water-surface changes. The root mean square error (RMSE) between the 10-min multi-GNSS water level time series and water level gauge records is 7.0 cm for a one-year period, with a 0.999 correlation coefficient. Our results demonstrate that GNSS-R can be used as a new complementary approach to study hurricanes or storm surges that cause abnormal rises of water levels.


Introduction
Global navigation satellite systems (GNSS) signal have been used for unanticipated purposes, such as remote sensing of the environment [1]. In the last two decades, GNSS reflectometry (GNSS-R) has emerged as a reliable method to estimate water levels, with many studies highlighting the technique's potential [2][3][4][5][6][7]. Geodetic GNSS-R uses near-surface or ground-based signals of opportunity at the L-band frequency to retrieve environmental variables such as snow depth, soil moisture, and water level at intermediate spatial scales that are proportional to the height of the antenna above the surface being observed (i.e., 1000 m 2 for a GNSS antenna installed on 2-m-tall monument) [4,[8][9][10][11][12]. This technique is based on the simultaneous reception of the direct or line-of-sight transmissions, along with their coherently reflected signals.
GNSS-R can retrieve absolute water levels by taking the vertical land motion into account. It should be noted that the vertical land motion of the tide gauge benchmarks ought to be corrected to isolate tide gauges from the land motion [13]. Furthermore, the international workshop on sea-level measurements in hostile conditions that was held from 13 to 15 March 2018 and was cosponsored by the Intergovernmental Oceanographic Commission (IOC) of the United Nations Educational, Scientific and Cultural Organization (UNESCO) has recommended the investigation of GNSS reflectometry (GNSS-R) "if conditions are too difficult for conventional observing methods"; this was in addition to the general recommendation that "tide gauges should be connected to benchmarks and co-located with GNSS as Therefore, geodetic GNSS-R would be especially valuable in the harsh environments where conventional tide gauges are very hard to maintain, and therefore, sea level observations are very sparse (i.e., sea-ice-covered regions).
The potential of GPS L1 reflections from a ground-based GNSS-R experiment first was demonstrated for sea-state monitoring in [14]. They used an upright right-handed circularly polarized (RHCP) antenna and an upside-down left-handed circularly polarized (LHCP) antenna to generate a times series of the interferometric complex field. The authors of [15] further developed the concept of the GNSS-R using GPS L1 phase delays from two geodetic-quality GNSS antennas.
The excess delay of water-surface reflections with respect to the line-of-sight transmission that is recorded with a commercial off-the-shelf GNSS receiver can be used to estimate the water level. The signal-to-noise ratio (SNR) observations from a geodetic-quality receiver/upright antenna at Kachemak Bay, Alaska were first demonstrated by [16] to retrieve the water levels. GNSS-R water levels using carrier-phase and SNR observations of Global Positioning System (GPS) and Russian global navigation satellite system (GLObal'naya NAvigatsionnaya Sputnikovaya Sistema; GLONASS) from a GNSS station in Sweden were compared with those of a co-located tide gauge sensor; the results showed correlation coefficients of 0.86 to 0.97 [17]. Using the GNSS station on top of the Cordouan Lighthouse for a three-month period, [18] reported a root means square error (RMSE) of 0.87 m between GNSS-R water levels and tide gauge (TG) records that were located ~9 km from the lighthouse, taking the dynamic of the reflection surface into account.
In this contribution, we used SNR observations for about a one-year period (from 27 June 2017 to 28 July 2018) at the Vianden Pumped Storage Plant (VPSP) in Luxembourg (Figure 1a) to report the water-surface levels to investigate if geodetic GNSS-R can detect sudden and large water-surface changes caused by nonastronomical forces. We corrected the estimated GNSS-R heights for tropospheric propagation [19]. Then, we corrected the reflector heights for vertical velocity and acceleration of the water surface [6]. The water level of the upper reservoir can abruptly change up to about 17 m every 4-8 h in either generating or pumping modes. This process repeats daily but at different times of the day when the energy demand is high. Therefore, the dynamic satellite and surface corrections are very important, as the water surface at VPSP has a unique water-surface regime. In the last decade, the velocity and vertical velocity corrections were estimated by fitting a daily sinusoidal fit based on the mean frequencies of the dominant tides [16]. However, the traditional method cannot be applied, as the water-surface changes are not caused by diurnal and semidiurnal tides at VPSP. Here, the estimated GNSS-R water levels are corrected using centered finite differences. Finally, the 10-min GNSS-R water levels with a two-hour window width are compared with those of a co-located water level gauge to assess the performance of the technique.

Experiment Setup
The Vianden power station that is located north of Vianden in Diekirch District, Luxembourg is a pumped storage power station that is used to store excess energy and generate peak current. The first four sets of pump generators in the Vianden pumped storage plant were commissioned in the winter of 1962/1963 [20]. The second phase of work was completed in 1964 by installing another five reversible turbine-pump generators. At times when the power consumption is low, water is pumped into the upper reservoir from the lower reservoir dam at VPSP. When the energy demand increases, the stored water is used to generate a peak current. The elevation difference between the upper and lower reservoirs is 282.8 m [20]. The VPSP hydroelectric plant generates an average of 1650 gigawatthours electricity every year. The upper reservoir with the total capacity of 7,340,000 m is equipped with ten pump-turbine generators. The 1st to 9th reversible turbine-pump generators with the turbine-runner diameter of 2.45 m can discharge 39.5 m /s of water when electricity is generated. The 10 th machine with a diameter of 4.40 m has a water flow of 77 m /s and 74 m /s in the turbine operation and pumping mode, respectively [20]. The 11 th pump generator with the turbine capacity of 200 megawatts was commissioned in December 2014.
The GNSS station VPSP (latitude: 49.94°, longitude: 6.18°, altitude: 561.96 m) was installed in late-June 2017 at VPSP. It was equipped with a Septentrio PolaRx5 receiver and a Trimble TRM59800.00 antenna with no external radome ( Figure 1b). The antenna is installed ~11.0 m above the mean water level, and the receiver was set at a 1-Hz sampling rate. As the station VPSP was installed specifically for GNSS-R, it is positioned towards the south direction in order to maximize the number of satellite tracks (Figure 1c).
VPSP is also equipped with Rittmeyer high-precision pressure gauge MPW2Q sensors. MPW2Q provides hydrostatic or pneumatic pressure measurements every 1 min. MPW2Q reports highprecision water levels (WL) with a resolution better than 1 mm [21].
An example of SNR observations of a rising arc at VPSP station is shown in Figure 2. The protected/encrypted signals (P(Y)) on both the L1 and L2 bands are ignored, as they produce nosier SNR fringes [12]. It can be seen that SNR fringes are very clear for all the seven signals, especially for low-elevation angle observations, i.e., observations below a 20° elevation angle. Second, GS5Q have the highest direct power, as no military signal is planned for L5 [11]. Third, there are up to 10-dB different direct power levels between the strongest (GS5Q) and the weakest (RS2P) signals.

Experiment Setup
The Vianden power station that is located north of Vianden in Diekirch District, Luxembourg is a pumped storage power station that is used to store excess energy and generate peak current. The first four sets of pump generators in the Vianden pumped storage plant were commissioned in the winter of 1962/1963 [20]. The second phase of work was completed in 1964 by installing another five reversible turbine-pump generators. At times when the power consumption is low, water is pumped into the upper reservoir from the lower reservoir dam at VPSP. When the energy demand increases, the stored water is used to generate a peak current. The elevation difference between the upper and lower reservoirs is 282.8 m [20]. The VPSP hydroelectric plant generates an average of 1650 gigawatt-hours electricity every year. The upper reservoir with the total capacity of 7,340,000 m 3 is equipped with ten pump-turbine generators. The 1st to 9th reversible turbine-pump generators with the turbine-runner diameter of 2.45 m can discharge 39.5 m 3 /s of water when electricity is generated. The 10th machine with a diameter of 4.40 m has a water flow of 77 m 3 /s and 74 m 3 /s in the turbine operation and pumping mode, respectively [20]. The 11th pump generator with the turbine capacity of 200 megawatts was commissioned in December 2014.
The GNSS station VPSP (latitude: 49.94 • , longitude: 6.18 • , altitude: 561.96 m) was installed in late-June 2017 at VPSP. It was equipped with a Septentrio PolaRx5 receiver and a Trimble TRM59800.00 antenna with no external radome ( Figure 1b). The antenna is installed~11.0 m above the mean water level, and the receiver was set at a 1-Hz sampling rate. As the station VPSP was installed specifically for GNSS-R, it is positioned towards the south direction in order to maximize the number of satellite tracks ( Figure 1c).
VPSP is also equipped with Rittmeyer high-precision pressure gauge MPW2Q sensors. MPW2Q provides hydrostatic or pneumatic pressure measurements every 1 min. MPW2Q reports high-precision water levels (WL) with a resolution better than 1 mm [21].
An example of SNR observations of a rising arc at VPSP station is shown in Figure 2. The protected/encrypted signals (P(Y)) on both the L1 and L2 bands are ignored, as they produce nosier SNR fringes [12]. It can be seen that SNR fringes are very clear for all the seven signals, especially for low-elevation angle observations, i.e., observations below a 20 • elevation angle. Second, GS5Q have the highest direct power, as no military signal is planned for L5 [11]. Third, there are up to 10-dB different direct power levels between the strongest (GS5Q) and the weakest (RS2P) signals.

GNSS-R Data Processing
SNR observations can be decomposed into a trend ( ) that is a monotonically increasing function of the elevation angle ( ) and detrended interference fringes ( ) . The detrended interference fringes that are superimposed on the trend can be expressed in terms of the detrended interference fringes' amplitude ( ) and reflection excess phase ( ) with respect to the direct phase [22]: The interferometric phase is a function of the coupled surface-antenna properties ( ) , geometric interference component ( ) , and right-handed circularly polarized (RHCP) antenna phase pattern evaluated in the direct path ( ) [23]: The compositional component of the interferometric phase can be expressed as a function of the complex-valued Fresnel reflection coefficients for same-( ) and cross-sense circular polarization ( ) and antenna complex vector effective length ( ) [11,22]: The geometric component of the interferometric phase can be modeled as a function of the wavenumber ( ): Assuming a flat and smooth surface, the interferometric propagation delay ( ) can be modeled as a function of the satellite elevation angle and an a priori reflector height ( )-the vertical distance between the GNSS antenna phase center and the reflecting surface [11]:

GNSS-R Data Processing
SNR observations can be decomposed into a trend (tSNR) that is a monotonically increasing function of the elevation angle (e) and detrended interference fringes (dSNR). The detrended interference fringes that are superimposed on the trend can be expressed in terms of the detrended interference fringes' amplitude (aSNR) and reflection excess phase (φ i ) with respect to the direct phase [22]: The interferometric phase is a function of the coupled surface-antenna properties (φ X ), geometric interference component (φ I ), and right-handed circularly polarized (RHCP) antenna phase pattern evaluated in the direct path (Φ R d ) [23]: The compositional component of the interferometric phase can be expressed as a function of the complex-valued Fresnel reflection coefficients for same-R S and cross-sense circular polarization R X and antenna complex vector effective length (L) [11,22]: The geometric component of the interferometric phase can be modeled as a function of the wavenumber (k): Assuming a flat and smooth surface, the interferometric propagation delay (τ i ) can be modeled as a function of the satellite elevation angle and an a priori reflector height (H)-the vertical distance between the GNSS antenna phase center and the reflecting surface [11]: Remote Sens. 2020, 12, 3614

of 12
As the GNSS antenna surroundings (i.e., water level) are an unknown condition, the GNSS-R inverse model is parameterized by a few biases assuming a priori values for the physical parameters. Therefore, the tSNR and dSNR in the electromagnetic forward model are defined as tSNR K −2 + K −2 B 2 and dSNR K −2 |B| −1 , respectively. The real-valued noise power K and complex-valued interferometric bias B with magnitude |B| 2 = 10 B dB /10 and phase φ B are introduced for the interferometric model as [12]: where the superscript in parentheses is an index; otherwise, it is a power exponent. While the direct and interferometric power biases are introduced only to improve the sinusoid fit to the SNR observations, the interferometric phase bias is used to retrieve the reflector height. The inverse model can use the forward model internally to simulate SNR observations to estimate parameter corrections by fixing them to their a priori values and corresponding covariance matrices. For more details, the reader is referred to [12] and [23]. Therefore, the interferometric phase can be updated in the forward model as: Finally, the linear phase bias can be used to estimate the total reflector height: where λ is the carrier wavelength. After inversion, a number of post-processing measures are applied to detect and reject outliers resulting from unexpected situations (i.e., secondary reflections, diminished interferometric powers) of each individual GNSS signal. However, SNR observations should be first clustered so retrieved biases will be compared under comparable environmental conditions. To do that, SNR observations are partitioned into satellite tracks, between rising/setting near the horizon, and culmination near the zenith. The satellite tracks are then clustered by taking advantage of GNSS ground tracks repeatability. Therefore, satellite tracks are discriminated by whether the satellite is ascending or descending near the same azimuth. Intra-signal and intra-cluster quality control (QC) compare the estimated linear phase biases within each individual cluster using robust estimators and statistical intervals and reject outliers resulting from anomalous conditions [9]. In the next step, as phase biases are highly correlated to one another, the constrained least-squares (LS) adjustment is implemented to retrieve constrained constant phase biases to improve the precision of the water level estimates [6].
The water level at VPSP abruptly changes with time. Therefore, the estimated water levels must be corrected for the dynamic water surface level [16]. The interferometric phase rate can be defined as a function of the stationary surface height (H) and height rate correction (∂H/∂k z ): where k z = 2k sin e is the vertical wavenumber [6]. The vertical velocity of the water surface level is an important correction, as the water level at VPSP varies up to 17 m every 4-8 h. Therefore, the height-rate correction can be modeled as: H tan e/ . e + ..
H∆t tan e/ . e The first term on the right-hand side represents the angular delay rate, and the second is the integrated delay rate. The height rate and velocity rate can be solved iteratively using the centered finite differences [6]. The elevation angel rate (∂e/∂t) is derived from the GNSS precise ephemerides. The authors of [6] reported a height rate and velocity rate of about 40 cm/h and 20 cm/h 2 , respectively, for a GNSS station with about 3-m daily tidal envelopes. Therefore, the height-rate correction cannot be neglected especially at VPSP; otherwise, GNSS-R water levels would be biased. It should be noted that vertical acceleration would be also beneficial, especially when sub-hourly water levels (10 min) with such a window width (2 h) are reported [6].
Tropospheric delays as a function of the satellite elevation angle and the reflector height were found as an important source of error in GNSS-R studies [19,24]. Therefore, here, the estimated water levels are corrected for atmospheric refraction after intra-signal and intra-cluster QC. The Vienna Mapping Function 3 (VMF3) that is based on the ray-traced delays of the ray-tracer RADIATE, together with the Global Pressure and Temperature 3 model (GPT3) on a 1 • × 1 • global grid that represents a comprehensive troposphere model [25], is used to estimate the tropospheric delay corrections. Using the minimum and maximum elevation angle observations of each individual satellite track, the along-path tropospheric delays can be estimated as: where n is the index of refraction along the instantaneous excess delay path [19,24]. All cluster-by-cluster tracks are combined to produce a regularly spaced signal-specific epoch-by-epoch water level time series together with their respective uncertainties based on a weighted moving average. The moving average returns water levels every 10 min using a 2-h window width [6]. It should be noted that the moving regression is not a simple moving average, as it also estimates the height and velocity rates. Finally, sub-hourly water level time series are reported using the signal-specific variance factors. The statistical inter-signal combination that offers a single multi-GNSS time series was found to improve the precision of geodetic GNSS-R retrievals [12]. The post-inversion algorithm of the water level retrieval at VPSP is summarized in Figure 3.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 13 The first term on the right-hand side represents the angular delay rate, and the second is the integrated delay rate. The height rate and velocity rate can be solved iteratively using the centered finite differences [6]. The elevation angel rate ( ⁄ ) is derived from the GNSS precise ephemerides. The authors of [6] reported a height rate and velocity rate of about 40 cm h ⁄ and 20 cm h ⁄ , respectively, for a GNSS station with about 3-m daily tidal envelopes. Therefore, the height-rate correction cannot be neglected especially at VPSP; otherwise, GNSS-R water levels would be biased. It should be noted that vertical acceleration would be also beneficial, especially when sub-hourly water levels (10 min) with such a window width (2 h) are reported [6]. Tropospheric delays as a function of the satellite elevation angle and the reflector height were found as an important source of error in GNSS-R studies [19,24]. Therefore, here, the estimated water levels are corrected for atmospheric refraction after intra-signal and intra-cluster QC. The Vienna Mapping Function 3 (VMF3) that is based on the ray-traced delays of the ray-tracer RADIATE, together with the Global Pressure and Temperature 3 model (GPT3) on a 1°× 1° global grid that represents a comprehensive troposphere model [25], is used to estimate the tropospheric delay corrections. Using the minimum and maximum elevation angle observations of each individual satellite track, the along-path tropospheric delays can be estimated as: where is the index of refraction along the instantaneous excess delay path [19,24]. All cluster-by-cluster tracks are combined to produce a regularly spaced signal-specific epochby-epoch water level time series together with their respective uncertainties based on a weighted moving average. The moving average returns water levels every 10 min using a 2-h window width [6]. It should be noted that the moving regression is not a simple moving average, as it also estimates the height and velocity rates. Finally, sub-hourly water level time series are reported using the signalspecific variance factors. The statistical inter-signal combination that offers a single multi-GNSS time series was found to improve the precision of geodetic GNSS-R retrievals [12]. The post-inversion algorithm of the water level retrieval at VPSP is summarized in Figure 3.  Before we combined all the independent GNSS-R water levels to provide a single multi-GNSS water level time series, first we compared all the signal-dependent water levels together. Here, the GS1C water level time series is used to evaluate the inter-signal consistency among all the GNSS-R retrievals. A systematic error at a decimeter level between the GPS-L1 and GPS-L2 signal water-surface level retrievals was reported in [26]. In this study, we found linear relationships between different GNSS-R water level estimates and GS1C, with correlations that exceeded 0.98 (Figure 4).

Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 13
Before we combined all the independent GNSS-R water levels to provide a single multi-GNSS water level time series, first we compared all the signal-dependent water levels together. Here, the GS1C water level time series is used to evaluate the inter-signal consistency among all the GNSS-R retrievals. A systematic error at a decimeter level between the GPS-L1 and GPS-L2 signal watersurface level retrievals was reported in [26]. In this study, we found linear relationships between different GNSS-R water level estimates and GS1C, with correlations that exceeded 0.98 (Figure 4).  Figure 5 shows the tropospheric delays at station VPSP. Due to the variations of the watersurface level, the atmospheric delay changes with respect to the satellite elevation angle. The tropospheric delay difference between the two modes of VPSP when the water-surface level is at its minimum and its maximum is 29.40 cm. The regression slope 0.984 m m ⁄ ± 0.001 m m ⁄ between the GNSS-R water level time series and WL records are statistically significant vis-à-vis its 95% confidence interval if the tropospheric delay corrections are neglected.  Figure 5 shows the tropospheric delays at station VPSP. Due to the variations of the water-surface level, the atmospheric delay changes with respect to the satellite elevation angle. The tropospheric delay difference between the two modes of VPSP when the water-surface level is at its minimum and its maximum is 29.40 cm. The regression slope 0.984 m/m ± 0.001 m/m between the GNSS-R water level time series and WL records are statistically significant vis-à-vis its 95% confidence interval if the tropospheric delay corrections are neglected.

Results
SNR observations from a GNSS station in Sweden, [17] demonstrated that GS2W and RS2P should be avoided for GNSS-R water level studies. The authors of [12] found that legacy GPS signals are worse in terms of accuracy, while modernized GPS and GLONASS signals are comparable. Therefore, signal-specific variance factors are used to report multi-GNSS GNSS-R time series, together with their formal errors ( Figure 6). The GNSS-R water level time series are compared with those of the WL records. First, the WL records are interpolated at the GNSS sampling epochs and vice versa by disabling the interpolations across the data outage periods. The RMSE error between the sub-hourly GNSS-R time series and WL records is 7.0 cm. The agreement between the WL records and GNSS-R water level retrievals is excellent, with a correlation coefficient that exceeds 0.999 (Figure 7). SNR observations from a GNSS station in Sweden, [17] demonstrated that GS2W and RS2P should be avoided for GNSS-R water level studies. The authors of [12] found that legacy GPS signals are worse in terms of accuracy, while modernized GPS and GLONASS signals are comparable. Therefore, signal-specific variance factors are used to report multi-GNSS GNSS-R time series, together with their formal errors ( Figure 6). The GNSS-R water level time series are compared with those of the WL records. First, the WL records are interpolated at the GNSS sampling epochs and vice versa by disabling the interpolations across the data outage periods. The RMSE error between the sub-hourly GNSS-R time series and WL records is 7.0 cm. The agreement between the WL records and GNSS-R water level retrievals is excellent, with a correlation coefficient that exceeds 0.999 (Figure 7). There is no scale error between the GNSS-R water level time series and WL records that can be seen in the van de Casteele diagram in Figure 8. The x-axis presents the difference between geodetic GNSS-R water levels as the test sensor, and the y-axis represents the WL records as the reference sensor. The two water levels deviate by −0.27% from a one-to one relationship. The authors of [4] reported a scale error of 0.84% for a GNSS site with ~2.5-m tidal amplitudes. We found that the RMS reduction at VPSP is 49% when the vertical acceleration is applied. In addition, the GNSS-R water level time series would be deviated by 1.97% from the WL records if they were only corrected for the height rate. There is no scale error between the GNSS-R water level time series and WL records that can be seen in the van de Casteele diagram in Figure 8. The x-axis presents the difference between geodetic GNSS-R water levels as the test sensor, and the y-axis represents the WL records as the reference sensor. The two water levels deviate by −0.27% from a one-to one relationship. The authors of [4] reported a scale error of 0.84% for a GNSS site with~2.5-m tidal amplitudes. We found that the RMS reduction at VPSP is 49% when the vertical acceleration is applied. In addition, the GNSS-R water level time series would be deviated by 1.97% from the WL records if they were only corrected for the height rate. We also compared the sub-hourly maximum and minimum GNSS-R water level retrievals and water level gauge records ( Figure 9). The RMSE between the maximum GNSS-R retrievals and water level records is 6.5 cm, but it is 11.1 cm between the minimum GNSS-R estimates and water level measurements. We also compared the sub-hourly maximum and minimum GNSS-R water level retrievals and water level gauge records ( Figure 9). The RMSE between the maximum GNSS-R retrievals and water level records is 6.5 cm, but it is 11.1 cm between the minimum GNSS-R estimates and water level measurements. the median offset between the maximum measurements is 0.5 cm, it is 4.7 cm for the minimum retrievals. On one hand, GNSS-R retrievals are less precise than those of the water level. On the other hand, the resolution of the Rittmeyer high-precision pressure gauge sensor is 10 times higher than the GNSS-R, and GNSS-R retrievals are averaged over the FFZs over 10 min. In addition, the offset error between minimum retrievals can be the result of the ramp slope of the dam. When electricity is generated, the water level goes down. Therefore, the vertical and horizonal distances between the GNSS antenna and water surface increase. There would be reflections from the ramp slope for higher elevation angle observations where FFZs are close to the antenna.

Conclusions
The geodetic GNSS-R is used to estimate water levels at the Vianden pumped storage plant (VPSP), where the water-surface level changes vertically ~17 m several times a day to generate a peak current. The retrieved heights are corrected for tropospheric delay using VMF3/GPT3 models before signal-specific water levels are estimated based on a weighted moving average using a 10-min postspacing and a 2-h averaging window width. The inter-cluster combination is implemented in a way to take the dynamic state of the sea surface into account by correcting the retrieved heights for the velocity and vertical velocity corrections. Finally, GNSS-R sub-hourly time series are reported using the signal-specific variance factors. The standard deviation between the sub-hourly GNSS-R water levels and co-located tide gauge records for about one year at the VPSP is 7.0 cm, with a correlation coefficient exceeding 0.999. Carrying out a regression between maxima measurements from the two sensors, the intercept −0.03 m ± 0.16 m and regression slope 0.99 m/m ± 0.03 m/m are not statistically significant vis-à-vis their 95% confidence intervals. The regression intercept of the minima GNSS-R retrievals and water level records is −0.02 m ± 0.19 m, and the regression slope is 1.01 m/m ± 0.02 m/m. While the median offset between the maximum measurements is 0.5 cm, it is 4.7 cm for the minimum retrievals. On one hand, GNSS-R retrievals are less precise than those of the water level. On the other hand, the resolution of the Rittmeyer high-precision pressure gauge sensor is 10 times higher than the GNSS-R, and GNSS-R retrievals are averaged over the FFZs over 10 min. In addition, the offset error between minimum retrievals can be the result of the ramp slope of the dam. When electricity is generated, the water level goes down. Therefore, the vertical and horizonal distances between the GNSS antenna and water surface increase. There would be reflections from the ramp slope for higher elevation angle observations where FFZs are close to the antenna.

Conclusions
The geodetic GNSS-R is used to estimate water levels at the Vianden pumped storage plant (VPSP), where the water-surface level changes vertically~17 m several times a day to generate a peak current. The retrieved heights are corrected for tropospheric delay using VMF3/GPT3 models before signal-specific water levels are estimated based on a weighted moving average using a 10-min post-spacing and a 2-h averaging window width. The inter-cluster combination is implemented in a way to take the dynamic state of the sea surface into account by correcting the retrieved heights for the velocity and vertical velocity corrections. Finally, GNSS-R sub-hourly time series are reported using the signal-specific variance factors. The standard deviation between the sub-hourly GNSS-R water levels and co-located tide gauge records for about one year at the VPSP is 7.0 cm, with a correlation coefficient exceeding 0.999.
The results confirm that GNSS-R can detect sudden and large water surface changes. Therefore, geodetic GNSS-R can be used to study water-surface changes caused by nonastronomical forces such as hurricanes or storm surges. The performance of GNSS-R water levels can be improved in terms of