Identifying 2010 Xynthia Storm Signature in GNSS-R-Based Tide Records

Abstract

In this study, three months of records (January–March 2010) that were acquired by a geodetic Global Navigation Satellite Systems (GNSS) station from the permanent network of RGP (Réseau GNSS Permanent), which was deployed by the French Geographic Institute (IGNF), located in Socoa, in the south of the Bay of Biscay, were used to determine the tide components and identify the signature of storms on the signal to noise ratio (SNR) during winter 2010. The Xynthia storm hit the French Atlantic coast on the 28th of February 2010, causing large floods and damages from the Gironde to the Loire estuaries. Blind separation of the tide components and of the storm signature was achieved while using both a singular spectrum analysis (SSA) and a continuous wavelet transform (CWT). A correlation of 0.98/0.97 and root mean square error (RMSE) of 0.21/0.28 m between the tide gauge records of Socoa and our estimates of the sea surface height (SSH) using the SSA and the CWT, respectively, were found. Correlations of 0.76 and 0.7 were also obtained between one of the modes from the SSA and atmospheric pressure from a meteorological station and a mode of the SSA. Particularly, a correlation reaches to 0.76 when using both the tide residual that is associated to surges and atmospheric pressure variation.

1. Introduction

The sea level rise (SLR), which will be caused by global warming, will impact low-lying coastal areas, not only through inundation, but also as an increase of storm surges repeatability and extreme astronomic tides will reach higher water levels [1,2,3]. Changes in climatic conditions are also likely to increase the elevation and the frequency of storm surges in various areas [4,5,6], causing a major threat to the increasing part of the world population that live in coastal regions within a few meters above sea level (SL) [7].

A surge is defined as the difference between the observed and astronomical modeled tides at a specific location. This difference is due to meteorological phenomena, like storms [8,9]. When there is a combination of high tide and storm events that are characterized by large rainfall and strong waves with important setup, surges can be responsible for the inundation of coastal areas and damage of coastal facilities [10,11].

If we consider the French case, this country has experienced many strong storms in recent years, such as Oritia storm in 2000, Cyclone Kyrill in 2007, European tornado outbreak in 2008, Xynthia storm in 2010, Joachim winter storm in 2011, tidal surge Xavier in 2013, Darwin storm in 2014, extratropical cyclone Zeus in 2017, and David storm in 2018. This series of extreme events put forward the urgent need for storm monitoring and an early warning system along the French/European coasts. Geodetic Global Navigation Satellite Systems (GNSS) stations that were located on top of cliffs can be used for this purpose. During the night of 27–28 February 2010, the Xynthia storm affected the south western part of Europe, mainly the Spanish and French coastal regions of Bay of Biscay, causing important damages, 59 casualties in Europe (47 in France), and a large marine submersion [12,13,14,15,16], as in the La Faute sur Mer area with 29 causalities. Winds up to 160 km.h⁻¹ and pressure down from ~1000 to 977 hPa that were recorded by meteorological stations during high spring tide were responsible for a huge storm surge along the coast of the Bay of Biscay [14,17]. For example a surge of 1.53 m was recorded at La Rochelle tide gauge (8.01 m above the hydrographic zero). Such a high tide level was never recorded since the set-up of this tide gauge in 1997. The surge was also greater than the largest recorded at Brest tide gauge (1.42 m) during the last 150 years [13]. For this reason, the densification of sensors and observations is crucial in establishing a well-structured surveillance and warning system to ensure the safety of populations. Currently, tide gauges that now used radar technique to measure the tides along the French coast ensured long-term monitoring. They are more and more co-located with either GNSS geodetic receivers or Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) space geodesy system to separate the changes in Sea Surface Height (SSH) from crustal motions (e.g., Wöppelmann et al., 2006 [18]).

The Global Navigation Satellite Systems-reflectometry (GNSS-R) technique has demonstrated a strong potential for the monitoring of sea surface variations since the mid-90s based on waveform analysis [19,20,21]. Nowadays, for monthly to interannual monitoring, SSH variations are mostly derived from the processing of the Signal-to-Noise Ratio (SNR) that was acquired by a geodetic GNSS receiver [22,23,24,25]. GNSS geodetic stations that are part of the national and international permanent GNSS networks and located in coastal areas can be used as “opportunistic” tide gauge, according to Larson et al. (2013), when they record the SNR data. Contrary to classical tide gauges that are located in protected areas to dampen the effect of the waves on the signal, the GNSS-R based time-series of Sea Surface Height (SSH) contains both of the signatures of the tide and other information, like the waves [26], as well as the surges signatures.

In this study, the SNR data records from the Socoa geodetic station (south west of France) were used to determine the SSH variations from January to March 2010 in the Saint Jean de Luz Bay. Two methods (Singular Spectrum Analysis—SSA and Continuous Wavelet Transform—CWT) were applied to the SSH derived from SNR records to separate the tide signature from other geophysical signals. Comparisons were performed against the tide gauge records from the Socoa station. Non-tidal signals were compared to surge estimates, pressure variations, and waves to determine their nature. The potential of this technique, the choice of the location of the geodetic GNSS station, and the choice of the statistical techniques are then discussed.

2. Study Area

The analyses that are presented below were performed in the Bay of Saint de Luz, which is located in the east of the Bay of Biscay, along the French Atlantic coast, a few kilometers north from the Spanish border (Figure 1a). For several centuries, powerful storms strongly affect this place and it has been regularly flooded. In the middle of the 19th century, a breakwater was constructed in order to protect the area from the ocean’s tides. The study area is located between the breakwater and the coastline of the bay, called hereafter inside part, where the wave effects are attenuated. The annual mean significant wave height in the inside part is 1.6 m with a maximum of 2.5 m [27,28].

The inside part of the bay is approximately 2 km long by 1 km wide. Owing to the protection of the breakwaters, this part is almost unaffected by the coastal currents and waves [29], with a tidal range of about 4.5 m at spring tides. The inside part also receives freshwater inflows from two small rivers (Untxin and Nivelle rivers in Figure 1a). The study area is in a semi-diurnal macro-tidal environment [30,31].

3. Datasets

3.1. GNSS Data

GNSS data, including SNR, have been continuously acquired since the 1st of April 2006 at a frequency of 1 Hz while using a Trimble NETR5 receiver (from 2009 until 2010) and a Trimble Zephyr model 2 (TRM55971.00) geodetic antenna in Socoa (43°23′42.83″N, 01°40′54.05″W) at 10.66 m above the mean sea surface (Figure 1b). This station (SCOA) is part of RGP (Réseau GNSS Permanent—http://rgp.ign.fr/) from IGN (Institut National de l’Information Géographique et Forestière). In this study, the data used were acquired at L1 frequency from GPS and GLONASS constellations from the 1st of January to 31 March 2010.

3.2. Tide Gauge Data

SSH records from the Socoa gauge station that was operated by the French REFMAR network (http://data.shom.fr/donnes/refmar/) were used for comparison with SSH estimates, derived from the SNR technique, from January to March 2010. SSH measurements are available with a temporal resolution of 10 min. The instrumental zero (of this tide gauge) is set to the hydrographic zero. The Socoa tide gauge is located (43°23′42.86″N; 01°40′53.83″W) at the end of Socoa’s breakwaters, protecting it against wave effects, 5 m away from the GNSS station. The Socoa station is equipped with a digital coastal tide gauge that is composed of an Optiwave 7300C sensor and a MARELTA acquisition unit.

3.3. Meteorological Data

The Socoa meteorological station provides hourly estimates of rainfall, wind speed, and atmospheric pressure. Records from the period January–March 2010 were used in this study. Infoclimat (https://www.infoclimat.fr) made them available. The atmospheric pressure was converted into inverted barometer effect ($h_{atmos})$ using the following equation:

$$h_{atmos}= -\frac{P_{atm}- \overline{P_{atm}}}{\mathsf{\rho }\mathrm{g}}$$

(1)

where $P_{atm}$ is the time varying mean of the global surface atmospheric pressure over the oceans; ρ = 1,020 gcm⁻³ is density of seawater; and, g = 9,81 cms⁻² is the mean acceleration of gravity [32].

3.4. Significant Wave Height Data

The directional wave buoy located at (43°31′56″N; 1°36′54″W) acquired significant wave height (SWH) data over the study period. SWH data are provided at the temporal resolution of 30 min by CANDHIS service from National swell off database from the Institute for maritime and inland waterways (CETMEF-Brest).

4. Methods

SNR GNSS data were first inverted to derive the water levels in the SCOA station using the methodology that was developed by [22]. The resulting water levels were then analyzed using two spectral methods i.e., SSA and CWT.

4.1. Inversion of the SNR Data

GNSS-R is an opportunistic technique that provides information regarding the properties of the Earth surface [33]. The inversion of the GNSS SNR can be used for the monitoring of SSH [24] while using the Interference Pattern Technique (IPT).

Following [34], the SNR quantity at any instant is described as:

$$\mathrm{SNR} = A_d^2+A_m^2+2A_d{A}_m\cos \psi$$

(2)

where $A_d$ and $A_m$ are the amplitudes of the direct and reflected signals, respectively, and φ is the phase difference between these two signals.

Assuming that $A_m$ ≪ $A_d$, thus SNR can be approximated by:

$$\mathrm{SNR} \approx A_d^2 + 2A_d{A}_m\cos \psi$$

(3)

The reflected signal perturbations are mainly visible for low satellite elevation angles [35], but recent development demonstrate the possibility to also use high satellite elevation [26]. Assuming a planar reflector that corresponds to the sea surface, the relative phase angle can be geometrically derived from the path delay δ of the reflected signal [36]:

$$\psi = \frac{2\pi }{\lambda }\mathsf{\delta } = \frac{4\pi h}{\lambda }\sin \left(\epsilon \right)$$

(4)

where λ the signal wavelength, ε the satellite elevation, and h the vertical distance between the antenna phase center and the reflecting surface. The frequency of the multipath oscillations can be derived from (3), as proposed by [37] and improved by [22]:

$$f_{\psi } = \frac{d\psi }{dt} = \frac{4\pi \dot{h}}{\lambda }\sin \left(\epsilon \right) + \frac{4\pi h}{\lambda }\cos \left(\epsilon \right)\dot{\epsilon }$$

(5)

$\dot{h}$ (=$\frac{dh}{dt}$) defines the vertical velocity and $\dot{\epsilon }$ (=$\frac{d\epsilon }{dt}$) defines the elevation angle velocity. Equation (4) can be simplified by making the following simple change of variable—x = $\sin \left(\epsilon \right)$:

$$\tilde{f} = \frac{d\psi }{dx} = \frac{4\pi }{\lambda }\left(\dot{h} \frac{\tan \left(\epsilon \right)}{\dot{\epsilon }}+h\right)$$

(6)

where $\tilde{f}$ is the frequency of the multipath oscillation.

For the static case that was proposed by [24]: $\dot{h}$ ≈ 0, thus $\tilde{f}$ is constant and is directly proportional to the receiver height above the reflecting surface.

In the dynamic case, $\dot{h}$ ≠ 0, so the time series and $\tilde{f}$(t) is proportional to the vertical antenna height velocity $\dot{h}$ and to the satellite elevation ε and its angular velocity $\dot{\epsilon }$. When considering Equation 5, only two unknown quantities (h and $\dot{h}$) have to be determined while using the Least Square Method (LSM) [22].

In this study, we applied the approach that was proposed by [22] to the 1 Hz L2C SNR from GPS and GLONASS satellites. The update interval was fixed to 10 min for height retrievals. A maximum range of variations of ±4 m from the elevation above mean sea surface (10.66 m) was considered for the inversion algorithm. Due to the configuration of the site, all reflections from azimuth angles between 60 and 210° were not taken into account to only keep the reflections that occur inside of the bay. Besides, due to the presence of mask and multipath due to buildings and dikes, only reflections from elevation angles ranging from 1° to 25° were considered.

4.2. Analysis of the GNSS-R-Based Water Levels

4.2.1. Singular Spectrum Analysis

SSA is a principal component analysis in the time domain, which is used to extract statistical patterns from short and noisy time-series without any a priori on the dynamics of the signal [38,39]. When considering a signal $x\left(t\right)$ composed of N samples taken in a window of size M, a lagged autocorrelation matrix is obtained for a maximum lag M equals the window size. The elements ($c_j$) of this matrix are defined, as follows:

$$c_j=\frac{1}{N-j}{{\displaystyle \sum }}_{i=1}^{N-j}x(t_i)x(t_{i+j}) \mathrm{with} 0 \le \mathrm{j} \le \mathrm{M}-1$$

(7)

where $t_i$ and $t_{i+1}$ are the sampling times.

Eigenvalues ${\lambda }_{\mathrm{k}}$ and eigenvectors $E_j^k$ (also known as empirical orthogonal functions) of the lagged autocorrelation matrix are determined and then sorted in descending order of the eigenvalues. Following [39,40], the $k^{th}$ principal components among M is given by:

$$a^k(t_i)={{\displaystyle \sum }}_{j=1}^{M}x(t_{i+j})E_j^k \mathrm{with} 0 \le \mathrm{i} \le \mathrm{N}-\mathrm{M}$$

(8)

The components of the original time-series analyzed using SSA can be reconstructed, as follows:

$$x^k(t_i)=\frac{1}{M}{{\displaystyle \sum }}_{j=1}^{M}a(t_{i-j})E_j^k \mathrm{with} \mathrm{M} \le \mathrm{i} \le \mathrm{N}-\mathrm{M}+1$$

(9)

where $x^k\left(t_i\right)$ is the $k^{th}$ reconstructed components (RC) of the original times series.

In this study, the time series are decomposed into principal components (PCs) and then reconstructed components (RCs) while using a modified version of the MATLAB algorithm [41].

4.2.2. Continuous Wavelet Transform

The CWT of a signal $x\left(t\right)$, is defined as the sum over all time of the signal multiplied by scaled, shifted versions of the wavelet function $\mathsf{\psi }\left(\mathrm{t}\right)$ [42]:

$$C_x\left(a,\tau \right)=\underset{-\infty }{\overset{+\infty }{{\displaystyle \int }}}x\left(t\right){\psi }_{a,\tau }^{*}dt ; {\psi }_{a,\tau }\left(t\right)=\frac{1}{a}\psi \left(\frac{t-\tau }{a}\right)$$

(10)

The results of the CWT are many wavelet coefficients $C_x\left(a,\tau \right)$, where * corresponds to the conjugated complex. The two parameters a and τ correspond, respectively, to the scale factor and the temporal translation (shift). Afterward, a wavelet spectrum $W_X\left(a,\tau \right)$ is constructed according to the relation below, defined as the modulus of the wavelet coefficients [43]:

$$W_X\left(a,\tau \right)=C_X\left(a,\tau \right)C_X^{*}\left(a,\tau \right)={\left|C_X\left(a,\tau \right)\right|}^2$$

(11)

In order to investigate the relationship in time-scale space between two given signals, the cross wavelet transform (XWT) is used. The XWT that was constructed from two CWTs exposes their common power and it provides further information concerning their relative phase. The wavelet cross spectrum $W_{XY}\left(a,\tau \right)$ between two signals $x\left(t\right)$ and $y\left(t\right)$ is given by:

$$W_{XY}\left(a,\tau \right)=C_X\left(a,\tau \right)C_Y^{*}\left(a,\tau \right)$$

(12)

where $C_X\left(a,\tau \right)$ and $C_Y^{*}\left(a,\tau \right)$ are, respectively, the wavelet coefficient of $x\left(t\right)$ and the conjugate of the wavelet coefficient of $y\left(t\right)$.

In this study, wavelet analyzes were performed while using a modified version of the cross wavelet and wavelet coherence toolbox from [44].

5. Results

In this study, two statistical approaches were applied to blindly decompose GNSS-R-based time series of SSH from Socoa: SSA and CWT. Our goal is to separate the major tide components from other geophysical signals, including the signature of the extreme storm event, as the Xynthia storm occurred during the study period.

5.1. GNSS-R SSH Time Series Analysis Using the SSA Method

The four first components (RC) of the SSA account for more than 95% of the explained variance. Figure 2 presents GNSS-R SSH and the fourth first components of the SSA.

RC1 and RC2 accounted for 47.2% and 44.5% of the explained variance, respectively, whereas RC3 and RC4 accounted for 2.7% and 1.0%, respectively. RC1 and RC2 exhibit oscillations of large amplitude (between 2 and 4 m) of frequency lower, ranging from less than a day to a quasi-monthly modulation, whereas RC3 and RC4 exhibit a less regular high frequency pattern of lower amplitude. As the study area is a semi-diurnal macro-tidal environment [31,45], the SSH is dominated by the tides signal; the major RC modes represent the tides. Comparisons between the original GNSS-R SSH, RC1, RC2, and sum of RC1 and RC2 modes and the tide gauge records at Socoa station were performed. Their results are presented in terms of root mean square error (RMSE) (Figure 3a) and correlation coefficient (R) (Figure 3b) for different window sizes (M). A good agreement is generally obtained for the original GNSS-R SSH, RC1, RC2, and sum of RC1 and RC2 modes and all window sizes. Nevertheless, it can be noticed that the agreement is lower with the tide gauge record when separately considering RC1 and RC2, being especially visible for the RMSE (Figure 3a) than considering the GNSS-R SSH and the sum of RC1 and RC2 modes. However, higher correlations and lower RMSE are observed for the sum of RC1 and RC2 modes than for GNSS-R SSH for all window sizes with an optimum at 6 h (R = 0.99 and RMSE = 0.16 m). Due to this very good agreement, the RC1 and RC2 modes correspond to the tide components of the GNSS-R SSH.

During storms, SSH is significantly affected from other effects, like atmospheric pressure (called inverted barometer—IB), wind [32,46], waves, which cause surges. The surge, as determined from the tide gauge records using T-Tide, the atmospheric pressure, converted into SSH ($h_{atmos})$ using the inverted barometer (IB) equation, and the wind speed (i.e., the latter two measured at the Socoa weather station) were compared to the non-tidal component i.e., the RC3 mode of the SSA. The result of the comparison between RC3 components with the wind speed is shown in Figure 4a, $h_{atmos}$ in Figure 4b, the surge in Figure 4c, and the SWH (measured at the Anglet buoy about 20 km away from the Socoa station) in Figure 4d. A good agreement was found between RC3 and both the atmospheric pressure (R = 0.70) and the surge (R = 0.72) over the whole study period. Low correlations with the wind speed and SWH were obtained (R equals to 0.53 and 0.53, respectively). These low correlations with the wind speed and SWH could appear surprising. However, in the inner part of the bay, the presence of the dike is responsible for a large decrease of the wave amplitude that is caused by the wind. If we focus our analysis on the four days of the Xynthia storm (from 27 February to 2 March 2010), the correlation increases with three environmental variables considered here: R equals 0.73, 0.77, and 0.65 for $h_{atmos}$, the surge, and the wind speed, respectively, as compared with periods before (R = 0.62, 0.70, 0.40, respectively) and after the storm (R = 0.74, 0.71, and 0.55, respectively). On the contrary, a better correlation between RC3 and SWH is found after (R = 0.64) than during (R = 0.54) Xynthia. The larger correlation and more constant over time correlation obtained between the surge and the third reconstructed mode (RC3) of the SSA showed that RC3 and surge are strongly related. The increase in correlation with the two other variables during the Xynthia storm event and also after can be accounted for the large and long-term impact of the surge on the SSH.

5.2. GNSS-R SSH Time Series Analysis Using CWT Method

A continuous wavelet analysis was also applied to identify the period and temporalities of the most energetic signals that are present in the time series of GNSS-R SSH from January to March 2010. The time-frequency diagram of the CWT clearly exhibits a large magnitude that is centered at the period of 12 h during the whole study period (Figure 5). As for the SSA, the CWT provides a clear identification of the main tide harmonics (semi-diurnal) that are present along the French Atlantic coast. A secondary maximum, starting around the 21st of February and finishing around the 10th of March, is also clearly visible at a period of 16 days. This maximum corresponds to the low frequency of the Xynthia storm. Another secondary maximum, which was well centered on the Xynthia storm, from the 26th of February to the 1st of March, had a period of 4 h–8 h.

6. Discussion

6.1. The Accidental Tide Gauge and More

Our results demonstrate the strong potential of the GNSS-R technique for the monitoring of the coastal ocean. Larson et al., 2013 [37] earlier demonstrated that a GNSS geodetic station that was located close to the shore could be used as an “accidental tide gauge”. In this study, we selected one of the stations of the RGP located on the shore and it obtained very accurate estimates of tides but also surge in spite of the complex environment surrounding the station, which caused perturbation of the signal due to:

-

the presence of several buildings, dikes, which mask part of the GNSS satellite and are likely to cause parasite multi-paths (more than one reflection), and

-

the presence of boats in the bay that are likely to also cause other parasite multi-paths.

This means that:

-

many other GNSS geodetic stations from permanent networks around the world can also be used as accidental tide gauges to complement or improve the existing tide gauge networks,

-

GNSS geodetic stations offer the opportunity to record other geophysical phenomena such SWH (e.g., Roussel et al., (2015) [22]) or surge or inverted barometer (this study), and

-

if located on a better environment as the top of a hill or a mast, better accuracy can be reached (e.g., 0.12 m see [47]). The choice of the location and the altitude of deployment can be facilitated by the use of dedicated softwares, such as GNSS Reflected Signals Simulations (GRESS) [48] or the GPS tool box [49], which provide a simulation of the position of the reflection points depending on the location of the GNSS geodetic station.

6.2. The Choice of the SSA and CWT for Separating Tides from Other Geophysical Parameters

One can wonder why not simply applying a harmonic analysis to the GNSS-R SSH, such as T-Tide harmonic analysis software [50], as performed with the tide gauge records to separate the tide tides from the residual (Figure 6a). When we applied T-Tide to the GNSS-R SSH, the signature of the surge during a Xynthia storm event is clearly visible. However, the identification of the tides is less satisfactory (presence of dissymetries) than when it is applied on the “wave protected” tide gauge records (Figure 6b). This is likely to be due to the larger complexity of the GNSS-R SSH that contains the signature of different geophysical phenomena. Such dissymetries could be due to: (i) the larger SWH amplitudes at high tides than at low tides, even if attenuated inside the bay; and, (ii) the impact of the Nivelle and Untxin river discharges on GNSS-R SSH that is composed of the reflections inside bay and not punctual as the tide gauge. Other point, the tide gauge is far from these river outlets/mouths and it is less impacted by the river flows. On the contrary, the use of the SSA or the CWT allows for a blind separation of the signal in different modes that we were able to relate to geophysical signals; contrary to the use of the harmonics analysis that allows for separation between the tides and a residual.

6.3. The Complementarity Between SSA and Inverse CWT (iCWT) to Separate Tides from Other Geophysical Signals

The GNSS-R SSH not only contains the tide components, but also other geophysical signals. Direct comparison with the tide gauge records provides quite good statistical results (R = 0.96 and RMSE = 0.30 m, see

Table 1. Comparisons between in situ tide gauge records and GNSS-R SSH based data, sum of (RC1 + RC2) using the SSA method, inverse CWT at 12 h frequency, inverse CWT from 6 h to 12 h frequencies (bias, RMSE and R).

	Bias (m)	RMSE (m)	R
SSH GNSS-R	0.001	0.30	0.96
RC1 + RC2	0.003	0.16	0.99
iCWT at 12 h	0.005	0.26	0.99
iCWT from 6 h to 12 h	0.005	0.25	0.97

and Figure 7a). Much better results are obtained using the two first reconstructed modes of the SSA (R = 0.99 and RMSE = 0.16 m, see Table 1 and Figure 7b). The use of the iCWT is more complex. Only considering the main tide period (12 h) of the iCWT, the reconstructed GNSS-R SSH is less accurate (R = 0.99, but RMSE = 0.26 m) than the GNSS-R SSH from RC1 and RC2 and it also exhibits a slightly higher bias against the tide gauge records (bias = −0.005 m, see Table 1 and Figure 7c). Considering a larger range of periods (6 to 12 h) for the iCWT reconstructed GNSS-R SSH, improving the slope of the regression ~1 between the tide gauge and GNSS-R SSH for similar accuracy (R = 0.97 and RMSE = 0.25 m, see Table 1 and Figure 7d). SSA and iCWT are complementary tools for analyzing the GNSS-R SSH estimates: the first one allows for a good estimate of the tide and surge (Figure 3, Figure 4, and Figure 7a) and the second is very useful in determining the typical time periods of the different phenomena (Figure 5) to relating them with their physical properties.

6.4. Is Each SSA Mode Related to a Single Geophysical Phenomena?

One can wonder whether each SSA reconstructed mode can be attributed to a single and specific geophysical phenomenon. The first two reconstructed modes were identified as the signature of different tidal harmonics, and their sum (RC1 + RC2) was found to be in very good agreement with the tide gauge records (R = 0.99, Figure 3 and Figure 7b). However, RC3 was found to be in very good agreement with the surge (R = 0.72, Figure 5c), but also with IB or h_atmos (R = 0.70, Figure 5b). A wavelet cross-correlation (XWT) and a liner correlation between RC3 and the sum of the surge and h_atmos were performed (Figure 8a,b), respectively. Higher correlations were found over a long time period. Very high correlations are observed between the two variables for periods that are higher than two weeks. High correlations were found on the periods 4–8 days before and after the Xynthia storm and 2–4 days during the Xynthia storm. High correlations are also found on smaller periods (from 4 h to two days) during storm occurrences. The XWT confirms that the Xynthia storm was the more energetic event that occurred during the observation period. Temporal correlation between RC3 and the sum of the surge and of the IB exhibits larger values than the correlation with either the surge or the IB at all of the time scales (0.77 against 0.72 and 0.70 between January and March 2010, 0.75 against 0.7 and 0.62 before the storm, 0.82 against 0.77 and 0.73 during Xynthia, 0.73 against 0.71 and 0.74 after the storm for the sum of the surge and IB, the surge, IB, respectively). As the surge and IB are strongly related, RC3 can be seen as the signature of the storm on the GNSS-R SSH signal.

7. Conclusions

This study is the first convincing example of the use of GNSS-R technique to detect the storm signature on SSH through blind signal decomposition techniques, such as SSA and CWT. One of the modes of the SSA decomposition was related to the temporal variations of the surge and atmospheric fluctuations through IB (R = 0.77 for the study period when combining both effects). CWT allows for identifying the main periods of the different geophysical signals that are present in the GNSS-R SSH. The high (4 h to 8 h of period) and low (4 and 16 days of period) frequencies were identified in the signature of the Xynthia event.

Our study confirms that the GNSS-R approach is able to estimate SSH with an almost similar accuracy as the tide gauges (R = 0.99 and RMSE = 0.16 m on the three months of our study period) while using the Socoa GNSS geodetic station from the RGP, which is not located at any ideal location for SNR inversion (masks from surrounding buildings, presence of dikes and boats causing multi-paths, proximity to two river mouths). Its main additional value was to demonstrate the ability of the GNSS-R technique that is based on SNR inversion for storm detection. Our study shows that SSH monitoring would benefit from the use of the SNR records that were acquired by geodetic stations located along the coasts, and this secondary goal could be taken into consideration when permanent GNSS stations are installed close the shore.

In France, the RGP is one of the oldest permanent GNSS network in the world (late 90s) and it has a dense coverage of the French metropolitan territory, including the coastal areas. These long-term records could be used for long-term monitoring of the tides, but also of the wave setup and surges, and perhaps for studying sea level rise due to global warming. These “accidental” records could also present a growing interest in the context of global warming with the possible intensification of storm frequencies. Besides, as the newly installed geodetic stations now includes the acquisitions from the new constellations, like GALILEO and BeiDou, but also new signals as L5, E5, and E6 will increase the accuracy in the GNSS-R SSH retrievals.

If we only consider the very low incidence angles (0°–5°), the distance between the GNSS geodetic station and the farthest reflection points can reach several kilometers: 8 km for an antenna height of 60 m above the sea level and to ~28 km for an antenna height of 250 m [48]. GNSS-R SSH could also be used for early warning in the case of a huge storm and tsunami.