First Evaluation of Topography on GNSS-R: An Empirical Study Based on a Digital Elevation Model

Abstract

Understanding the effects of Earth’s surface topography on Global Navigation Satellite Systems Reflectometry (GNSS-R) space-borne data is important to calibrate experimental measurements, so as to provide accurate soil moisture content (SMC) retrievals. In this study, several scientific observables obtained from delay-Doppler maps (DDMs) $⟨{{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{o}}(\mathsf{\tau },\mathrm{f})\right|}}^2⟩$ generated on board the Cyclone Global Navigation Satellite System (CyGNSS) mission were evaluated as a function of several topographic parameters derived from a digital elevation model (DEM). This assessment was performed as a function of Soil Moisture Active Passive (SMAP)-derived SMC at grazing angles ${\mathsf{\theta }}_{\mathrm{e}}$~ [20,30] ° and in a nadir-looking configuration ${\mathsf{\theta }}_{\mathrm{e}}$~ [80,90] °. Global scale results showed that the width of the trailing edge (TE) was small $\mathrm{T}\mathrm{E}$~ [100, 250] m and the reflectivity was high $\mathsf{\Gamma }$~ [–10, –3] dB over flat areas with low topographic heterogeneity, because of an increasing coherence of Earth-reflected Global Positioning System (GPS) signals. However, the strong impact of several topographic features over areas with rough topography provided motivation to perform a parametric analysis. A specific target area with little vegetation, low small-scale surface roughness, and a wide variety of terrains in South Asia was selected. A significant influence of several topographic parameters i.e., surface slopes and curvatures was observed. This triggered our study of the sensitivity of $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ to SMC and topographic wetness index ($\mathrm{T}\mathrm{W}\mathrm{I}$). Regional scale results showed that $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ are strongly correlated with the $\mathrm{T}\mathrm{W}\mathrm{I}$, while the sensitivity to SMC was almost negligible. The Pearson correlation coefficients of $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ with $\mathrm{T}\mathrm{W}\mathrm{I}$ are ${\mathrm{r}}_{\mathsf{\Gamma }}$~ 0.59 and ${\mathrm{r}}_{\mathrm{T}\mathrm{E}}$~−0.63 at ${\mathsf{\theta }}_{\mathrm{e}}$~ [20, 30] ° and ${\mathrm{r}}_{\mathsf{\Gamma }}$~ 0.48 and ${\mathrm{r}}_{\mathrm{T}\mathrm{E}}$~ −0.50 at ${\mathsf{\theta }}_{\mathrm{e}}$~ [80, 90] °, respectively.

1. Introduction

Soil moisture content (SMC) is an important component in the global water and energy balance because it determines the flux of water that infiltrates to groundwater or drains via surface water [1]. SMC determines the water and energy available for evapotranspiration. The dependence of SMC on topographic convergence suggests that topography is an important element in the spatial distribution of SMC [2]. The state and pattern of SMC affect the hydrological response of a catchment area [1]. Temporal dynamics in the state of SMC depend on rates of outgoing horizontal and vertical fluxes, e.g., drainage and evaporation. On the other hand, the dependence of SMC on topographic convergence [2] suggests that topography is an important factor that partially determines the spatial distribution of SMC. As such, processes that lead to spatial organization in SMC could be represented by topographic parameters derived from a digital elevation model (DEM). Indeed, topography-based wetness indices can be used to predict the spatial distribution of wetlands because topography is one of the most important elements that influences the spatial patterns of saturated areas. Alternatively, wetlands mapping and spatial patterns of wetness can be generated in catchments under the assumption that groundwater tables follow topography holds.

The future of Earth observation in SMC determination is strongly based on the ability of the next generation of space-borne sensors to improve the Earth’s surface spatio-temporal sampling properties. Global Navigation Satellite Systems Reflectometry (GNSS-R) appears in this framework, with an intrinsic capability to advance over more classical remote sensing approaches. GNSS-R is a sort of L-band passive bistatic radar (as many transmitters as navigation satellites are in view) that provides a wide swath up to ~ 1500 km. It was first conceived as an interferometric technique or iGNSS-R to provide higher precision in ocean altimetry measurements [3,4]. However, it can also be used for land surfaces applications. In this scenario, the conventional GNSS-R technique or cGNSS-R is more appropriate because lower coherent and incoherent integration times are required, so that the associated spatial resolution is improved [5]. At present, several works have demonstrated the feasibility of SMC estimation using space-borne cGNSS-R data from the United Kingdom (U.K.) TechDemoSat-1 [6,7,8], Soil Moisture Active Passive (SMAP) [9], and Cyclone Global Navigation Satellite System (CyGNSS) [10,11] missions.

CyGNSS was selected by the National Aeronautics and Space Administration (NASA) in 2012 as a highly innovative Earth Venture space system, and it was launched into space on 15 December 2016 [12]. CyGNSS is an 8-microsatellite constellation. Each single CyGNSS satellite carries on board two down-looking antennas with left hand circular polarization (LHCP). Both antennas point to the Earth’s surface with an elevation angle ${\mathsf{\theta }}_{\mathrm{e}}$~ 62 ° on each side of the satellite ground track (${\mathsf{\theta }}_{\mathrm{e}}$ ~ 90 °at the normal direction to the Earth’s surface). The gain of the antennas is ~ 14.7 dB (antenna boresight).

Topographic inhomogeneities can distort GNSS-R measurements [13,14]. They have to be compensated to improve the accuracy in the retrieval of SMC. In a bistatic radar geometry, GNSS scattered signals can be collected along the specular direction (elevation angle of the incident wave equals elevation angle of the scattered wave i.e., ${\mathsf{\theta }}_{\mathrm{e},\mathrm{i}}$=${\mathsf{\theta }}_{\mathrm{e},\mathrm{s}}$=${\mathsf{\theta }}_{\mathrm{e}}$). The scattering is only strong over an area around the nominal specular point. On the other hand, the scattering area depends on the local topography as well as on the GNSS satellites’ elevation angles. Local slopes modify the scattering area as compared to that corresponding to a flat Earth assumption. The scattering area may be a complex and sometimes non-contiguous area on ground. However, it remains uncertain to what extent topography influences GNSS-R data and how these effects can be corrected. The main purpose of the present investigation was to further understand the effects of rough topography using data collected by CyGNSS. The counterpart of the radar pulse width in CyGNSS is ~ 300 m (1 C/A code chip). As such, reflected delay Doppler maps (DDMs) are potentially sensitive to the relative orientation of the scattering surfaces, but also to the vertical relief distribution. In the first part of this study, the influence of six different topographic parameters derived from the Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010) model on DDMs derived from CyGNSS Level 1 Science Data Record Version 2.1 [12,15,16] was evaluated as a function of SMC derived from SMAP-enhanced L3 Radiometer Global Daily 9 km Level L3 SPL3SMP_E Version 1.0 product. In the last part, the following question was addressed: Can topographic features determine the main spatial pattern of Earth’s surface reflectivity as it is measured by CyGNSS?

The remainder of this paper is organized as follows. Section 2 provides a theoretical background on the effects of topography on GNSS-R. Section 3 describes the methodology. Section 4 describes the effects of topography on CyGNSS-derived DDMs based on several parameters derived form a DEM. Additionally, it investigates the relationship between DDMs and a topography-based surface wetness index. Section 5 includes final discussions, and the main conclusions are discussed in Section 6.

2. Theoretical Background

CyGNSS allows sampling of the Earth’s surface along 32 tracks simultaneously, within a wide range of satellites’ elevation angles ${\mathsf{\theta }}_{\mathrm{e}}$~ [20, 90] ° over the tropical latitudinal band Latitude ~ [−40, 40]°. Geophysical parameters retrieval from CyGNSS requires that the GNSS-R payload cross-correlates directly the collected signals against a known replica of the GPS L1 C/A code, after compensating for the Doppler shift. Consequently, 32 complex reflected DDMs,${\mathrm{Y}}_{\mathrm{r}}$, are generated on board at a time, using a short coherent integration time,${\mathsf{{\rm T}}}_{\mathrm{c}}$~ 1 ms. Complex DDMs, ${\mathrm{Y}}_{\mathrm{r}}$, are incoherently averaged along ~ 1 s to finally generate the so-called power DDMs.

Power DDMs, $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$, can be modelled using geometrical and scattering related parameters as follows [17]:

$$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle =\frac{{{\displaystyle {{\displaystyle \mathrm{T}}}_{\mathrm{c}}^2\mathrm{P}}}_{\mathrm{T}}{{\displaystyle \mathsf{\lambda }}}^2}{{{\displaystyle (4\mathsf{\pi })}}^3}{\displaystyle \iint \frac{{{\displaystyle \mathrm{G}}}_{\mathrm{T}}{{\displaystyle \mathrm{G}}}_{\mathrm{R}}{{\displaystyle \left|\mathsf{\chi }(\mathsf{\tau }-({\mathrm{R}}_{\mathrm{T}}+{\mathrm{R}}_{\mathrm{R}})/\mathrm{c},\mathrm{f}-{\mathrm{f}}_{\mathrm{c}})\right|}}^2{\mathsf{\sigma }}^0}{{{\displaystyle \mathrm{R}}}_{\mathrm{T}}^2{{\displaystyle \mathrm{R}}}_{\mathrm{R}}^2}}{{\displaystyle \mathrm{d}}}^2\overline{\mathsf{\rho }},$$

(1)

where ${{\displaystyle \mathrm{P}}}_{\mathrm{T}}$ is the power of the right HCP(RHCP) transmitted signals, $\mathsf{\lambda }$~ 19 cm is the wavelength of the signals at L1, ${{\displaystyle \mathrm{G}}}_{\mathrm{T}}$ and ${{\displaystyle \mathrm{G}}}_{\mathrm{R}}$ are the transmitting and receiving antenna gains, ${{\displaystyle \mathrm{R}}}_{\mathrm{T}}$ and ${{\displaystyle \mathrm{R}}}_{\mathrm{R}}$ are the ranges from the transmitter and the receiver to the specular point, respectively, $\mathsf{\chi }$ is the Woodward ambiguity function (WAF), $\mathsf{\tau }$ is the delay of the signal from the transmitter to the receiver, $\mathrm{f}$ is the Doppler shift of the reflected signal, ${\mathrm{f}}_{\mathrm{c}}$ is aimed to compensate the Doppler shift of the signal, ${\mathsf{\sigma }}^0$ is the bistatic scattering coefficient at LHCP, and $\overline{\mathsf{\rho }}$ is the positioning vector of the scattering point.

The bistatic scattering coefficient can be defined as follows [18]:

$${\mathsf{\sigma }}^0={\mathsf{\sigma }}^{\mathrm{c}\mathrm{o}\mathrm{h},0}+{\mathsf{\sigma }}^{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h},0},$$

(2)

where ${\mathsf{\sigma }}^{\mathrm{c}\mathrm{o}\mathrm{h},0}$ and ${\mathsf{\sigma }}^{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h},0}$ are the coherent and the incoherent scattering terms respectively. Consequently, DDMs consist on a sum of two terms as follows [19,20,21]:

$$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle =\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{c}\mathrm{o}\mathrm{h}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle +\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle .$$

(3)

$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{c}\mathrm{o}\mathrm{h}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ accounts for coherent reflections from the surface, while $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ is responsible for the diffuse scattering.

$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{c}\mathrm{o}\mathrm{h}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ can be expressed analytically as follows [10]:

$$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{c}\mathrm{o}\mathrm{h}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle =\frac{{{\displaystyle {{\displaystyle \mathrm{T}}}_{\mathrm{c}}^2\mathrm{P}}}_{\mathrm{T}}{{\displaystyle \mathsf{\lambda }}}^2{{\displaystyle \mathrm{G}}}_{\mathrm{T}}{{\displaystyle \mathrm{G}}}_{\mathrm{R}}{{\displaystyle \left|\mathsf{\chi }(\mathsf{\tau },\mathrm{f})\right|}}^2}{{{\displaystyle (4\mathsf{\pi })}}^2{({\mathrm{R}}_{\mathrm{T}}+{\mathrm{R}}_{\mathrm{R}})}^2}{\left|\mathrm{R}({\mathsf{\theta }}_{\mathrm{e}})\right|}^2\mathsf{\gamma }\exp (-{(2\mathrm{k}\mathsf{\sigma }\sin {\mathsf{\theta }}_{\mathrm{e}})}^2),$$

(4)

where $\mathrm{R}$ is the LHCP Fresnel reflection coefficient, $\mathrm{k}$ is the signal angular wavenumber, $\mathsf{\sigma }$ is the surface height standard deviation (related to small-scale surface roughness), and $\mathsf{\gamma }$ is the transmissivity of the vegetation.

$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ includes contributions from vegetation, small-scale surface roughness, and topography. Further details on the incoherent scattering over rough surfaces can be found in Reference [22]. The impact of vegetation should be also considered to establish a model of the incoherent scattered field. A comprehensive method can be found in References [23,24].

Over regions without rough topography, DDMs can be expressed as follows [25]:

$$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle ={\left|\mathsf{\chi }(\mathsf{\tau },\mathrm{f})\right|}^2\ast \ast {\left|{\mathsf{\sigma }}^0(\mathsf{\tau },\mathrm{f})\right|}^2,$$

(5)

where $\ast \ast$ is the 2-D convolution in both domains, delay and Doppler. The WAF can be defined as follows [17]:

$$\left|\mathsf{\chi }(\mathsf{\tau },\mathrm{f})\right|=\mathsf{\Lambda }(\mathsf{\tau })\mathrm{S}(\mathrm{f}),$$

(6)

where $\mathsf{\Lambda }$ is the autocorrelation function of the pseudo-random noise (PRN) codes and $\mathrm{S}$ is the system impulse response in the frequency domain.

Regions with rough topography are characterized with a root mean square error (RMSE) of the surface height variation higher than the length of the GNSS code under study. The effects of rough topography in the DDMs can be understood as follows:

$$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{o}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle =\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle \ast \ast \xi (\mathsf{\tau },\mathrm{f}),$$

(7)

where the function $\xi$ represents the mean density of scattering points as a function of the delay, and the Doppler shift.

3. Data and Methodology

3.1. Topographic Parameters

Several topographic parameters were derived from the ~250 m GMTED2010 model based on a 3 × 3 moving window (Figure 1). This window was composed of the focal cell and its eight surrounding cells. The DEM surface was approximated by a bivariate quadratic equation [26,27]:

$$\mathrm{Z}=\mathrm{a}{\mathrm{X}}^2+\mathrm{b}{\mathrm{Y}}^2+\mathrm{c}\mathrm{X}\mathrm{Y}+\mathrm{d}\mathrm{X}+\mathrm{e}\mathrm{Y}+\mathrm{f},$$

(8)

where $\mathrm{Z}$ represents the height of the DEM surface and $\mathrm{X}$ and $\mathrm{Y}$ are the horizontal coordinates. The coefficients in Equation (8) were solved within the 3 × 3 window using simple combinations of neighboring cells, based on the so-called Wood’s approach [28].

The topographic parameters were used to investigate the impact of different topographic features on $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{o}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$. In particular, three heterogeneity variables (topographic position index or $\mathrm{T}\mathrm{P}\mathrm{I}$, terrain ruggedness index or $\mathrm{T}\mathrm{R}\mathrm{I}$, vector ruggedness measure or $\mathrm{V}\mathrm{R}\mathrm{M}$), slope ($\mathsf{\beta }$), and two curvature variables (tangential ($\mathrm{T}\mathrm{C}$) and profile curvature ($\mathrm{P}\mathrm{C}$) were used in the parametric study. They have the potential to enable a better understanding on the properties of the Earth’s surface than raw DEM data. Topographic heterogeneity is described as the variability of surface elevations within an area. On the other hand, curvature attributes are based on the change of slope in a particular direction. More specifically, the selected topographic parameters (Figure 2) were defined as follows [26,27,28]:

• $\mathrm{T}\mathrm{P}\mathrm{I}$ is the difference between the elevation of a focal cell and the mean of its eight surrounding cells. Positive and negative values correspond to ridges and valleys, respectively. Zero values correspond to flat areas. It was computed as follows:

  $$\mathrm{T}\mathrm{P}\mathrm{I}={\mathrm{Z}}_{(0,0)}-{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{Z}}_{\mathrm{i}}/\mathrm{m}},$$

  (9)

  where $\mathrm{m}$ is total number of surrounding points employed in the evaluation.
• $\mathrm{T}\mathrm{R}\mathrm{I}$ is the mean of the absolute differences in elevation between a focal cell and its eight surrounding cells. It quantifies the total elevation change across the 3 × 3 moving window. Flat areas have a value of zero and mountain areas with steep ridges have positive values. It was computed as follows:

  $$\mathrm{T}\mathrm{R}\mathrm{I}=\frac{{\displaystyle {\sum }_{\mathrm{i}=-\mathrm{N}}^{\mathrm{N}}{\displaystyle {\sum }_{\mathrm{j}=-\mathrm{N}}^{\mathrm{N}}\left|{\mathrm{Z}}_{(\mathrm{i},\mathrm{j})}-{\mathrm{Z}}_{(0,0)}\right|}}}{({\mathrm{n}}^2-1)},$$

  (10)

  where $\mathrm{N}=(\mathrm{n}-1)/2$ and $\mathrm{n}$ may be any odd integer smaller than the number of cells in the shortest side of the raster.
• $\mathrm{V}\mathrm{R}\mathrm{M}$ quantifies the terrain ruggedness based on the measurements of the dispersion of vectors orthogonal to the surface. VRM quantifies the local variation of slope in the terrain more independently than TPI and TRI. VRM values range from 0 over flat areas to 1 over rough areas. It was computed using vector analysis [26].
• $\mathsf{\beta }$ is defined as the rate of change of elevation in magnitude for the steepest descent vector. It influences hydraulic gradients driving any surface flows and also sub-surface flows when the water table has a similar slope to the ground surface. Slope and curvature serve as useful input variables in erosion and hydrological models. It was computed as follows:

  $$\mathsf{\beta }=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}(\sqrt{{\mathrm{d}}^2+{\mathrm{e}}^2}).$$

  (11)
• $\mathrm{T}\mathrm{C}$ measures the rate of change perpendicular to the slope gradient. It is related to the convergence and divergence of lateral flows across a surface. It was computed as follows:

  $$\mathrm{T}\mathrm{C}=\frac{200(\mathrm{b}{\mathrm{d}}^2+\mathrm{a}{\mathrm{e}}^2-\mathrm{c}\mathrm{d}\mathrm{e})}{{({\mathrm{e}}^2+{\mathrm{d}}^2)}^{1.5}}.$$

  (12)
• $\mathrm{P}\mathrm{C}$ measures the rate of change of slope along a flow line. It affects the acceleration or deceleration of surface flows along the surface and thus it influences erosion and deposition of soils. Positive and negative values indicate convex and concave surfaces, respectively. It was computed as follows:

  $$\mathrm{P}\mathrm{C}=\frac{-200(\mathrm{a}{\mathrm{d}}^2+\mathrm{b}{\mathrm{e}}^2+\mathrm{c}\mathrm{d}\mathrm{e})}{({\mathrm{e}}^2+{\mathrm{d}}^2){(1+{\mathrm{e}}^2+{\mathrm{d}}^2)}^{1.5}}.$$

  (13)

3.2. Topographic Wetness Index

The topographic wetness index ($\mathrm{T}\mathrm{W}\mathrm{I}$) is a parameter used to match runoff- producing elements in the landscape [29]. Different target areas with similar $\mathrm{T}\mathrm{W}\mathrm{I}$ values should have similar hydrological dynamics [30]. As such, $\mathrm{T}\mathrm{W}\mathrm{I}$ is an indicator of the relative propensity of soils to become saturated to the surface because of the local topography. In this work, $\mathrm{T}\mathrm{W}\mathrm{I}$ values were calculated based on the GA2 algorithm [31]. GA2 calculates the outflow gradient of each target area and uses precalculated up-land values from HydroSHEDS for the catchment area [32]. It can be expressed as follows [31]:

$$\mathrm{T}\mathrm{W}\mathrm{I}=\ln \left(\frac{\mathrm{A}}{\tan \mathsf{\beta }}\right),$$

(14)

where $\mathrm{A}$ is the specific catchment area [33,34]. $\mathrm{T}\mathrm{W}\mathrm{I}$ values are low at ridges and high in valleys. Generally, humid areas generate higher $\mathrm{T}\mathrm{W}\mathrm{I}$ values, although there are some exceptions such as desert areas where high $\mathrm{T}\mathrm{W}\mathrm{I}$ values do not correlate with high flow accumulation.

3.3. Delay Doppler Maps Parameters

CyGNSS Level 1 Science Data Record Version 2.1 [12,15,16] was selected for this study. Original data were first filtered out using an equivalent “CyGNSS overall quality flag” over land surfaces [35]. Reflected delay waveforms $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{r}\mathrm{a}\mathrm{w}}$ were obtained from the original DDMs $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{o}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ at zero Doppler frequency:

$$\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{r}\mathrm{a}\mathrm{w}}=\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{o}}(\mathsf{\tau },\mathrm{f}=0)\right|}}^2\rangle .$$

(15)

The delay bin resolution of the original 17 lag waveforms $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{r}\mathrm{a}\mathrm{w}}$ was ~ 0.2552 GPS C/A chips. A re-sampling and interpolation of the resulting 1700 lag waveforms was performed using a spline method to increase the accuracy of the waveforms, before applying the algorithms to extract the observables.

The trailing edge (TE) width was defined as the lag difference between the 70% power threshold of the high resolution waveforms $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{t}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}}$, and that corresponding to the maximum power of the waveforms $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}$ [9]:

$$\mathrm{T}\mathrm{E}={\mathsf{\tau }}_{\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{t}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}}}-{\mathsf{\tau }}_{\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}}.$$

(16)

Different power thresholds were tested. The 50% power threshold was found to cut off some lags, when this threshold was out of the original 17 lag waveforms. On the other hand, the 90% power threshold provided a lower dynamic range as compared to the 70%. The incoherent scattering term ${\mathsf{\sigma }}^{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h},0}$ was the main contribution to $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{t}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}}$, while both the coherent ${\mathsf{\sigma }}^{\mathrm{c}\mathrm{o}\mathrm{h},0}$ and the incoherent ${\mathsf{\sigma }}^{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h},0}$ scattering terms significantly contributed to the peak power of the waveforms $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}$.

The reflectivity $\mathsf{\Gamma }$ was estimated as the ratio of the reflected $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}$ and the direct $\mathrm{W}{\mathrm{F}}_{\mathrm{d},\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}$ power waveforms’ peaks, after compensating for the noise power floor and the antennas’ gain patterns as a function of the elevation angle ${\mathsf{\theta }}_{\mathrm{e}}$ as in [10]:

$$\mathsf{\Gamma }=\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}/\mathrm{W}{\mathrm{F}}_{\mathrm{d},\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}.$$

(17)

The antennas’ gain patterns were compensated versus the gain at the corresponding boresight direction, and the difference of both antennas at boresight [10]. The nominal mission lifetime was 2 years. This milestone was achieved in March 2019. At present, the mission is operating nominally, and the automatic gain control (AGC) was disabled from August 2018 [12]. This fact improves the accuracy in the estimation of $\mathsf{\Gamma }$ using Equation (17). On the other hand, $\mathsf{\Gamma }$ could also be inverted from scattering models. This approach relies on the assumption that the scattering is totally coherent [35] or totally incoherent [36]. The assumption in [35] is valid over smooth surface areas with little vegetation. The assumption in [36] is valid if inland water bodies are removed. However, in a more general scenario, the total scattered electromagnetic field is composed of both a coherent and an incoherent contribution [18] in different proportions depending on the dielectric and geometrical properties of the scattering medium, and the directions of incoming and outgoing electromagnetic waves. Both, the coherent ${\mathsf{\sigma }}^{\mathrm{c}\mathrm{o}\mathrm{h},0}$ and the incoherent ${\mathsf{\sigma }}^{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h},0}$ scattering terms contribute to the peak power of the waveforms. This was the main motivation to use Equation (17) in this work.

4. Understanding the Effects of Rough Topography on GNSS-R Space-Borne Data

4.1. Global Scale

The main goal of this work was to study the impact of rough topography on CyGNSS-derived $\mathrm{T}\mathrm{E}$ & $\mathsf{\Gamma }$ as a function of SMAP-derived SMC. As a first step, the selected topographic parameters (Figure 2), $\mathrm{T}\mathrm{W}\mathrm{I}$ (Figure 3a), CyGNSS-derived $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ (Figure 3b,c), and SMAP-derived SMC (Figure 3d) were displayed over the complete coverage of the Earth’s surface enabled by CyGNSS. SMC data were derived based on the V-pol single channel algorithm [37,38]. Corrections for roughness, effective surface temperature, and vegetation water content (VWC) were also applied. Six months of CyGNSS and SMAP data were selected from 1 August 2018 to 31 January 2019. This specific time period was selected because the AGC was disabled in July 2018, so as to improve the performance in the estimation of $\mathsf{\Gamma }$(Equation (17)). The selected temporal length was as large as possible at the time of starting this work. All the parameters were averaged using a 0.1° × 0.1° latitude/longitude grid with a moving window of 0.2° in steps of 0.1°. This gridding strategy enabled the analysis using auxiliary data and different sensors. The resulting spatial resolution was ~ 20 km × ~ 20 km at equatorial latitudes, which was wide enough to account for the across-track spreading ~ 20 km of the DDMs due to $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$.

Overall, it was found that the spreading of $\mathrm{T}\mathrm{E}$ was low ~ [100, 250] m and $\mathsf{\Gamma }$ is high ~ [−10, −3] dB over flat areas with low topographic heterogeneity and small surface elevation gradients (Figure 2 and Figure 3). Over these areas, the scattering was quite specular and the coherence of the Earth-reflected GPS signals could be assumed to be high [35]. Under these circumstances,$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{c}\mathrm{o}\mathrm{h}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle \gg \langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$. Thus, the power of the reflected signals could be assumed to be roughly independent of the platform’s height (Equation.4). This point justifies the strong values of $\mathsf{\Gamma }$ and the reduced spreading of the waveforms. SMC was high over regions that do not always correspond to flat areas (Figure 3d). On the other hand, the spatial patterns of $\mathrm{T}\mathrm{E}$ & $\mathsf{\Gamma }$ showed a certain agreement with $\mathrm{T}\mathrm{W}\mathrm{I}$. An exception is found over e.g., tropical rainforests. Here, the spreading of $\mathrm{T}\mathrm{E}$ was high and $\mathsf{\Gamma }$ was low because specular reflection is weakened by vegetation attenuation, $\mathsf{\gamma }$, and vegetation structure contributes to diffuse scattering, both of which significantly modify the shape of $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{o}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ [39]. A more quantitative analysis is performed at a regional scale over an area with little vegetation and low small-scale surface roughness. Over desert areas e.g., Sahara, Kalahari, and Australia it was found that $\mathsf{\Gamma }$ was also quite high, despite SMC being almost negligible. This aspect was previously attributed to the effect of sub-surface scattering over rich sand content areas [6,9]. In very dry conditions, GNSS-R observables have shown a certain sensitivity to the soil moisture at deeper levels [8].

In addition, specular scattering over areas with high $\mathrm{T}\mathrm{W}\mathrm{I}$ could significantly increase $\mathsf{\Gamma }$, despite SMC being low (Figure 3). Indeed, $\mathrm{T}\mathrm{E}$ was rather low over areas with high $\mathrm{T}\mathrm{W}\mathrm{I}$, which is an additional indication of the specular nature of the scattering. $\mathrm{T}\mathrm{W}\mathrm{I}$ accounted for saturation of SMC in local topographic converging areas and was used to predict zones of surface saturation and therefore organized spatial fields of SMC. It could be expected that inland water bodies, i.e., rivers, are predominant in regions with high $\mathrm{T}\mathrm{W}\mathrm{I}$, which could justify the specular nature of the scattering [29]. On the other hand, topographic heterogeneity ($\mathrm{T}\mathrm{P}\mathrm{I}$, $\mathrm{T}\mathrm{R}\mathrm{I}$, and $\mathrm{V}\mathrm{R}\mathrm{M}$), $\mathsf{\beta }$, and curvature parameters ($\mathrm{T}\mathrm{P}\mathrm{I}$, $\mathrm{T}\mathrm{R}\mathrm{I}$) are high over the Andes and Himalayan mountains, where $\mathrm{T}\mathrm{E}$ is high because of the effect of diffuse scattering $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$.

4.2. Regional Scale

4.2.1. Introduction

A specific target area over South Asia was selected to further investigate topographic effects in $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{o}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$, without the influence of high levels of $\mathsf{\sigma }$, above ground biomass (AGB), and canopy height (CH). The coordinates of the selected target area were as follows: latitude = [5, 36] ° and longitude = [64, 94] °. The characteristics of this region enabled this study to be performed over a rich variety of terrains, including the Tibetan Plateau, Great Himalaya, Indo-Gangetic Plains, and Hindustan. As such, the impact of different topographic features in $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ can be properly evaluated (Figure 4, Figure 5 and Figure 6). In so doing, $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{o}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ were classified in two groups to evaluate the effects of topography over two different scenarios. At ${\mathsf{\theta }}_{\mathrm{e}}$~ [20, 30] °, the main contribution to the generation of the reflected signal in a specular manner (elevation angle of the incident wave equals elevation angle of the scattered wave i.e., ${\mathsf{\theta }}_{\mathrm{e},\mathrm{i}}$= ${\mathsf{\theta }}_{\mathrm{e},\mathrm{s}}$= ${\mathsf{\theta }}_{\mathrm{e}}$) came from surface facets over high elevation mountains. Some low elevation areas could be occulted by the top of the mountains. On the other hand, at ${\mathsf{\theta }}_{\mathrm{e}}$~ [80, 90] °, both valleys and mountain peaks could reflect signals in a specular manner towards the CyGNSS. This study was additionally complemented with SMAP-derived SMC and its coefficient of variation (CV). A wide range of SMC levels could be observed. The highest SMC values were associated with both the Indo-Gangetic Plains and mountains regions in the West of Hindustan, while the lowest SMC levels were found over the Tibetan Plateau and the Thar Desert. The spatial patterns of SMC did not appear correlated with the topographic parameters. Significant SMC levels could be found over mountains but also over rivers basins.

SMAP SMC data were used as an auxiliary product to help in the interpretation of the experimental results. The SMAP baseline accuracy requirement of 4% is not guaranteed over areas with very rough topography. SMC data are here provided with a lower accuracy. Nonetheless, 6 months of SMC data were averaged using a 0.1 ° x 0.1 ° latitude/longitude grid with a moving window of 0.2 ° in steps of 0.1 °. This strategy mitigated the effects of potential artefacts in SMAP-derived SMC. Additionally, this moving filter smoothed potential residual errors in the GMTED2010. Some errors have been reported over areas above Latitude ~ 60 ° [26]. However, this study was focused on areas within CyGNSS coverage Latitude ~ [−40, 40] °.

4.2.2. Interpretation of the Empirical Results

Figure 4 shows that $\mathrm{T}\mathrm{E}$ increases and $\mathsf{\Gamma }$ decreases with increasing values of the heterogeneity variables i.e.,$\mathrm{T}\mathrm{P}\mathrm{I}$ (Figure 4a), $\mathrm{T}\mathrm{R}\mathrm{I}$ (Figure 4b), and $\mathrm{V}\mathrm{R}\mathrm{M}$ (Figure 4c). Additionally, the measurements of $\mathsf{\beta }$ (Figure 4d), $\mathrm{T}\mathrm{C}$ (Figure 4e), and $\mathrm{P}\mathrm{C}$ (Figure 4f) have a clear impact on both observables. The ranges of $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ are the same in all the plots. This strategy was assumed to provide intercomparable plots, and at the same time show full variability. Over the Himalayan mountains, the spreading of $\mathrm{T}\mathrm{E}$ was very high (up to ~ 600 m) and $\mathsf{\Gamma }$ was very low (down to ~ −25 dB). On the other hand, an inverse behavior appeared over the Indo-Gangetic Plains. Here, the scattered electromagnetic field became mostly coherent at grazing angles because the exponential factor (Equation (4)) increased significantly with decreasing angles [40]. As a consequence, there was an increment of $\mathsf{\Gamma }$, which was distributed in a rather homogeneous manner (Figure 4j), despite the strong SMC gradient from east to west (Figure 4i). In Reference [40] it was previously found that the $\mathsf{\Gamma }$ dynamic range decreased at grazing angles. In other words, the sensitivity to SMC reduced at this geometry, which helps to justify the homogeneous distribution of $\mathsf{\Gamma }$. On the other hand, $\mathsf{\Gamma }$ appeared more associated with the variability of SMC in a nadir-looking configuration, despite being lower than at grazing angles [8,40].

Figure 5 and Figure 6 show the relationship of $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ with the selected topographic parameters as a function of SMAP-derived SMC, over the two angular ranges used in this study ${\mathsf{\theta }}_{\mathrm{e}}$ ~ [20, 30] ° and ${\mathsf{\theta }}_{\mathrm{e}}$ ~ [80, 90] °: $\mathrm{T}\mathrm{P}\mathrm{I}$ (Figure 5a,d,g,j), $\mathrm{T}\mathrm{R}\mathrm{I}$ (Figure 5b,e,h,k), $\mathrm{V}\mathrm{R}\mathrm{M}$ (Figure 5c,f,i,l), $\mathsf{\beta }$ (Figure 6a,d,g,j), $\mathrm{T}\mathrm{C}$ (Figure 6b,e,h,k), and $\mathrm{P}\mathrm{C}$ (Figure 6c,f,i,l). The influence of SMC was evaluated at four differentiated ranges from ~0 m³/m³ to ~0.4 m³/m³ in steps of ~0.1 m³/m³.

For low SMC ~ [0, 0.2] m³/m³, the impact of the heterogeneity variables on both $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ was quite strong in the ranges $\mathrm{T}\mathrm{P}\mathrm{I}$~ [−0.01, 0.01] m, $\mathrm{T}\mathrm{R}\mathrm{I}$~ [0, 20] m, and $\mathrm{V}\mathrm{R}\mathrm{M}$~ [0, 0.005] °. This point could be just justified because for low SMC levels: (a) the geometrical properties of the scattering medium play a stronger role than the dielectric properties and (b) the GNSS-R observables become saturated for a smaller topographic roughness.

Increasing SMC levels reduced the gradients because of the waveforms’ peak power increment. Consequently, the 70% power threshold appeared at smaller $\mathrm{T}\mathrm{E}$ values as compared to dry soils. The variability of $\mathrm{T}\mathrm{E}$ with SMC was higher in a nadir-looking configuration than at grazing angles because of the lower topographic influence. On the other hand, $\mathsf{\Gamma }$ increased with increasing SMC in a nadir-looking configuration. However, it was roughly independent of SMC at grazing angles [8]. As such, the spreading of $\mathrm{T}\mathrm{E}$ and the decrease of $\mathsf{\Gamma }$ showed a much more gradual evolution for increasing heterogeneity levels up to $\mathrm{T}\mathrm{E}$ ~ 600 m and $\mathsf{\Gamma }$ ~ – 22 dB at ${\mathsf{\theta }}_{\mathrm{e}}$ ~ [80, 90] ° and $\mathrm{T}\mathrm{E}$ ~ 600 m and $\mathsf{\Gamma }$~ - 15 dB at ${\mathsf{\theta }}_{\mathrm{e}}$~ [20, 30] °.

Finally, it is pointed out that the topographic dynamic range of both GNSS-R parameters $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ was higher for SMC ~ [0, 0.2] m³/m³ than for SMC ~ [0.2, 0.4] m³/m³. This observation is also attributed to the stronger influence of the geometrical properties for decreasing SMC levels. A similar behavior was also pointed out in the case of the small-scale surface roughness [40].

The impact of surface slopes and curvature variables on both observables, $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$, was quite strong at low SMC ~ [0, 0.2] m³/m³, over the following ranges: $\mathsf{\beta }$~ [0, 5] °, $\mathrm{T}\mathrm{C}$~ [-1, 2] rad/m, and $\mathrm{P}\mathrm{C}$~ [−2, 1] rad/m. Increasing SMC levels reduced the gradients in agreement with the findings using the heterogeneity variables (Figure 5). The spreading of $\mathrm{T}\mathrm{E}$ and the decrease of $\mathsf{\Gamma }$ showed a much more gradual evolution for increasing curvature levels up to (a) ${\mathsf{\theta }}_{\mathrm{e}}$~ [20, 30] °: $\mathrm{T}\mathrm{E}$~ 600 m and $\mathsf{\Gamma }$ ~ – 13 dB for SMC ~ [0.3, 0.4] m³/m³, $\mathrm{T}\mathrm{E}$~ 600 m and $\mathsf{\Gamma }$ ~ – 13 dB for SMC ~ [0.2, 0.3] m³/m³, $\mathrm{T}\mathrm{E}$ ~ 400 m and $\mathsf{\Gamma }$~ – 13 dB for SMC ~ [0.1, 0.2] m³/m³, and $\mathrm{T}\mathrm{E}$~ 300 m and $\mathsf{\Gamma }$~ – 13 for SMC ~ [0, 0.1] m³/m³ and (b) ${\mathsf{\theta }}_{\mathrm{e}}$~ [80, 90] °: $\mathrm{T}\mathrm{E}$~ 600 m and $\mathsf{\Gamma }$~ – 16 dB for SMC ~ [0.3, 0.4] m³/m³, $\mathrm{T}\mathrm{E}$~ 600 m and $\mathsf{\Gamma }$~ – 18 dB for SMC ~ [0.2, 0.3] m³/m³, $\mathrm{T}\mathrm{E}$~ 375 m and $\mathsf{\Gamma }$~ – 22 dB for SMC ~ [0.1, 0.2] m³/m³, and $\mathrm{T}\mathrm{E}$~ 200 m and $\mathsf{\Gamma }$~ – 23 dB for SMC ~ [0, 0.1] m³/m³.

Surface slopes modify the scattering area with respect to that corresponding to a flat Earth assumption. They are a key parameter in GNSS-R studies over land surfaces, which have a certain degree of correlation with $\mathrm{T}\mathrm{R}\mathrm{I}$ and $\mathrm{V}\mathrm{R}\mathrm{M}$ because increasing topographic heterogeneity belongs to increasing surface slopes. On the other hand, the $\mathrm{P}\mathrm{C}$ and $\mathrm{T}\mathrm{C}$ have a much clearer differentiated impact on $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$. These topographic variables are directly related to water accumulation over the surface. Over flat surfaces, $\mathsf{\Gamma }$ is high and $\mathrm{T}\mathrm{E}$ is low because the curvature is negligible, and the scattering is quite specular. However, the impact of SMC is lower compared to areas with higher curvatures. In these areas, $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ have a clearer differentiated response to different SMC levels. Finally, it was found that the variability with SMC is higher in a nadir-looking configuration than at grazing angles. It increases for negative $\mathrm{P}\mathrm{C}$ levels (concave), while for positive curvature levels (convex) there is a negligible of influence of SMC.

4.2.3. Sensitivity of CyGNSS to TWI vs. SMC

The strong influence of primary terrain attributes on GNSS-R observables (Figure 5 and Figure 6) triggered the sensitivity study of $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ with $\mathrm{T}\mathrm{W}\mathrm{I}$ vs. SMC. To do so, a specific target area (Latitude ~ [20, 27] ° and Longitude ~ [72, 84] °) was selected with a strong gradient of SMC, and a wide topographic heterogeneity (Figure 7). Figure 8 shows the density scatter plots of $\mathsf{\Gamma }$ vs. $\mathrm{T}\mathrm{W}\mathrm{I}$ at ${\mathsf{\theta }}_{\mathrm{e}}$~ [20, 30] ° (Figure 8a) and ${\mathsf{\theta }}_{\mathrm{e}}$~ [80, 90] ° (Figure 8b), $\mathsf{\Gamma }$ vs. SMC at ${\mathsf{\theta }}_{\mathrm{e}}$~ [20, 30] ° (Figure 8c) and ${\mathsf{\theta }}_{\mathrm{e}}$~ [80, 90] ° (Figure 8d), $\mathrm{T}\mathrm{E}$ vs. $\mathrm{T}\mathrm{W}\mathrm{I}$ at ${\mathsf{\theta }}_{\mathrm{e}}$~ [20, 30] ° (Figure 8e) and ${\mathsf{\theta }}_{\mathrm{e}}$~ [80, 90] ° (Figure 8f), and $\mathrm{T}\mathrm{E}$ vs. SMC at ${\mathsf{\theta }}_{\mathrm{e}}$~ [20, 30] ° (Figure 8g) and ${\mathsf{\theta }}_{\mathrm{e}}$~ [80, 90] ° (Figure 8h).

Both observables $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ are more correlated with $\mathrm{T}\mathrm{W}\mathrm{I}$ than with SMC (

Table 1. Statistics (Pearson correlation coefficient $\mathrm{r}$, RMSE, and slope of the linear fits) of the scatter plots of CyGNSS-derived $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ with $\mathrm{T}\mathrm{W}\mathrm{I}$ and SMAP-derived SMC over the selected target area over North India.

	$\mathbf{Topographic} \mathbf{Wetness} \mathbf{Index} \left(\mathbf{T}\mathbf{W}\mathbf{I}\right)$				SMC [m3/m3]
	$\mathbf{T}\mathbf{E}$ [m]	$\mathsf{\Gamma }$ [dB]	$\mathbf{T}\mathbf{E}$ [m]	$\mathsf{\Gamma }$ [dB]	$\mathbf{T}\mathbf{E}$ [m]	$\mathsf{\Gamma }$ [dB]	$\mathbf{T}\mathbf{E}$ [m]	$\mathsf{\Gamma }$ [dB]
	[20, 30] °	[20, 30] °	[80, 90] °	[80, 90] °	[20, 30] °	[20, 30] °	[80, 90] °	[80, 90] °
$\mathrm{r}$	−0.63	0.59	−0.50	0.48	−0.08	0.27	−0.06	0.30
RMSE	95.3	1.7	68.5	2.7	121.9	2.1	80.7	2.8
Slope	−75.13	1.40	−18.65	1.57	−175.21	8.4078	−63.47	13.55

). This point could be justified because the potential presence of inland water bodies over areas with high $\mathrm{T}\mathrm{W}\mathrm{I}$ have a stronger impact on GNSS-R observables than SMC. In this analysis, areas with high elevation and very rough topography (i.e., Himalayan mountains) were excluded (Figure 7) because SMAP-derived SMC could have some uncertainties over areas with very rough topography. As such, SMC is delivered within the baseline accuracy in this target area.

The Pearson correlation coefficients increased at grazing angles because of the stronger impact of topography at this geometry as compared to a nadir-looking configuration. On the other hand, the correlation of $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ with SMC is almost negligible despite the moderate-to-high SMC dynamic range ~ [0, 0.4] m³/m³. The slopes of the linear regressions were also included in Table 1, for both angular ranges. It is found that the sensitivity of $\mathrm{T}\mathrm{E}$ to $\mathrm{T}\mathrm{W}\mathrm{I}$ and SMC increased at grazing angles, while that of $\mathsf{\Gamma }$ was stronger in a nadir-looking configuration in agreement with the results in Reference [8]. This empirical evidence suggests that topographic features can determine the main spatial pattern of Earth’s surface reflectivity, at least over some specific target areas.

Finally, it is worth pointing out that these results at a regional scale support qualitative results obtained at a global scale (Figure 1). As such, attention should be paid when analyzing GNSS-R data because soils saturated by the effect of topography could introduce confounding effects with respect to the impact of SMC.

5. Final Discussion

CyGNSS space-borne GNSS-R data appear to have a strong dependence on topography. This dependence was found to be different from that corresponding to small-scale surface roughness. The impact of topographic roughness increased with decreasing angles ${\mathsf{\theta }}_{\mathrm{e}}$, while the impact of surface roughness increased with increasing ${\mathsf{\theta }}_{\mathrm{e}}$. This latter point can be justified based on the Rayleigh criterion that establishes that a surface (without vegetation and rough topography) could be considered smooth if the phase difference between two reflected electromagnetic waves is lower than π/2 rad. On the other hand, it was found that topographic roughness was a much more complex phenomenon, which depended on several topographic features in a differentiated manner. Some Earth observation missions such as the Soil Moisture Ocean Salinity (SMOS) mission used a method to flag the pixels according to the relative impact of topography on the brightness temperature [41]. The GNSS-R scenario appears more complex because the bistatic scattering includes both a coherent ${\mathsf{\sigma }}^{\mathrm{c}\mathrm{o}\mathrm{h},0}$ and an incoherent ${\mathsf{\sigma }}^{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{h},0}$ component, in different proportions depending on the dielectric and geometrical properties of the scattering medium, and the directions of incoming and outgoing electromagnetic waves. On the other hand, the spatial resolution depends on the coherent-to-incoherent scattering ratio. This point complicates the topographic masking because of the impact of ${\mathsf{\theta }}_{\mathrm{e}}$.

The scattering over surfaces with rough topography could be assumed to be totally incoherent. However, experimental results showed that $\mathsf{\Gamma }$ increased with decreasing ${\mathsf{\theta }}_{\mathrm{e}}$ even over regions with very rough topography. This angular behavior is classically associated with the coherent scattering term ${\mathsf{\sigma }}^{\mathrm{c}\mathrm{o}\mathrm{h},0}$ because the effective surface roughness is much lower at grazing angles [40]. At this point, it could be hypothesized that the scattering could be locally coherent over specific mountains facets. The total scattered electromagnetic field could be coherent if only a single mountain facet contributes in a specular manner to the reflected signal as collected by CyGNSS.

6. Conclusions

Several topographic parameters, $\mathrm{T}\mathrm{P}\mathrm{I}$, $\mathrm{T}\mathrm{R}\mathrm{I}$, $\mathrm{V}\mathrm{R}\mathrm{M}$, $\mathsf{\beta }$, $\mathrm{T}\mathrm{C}$, and $\mathrm{P}\mathrm{C}$ were derived from the ~250 m GMTED2010 DEM. The impact of these parameters on different GNSS-R observables $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ was assessed over land surfaces as a function of SMAP-derived SMC and ${\mathsf{\theta }}_{\mathrm{e}}$. A parametric study was then performed over a target area in South Asia. For low SMC levels, $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ presented a strong fluctuation from flat surfaces to areas with moderate topographic roughness. On the other hand, increasing SMC levels reduced the gradients. Additionally, it was found that the sensitivity of both observables to SMC increased in a nadir-looking configuration as compared to grazing angles. As a final remark, $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ converged with increasing levels of $\mathrm{T}\mathrm{P}\mathrm{I}$, $\mathrm{T}\mathrm{R}\mathrm{I}$, $\mathrm{V}\mathrm{R}\mathrm{M}$, and $\mathsf{\beta }$; while they diverge with increasing levels of $\mathrm{T}\mathrm{C}$ and $\mathrm{P}\mathrm{C}$ for different SMC values.

In the last part of this study, the sensitivity of $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ to $\mathrm{T}\mathrm{W}\mathrm{I}$ and SMC was evaluated over a specific target region in North India, with a strong SMC gradient and a wide topography heterogeneity. Areas with very rough topography were excluded. The Pearson correlation coefficients of $\mathrm{T}\mathrm{E}$ and $\mathsf{\Gamma }$ with $\mathrm{T}\mathrm{W}\mathrm{I}$ were rather high. On the other hand, the correlation with SMC was almost negligible. This empirical evidence suggests that topographic features could determine the main spatial patterns of $\mathsf{\Gamma }$, at least over some areas. $\mathrm{T}\mathrm{W}\mathrm{I}$ provided an indication of the relative propensity of soils to become saturated to the surface due to the effect of the local topography.

Present and future applications of GNSS-R for Earth Observation over land surfaces could benefit from the results of this work. Accuracy on geophysical parameters retrieval e.g., SMC, could be improved after properly calibrated reflected DDMs. The estimation of SMC based on GNSS-R observables relies on their ability to represent the dielectric constant of the soil. However, topography significantly distorts these observables. The empirical relationships that were found between the selected topographic parameters and the observables pave the way to compensate underestimated SMC levels due to topographic roughness. Future work should include ground truth data to further validate these results.