Single-Pass Soil Moisture Retrievals Using GNSS-R: Lessons Learned

Abstract

In this paper, an algorithm to retrieve surface soil moisture from GNSS-R (Global Navigaton Satellite System Reflectometry) observations is presented. Surface roughness and vegetation effects are found to be the most critical ones to be corrected. On one side, the NASA SMAP (Soil Moisture Active and Passive) correction for vegetation opacity (multiplied by two to account for the descending and ascending passes) seems too high. Surface roughness effects cannot be compensated using in situ measurements, as they are not representative. An ad hoc correction for surface roughness, including the dependence with the incidence angle, and the actual reflectivity value is needed. With this correction, reasonable surface soil moisture values are obtained up to approximately a 30° incidence angle, beyond which the GNSS-R retrieved surface soil moisture spreads significantly.

1. Introduction

The first evidence of soil moisture (SM) signatures in GPS reflected signals date back to 1996 [1]. Exploiting this technique, accurate soil moisture retrievals were successfully conducted using ground-based instruments using the so-called “interference pattern technique” (IPT) in which the direct and reflected signals are collected simultaneously by the receiving antenna, creating fringes in the received power as navigation satellites elevation varies over time. This technique can be applied using geodetic zenith-looking right-hand circularly polarized (RHCP) antennas [2], and the reflections are collected at low elevation angles for which the reflected signal is also RHCP; or it can be applied using linearly polarized (vertical, or vertical and horizontal) antennas pointing to the horizon [3,4] which allows collecting reflected signals over a much wider range of elevation angles. The latter technique also allows inferring vegetation height, and with two linear polarizations, surface roughness effects can be corrected. However, none of these techniques can be applied to airborne or spaceborne instruments.

Retrieving soil moisture from airborne and spaceborne instruments can only be performed using the reflected signal power, which requires on the instrument side the compensation of platform attitude impact on the antenna pattern, receiver’s noise… and on the target side the compensation of surface roughness (small scale as opposed to topography or large scale), and vegetation effects (both attenuation and possibly scattering). So far, two main approaches have been presented in the literature: using a change detection approach by looking to a time-series of GNSS-R data, or attempting a single-pass retrieval.

The time-series approach [5,6,7] has been able to produce operationally a soil moisture product using NASA CYGNSS (Cyclone GNSS) mission data [8]. The main advantages of the time-series approach are that assuming surface roughness and vegetation effects are constant within a grid cell (36 km in [8]), and the change in the observed reflectivity is only due to the change in the soil moisture. The main drawbacks are the difficulty to really achieve collocated GNSS-R observations over time as in the case of land, GNSS reflections exhibit a very high spatial resolution (several hundreds of meters even from space, since the spatial resolution is given by the size of the first Fresnel zone [9]), the implicit assumption of the temporal coherence of data (i.e., surface roughness or vegetation can change), and the strong non-linearities of the sensitivity of GNSS-R observables to soil moisture, which is very noticeable at low soil moisture values (Figure 1). To overcome this, in [8], an ad hoc calibration had to be performed as “the product was developed by calibrating CYGNSS reflectivity observations to soil moisture retrievals from NASA’s Soil Moisture Active Passive (SMAP) mission”.

On the other hand, single pass retrievals [10,11,12,13,14] could provide a faster revisit time, as SM maps could be generated in just one pass, provided that the instrument is accurately calibrated, and suitable ancillary data are used to compensate for surface roughness and vegetation effects without tuning the algorithm for each pixel. In [10,11], the sensitivity to soil moisture was investigated, together with a first attempt to compensate for vegetation effects using the Normalized Differential Vegetation Index (NDVI). Surface topography, spatial scales, and some subsurface “artifacts” were already reported using spaceborne data from the UK TechDemoSat-1 (TDS-1) data. Using airborne data, in [12], the sensitivity to soil moisture was studied, as well as vegetation effects, which were parameterized in terms of the Leaf Area Index (LAI) and a single scattering albedo. More recently [13], IceSat-2 data have been proposed to correct for surface roughness effects and NASA SMAP data to correct for surface roughness and vegetation optical depth (VOD) effects in the Fresnel reflection coefficients at incidence angles smaller than 35°. Even, more recently, Machine Learning techniques have been applied [14].

Despite the encouraging results of these previous works, to the authors’ knowledge, the correction of surface roughness and vegetation effects are not yet properly understood so as to perform an accurate soil moisture retrieval. In this work, the goodness of these corrections are analyzed using Sentinel-2 satellite data; GNSS-R, L-band microwave radiometry, VNIR (visible and near-infrared) hyperspectral and TIR (thermal infrared) airborne data, as well as in situ surface roughness, soil type, and soil moisture samples aiming at performing single-pass surface soil moisture (SM) retrievals using GNSS-R observations.

In order to address the above questions, an incremental approach was followed. First, based on our previous experience, very different sensitivities to soil moisture were encountered, from as high as 38 dB/100% for bare dry soils [10,15] to a more moderate approximately 9 dB/100% at global scale [11]. This large variability has to be understood first in order to properly model the sensitivity to soil moisture in the retrieval algorithm.

As it can be appreciated in Figure 1, computed using a comprehensive GNSS-R simulator [16], assuming a flat soil surface and no vegetation, the variation of the peak of the DDM (Delay and Doppler Map) with respect to SM is highly non-linear. Depending on the range of soil moisture values encountered in the field experiments, different slopes can be found if a linear approximation is assumed to compute the sensitivity. For example, in the range 0 to 0.15 m³/m³, the sensitivity is as high as 43 dB/100%, while in the range 0.15 to 0.45 m³/m³, the sensitivity is much lower 8.6 dB/100%, and if the whole dynamic range is considered 0 to 0.45 m³/m³, the sensitivity is 17 dB/100%.

Actually, these values demonstrate that a linearized (in dB) dependence ΔΓ/ΔSM is not valid for the whole range of SM values, and it cannot be used in SM retrieval algorithms. In addition, time-series algorithms based on the detection of reflectivity changes must also include this non-linearity as depending on the actual SM value, a given reflectivity change will correspond to different soil moisture change.

Second, the retrieval algorithm from GNSS-R observations (airborne and/or spaceborne) proceeds as follows:

• From the peak of the DDM (direct and reflected signals), and assuming that the scattered signal has a dominant coherent component, the calibrated reflection coefficient is computed taking into account the antenna pattern gains in the directions where the signals are collected, and the distances from the transmitter to the receiver, and from the transmitter to the specular reflection point plus from the specular reflection point to the receiver.
• Then, following an approach similar to the Single Channel Algorithm in SMAP [17], or in [10,12], but neglecting the single scattering albedo (ω = 0), the computed reflectivity is compensated for surface roughness (h) and vegetation effects (τ: vegetation optical depth, or VOD), in order to derive the flat surface reflectivity:

$${\Gamma }_{bare \left[dB\right]}={\Gamma }_{\left[dB\right]}+10log(e^{h·co{s}^{\mathrm{n}}\left({\theta }_i\right) })+10\log \left(e^{2·\tau /cos \left({\theta }_i\right) }\right),$$

(1)

where θ_i is the incidence angle, and n is an empirical factor to account for the angular dependence of the surface roughness. In theory [18,19], it should be n = 2. However, in SMAP, the value used is n = 0. This effect is discussed afterwards. The factor “2” in the second exponential term of the VOD accounts for the two-way attenuation through the vegetation canopy. Considering the exponential dependence with h, this term is sometimes expressed in dB. Note that at circular polarization (LHCP: left hand circular polarization), at a constant frequency (f = 1575.42 MHz), the “flat surface reflectivity” (${\Gamma }_{\left[dB\right]}$) has a marginal variation with the incidence angle up to approximately 30° to 40°, and the largest variability is linked to the dielectric constant.

3.

Then, considering the variation of the SM in the dielectric constant model for soil (among other variables, such as the clay fraction and physical temperature…), the soil moisture value can be in principle estimated through a minimization process.

The volumetric SM retrieval process is sketched in Figure 2.

At this stage, it is important to remind that in SMAP, there are three main soil moisture retrieval algorithms using microwave radiometer data: the Single Channel Algorithm (SCA), which uses the brightness temperature observations at either horizontal or vertical polarizations, and a vegetation climatology to retrieve soil moisture; and the Dual Channel Algorithm (DCA), which uses the brightness temperature observations at both horizontal and vertical polarizations, while inferring vegetation parameters as well. SMAP also uses static (climatologic) ancillary data such as soil texture and land cover class. The SCA using the brightness temperature at vertical polarization is the current baseline algorithm for SMAP level 2 soil moisture products. A recent study [20] shows that “while the DCA has its lowest unbiased root man squared error (ubRMSE) and highest R² when ω is non-zero, the SCA have their lowest ubRMSE and highest R² when ω = 0, and the dry bias of all algorithms increases as ω increases. Errors in soil texture are not significant, but soil surface roughness should not be static and have a higher overall value.”

When in point 2 above, it is stated that an algorithm similar to the SMAP SCA is used (with ω = 0 as found in [20]), the authors refer to the functional relationships used to account for surface roughness and vegetation effects (Equation (1)), but not to the use of SMAP climatology and ancillary data. Actually, in situ, aiborne and satellite data were acquired to parameterize the vegetation and surface roughness effects, as described in Section 2.

The rest of this work is organized as follows: Section 2 describes an airborne field experiment that was conducted to assess the SM retrieval capabilities of GNSS-R. Section 3 shows the results of the SM retrieval algorithm performance and the limitations of the surface roughness and vegetation corrections. In Section 4, additional considerations for SM retival approaches are discussed. Finally, Section 5 summarizes the main conclusions of this study.

2. Methodology: Field Experiment Description

In 2017, a first airborne flight was conducted in the Tremp area (42.2°N, 0.89°E) of Lleida, Spain. The area was selected because of the availability of soil moisture stations for validation purposes. However, topography effects and the strong wind gusts affected the aircraft trajectory and attitude and did not allow for conducting the SM retrievals successfully. A second airborne campaign was conducted on October 22nd, 2018, in a flatter agricultural region around Balaguer (41.9°N, 0.80°E), Lleida, Spain (Figure 3). In situ and satellite data were acquired to complement the airborne data and allow the estimation of the surface roughness and vegetation optical depth, necessary to perform the vegetation corrections. The acquired data and the sensors are presented below.

2.1. Airborne Instrumentation and Configuration

The plane used was a CESSNA Caravan (Figure 4a) from the Institut Cartogràfic i Geològic de Catalunya, with an autonomy of 4 h. The flight height was 1200 m. The instrumentation used was the following:

• AISA [21]: VNIR (406–993 nm) pushbroom hyperspectral sensor with a maximum of 254 bands and a spectral width of 2.4 nm (Figure 4b). The sensor was configured with 126 bands, combining together every two bands leading to an average spectral width of approximately 4.8 nm.
• TASI [22]: TIR (8050–11450 nm) hyperspectral with 32 bands (Figure 4b), and an average spectral width of 109.5 nm.
• ARIEL [23]: L-band radiometer with central frequency at 1.413 GHz and bandwidth of 20 MHz (Figure 4c), and a radiometric accuracy of 0.7 K at 1 Hz. ARIEL was designed by the Universitat Politècnica de Catalunya and is commercialized by BALAMIS [24].
• CORTO (COmpact Reflectometer for Terrain Observations) is a miniature version of the LARGO [25] reflectometer at L1 GPS bands developed by the Universitat Politècnica de Catalunya. It performs the correlations between the direct and the reflected GPS signals and outputs reflectivity.

The APPLANIX [26] Inertial Navigation System (INS) on-board is a precise GNSS receiver + IMU (Inertial Measurement Unit) System that was used to geo-locate the data acquired and to perform the geometric correction on the images.

In order to receive the GPS signal for the GNSS-R instrument, there are two antennas: an RHCP up-looking GPS antenna already installed on the aircraft for navigation purposes and a LHCP down-looking one that covered the L1/L2 GPS bands. To avoid picking the direct signal through a surface wave propagating around the plane fuselage, a choke ring antenna with a carefully designed mounting was employed. The antenna radiation pattern in that configuration was measured in the UPC anechoic chamber [27] using a dummy of the bottom part of the fuselage.

2.2. Vegetation Optical Depth (VOD) Estimation from Normalized Differential Vegetation Index (NDVI)

VOD is needed to compensate for the vegetation effects (Equation (1)). At the L-band, and taking into account that it is an agricultural area, with some trees, but no dense forests, the SMAP approach can be in principle used [10,12,13]. High-resolution (0.8 m) NDVI was computed from AISA data. The Sentinel-2 L1C NVDI product corresponding to Octobr 23rd, 2018 was also downloaded from the Sentinel Hub [28] (Figure 5). The VOD (τ) can be computed as:

$$\tau =b·VWC,$$

(2)

where b is a proportionality constant that depends on the vegetation structure and the frequency, and VWC is the vegetation water content, which can be calculated using the following relationship:

$$VWC=b·1.9134·NDV{I}^2-0.3215·NDVI+\frac{StemFactor\left(NDVI-0.1\right)}{0.9},$$

(3)

where the NDVI is given by the ratio of the difference and the sum of the reflectivities in the near infrared (NIR) and the red:

$$NDVI=\frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+{\rho }_{Red}}.$$

(4)

The Stem Factor in Equation (3) is the product of the average vegetation height of a land cover class and the ratio of sapwood area to leaf area [17]. Sample values of b, StemFactor, and h are given in Table 3 of [17] for the land use of each observation according to the cadaster, and these are confirmed by in situ observations.

At this stage, it is important to highlight two potential error sources in these parameterizations:

• that the value of the b parameter directly determined at the satellite scale “is nearly identical to what was proposed for crops in the ESA SMOS (Soil Moisture and Ocean Salinity mission) algorithm, but half as large as what is currently used by SMAP” [29], and
• the current values of the h parameter may be too smooth (low) for crop regions, as reported in [20] for the South Fork SMAP Core Validation Site in the Corn Belt state of Iowa.

AISA NDVI was calculated using a narrowband (central wavelength) and a wideband approach, with a parametric adjustment to Sentinel-2 spectral sensitivity with respect to the wavelength.

Figure 6 shows the difference between Sentinel 2 NDVI and AISA NDVI computed either with the narrowband or wideband approaches. As it can be appreciated, the narrowband approach introduces larger errors (overestimation/underestimation) when compared to Sentinel 2 NDVI, while results with the wideband approach match much better Sentinel-2 results. As a result of the higher spatial resolution, the AISA NDVI using wideband data is used in Equations (3) and (4) to try to compensate for the vegetation effects in Equation (1). The bottom panel in Figure 6 shows the TASI thermal image in [°C]. The calibrated physical land surface temperature (LST) is lower in the bottom right corner due to the presence of clouds.

2.3. Soil Surface Roughness and Soil Moisture

A field campaign to validate the land use and to perform the surface roughness, soil sampling, and in situ soil moisture measurements was performed the same day as the airborne campaign. Before the experiment, the Hydrogeology and Soils Unit of the Catalan Institute of Cartography and Geology (ICGC) studied the region and selected 5 different areas representative of the different agricultural activities encountered. The locations of the 29 sampling points are indicated in Figure 7.

Soil surface roughness was measured with a laser profiler based on the Leica Disto Pro4a laser distance meter, with millimetric precision, in 10 mm steps. Soil Moisture samples were acquired with a Theta Probe sensor, acquiring at least 3 samples per “sampling point” and computing the median value. Then, soil samples were analyzed in the laboratories of the IRTA (Instituto de Investigación y Tecnología Agroalimentaria)/Lleida University facilities and used for calibration of the Theta Probe sensor.

The largest problem with the “effective” surface soil roughness is its difficult estimation, as it depends on numerous additional factors such as the soil composition, the agricultural activity, the state of the crops, etc. [30]. Moreover, this seasonal dependency only correlates with similar fields.

In the aircraft, the ARIEL L-band radiometer was flown to provide a complete 0.8 m resolution SM map (Figure 8) using a downscaling algorithm that combines TASI and AISA data [31,32]. This algorithm has been validated [33] extensively, and in this experiment, the root mean square error (rmse) and bias were approximately 0.06 m³/m³, and −0.01 m³/m³ with respect to the in situ soil moisture measurements. As a qualitative validation, Figure 9 shows the Sentinel 2 Normalized Difference Water Index (NDWI) for October 23rd, 2018. Similarly to the NDVI, the NDWI index uses the green and near infrared bands (bands 3 and 8), and because of the strong absorbability and low radiation in the range from visible to infrared wavelengths, it is most appropriate for water body mapping, to assess the surface soil moisture and vegetation hydric stress. Figure 10 shows the GNSS-Reflectivities overlaid with the Sentinel-2 L1C NDWI.

3. GNSS-R Soil Moisture Retrieval: Results

In this section, the resulting SM maps are assessed for different corrections applied. The results are shown over an ortho-photo of the Region Of Interest (ROI). There are two layers of data: ARIEL downscaled SM (0.8 m spatial resolution), and overlapped samples of the GNSS-R retrieved SM at the same resolution. The ARIEL SM retrievals are used to validate the estimations of GNSS-R SM retrievals.

3.1. In Situ Surface Roughness Correction

In this first attempt, the measured surface roughness with the laser profiler was used to compute the roughness parameter required in Equations (1) and (11) of [18] from the soil surface roughness (σ):

$$h=4 k^2{\mathsf{\sigma }}^2.$$

(5)

Due to the limited accuracy and coarse sampling of the soil surface roughness, the results show a large disparity between the ARIEL and the GNSS-R SM retrievals. As it can be appreciated in Figure 11a, for a single track of GNSS-R reflections over a relatively flat area (no topography), the retrieved SM values are highly variable with the incidence angle (see alternating SM retrievals between approximately 0.14 and 0.27), and these are quite different than the ARIEL downscaled SM ~0.21, although on average, it is relatively close (e.g., 0.21 ≈ (0.14 + 0.27)/2). Apart from the different scales (0.8 m versus 10 m), note the high visual correlation between Figure 11a, ARIEL downscaled soil moisture, and Figure 11b Sentinel 2 L1C NDWI.

3.2. Ad Hoc Soil Surface Roughness Correction

The lack of success in this first SM retrieval approach with the in situ measured soil surface roughness led us to consider that the surface roughness correction applied was not appropriate, either in the value of the h parameter, and/or the power of the cosine function (n = 0, 1, or 2) in Equation (1). As discussed in page 35 of [17], the cos²θ term is sometimes changed to cosθ or even dropped completely to avoid overcorrecting for roughness.

To derive the “correct” surface roughness correction, the same procedure was applied, but adding an ad hoc surface roughness correction term h, as in $e^{h·co{s}^{\mathrm{n}}\left({\theta }_i\right) }$, so that the retrieved GNSS-R SM values match as close as possible the ARIEL-derived ones. The step-by-step procedure is described below:

• Compensation of vegetation effects (Equations (1)–(4)).
• Compensation of surface roughness effects (as in Equation (1) and (2): $10·log(e^{h·co{s}^{\mathrm{n}}\left(\theta i\right) })$, for n = 0, 1, and 2.
• Apply SM retrieval algorithm by minimizing the difference between the computed flat bare soil reflectivity using ARIEL-derived SM, and GNSS-R calibrated reflectivity (${\Gamma }_{\left[dB\right]}$) plus the ad hoc surface roughness compensation term (in [dB]).

It was found that for n = 0, the scatter in the ad hoc surface roughness corrections was the smallest one obtained. For n = 2, the corrections were more exaggerated, indicating they were artificially high. Figure 12 shows the map of the computed surface roughness parameters (h). Apparently, there is no evident correlation neither with the different fields nor with the soil moisture values, although some of the largest values occur in the plots with higher humidity, but not all. In addition, as expected, by construction, the matched histograms of the SM values retrieved by downscaling ARIEL’s measurements and from CORTO’s reflectivity with the ad hoc surface roughness correction are basically the same (Figure 13), except for numerical round off errors in the surface roughness correction (0.1 dB) during the optimization process (third bullet in the above list, or the right box in the diagram in Figure 2).

Figure 14 shows the estimated values of h (as in $e^{-h}$). A more detailed analysis of the ad hoc surface roughness correction reveals that it increases for lower reflectivity values. This apparently has no physical meaning, and it is difficult to interpret as an instrumental error. It was also considered that it could be an imperfect estimation of the noise floor of the DDM background, which is needed in the calibration, but this does not explain why, for the same reflectivity values, the correction varies with the incidence angle, increasing for the larger ones. Additionally, negative values in h (roughness correction larger than 0 dB at low reflectivities) do not have a physical meaning, but they led us to think that actually the compensation of the vegetation effects as in Equations (1)–(4) was too high, as already shown in Figure 13 of [11]. This effect is examined in more detail in the next section.

3.3. Vegetation and Roughness Compensation for Different Incidence Angles

Figure 15 shows the scatter plot of the vegetation compensation [dB] (Equation (1)) with respect to the incidence angle. The compensation of vegetation effects increases with incidence angle, as it has a functional dependence with $1/cos\left({\theta }_i\right)$ (see the variation of peak values wrt. incidence angle). The compensation of the roughness effects has also a dependence on the incidence angle, but it can be modeled neither with n = 1, nor with n = 2.

Since the compensation of vegetation effects is fixed and depends on the VOD (which is calculated using the NDVI) and the incidence angle, an overcompensation of the vegetation attenuation forces the algorithm to reduce the compensation of the surface roughness component, and that is why the algorithm even finds negative surface roughness coefficients. This effect is more important for large incidence angles, and it suggests that the vegetation compensation is too high at large incidence angles, which is in agreement with [13], who found that corrections were only applicable up to approximately a 35° incidence angle.

Other sources of variability are:

• high reflectivity values encountered in vegetated areas (trees, higher soil moisture), but from the processing point of view, require a too high attenuation compensation, leading to “corrected” flat bare reflectivity values larger than one (0 dB) in Figure 14. This is illustrated in Figure 16 by the dark blue dot over a forested area in the middle of the image, pointed by the triangle, and legend with the parameters.
• there is likely a subsurface variability of the soil moisture that is detected by GNSS-R, as there are regions in the fields where the reflectivity suddenly increases (see lines of light blue or green values in Figure 10, with some dark blue—high reflectivity—observations in the middle), while all other parameters are nearly constant.

3.4. Empirical Roughness Compensation

Not being satisfied with the result that accurate single-pass soil moisture retrievals were only feasible if a pixel-based ad hoc correction of the surface roughness was used, an empirical correction (1st order polynomial) for all pixels within a given range of incidence angles is attempted. The procedure is the same as in the pixel-by-pixel approach described before, but now computing the minimum root mean squared error (RMSE) for all pixels simultaneously. Figure 17 shows the histograms of the ARIEL-retrieved SM values, and the GNSS-R retrieved ones using a generic empirical 1st order polynomial correction for all pixels within a given range of incidence angles. As it can be appreciated, except for a few outliers, the matching is quite reasonable up to approximately a 30° incidence angle, and beyond that, the retrieved histograms widen and GNSS-R derived SM values tend to widen and scatter significantly (sensitivity to SM decreases, while the effect of roughness is still non-negligible). This result is also in agreement with [13] where corrections were only applicable up to an approximately 35° incidence angle.

Figure 18 shows the soil moisture maps for two regions (left and right column, corresponding to regions A and B, respectively in Figure 8, Figure 9 and Figure 10) during the flight for all ranges of incidence angles in steps of 15°. Squares correspond to GNSS-R-derived soil moisture using the generic angular-dependent roughness correction.

Finally, Figure 19 presents the downscaled ARIEL and GNSS-R-retrieved SM overlaid maps for all incidence angles included, showing a reasonable agreement, which degrades for increasing incidence angles.

4. Discussion

In the previous sections, it was shown that the sensitivity of the calibrated reflectivity to Soil Moisture is non-linear and varies with the Soil Moisture value (Figure 1) [16]. Linearized retrieval approaches must be cautious in order not to introduce biases and other artifacts. Soil moisture retrieval algorithms require that vegetation and surface roughness effects must be properly compensated first.

Current models adopted from L-band microwave radiometry (1.413 GHz) overcorrect vegetation effects [11], which leads to an overcorrection of the surface roughness to try to compensate for the vegetation correction. Therefore, techniques to better predict the vegetation corrections are needed, as current models adopted from L-band microwave radiometry (1.413 GHz) using a VOD estimate from the NDVI [17] seem to overestimate it, which leads to an overcompensation of the surface roughness effects to counteract them. A potential source of this discrepancy may be in the values of the b parameter, which are too high in crop regions, as recently reported in [29].

In situ measurements of the surface roughness are not suitable to predict the surface roughness parameter to be used in the compensation, as different roughness values are obtained with different techniques [34], i.e., a laser profiler with respect to mechanical profiler, and that there is an impact of the profiler length and number of independent measurements. It is also possible that the reflection is taking place effectively in the water table underground, behaving then as a much “flatter” surface than the air–soil interface. Volume scattering effects in the soil cannot be neglected, either. The use of a constant value of “h” (e.g., as in SMAP) is not suitable either, even if it is different for each field. In addition, it has been recently reported that the current values of the h parameter may be too smooth (low) for crop regions [20]. Generic surface roughness empirical corrections as a function of the range of incidence angles are also possible, leading to acceptable results: the root mean squared error between downscaled ARIEL and GNSS_R SM map is < 5%, while today’s UCAR/CU soil moisture product developed by calibrating CYGNSS reflectivity observations to soil moisture retrievals from NASA SMAP mission exhibits a median unbiased root-mean-square error (ubRMSE) of 0.049 m³/m³ [8].

5. Conclusions

Recent results [13] claim that IceSat-2 surface roughness data can be used to correct for these effects at a global scale. However, the authors of this study believe that further improvements using auxiliary data from other sensors are unlikely, as it is very difficult to match all the observations in time and space, and that they can vary over time, and techniques to infer the effective soil surface roughness parameter, if possible from the data themselves, are needed. A clear example is the estimation of the surface roughness from the GNSS-R observables [35].

Further research is also needed to explore the simultaneous use of several GNSS bands and other GNSS-R observables with variable coherent and incoherent integration times, and to include the non-linear dependence of the VOD with the above ground biomass (ABG), as shown in Supplementary Figure 4 from [36,37]. Vegetation scattering and wave depolarization have to be accounted for as well.

In the meantime, while these corrections are not improved for single pass retrievals, time-series approaches based on temporal changes of the reflectivity (e.g., [5,6,7,8]) can be applied, provided that the non-linearities of the dependence of the GNSS-R observables and the SM are properly taken into account, and self-consistency of the retrieved products is demonstrated.

The most relevant result of this study is that despite all the above limitations in the surface roughness and vegetation effects corrections, single-pass GNSS-R soil moisture values seem feasible with a reasonable error (<5%) for incidence angles up to 30°, provided that an incidence angle dependent soil surface roughness correction is applied.