GNSS-R monitoring of soil moisture dynamics in areas of severe drought: example of Dahra in the Sahelian climatic zone (Senegal)

ABSTRACT

With population growth, water will increase in the following decades tremendously. The optimization of water allocation for agriculture requires accurate soil moisture (SM) monitoring. Recent Global Navigation Satellite System Reflectometry (GNSS-R) studies take advantage of continuously emitted navigation signals by the Global Navigation Satellite System (GNSS) constellations to retrieve spatiotemporal soil moisture changes for soil with high clay content. It presents the advantage of sensing a whole surface around a reference GNSS antenna. This article focuses on sandy SM monitoring in the driest condition observed in the study field of Dahra, (Senegal). The area consists of 95% sand and in situ volumetric soil moisture (VSM) range from ~3% to ~5% durinf the dry to the rainy season. Unfortunately, the GNSS signals’ waves penetrated deep into the soil during the dry period and strongly reduced the accuracy of GNSS reflectometry (GNSS-R) surface moisture measurements. However, we obtain VSM estimate at low/medium penetration depth. The correlation reaches 0.9 with VSM error lower than 0.16% for the 5–10-cm-depth probes and achieves excellent temporal monitoring to benefit from the antenna heights directly correlated to spatial resolution. The SM measurement models in our research are potentially valuable tools that contribute to the planning of sustainable agriculture, especially in countries often affected by drought.

KEYWORDS:

• GNSS-R
• signal-to-noise ratio (SNR)
• Interference Pattern Technique (IPT)
• phase unwrapping
• soil moisture

1. Introduction

Water resource management is considered one of the critical problems in the world, especially in the context of the growth of the population, the development of the industry, the need for increased food production and the pollution of the water resource due to various human activities. In addition, it is also necessary to ensure sustainable agricultural development (Zhang et al., 2017) and flood control (Rodriguez-Iturbe & Porporato, 2004).

A vital parameter of the water resource is soil moisture (SM). Indeed, its impact on the soil-vegetation-atmospheric exchange is exceptionally high in three regions: the Great Plains of North America (Xu et al., 2018), India and the Sahel, where rainfall is strongly influenced by SM (Koster et al., 2004). In many countries like Mali, or Senegal, water resources basically come from groundwater reservoirs that are declining and not renewable. There are no clear conclusions in previous studies whether or not this groundwater depends on SM variations due to the limitation of SM data and the complexity of hydrological phenomena (Fatras et al., 2012; Lebel et al., 2009; Lopez et al., 2016; Mougin et al., 2009). Knowledge of these phenomena is fundamental for a proper understanding (Rascon et al., 2021; Steelman et al., 2012) of water resources management, in ecology, in the development of agricultural strategies, as well as in case of extreme events (He et al., 2012; Martínez-Fernández et al., 2021; Sun et al., 2015).

As a result of its importance, various methods have been developed to measure SM indices, with three main approaches. The most direct one consists of using gravity surveys to carry out measurements in the field (Gehman et al., 2009), and it enables calibrating the measurement probes. These probes are today the most widely used field measurement techniques, even with the information being scarce and very local. Modeling enables the provision of large-scale SM values but remains dependent on the availability of other variables such as meteorological data (Kumar et al., 2008; Mitchell et al., 2004), the selection of model parameters and the model itself (e.g. WGHM (Döll et al., 2003), Global Land Data Assimilation System models (Rodell et al., 2004)).

Remote sensing in the microwave domain is based on the gap presets and is advantageous in that it can work any time of day or night, regardless of the weather conditions at the time of acquisition. It records information directly related to SM using the soil dielectric properties and electromagnetic (EM) radiation in the microwave region 0.5–100 cm (Karthikeyan et al., 2017; Sun & Cui, 2021; Wang & Qu, 2009). Many research and experiments have used the SM inversion method from an active microwave, which assumes that heat exchange in the ground happens through thermal conduction (Han et al., 2018; Kerr, 2007; Xu et al., 2020). However, these spectral windows are also susceptible to cloud cover, limiting their use due especially to their low spatial and temporal coverage.

Meanwhile, passive microwaves use a sensor which emits an EM wave in a known frequency and polarization in the direction of a target, such as its surface (Xu et al., 2020). Passive microwaves are preferred for estimating SM (Chew & Small, 2018; Entekhabi et al., 2010; Peng & Loew, 2017). Passive microwave sensors including Scanning Multichannel Microwave Radiometer, Special Sensor Microwave/Imager sensors, Advanced Microwave Scanning Radiometer Earth Observing System and Soil Moisture Ocean Salinity Satellite/Soil Moisture Active Passive mission have been used successfully to measure SM (Ahmad et al., 2011; Jones et al., 2010; Kerr et al., 2001; Njoku & Chan, 2006; Paloscia et al., 2001; Wen et al., 2005; Wigneron et al., 2021). However, spatial and temporal resolution is significantly lower than the infrared range and is thus considered a significant limiting factor for SM measurements.

In addition, recent used EM waves come from Global Navigation Satellite System(GNSS) (GPS, GLONASS, Galileo, etc.), which extract the various geographic parameters of the earth surface with high spatiotemporal resolution (Darrozes et al., 2016). This opportunistic remote sensing technique is called GNSS reflectometry (GNSS-R), based on the evaluation of GNSS waves obtained by an antenna after reflection on the earth surface (Aranyossy & Ndiaye, 1993; Alonso-Arroyo et al., 2014; Darrozes et al., 2016; Roussel et al., 2014).

Most GNSS signals are received directly by a geodetic antenna’s upper hemisphere, as Right-Hand Circular Polarization (RHCP) waves. Meanwhile, a portion comes from below the horizon of the antenna after reflections from surrounding surfaces. Fortunately, the lower antenna gain pattern is mainly Left-Hand Circular Polarization (LHCP) and records LHCP signal corresponding to the reflected signal. EM interferences between direct (RHCP from upper part of the antenna) and reflected waves (LHCP from lower part) are especially obvious in the total signal-to-noise ratio (SNR) recorded by the receiver.

This SNR technique is applied in various studies of Earth Sciences, particularly SM analysis. The temporal resolution can be several minutes, hours or days, depending on the calculating algorithm applied to assess the SM variation (Chew et al., 2014; Larson et al., 2008).

The reflected GNSS signal changes are due to the surface variations (soil composition, moisture, roughness) and their characteristics. Therefore, the SNR is then evaluated to predict SM content. This SNR method is applied using any conventional GNSS receiver/antenna without any hardware modifications. The GNSS R technique has SM retrieval accurately at 2–5 cm depth over clayed soils (Rodriguez-Alvarez et al., 2011; Roussel, 2015). Nevertheless, if the SM content is less than 10% volumetric soil moisture (VSM), it causes cyclic wrapping of the phase.

The present study aims to develop a comprehensive approach for desert or semiarid climates where the SM levels are weak. According to the theory of radar waves, these weak SM levels induce a significant wave penetration (Ulaby & Moore, 1986). This penetration depends on the soil composition. For clay soils or rich clay soil content, even in case of extreme dryness, the waves are reflected in the first centimeters, and the SM retrieval on the surface is well correlated to the in situ SM probe at 5 cm depth (Vey et al., 2016; Zhang et al., 2018). For sandy soils, the penetration is more critical and exceeds the meter in conditions of intense drought. Our approach estimating SM is based on state-of-the-art GNSS-R for SM retrieval. It is combined with the phase unwrapping, generally used in SAR imagery (Benoit et al., 2020; Lu & Zebker, 1998) for very weak values of SM, which induces the intense penetration for the sandy soil of the study area. The method using SNR data incorporates both direct and reflected signals, continuously collected by the geodetic antenna. The variation in soil moisture is not the only parameter that will affect the reflected signal of radar waves; two other aspects will interfere with bare soil: (i) the surface roughness of the soil and (ii) variations in soil composition. However, of these three parameters, only variations of the surface soil moisture (SSM) has a cyclical response related to seasonal changes of this Sahelian area. SSM changes the EM properties of the reflected waves, i.e. the complex dielectric constant that depends on electrical conductivity and relative permittivity of the media. Two parameters, namely the phase and the frequency of the multipath contribution to SNR, are analyzed to connected with the antenna height above the reflecting surface (Lopez et al., 2016). Results obtained from the total of the different GNSS satellite constellations will be evaluated applying independent SM records measured in the field with Theta probe sensors (Roussel, 2015).

2. Methodology

2.1.1 Interference Pattern Technique

Interference Pattern Technique (IPT) is a method based on analyzing the multipath signature present in the GNSS SNR data recorded continuously by a GNSS receiver and its associated single geodetic antenna (Figure 1a). Multipath propagation affects the SNR recorded by a GNSS receiver by generating frequency oscillations with lower amplitudes than the direct signals. These frequency oscillations depend on the satellite elevation (SE) angle and vary from the intermediate frequency at Low Satellite Elevation (LSE) to the high frequency at High Satellite Elevation (HSE).

Figure 1. (a) Transmission of EM waves to the single geodetic antenna: direct GNSS signal path with RHCP and the reflected path with LHCP; (b) SNR in dB-Hz of GNSS signals and SNR multipath (SNRm) signals removed the direct contribution of signal (parabolic red line) from the initial SNR profile (Modified from Roussel et al., 2016 and Darrozes et al., 2016).

Because the antenna gain pattern prefers the direct signal over the reflected signal by construction, and during reflection will also reduce the signal energy; the primary influence of the structure of the SNR time series is the direct one (Figure 1b). The SNR total show a low frequency which corresponds to the main variation of the energy of the direct signal as a function of the SE angle. The multipath (Figure 1b, lower graph) under medium frequency oscillation is visible for low elevations while for higher elevations of the satellite these multipaths do not emerge or poorly emerge of noise level.

For multipath component analysis, the direct signal contribution must be eliminated by fitting and subtracting a simple parabolic polynomial from the raw time series of SNR to obtain the residual of SNR multipath (SNR_m) component (Larson et al., 2008; Figure 1b, red curve).

The influence of multipath interference appears in SNR measurements at different frequencies of the L-bands, such as L1, L2 and L5. According to Larson et al. (2010), for a static or quasi-static surface (vertical velocity of the surface is lower than 10^–6 m.s⁻¹). the oscillations frequency is approximated as a sinusoidal function of the SE:

(1) $\tilde{f}\approx \frac{4\pi }{\lambda }{h}_{eff}$(1)

where $h_{eff}$ is the effective antenna height measured by the antenna height H_o (Figure 1a) and λ is the signal wavelength. The multipath oscillation frequency $\tilde{f}$ and the effective antenna height $h_{eff}$ are directly proportional, as shown in equation above, and can be calculated for each satellite at each epoch and for each specular point of the time series where the moisture measurement occurs.

When the GNSS wave penetrates the reflecting surface, $h_{eff}$ varies depending on the characteristics of the non-moving or weakly moving surface (e.g. vegetation density, SM, soil type, Figure 2a).

Figure 2. (a) Theoretical penetration depth of the GNSS waves when the ground consists of 100% clay (red), 95% sand (blue, our study) and 100% sand (green); (b) theoretical penetration depth of the various GNSS wavelengths and for various SE for a sandy soil; (c) representation of the unwrapping phase process (Zebker & Yanping, 1998). Red arrow shows a positive phase slip while green arrow shows a negative one.

Frequency $\tilde{f}$ is estimated, based on harmonic Fourier decomposition of the Lomb–Scargle periodogram (Lomb, 1976; Scargle, 1982), using a moving window from the SNR_m time series. We also add, before Lomb–Scargle analysis, a step of phase unwrapping (Figure 2, see Section 2.2). The retained periodogram peaks must reach a significance level equal to 1% (Li et al., 2018) and defining a local maximum between ${\tilde{f}}_{max}$ and ${\tilde{f}}_{min}$ using user-defined h_min,h_max,; in our case, the penetration of the wave leads to h_max = H_o + P_d, where P_d is the penetration depth of the wave (see next section).

According to the results of Larson et al. (2010), and a given quite constant $h_{eff}$, SNR_m can be expressed as follows:

(2) $SN{R}_m=A_{m}cos\left(\frac{4\pi h_{eff}}{\lambda }sin\left(\theta \right)+{\mathrm{\Phi }}_m\right)$(2)

where $A_m$ is the amplitude of the multipath oscillations and ${ }_m$ is the phase of the multipath oscillations. $A_m$ denotes the backscatter intensity of ground reflections, including both the gain pattern and multipath power. According to field studies, $A_m$ and ${\mathrm{\Phi }}_m$ are strongly correlated with the level of SM (Larson et al., 2010), and SM is more sensitive to ${\mathrm{\Phi }}_m$ than $A_m$ (Chew et al., 2015; Roussel et al., 2016). Consequently, the three metrics of $A_m$, ${\mathrm{\Phi }}_m$ (above) and $h_{eff}$ (derived from $\tilde{f}$, see Equation (1)(1) $\tilde{f}\approx \frac{4\pi }{\lambda }{h}_{eff}$(1) ) can be used to estimate SM variations.

2.2 Phase unwrapping method

2.2.1 Unwrapping SNR phase method

The GNSS waves penetration depth depends on the soil dielectric properties (Njoku & Entekhabi, 1996), the SM content and the SE (Figure 2a, 2b). The zenith signal’s penetration depth $P_d$, according to Njoku and Entekhabi (1996) and Behari (2006), can be expressed as follows:

(3) $P_d=\frac{\lambda \sqrt{R_e\left({\epsilon }_r\right)}}{2\pi I_m\left({\epsilon }_r\right)}$(3)

where $R_e\left({\epsilon }_r\right)$ and $I_m\left({\epsilon }_r\right)$ are the real and imaginary parts of the reflecting surface relative permittivity, respectively.

Recent studies applying IPT on clay soils give excellent results with high accuracy SM estimation at 2–7 cm depth because the penetration depth of the wave is less than the wavelengths of the GNSS signals (Figure 2a, b). However, if we are in the conditions of sandy soils and VSM has low values, i.e. in the driest season, the wave penetration can exceed the wavelength $\lambda$ and causes the cyclic phase wrapping (Figure 2a). Theoretical, the penetration depth of GNSS signals shows a significant increase when the soil is 100% clay (red) or 100% sand (green), especially for sandy soil when VSM is very low (<10%) (Figure 2a). This penetration also depends on the SE and the wavelength of the GNSS signal (Figure 2b). When SE is low, the penetration of the wave will be lower, while for high SE, the energy of the EM signal being more important, the depth of penetration will be greater. Still, from a theoretical point of view, we see that when the wavelength of the signal increases, its penetration depth also increases (Figure 2b). Thus, when the penetration depth is greater than the wavelength $\lambda$ used, the phase undergoes one or more cycle slips. As a result, the three metrics estimated from SNR_m does not correlate with the SM. Therefore, in the sandy field, we propose using an additional phase unwrapping algorithm to solve this problem by reconstituting the phase cycle slip that appears in the estimated ${\mathrm{\Phi }}_m$, thus obtaining continuous phase variations (Figure 2c).

Figure 2c shows the example of a discrete noisy phase signal $\phi \left(x\right)$ whose amplitude exceeds the range [-π, π] by processing phase wrapping operation $W\left(\phi \left(x\right)\right)$ with the original continuous noisy phase signal. The wrapped phase $\psi =W\left(\phi \left(x\right)\right)$ is calculated using the four-quadrant arctangent (atan2) function of $sin\left(\phi \left(x\right)\right)$ and $cos\left(\phi \left(x\right)\right)$. We must remove the 2π slips that are present in the wrapped phase signal (Figure 2c) to eliminate the phase cycle slip and restore the continuous form, allowing the phase variation $\phi \left(x\right)$ usable in analysis or further processing (Figure 2d). This process is called phase unwrapping and can be represented as follows:

(4) ${\phi }_{unw}\left(x\right)=U\left[\psi \right]=\psi +2k\pi$(4)

where ${\mathrm{\Phi }}_{unw}\left(x\right)$ is the unwrapped phase signal, $U\left[\psi \right]$ is the operation of phase unwrapping and $k$ is an integer relating to the cycle number of phase slips. To apply the phase unwrapping procedure for the estimated phase ${\mathrm{\Phi }}_m$ from the SNRm time series, we must calculate the difference $\mathrm{\Delta }\psi$ between two adjacent successive samples on its left. Once a phase wrap is detected, $\mathrm{\Delta }\psi$ is less than –π or greater than +π, 2π is either subtracted (negative phase slip, k < 0) or added (positive phase slip, k > 0), from this sample to the right-hand side. In addition, this phase unwrapping method is ideally suited to noisy data, such as GNSS signals (Gdeisat, 2018) or InSAR images (Hooper & Zebker, 2007; Smith, 2002; Zebker & Yanping, 1998).

2.2.2 New metrics

Following Equation (2)(2) $SN{R}_m=A_{m}cos\left(\frac{4\pi h_{eff}}{\lambda }sin\left(\theta \right)+{\mathrm{\Phi }}_m\right)$(2) , we can obtain the three metrics $h_{eff}$, A_m and Φ_m from SNR_m for each GNSS constellation (see IPT). Considering the A_m depending on the antenna gain pattern and multipath intensity, it is constant. However, the phase unwrapping algorithm processing changes the phase ${\mathrm{\Phi }}_m$ to obtain a variation of continuous phase ${\mathrm{\Phi }}_{unw}$, the effective antenna height $h_{eff}$ (proportional to the frequency $\tilde{f}$, Equation (1)(1) $\tilde{f}\approx \frac{4\pi }{\lambda }{h}_{eff}$(1) ) must be re-calculated. The new height effective $h_{eff}^{unw}$ after using the unwrapping phase is determined as follows:

(5) $h_{eff}^{unw}=\left({cos}^{-1}\left(\frac{SN{R}_m}{A_m}\right)-{\mathrm{\Phi }}_{unw}\right)\times \frac{\lambda }{\left(4\pi \times sin\left(\theta \right)\right)}$(5)

2.3 Estimating SM from SNR signals using IPT and unwrapping phase method

Following the GNSS-R theory presented in the previous parts, a complex processing chain integrating IPT and phase unwrapping methods has been developed to determine SSM in sandy environments (Figure 3). In this processing chain, we employ a new time series of $h_{eff}^{unw-cal}$, calibrated with in situ VSM to remove negative correlation of $h_{eff}^{unw}$, to better retrieve variations of SM and to establish the SM maps.

Figure 3. The flow chart represents the complex processing chain for SSM detection/observation in sandy environments using an integration of IPT and phase unwrapping methods: the input, the various processing steps and the output.

As defined in Equation (3)(3) $P_d=\frac{\lambda \sqrt{R_e\left({\epsilon }_r\right)}}{2\pi I_m\left({\epsilon }_r\right)}$(3) , the new variable $P_{d\mathrm{\_}eff}$ is computed using the effective height $h_{eff}^{unw}$ derived from the unwrapping phase process:

(6) $P_{d\mathrm{\_}eff}=\left(h_{eff}^{unw}-H_o\right)\pm E_r$(6)

where $H_o$ is the antenna height from the phase center to the ground without any penetration and $E_r$ corresponds to the total error term due to the measurement error for each GNSS wavelength and errors associated with the phase unwrapping and the surface roughness, which will further increase the measurement noise.

A Savitzky–Golay (Sgoley) filtering method (Baba et al., 2014; Schafer, 2011) was applied for $h_{eff}^{unw}$ to minimize errors caused by interference values (noisy signal whose frequency span is large) to obtain the filtered effective height $h_{eff}^{unw-cal}$ (Equation (7))(7) $H_{Sgo}\to \widehat{H_{Sgo}}=\frac{1}{L}{\sum }_{n=-k}^k{H}_{Sgo+n}$(7) . This improvement in the height reduces the noise of the time series and increases the accuracy of SM estimations. The principle of this filter is to replace each $h_{eff}^{unw}$ data point with some combination of the values contained in the moving window centred at that point. Sgoley filtering process can be described as follows:

(7) $H_{Sgo}\to \widehat{H_{Sgo}}=\frac{1}{L}{\sum }_{n=-k}^k{H}_{Sgo+n}$(7)

where $\widehat{H_{Sgo}}$ is the result from the Sgoley filter, $H_{Sgo}$ is the input data and $L=2k+1$ is the moving window length.

New time series ${\mathrm{\Phi }}_{unw}$, $h_{eff}^{unw}$ and $h_{eff}^{unw-cal}$ were normalized to minimum and maximum of in situ SM variables recorded by a classical SM sensor (Theta probe). Pearson correlation coefficient R is quantified to determine the correlation between these normalized time series with VSM measured by the Theta probe. We can also describe the differences between them in terms of the totalized root mean squared error ($RMS{E}_{total}$) as follows:

(8) $RMS{E}_{total}=RMS{E}_{rp}+RMS{E}_{sp}+RMS{E}_{SE}$(8)

where $RMS{E}_{rp}$ is an error caused by the difference in position between the reflected point used to retrieve the SM and the Theta probe used to measure the reference moisture; $RMS{E}_{sp}$ is the error caused by the difference in spatial observation of the two moisture measurement methods (ML3 Theta probe is $\text{~}$0.03 m², GNSS-R IPT is from 0.8 to 137 m² for an antenna height H_o = 1.61 m with SE of 2–70°); and $RMS{E}_{SE}$ is the error caused by the GNSS instrument affected by SE, the roughness of the ground and the GNSS wave penetration depth.

Although the statistics parameters described above are helpful in research for linear problems, they cannot reflect all of the correlations seen in wave field studies. As a result, we propose an additional method for analyzing waveform data to consider the correlations known as Wavelet Transform Coherence (WTC) (Grinsted et al., 2004). The WTC is a method to analyze the phase lag and coherence between two-time series S (time series of ${\mathrm{\Phi }}_{unw}$/$h_{eff}^{unw}$/$h_{eff}^{unw-cal}$ and in situ measurements) as a function of time and frequency.

To compute the WTC analysis, we first have to transform our time series into a wavelet coefficient map. The wavelet coefficient power map $V_{\left(a,n\right)}^{Si}$ corresponds to the convolution in the frequency domain of the times series $S_{\mathrm{i}}$ and a set of wavelet functions ${ }_{\left(\mathrm{a},\mathrm{n}\right)}$ (Equation (9)(9) $V_{\left(a,n\right)}^{Si}={}_{\left(\mathrm{a},\mathrm{n}\right)}\ast S_{\mathrm{i}\left(\mathrm{n}\right)}$(9) ). This set is derived from a mother function well known as Morlet wavelet function, a Gaussian-windowed complex sinusoid. For each wavelet function, a coefficient is computed for the time series’ time index n, and a defined period a. This period a varying from the minimal period, i.e. 1/5 of day, to the maximal period ~1 month.

(9) $V_{\left(a,n\right)}^{Si}={}_{\left(\mathrm{a},\mathrm{n}\right)}\ast S_{\mathrm{i}\left(\mathrm{n}\right)}$(9)

WTC analysis performs detection of phase difference (= time delay), phase angle (cross magnitude), non-stationary and also wavelet coherence $M^2$ between the two-time series (Equation (10)(10) $M_n^2\left(s\right)=\frac{{\left|S\left(s^{-1}{V}_n^{XY}\left(s\right)\right)\right|}^2}{S\left(s^{-1}{\left|V_n^X\left(s\right)\right|}^2\right).S\left(s^{-1}{\left|V_n^Y\left(s\right)\right|}^2\right)}$(10) ). The phase difference between the two-time series is indicated by the arrows’ directions, and we can gain confidence in causal and geophysical test relationships using phase angle statistics. The two-time series are in-phase if the arrow is horizontal and pointing to the right, in-phase opposition if the arrow is horizontal but pointing to the left, and phase quadrature if the arrow is pointing up or down. According to Torrence and Webster (1999) and Grinsted et al. (2004), wavelet coherence $M^2$ can be calculated using the following formula:

(10) $M_n^2\left(s\right)=\frac{{\left|S\left(s^{-1}{V}_n^{XY}\left(s\right)\right)\right|}^2}{S\left(s^{-1}{\left|V_n^X\left(s\right)\right|}^2\right).S\left(s^{-1}{\left|V_n^Y\left(s\right)\right|}^2\right)}$(10)

3 Experimental setup

3.1 Measurement site

Senegal has a strong diversity in climate, and it crosses three different climate zones: in the north the desert, and then crosses the semiarid zone, which leaves the place, in its southward part, to tropical climate (Figure 4). Senegal is crossed by the West African river, which gave name to the country. Senegal River is 1800 km long. However, the catchment area is six times smaller than the two principal African basins, i.e. Congo and Nile; its area covers ~0.5 million km², with an average discharge of ~630 m³. s⁻¹.

Figure 4. Map of African Climate using Köppen climate classification, the studied watershed is delimited by a dark red line and associated hydrographic network corresponds to dark blue line for main rivers until white line for small tributaries or gully (Modified from Peel et al., 2007; Ramillien et al., 2021). The white rectangle corresponds to the study area presented in Figure 6.

The study area is located near Dahra, Senegal, where the climate is semiarid. The main rainfall occurs during the rainy season, extending from mid-June to mid-October (Figure 5) with relative humidity ranging from 60% to 75%. It is followed by a dry season that extends from late October to late May, with relative humidity ranging from 20% to ~0%. This season is marked by a rise in temperatures, reaching 33°C, and a sharp drop in rainfall. It is also a slightly windier period with winds up to 3.5 m. s⁻¹. A detailed description of the field site can be found in Tagesson et al. (2015). Agriculture is one of Senegal’s strengths, but while the south (tropical climate) has very favorable conditions, this is less true in the semiarid Sahelian zone. We then understand the importance of mapping SM to water the fields while saving this rare water commodity.

Figure 5. Dahra region charts of (a) monthly precipitation (mm), (b) relative humidity (%), (c) air temperature at 2 m (°C) and (d) wind speed at 2 m (m . s⁻¹) from 2015 to 2017. Black vertical dashes correspond to January for each year (Modified from Laouali et al., 2021).

To make this mapping, we have installed a Leica GR25 receiver with an AR10 antenna in a sparsely wooded grassland in Dahra experiment (Figure 6a, b). In terms of geology and soil, we have a homogeneous area of rubified sand from a continental dune (Faye et al., 2017) of further square kilometers. Therefore, this field is considered to have no significant topography changes (Figure 6c), and in the modeling of the specular point, we assume that it is a flat area, and we do not use digital elevation model corrections. GPS and GLONASS signals were continuously acquired at a 1 Hz sampling frequency from October to December 2016.

Figure 6. (a) The experiment site of Dahra (Senegal), (b) location of the AR10 antenna (15°24ʹ5.36” N; 15°26ʹ0.13” W); (c) photography of the AR10 antenna with a Leica GR25 receiver installed on sandy soil, this flat area is covered by grassland and sparce trees. The ring of thorny trees surrounding the antenna is clearly visible and serves to protect the equipment from animals.

3.2 Data used for validation

The SSM is monitored by an SM station located about 500 m away from the GNSS receiver (Figure 7a). Some authors (e.g. Baup et al., 2007; Rosnay et al., 2009) have shown that in semiarid Sahelian areas, with equivalent soil and vegetation, in situ SM measurements are representative of much larger area. This is why, we consider in this case of very homogeneous geological/soil (Faye et al., 2017), the moisture station represents a large area, including the GNSS-R station and can be used for calibration. Soil content of the Dahra experiment is ~95% of sand associated with few percent of organic matter in surface and few percent of clay. In order to avoid surface disturbances (presence of this small concentration of organic matter), we took soil samples at 10 cm depth. These samples, as shown in Figure 7a, are representative of the soil column down to at least 1 m depth (no visible soil change).

Figure 7. (a) Configurations of the SM station (15°43ʹ42.89” N, 15°00ʹ55.69” W). The station equipped with six Theta probes (green): (number 1, 3): 5 cm, (2, 4): 10 cm, (6): 100 cm; (b) VSM recorded by Theta probes during the study period in 2016 which corresponds to the transition between rain season and dry season. The gray rectangle corresponds to a non-acquisition period, light blue corresponds to the precipitation (mm of water) during the study period.

For the GNSS frequencies which are in L-band, we calculated the theoretical curves for a soil consisting of 100% sand and another consisting of 100% clay (Figure 2a). This figure shows that the depth of penetration becomes important as soon as the SM is weak and below 10% the penetration exceeds 7 cm to reach 40 cm for SM~1%. The response of our soil (blue curve) is close to the theoretical green curve (100% sand).

The SM station is equipped with two soil temperature sensors and 5 ML3 Theta probe sensors with a sampling frequency of 15 min (accuracy of ± 0.1%), installed at 5 cm (1, 3), 10 cm (2, 4) and 1 m depth (5 in Figure 7a). Figure 7b shows the variation of SM measured by the Theta probes. During the transition phase from the rainy season to the dry period, the VSM at 5 cm depth (P_5cm) varied from 2.23% to 4.39%, the VSM at 10 cm depth (P_10cm) varied from 2.87% to 4.76% and the VSM at 1 m depth (P_1m) varied from 5.03% to 5.63%. At P_5cm and P_10cm depth, VSM shows a daily sinusoidal pattern.

4 Results

4.1 Importance of the unwrapping in sandy soil

For sandy soils of Dahra, we have a very noisy phase dataset (Figure 8a) without any correlation with the in situ probes. So we applied the unwrapping technique to this phase dataset and we make an inspection of the unwrapped/wrapped signals to identify possible fake slips due to the noise. To calculate the number of fake slips, we compute the difference between two adjacent samples of the wrapped signal. If the computed difference is upper than twice the noise variance threshold criteria (TC), we assume we have a fake slip. In that case, the signal after the fake slip is affected by “miss wrap” due to the error propagation.

Figure 8. Impact of the phase unwrapping for sandy soils (a) before the unwrapping and (b) after the unwrapping. The correlation factor reaches 0.9.

This part of the time series is affected by a phase offset that can be positive or negative depending on the noise difference sign. We need to apply the unwrapping process to correct these harmful effects until all the computed differences are lower than TC. The resulting phase correlates with the in situ VSM measurements (P_5cm, (1) Figure 7a) with a calibration function in this type of soil of VSM~1 x unwrapped phase where the bias is close to 0 (Figure 8b) which shows that we do not need to calibrate the phase to obtain VSM because of VSM ~${\mathrm{\Phi }}_{unw}$. The rest of the study will use the ${\mathrm{\Phi }}_{unw}$ dataset to compute the antenna height parameters (see 4.2).

4.2 Comparisons between $\hspace{0.17em}{\mathrm{\Phi }}_{unw}$, $h_{eff}^{unw}$ and $h_{eff}^{unw-cal}\hspace{0.17em}$ parameters for estimating SM in Dahra experiment

Table 1 lists the results of root mean square error (RMSE) and correlation coefficient of Pearson (R) between ${\mathrm{\Phi }}_{unw}$,$h_{eff}^{unw}$, and $h_{eff}^{unw-cal}$ ($h_{eff}^{unw}$ after Sgoley filtering) with in situ VSM measured by the Theta probes at P_5cm and P_10cm for different constellations from October to December 2016 in the study area. We do not use P_1m due to the weak correlation with subsurface SNR_m measurement. This observation shows that the wave penetration depth is less than 1 m, which agrees with the theoretical model for VSM upper than 2.5% (Figure 2a). In the Darha experiment, where the soil consists of over 95% sand, retrieving SM from GNSS-R data became more complicated. The A_m values can be disturbed by vegetation on the surface around the antenna. Thus, we used $h_{eff}^{unw}$ $,h_{eff}^{unw-cal}$ and ${\mathrm{\Phi }}_{unw}$ to retrieve variations of the VSM. The three parameters show a good correlation with the in situ measurements of SM at 5 and 10 cm. Indeed, R is always higher than 0.73 in the rainy season and higher than 0.51 during the dry season. However, during the rainy season, ${\mathrm{\Phi }}_{unw}$ has a lower correlation than $h_{eff}^{unw-cal}$ and results quite similar to those of $h_{eff}^{unw}$ regardless of the probe depth, the SE or even the constellation used. Conversely, ${\mathrm{\Phi }}_{unw}$ has the slightest deviation from the in situ measurements during the dry period, with an RMSE varying between 0.18% and 0.08%. It also has the best correlation, with R ranging between 0.92 and 0.95 (Table 1).

Table 1. Comparison between (${\mathrm{\Phi }}_{unw}$, ${\mathrm{h}}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{u}\mathrm{n}\mathrm{w}}$, and ${\mathrm{h}}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{u}\mathrm{n}\mathrm{w}-\mathrm{c}\mathrm{a}\mathrm{l}}$) and VSM measured by P_5cm/P_10cm (RMSE_Total/R) for difference constellations, and for SNR from LSE (2–30°) and HSE (30–70°) from October to December 2016 in Dahra site.

GNSS	SNR	RMSETotal(%)/R
		Rainy period (01/10/2016–15/10/2016)				Dry period (16/10/2016–06/12/2016)
		SE2-30		SE30-70		SE2-30		SE30-70
		P5cm	P10cm	P5cm	P10cm	P5cm	P10cm	P5cm	P10cm
GPS	Φunw	0.20/0.77	0.12/ 0.80	0.21/0.78	0.10/0.85	0.14/0.86	0.08/0.96	0.18/0.76	0.09/ 0.95
	$h_{eff}^{unw}$	0.21/0.74	0.12/0.89	0.22/0.74	0.14/0.87	0.19/0.51	0.11/0.84	0.19/0.56	0.11/0.81
	$h_{eff}^{unw-cal}$	0.15/0.89	0.08/0.96	0.17/0.87	0.08/0.97	0.11/0.85	0.07/0.93	0.15/0.77	0.08/0.92
GLO	Φunw	0.19/0.77	0.12/ 0.77	0.19/0.78	0.10/0.83	0.16/0.81	0.09/0.95	0.16/ 0.80	0.09/ 0.94
	$h_{eff}^{unw}$	0.22/0.74	0.12/0.89	0.22/0.7	0.13/0.89	0.19/0.56	0.11/0.82	0.19/0.51	0.11/0.82
	$h_{eff}^{unw-cal}$	0.16/0.87	0.09/0.95	0.15/0.88	0.09/0.95	0.14/0.78	0.08/0.92	0.14/0.78	0.07/0.92
GPS+GLO	Φunw	0.21/0.73	0.11/ 0.82	0.19/0.79	0.10/0.83	0.17/0.79	0.09/0.95	0.16/ 0.80	0.09/ 0.95
	$h_{eff}^{unw}$	0.20/0.79	0.12/0.90	0.23/0.74	0.14/0.87	0.18/0.59	0.10/0.84	0.19/0.53	0.11/0.82
	$h_{eff}^{unw-cal}$	0.16/0.87	0.09/0.95	0.16/0.87	0.09/0.95	0.14/0.78	0.08/0.91	0.14/0.77	0.08/0.92

Finally, the parameter that gives the best results over the two seasons is clearly $h_{eff}^{unw-cal}$ with the best correlation, with the in situ measurements at 10 cm depth, which is always higher than 0.91 and an RMSE lower than 0.16%. This parameter is less influenced by the vegetation density and soil roughness. It is also less sensitive to the constellation used: GPS and GLONASS alone and GPS+GLONASS have very closed results. Nevertheless, the number of specular points increases when both constellations are used simultaneously, resulting in more accurate moisture maps (more points with a better spatial distribution). Furthermore, the elevation of GNSS satellites has little influence on this parameter, and the results are very favorable for low and high angles of incidence (Table 1). On the other hand, the parameter $h_{eff}^{unw}$, even if it presents globally good results, generally shows the weakest correlation of the three parameters studied here and the most significant differences with the in situ measurements of SM. In bold: results with a correlation greater than or equal to 0.9.

4.3 Wavelet analysis of the SNR time series

For WTC analysis, it is necessary to construct wavelet coefficient maps of two-time series, i.e. $V_{\left(a,n\right)}^{VSMin-situ}$ and $V_{\left(a,n\right)}^{VSM{ }_{unw}}$ derived from the in situ time series of the SM probe and the SM retrieved from ${\mathrm{\Phi }}_{unw}$. We used only the probe at 5 and 10 cm depth. The sensor probe at 1 m is too deep, and GNSS waves cannot reach this depth due to the moisture conditions (>2.4%) (Figure 2a).

Figure 9 shows the time series of VSM P_10cm and VSM derived from the $\mathrm{\Phi }{\$}_{unw}$ and $h_{eff}^{unw\mathrm{\_}cal}$ (Figure 9a) and the wavelet coefficient map of the three time series of VSM P_10cm (Figure 9b), $VSM{ }_{unw}$ (Figure 9c) and $VSM{h}_{eff}^{unw\mathrm{\_}cal}$ (Figure 9d), respectively. These time series have quite similar coefficient maps with the main period of more than 8 days. The results show a normalized coefficient ranging from 0.6 to 0.8, with a period of 1 day for the in situ time series and a period of 1–2 days for the GNSS-R time series derived from the ${\mathrm{\Phi }}_{unw}$ and $h_{eff}^{unw\mathrm{\_}cal}$. Two GNSS-R time series have a power noise close to 0.2/0.3 at the period varying from half a day to a quarter of the day. Figure 9c shows a good correlation between VSM P_10cm with $VSM{h}_{eff}^{unw\mathrm{\_}cal}$ in the rainy season.

Figure 9. Time series and magnitude maps calculated from Torrence and Compo (1998). (a) In situ P_10cm and ${ }_{unw}$, $h_{eff}^{unw\mathrm{\_}cal}$ time series and associated wavelet coefficient maps; (b) $V_{\left(a,n\right)}^{VSMin-situ}$; (c) $V_{\left(a,n\right)}^{VSM{ }_{unw}}$; (d) $V_{\left(a,n\right)}^{VSM{h}_{eff}^{unw\mathrm{\_}cal}}$.

4.4 Wavelet Coherency Analysis between VSM estimated from $\hspace{0.17em}{\mathrm{\Phi }}_{unw}$ and $h_{eff}^{unw-cal}$ with VSM measured by in situ probe

The correlation between VSM estimated from ${\mathrm{\Phi }}_{unw}$ and $h_{eff}^{unw-cal}$ with VSM measured by probes at P_5cm (Figure 10a, b) and P_10cm (Figure 11a, b) exhibited significant changes by period. Figures show that the VSM derived from ${\mathrm{\Phi }}_{unw}$ _/ $h_{eff}^{unw-cal}$ and in situ VSM measured has no significant correlation at the period of less than a half-day. During the dry season, a negligible correlation was found in the period of fewer than 4 days. This decorrelation is mainly due to the wave penetration. In this period GNSS wave reaches more than 10 cm and does not correspond to the surface condition.

Figure 10. The wavelet coherence (M²) of WTC analysis between VSM derived from GNSS-R parameters and in situ VSM at P_5cm (a) ${ }_{unw}parameter$, (b) $h_{eff}^{unw-cal}parameter$.

Figure 11. The wavelet coherence (M²) of WTC analysis between VSM derived from GNSS-R parameters and in situ VSM at P_10cm. (a) ${ }_{unw}parameter$, (b) $h_{eff}^{unw-cal}parameter$.

Nevertheless, the results in Figure 10a show a phase-shift (arrow field oriented toward the left) and ranging from π to π/2 between ${\mathrm{\Phi }}_{unw}$ and SM measured by the probe at more or less than 1-day period. While the VSM estimated from $h_{eff}^{unw-cal}$ has a weak phase shift ranging from 0 to π/4 (arrow oriented ~horizontally toward the right) with VSM in situ. Low frequencies show high homogeneity with a high value of M² and no phase delay (horizontal vector field) between in situ and GNSS-R data for periods ranging from 8 days to a month.

Inflexion of the trend related to the transition between the rainy and the dry seasons is observed during the period ranging from 29 October 2016 to 11 November 2016 over 2-/4-day period with the probe at 10 cm (Figure 11a, b) and to a small extent with the one at 5 cm over a very temporally localized to 2-day period and only for the parameter ${ }_{unw}$ (Figure 10a). Figure 11 shows a significant phase shift ranging from π/2 and π/4 between in situ data and GNSS-R ones during this transition period. It means that GNSS-R detects SM variations a little bit before the probe at 10 cm records SM. Figure 11b shows a better correlation between $h_{eff}^{unw-cal}$ and the in situ P_10cm data during the transition and dry periods. In general, the $h_{eff}^{unw-cal}$ parameter presents a better correlation than ${\mathrm{\Phi }}_{unw}$. Globally, $h_{eff}^{unw-cal}$ parameter is more homogeneous in terms of phase shift and M² over all periods and for both in situ P_5cm and P_10cm measurements.

4.5 SM mapping based on GNSS-R data

As final results, Figure 12 shows SM time series obtained combining GLONASS and GPS and using of $h_{eff}^{unw-cal}$ parameter. These time series range between October and December 2016, i.e. from rainy to dry season. However, the VSM at P_5cm and P_10cm depth only changed 2% during the entire study period. At the end of the rainy season of 2016, we see that the rains are already almost non-existent. Figure 12 shows that less than 3cm of water fell during October to decrease to less than 1 mm in November and disappear entirely in early December. The SM is following the same trend with an apparent decrease from ~4.5% of volume moisture to almost 3% for P_10cm and even 2.5% for P_5cm. The transition period is marked by quasi-constant SM of ~3.7%.

Figure 12. Evolution of the final results for all the time series $h_{eff}^{unw-cal}$ – during period from October to December 2016. Results are obtained using Unwrapping phase algorithm for GPS+GLO satellites (dark blue dot) and compared to the in situ VSM measured by P_10cm (a) and P_5cm (b); $bluedotcorrespondto{h}_{eff}^{unw-cal}$ in m and the curve in red is a fit of $h_{eff}^{unw-cal}$ using a cubic spline. Rectangle blue correspond to the precipitation during this period.

The time series of $h_{eff}^{unw-cal}$ shows a very clear anti-correlation with SM. Indeed, when the SM decreases, the GNSS waves penetrate more strongly into these sandy soils, resulting in an increase of $h_{eff}^{unw-cal}$. We also notice that $h_{eff}^{unw-cal}$ has a minimum height that can be lower than the measured height of the ground, which is 1.61 m. This result is due to different factors. The first one is the multi-centimetric measurement noise which explains this phenomenon partly. Another one is the existence of some periodogram of a double peak indicating two possible surfaces of reflection. Moreover, a significant decrease in satellites visibility causes an increase in measurement error, as demonstrated by Zeiger et al. (2021).

Figure 13 shows SM maps created from the results of $h_{eff}^{unw-cal}$ in the study area. The specular points used to retrieve the SM maps have been computed thanks to the simulator developed by Roussel et al. (2014) which allows to do the accurate ray tracing for each satellite according to the following parameters: orientation and elevation of the satellite, tropospheric refractive indices, local topography around the antenna, antenna height and the altimetric reference surface chosen – here the ellipsoid of the WGS84 datum. The time step for calculating spatial variability is set to 3 days. SM content changed from 2.87% to 4.76% during the rainy, transition and dry periods from October to December 2016. The temporal variation of the SM maps shows that the SM content gradually decreases as it transitions to the dry period. This soil drying increases toward the southeast with a spatial variability lower than 0.01%. The results demonstrate that the unwrapping phase technique improves the accuracy of SM retrievals compared to phase measurements.

Figure 13. Dahra SM maps generated from the results of 3 days $h_{eff}^{unw-cal}$ during the rainy, transition and dry periods in 2016 by interpolating the specular reflection points (red dots) around the GNSS station (the ground surface is considered flat).

5 Discussion

In the study area, the GNSS-R signal used to retrieve soil geophysical parameters of the soil represents a more complex observation of soil surface moisture. In areas with very low SM, $h_{eff}^{unw}$ still has a good correlation but is negatively correlated due to wave penetration. The effective height time series shows the best correlation with SM after it is corrected by the unwrapping SNR method and filtered by the Sgoley method to reduce the influence of vegetation and surface roughness. Equation (5)(5) $h_{eff}^{unw}=\left({cos}^{-1}\left(\frac{SN{R}_m}{A_m}\right)-{\mathrm{\Phi }}_{unw}\right)\times \frac{\lambda }{\left(4\pi \times sin\left(\theta \right)\right)}$(5) shows that the $h_{eff}^{unw}$ has an anti-correlation with ${\mathrm{\Phi }}_{unw}$, and it could be explained as the increases of SM in any soil leads to the penetrations depth of GNSS signals decreases. Hence, as the effective height $h_{eff}^{unw}$ also decreases as demonstrated in previous studies using ${\mathrm{\Phi }}_m$ to estimate SSM and positively correlates with the in situ measurements (Roussel et al., 2016).

Another advantage of this method is the ability to provide spatial and temporal variability of the SM (Figure 13). The illuminated spatial area depends on a critical parameter: the antenna height. In this case, 1.61 m antenna height corresponds to a radius of the reflection area of about 12 m, an average area of 452 m² and a mean volume of ~100 m³. The comparison between GPS and GLONASS constellations shows that GPS satellites have a better correlation with in situ probes than GLONASS satellites. Here we show the effect of the signal noise of GLONASS constellations (lower SNR) which is more important than that of the GPS constellation. We choose to use both constellations to increase the number of specular points, and consequently, it improves over time to obtain long time series with a proper time step. It is possible to obtain two decades of data with time steps close to 30 s or less using the international GNSS network, which is essential to look at the effects of climate change, for example, and to highlight local and rapid phenomena like agricultural irrigation.

Another important aspect is that the antenna height adjustment can calibrate many space missions, e.g. a 2 m antenna height, with a coverage radius of few ten meters, could be useful to calibrate the Sentinel I, RADARSAT-2, ALOS-2 missions. On the other hand, a 60-m antenna with a coverage radius of few kilometers could be used to calibrate CYGNSS or SMAP data worldwide using GNSS stations of national (RGP for France; NGS for USA, CMONOC GNSS network for China, etc.) or international networks (IGS, EUREF, etc.). These stations are selected to identify those that cover either bare ground, agricultural fields or even forests. From the investigated surfaces, the results pointed out that the total of the first Fresnel surface depends on three main factors: the height of the antenna; the height of satellite ε, which can receive from broadcast ephemeris or precise one; and the many satellites visible at time t. So the different scales can be determined by dividing the visible satellites into low elevations and high elevations. Then, surfaces of a few square meters can be evaluated using the antenna approximately 5 m with an ε range of 35–70°. While surfaces from 100 to 10,000 m² can be analyzed using the 100-m antenna with ε < 35 °. However, this method also has disadvantages due to the maximum height of the antenna, which must be less than the wavelength of the C/A code (293 m); otherwise, there is the problem of the construction of the interference pattern (Ribot et al., 2016).

GNSS-R with the single antenna is considered an excellent alternative solution to complement the in situ technique in order to monitor SM for various types of soil: bare soil (Roussel et al., 2016), bare and vegetated soil (Chew et al., 2016), cobbly clay loam (Larson et al., 2010), etc., and even areas where the VSM is too low or rapidly changes. The in situ technique requires calibrated reference data by applying the reference technique in the laboratory condition (Walker et al., 2004). Various sensors depend on daily temperature variation like the Theta probe (Holzman et al., 2017). Other sensors like time domain reflectometry are less sensitive to thermal artifacts (Walker et al., 2004), while GNSS-R is more sensitive to thermal artifacts. Besides, GNSS-R integrates a larger volume and deeper penetration than the Theta probe method, and its signal is more correlated to SM at ground level as humidity increases. Roxy et al. (2010) indicate that daily moisture variations correspond to daily variations of surface albedo. These variations decrease from the dew point until the maximum sun illumination and reversely.

One of the challenges in improving the spatiotemporal resolution is that they are inversely proportional to each other. Various researches have justified the good results of the humidity calibration method using in situ probes in homogeneous soil (Holzman et al., 2017; Koyama et al., 2017). However for a non-homogeneous soil with strong variations in humidity, several in situ sensors have needed to monitor SM in a large area. GNSS-R allows a sensor to measure SM over a wide area to limit these gaps and reduces monitoring and human resources costs.

6 Conclusions

Long-term continuous SM data are essential in climatic research (Longobardi, 2008) and describe the region characteristics (Cleverly et al., 2016; Kirkby, 2016). Various automated techniques, namely in situ measurements and remote sensing, have been usable in the last years. The IPT is presented by Larson et al. (2008) justified an excellent retrieval of SM for the soil rich in clay (Chew et al., 2014; Vey et al., 2016; Zhang et al., 2017). Our study developed a novel combination between technical IPT and phase unwrapping method to obtain a continuous phase measurement of the SNR multipath, used to monitor SM even for dry to very dry sandy soils (2.5% <VSM <8%). In this case, the correlation coefficient between SM variations in our technique and VSM observed in situ at ~10 cm depth is 0.97 over the study area (95% of sand). Furthermore, the combination of all GPS/GLONASS/GALILEO/BEIDOU satellites in the SNR data treatment can be applied to generate a complete series of time precision, which considerably increase the spatial and temporal resolution from 1 day to 10 min.

Currently, GNSS networks exist in various countries and can be applied GNSS-R tool for free. Another critical point is the number of GNSS satellites which increase continuously from 1990 (~40 compared to more than 120 in 2030 (Gao & Enge, 2012)). This factor increases the spatiotemporal monitoring of the Earth, except in the pole areas due to the structures of GNSS orbital planes. GNSS-R SM maps can be established continuously worldwide and can provide powerful time series. Monitoring SM by GNSS-R will be a crucial parameter to improve the understanding of spatial variability of droughts at all spatial and temporal scales in the context of global climate change.

Acknowledgments

This research was funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 105.08-2020.17, by the University of Toulouse, and by the foundation of Thematic Network for Advanced Research (RTRA) in Science and Technology for Aeronautics and space (STAE) under the MOMMHYS (Mesures d’Opportunités pour la Modélisation Météo-Hydro-Surfaces continentals) project.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author, upon reasonable request.

Additional information

Funding

This work was supported by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) [105.08-2020.17].