Retrieval of Subsurface Soil Moisture and Vegetation Water Content From Multifrequency SoOp Reflectometry: Sensitivity Analysis

Signals of opportunity reflectometry (SoOp-R), the reutilization of noncooperative satellite transmissions for communication and navigation, is a promising approach to remote sensing of root-zone soil moisture (RZSM). Satellite transmissions in the frequency ranges of 137–138, 240–270, and 360–380 MHz are of interest due to the increased penetration depth. These can be combined with global navigation satellite system reflectometry (GNSS-R) in L-band (1575.42 MHz) to estimate the subsurface SM profile. The objective is to define requirements (e.g., frequency and polarization combinations, observation error, and temporal coincidence of multisource observations) for satellite-based remote sensing of RZSM. Our approach is to use synthetic observations generated from multiyear time series of in situ SM measurements from seven U.S. climate reference network (USCRN) sites and dynamic vegetation structure based on a simple scaling method. A multifrequency/polarimetric retrieval algorithm is developed and applied to these synthetic observations and used to predict retrieval errors for a range of changes in system parameters. We found that the use of both high and low frequencies improves retrieval accuracy by limiting uncertainties from vegetation and surface SM and providing sensitivity to deeper layers. Moreover, the retrieval errors were found to increase linearly with the reflectivity error and inter-frequency time delays. A bivariate model derived from this linear relationship will be useful for developing requirements on reflectivity precision based upon science requirements for SM/vegetation water content (VWC) retrievals. Although orbits of specific transmitter constellations were used to generate realistic distributions of incidence angle combinations, the method and results could be applied more generally.


I. INTRODUCTION
S OIL moisture (SM) is a key variable in the processes of evapotranspiration, infiltration, and runoff.Observing and modeling these processes are critical to understanding the water and energy cycle and exchanges between the land surface and atmosphere.Root-zone soil moisture (RZSM), in Manuscript submitted February 16, 2023.This work was sponsored in part by NASA under Grant 80NSSC21K1629, Grant 80NSSC18K0245, and Grant 80NSSC18K1524 and in part by Purdue University under Ross Fellowship F.00038215.02.046 (Corresponding author: Seho Kim).
Seho Kim and James L. particular, is a critical variable in plant-soil-atmosphere interactions and in monitoring and forecasting agricultural drought as it is closely related to plant growth.As the name implies, RZSM is defined as water content enclosed by the roots of vegetation within the top meter of the soil.This is based on an assumption that this depth range generally includes the highest density of plant roots.Despite its importance, RZSM is not globally measured using satellite remote sensing.
Several active and passive microwave sensors have monitored global surface SM over the last four decades [1].Two L-band passive microwave sensors, SMAP [2] and SMOS [3] are reliably observing surface SM in the top 5 cm soil, but cannot provide direct measurement of RZSM due to limited penetration depth of the signal at L-band.Although the assimilation of surface SM can estimate RZSM [4], [5], a disconnect between surface and subsurface SM dynamics often leads to large uncertainties in a hydrological model.
AirMOSS, an airborne P-band Synthetic Aperture Radar (SAR), was able to obtain RZSM observations with a retrieval error of 0.05 m 3 /m 3 for the top 25 cm of soil after bias removal [6].Extending P-band SAR to a space mission, however, faces technical challenges such as competitive spectrum allocation, large antenna size, and substantial radio frequency interference (RFI) [7].Signals of Opportunity Reflectometry (SoOp-R) has emerged in recent years as an alternative approach to observe RZSM by re-utilizing existing transmissions from communication or navigation satellites in a multistatic radar configuration [8], [9].For example, Yueh et al. proposed a mission concept based on the P-band SoOp-R to remotely sense RZSM and snow water equivalent [10].SoOp-R enables observations outside of the frequency allocations protected for science, using powerful anthropogenic sources that would be considered RFI in traditional active or passive systems.SoOp-R uses process gain to increase the signal to noise ratio.For coherent reflections, resolution is determined by the Fresnel zones, not the antenna beamwidth.Large antennas are therefore not required, even for observations at meter-scale wavelengths.
Recent experiments and simulation studies have shown the potential of multi-frequency SoOp-R to measure RZSM [10], [16]- [19].A single frequency SoOp-R observation will be sensitive to the entire SM profile within the sensing depth [18].It is therefore necessary to combine observations at multiple frequencies and solve an inverse problem to estimate this full profile and obtain SM retrievals at different depths.Observations at different frequencies are not available at the same time or with the same incidence angle, since the transmitters are in different orbits.It is expected that the accuracy of a multi-frequency retrieval would depend somewhat on the time delays between the different observations being combined.Designed for communications, all of the source signals are circularly polarized.Scattering from the surface would change the polarization, so additional information could be obtained from polarimetric observations.The vegetation canopy over a soil surface is one of the most important factors that obscure the SM observation, so the effect of vegetation must be included in any SM profile retrieval.Lastly, the accuracy of SM profile retrievals would depend upon the accuracy of the reflectivity observations made in each individual frequency.
The objective of this study is to determine the sensitivity of multi-frequency retrievals to the various parameters listed above.This will help advance the SoOp-R technique to a future Earth science mission as well as provide insights on potential improvement of the inversion algorithm.To achieve this objective, we developed a practical multi-frequency retrieval algorithm that integrates a hydrological SM profile model into a forward model and inversion scheme.This algorithm is evaluated using synthetic reflectivities generated from multiyear in-situ SM data and satellite-based NDVI data under different measurement configurations in order to evaluate the retrieval accuracy and define requirements for spaceborne remote sensing of RZSM using SoOp-R.
The remainder of this article is organized as follows.Section II describes our approach to the truth modeling, the retrieval algorithm, and accuracy assessment.SoOp-R system parameters and ancillary parameters used in our simulations are detailed in Section III.Section IV illustrates the simulation environment setup, procedure, and simulation results.Section V discusses significant findings and some limitations, followed by concluding remarks in Section VI.

II. APPROACH
A. Soil Moisture Data Source Evaluation of retrieval algorithms and their sensitivity to various parameters requires a long period of SoOp-R observations at the frequencies of interest, collected under a range of different conditions, along with in-situ SM profiles to be used as "truth".At the present time, such data do not exist.We will, therefore, utilize synthetic observations.It is important that such observations accurately represent the expected SM profiles and their evolution in time (to study the effect of noncoincident observations).We determined that the best approach to assure a distribution of observations with realistic properties was to generate the simulated reflectivity observations from actual in-situ SM measurements.We obtained data from seven (7) stations in the United States Climate Reference Network (USCRN) [20] operated by the National Oceanic and Atmospheric Administration (NOAA).These stations, listed in the ascending order of clay fractions in Table I, were selected based on their sensing depths, SM variability, data continuity, and containing vegetation classified as grasslands.The vegetation classification is based on International Geosphere-Biosphere Program (IGBP) obtained from Moderate Resolution Imaging Spectroradiometer (MODIS) land cover product (MCD12Q1) [21].Each station is equipped with a total of 15 probes (Hydra Probe II Soil Sensor Model SDI-12 of Stevens Water Monitoring Systems, Inc) installed in three plots at five depths of 5, 10, 20, 50, and 100 cm.SM and soil temperature are averaged from three independent measurements for each depth.Hourly in-situ SM measurements spanning 5 years from 2016 to 2020 from the 7 USCRN stations were used to create a sample set of SM profiles for input to the reflectivity simulation.Samples with air temperature or soil temperature of any layer below 0 °C were filtered out to eliminate frozen conditions.If a sample of any layer was missing, all samples at that time from that station were discarded.Soil fractions taken from Atmospheric Turbulence and Diffusion Division (ATDD-NOAA) ftp site [22] were used to classify soil texture for each station.

B. Truth Model for Synthetic Reflectivity Observation
A truth model is used to generate the time series of synthetic reflectivities from the in-situ SM profiles covered with vegetation.
1) Soil Moisture Profile Modeling: The electromagnetic model uses a volumetric SM profile, Θ v (z), at very small discrete steps in depth, z (1 cm in our application).Insitu data, however, is only available at the five (5) discrete depths.Various approaches, which can be categorized as traditional methods [23], have been used to generate a SM profile including regression, inversion, knowledge-based, and combinations of remotely sensed data and water balance models.Probabilistic methods include maximum entropy [24], [25], exponential filter [26], artificial neural network [27], [28], and wavelet transform [29], [30].The choice of a method is usually based on the amount of available a priori information and application scale.In the remote sensing community, the exponential filter-based method has received a lot of attention in recent years since it only requires surface SM data as input [26], [31], [32].However, it often fails to estimate RZSM properly under dry conditions due to its inherent assumption of no transpiration and constant hydraulic connectivity between surface and root-zone [33].
Another statistical model becoming popular in the hydrologic community is the Principle of Maximum Entropy (POME) model [24], [34], [35].A recent study has shown that the POME model distributes moisture content through a soil column more accurately than the exponential filter method [33].In principle, the maximization of entropy characterizes the diffusion of moisture movement through a soil column and guarantees the minimum variance unbiased profile given boundary conditions and mean moisture content.The POME uses (in "effective" values) bottom-most SM, Θ L and mean SM, Θ M , in addition to the surface SM, Θ 0 , to provide a better connection to the subsurface SM dynamics.The mean SM ultimately relates other hydrologic variables, such as evapotranspiration, runoff, and precipitation, as well as soil parameters in a field environment, such as soil porosity, to the SM distribution.
The POME model requires two constraints: (a) the total probability constraint, is a probability density function of the effective moisture content Θ.The second constraint relates the first moment in probability space to the mean SM of the soil column in physical space.Maximizing the Shannon entropy [36] given the constraints, a monotonic SM profile can be modeled by the main POME equation [24], in which Θ(z) is the effective moisture content at layer depth z, L is the total soil column depth which can be parameterized or act as a calibration factor, and λ 1 , λ 2 are Lagrange multipliers.The effective SM is defined by Θ where Θ v is volumetric SM as defined earlier, Θ r is residual SM, and n is soil porosity.The hydrologic soil properties corresponding to the soil texture of each USCRN site are taken from [37] to map effective SM (Θ) to volumetric SM (Θ v ) and vice versa.The Lagrange multipliers can be solved from application of the constraints and the boundary conditions.The model can represent any one of three distinct SM profiles: (1) dry case (+ sign in (1)), corresponding to a long time after a wetting event where the upper layer is drier than the lower, (2) wet case (− sign in (1)), during or soon after a wetting event where the upper layer is wetter than the lower, and (3) dynamic case (a combination of dry and wet cases), corresponding to a short time after a wetting event where one prominent inflection point exists.A complete model is derived in [24].
We apply the POME method to the time series of in-situ observations at the five discrete depths to model the vertical distribution of RZSM.Boundary conditions (Θ 0 , Θ L ) were set by the in-situ data at surface and bottom-most layers.Θ M was calculated from the SM data of all available depths in the soil column through numerical integration using the trapezoidal rule.We assume that the total soil column depth L is 100 cm and the inflection point is at 20 cm depth where it typically presents as suggested by [34].If an inflection point exists, Equation ( 1) is applied twice to model two profiles (upper and lower profiles with respect to the inflection point) while keeping the mass balance constraint.The lower boundary of the top profile and the upper boundary of the bottom profile will match.Following an idea in [33], the initial SM value at the inflection point is set as field capacity and decreased in the iterative operations until the mass balance constraint is satisfied.
The POME model is used both to generate synthetic reflectivities, and is included in our inversion algorithm for estimating the SM profile from the simulated SoOp-R measurements.
2) Dynamic Vegetation Structure Modeling: We represent vegetation as a layer of discrete scatterers located above the ground using a geometric modeling approach, in which the relative distribution and orientation of each component of the plant canopy, such as leaves and stalks, are modeled by their geometrical shapes.The advantage of this approach is that the attenuation and scattering properties of the vegetation canopy are expressed in terms of the biophysical quantities of individual plants including shape and orientation relative to the incident wave's electric field.Modeling seasonal variations of vegetation with this approach, however, requires detailed information about individual plant elements that change over time, which is difficult to obtain at a global scale.Thus, we approximate the geometrical and physical parameters of a vegetation component using the following linear relationship with the vegetation water content (VWC) data (W , in kg/m 2 ) [38]: where ρ is the number density of the scatterers (in #/m 3 ), V and m v are the volume (in m 3 ) and volumetric moisture content (in kg/m 3 ) of a single scatterer, respectively, and d is the canopy height in meters.A vegetation structure vector, a, contains the parameters a i which model the temporallyvarying geometrical and physical structure of vegetation.
Given a reference set of in-situ measurements of vegetation properties, a 0 , and a time series of VWC data, W (t), the element a i of the vector a at time t can be computed by The quantity a i,0 is the i-th element of the nominal parameter vector a 0 and W 0 is its corresponding VWC calculated by Equation (2).The quantity s i ∈ s is a scale parameter that can be selected empirically, taking into account the appropriate range of corresponding vegetation parameter values.This scaling method is intended to model a simple time-varying vegetation structure rather than accurate vegetative growth.
The vegetation parameters and corresponding scale parameters used in this study will be presented later in Section IV-A.
The VWC data used in this study is computed by using normalized difference vegetation index (NDVI) through the empirical derivation of the SMAP mission [39]; where NDVI max is the annual maximum NDVI and NDVI min is the annual minimum NDVI.SMAP's VWC data report suggests that the current NDVI substitutes for NDVI max in the case of grasslands and croplands and a global constant value of 0.1 is used for NDVI min for all vegetation types [39].The stem factor is a coefficient to estimate the stem water content whose values for different land cover is available from [39].We used satellite-based 16-day NDVI data (MOD13Q1) [40] at 250-m resolution from MODIS to produce daily time series of NDVI for each of the seven USCRN stations.The NDVI values were extracted for each location, filtered out by the quality flags, and linearly interpolated into daily NDVI from 2016 to 2020.
3) Forward Scattering Modeling: The volumetric SM profile produced from the in-situ SM data through the POME model is sampled at 1-cm depth increments.Given the frequency of the signal, the SM at each layer is converted into a complex dielectric constant using the Mironov model [41] to construct a multilayer dielectric structure.The soil clay content of each USCRN site listed in Table I is used as input to the Mironov model, assuming a homogeneous soil texture composition in all layers.Dielectric constants for vegetation components are also calculated.We utilize Ulaby and El-Rayes's model [42], one of the most widely used vegetation dielectric models applicable to all vegetation.Other input parameters of this model are volumetric moisture, which can be provided via Equation (3), and salinity, which is assumed to be fixed at 8.5 psu.The dielectric constant of air is unity and the air-vegetation interface is diffuse.
We employ the approach used in the Signals of opportunity Coherent Bistatic scattering model for Vegetated terrains (SCoBi-Veg) [43] to simulate polarimetric reflections on the dielectric structure through application of Maxwell's equations.The SCoBi-Veg model is a comprehensive electromagnetic forward scattering model that calculates power and complex field quantities of the direct, specular (coherent), and diffuse (incoherent) components of the signal received from the transmitter.It is capable of modeling the effects of a scattering environment including geometry, multilayer ground reflection, vegetation structures, and surface roughness, and of the transmitter and receiver properties such as antenna pattern, polarization mismatch, and frequencies through analytical wave theory and distorted Born approximation.
In general, the surface roughness relative to the wavelength of incident radiation determines the dominant scattering among those that are coherent and incoherent.Based on previous research findings [9], [44], it is reasonable to assume, for the purpose of this study, a dominantly coherent reflection.SCoBi-Veg models the specular reflectivity (coherent component), Γ s , as the matrix equation [18], [43], where î− s and ô+ s are the unit vectors indicating the incoming and outgoing directions of the wave propagation, G r and G t are the antenna normalized (power) pattern matrices for the receiving and transmitting antennas, respectively, U s→r and U t→s are the polarization basis rotation matrices between specular-point-to-receiver and transmitter-to-specularpoint, respectively, and E t is the coherency vector of nominal polarization state of the transmitting antenna.The term R s describes the propagation and reflection processes within the vegetation and subsurface soil.The specular reflectivity Γ s includes both co-and cross-polarized reflectivities in its vector form and can be obtained from the effective reflectivity (i.e.ratio of direct and reflected signal power) measured at the receiving antenna considering the observation geometry and measurement errors.
4) Simulated Observables: Equation ( 5) models the coherent reflectivity of soil+vegetation.The co-and cross-pol elements represent this reflectivity, Γ veg , as a function of the SM profile Θ v , a vegetation structure a, and a vector of known parameters, β.
, where f is the frequency of the transmitted signal, θ is the incidence angle of the signal at the specular point, and p and q are the polarization of the transmitter antenna and the receiver antenna, respectively.Given a "true" SM profile and VWC, Γ veg is simulated for one observation at a specific received frequency and polarization.This simulated reflectivity is then perturbed by adding measurement error, ν, which is assumed to be a zero-mean Gaussian random process.The standard deviation, σ, of ν is set to be a fraction of the reflectivity, σ = ηΓ veg .This will allow a better intercomparison of the effect of errors on different observations, given that the magnitude of reflectivity is dependent on the frequency and polarization.We use tilde ( ) to indicate the calculation+noise (simulated measurement) and hat ( ) to the estimation introduced in Section II-C2.A "true" value (the calculation without noise) is written without neither tilde nor hat.

C. Retrieval Algorithm
A practical retrieval algorithm was developed to estimate a SM profile and VWC given the synthetic reflectivity observations.
1) Forward Model: A forward model is used within the retrieval algorithm to demonstrate the retrieval of an estimated SM profile and VWC from the simulated noisy observations.Since the canopy structure of the observed scene is generally unknown, it is impractical to implement the discrete scatterer approach, used in the truth model to generate the synthetic observations, into the retrieval.Instead, we use the vegetation optical depth (VOD) to model the vegetation effect on the reflectivity.The decrease in reflectivity due to the soil layer is approximately an exponential function of the vertical (i.e., θ = 0) VOD at receiver polarization state q and frequency f , τ f,q [45]; where Γ soil = g(Θ v , 0; β) is the coherent reflectivity of bare soil estimated from the forward scattering model including the effect of surface roughness.The factor of 2 in the exponential in ( 8) accounts for the two-way attenuation along the incident and reflected paths through the vegetation canopy.
VOD can be linearly related to VWC for most crops and low vegetation, τ f,q = b f,q W . Parameter b f,q is dependent on the vegetation structure, microwave frequency, polarization, and incidence angle [45]- [47].This allows substitution of b f,q W into Equation ( 8) to include a single variable, W (VWC), to represent the net effect of the vegetation.The SM profile Θ v , generated by the POME, is determined by three variables defining the surface, bottom-most, and mean SM in effective values.A practical forward model, as a function of four (4) variables, x = [Θ 0 , Θ L , Θ M , W ], is therefore To our knowledge, numerical values for b f,q for grasslands with the I-band and P-band frequencies of interest have not been presented in the literature.We obtained these numerical values by running the SCoBi-Veg model for 10% of the data set that were randomly selected and then ran the model again using bare soil, and otherwise identical conditions.The VOD was then calculated by rearranging Equation ( 8), b f,q was then found from a least squares fit of Equation (10) to the the simulated τ f,q and W . Table II presents numerical values for b f,q obtained through this approach.b f,q at Lband was found to be somewhat larger than that used by the CYGNSS mission for grasslands [47], but is still within the range of experimental data [45].Generally, these are in reasonable agreement with the expectation that b f,q will be linearly related to the inverse of wavelength [45], [46].Note that the value for b f,q cannot be directly compared to that used in SMAP because it was defined for a different set of polarizations.2) Inversion Method: A number of geophysical parameters, including vegetation canopy, surface roughness, and soil layers, have complicated the SM retrieval of microwave sensors and often led to an ill-posed problem due to the lack of available observations to resolve the unknowns.Several inversion methods such as sequential assimilation of surface SM [48], [49], global and local optimization techniques [50], [51], and artificial neural networks [52], [53] have been utilized in the retrieval of SM profiles to overcome the nonlinearity and modeling complexity.In this study, we use a global optimization scheme known as adaptive simulated annealing (ASA) [54] for a generic multi-frequency SoOp-R system to directly retrieve a SM profile (i.e., surface, mean, and bottommost SM values) and VWC.ASA is a modified version of the classical simulated annealing method that statistically ensures faster convergence to a global optimal solution by adopting a new exponential annealing schedule with new distributions for the state and acceptance functions.The term adaptive comes from the mechanism where the annealing temperature and its rate of change is adjusted independently for each state parameter throughout the convergence process.It is particularly useful for applications in which a cost function has different sensitivities along different dimensions of the state space.Detailed descriptions can be found in [54]- [56].
The inverse problem is to determine the best estimates of the state vector x, defining the SM profile, Θ v , and the VWC, W , given a vector of N obs observations, T , consisting of reflectivities measured at specific frequencies and polarizations, incorporating independent measurement errors, ν i , as in Equation (7).
To solve the inverse problem, we pose it as an optimization problem to minimize the following least squares cost function using ASA.Each difference (computed minus observed) in the numerator is normalized by the observation in the denominator to balance the contribution of each observation to the cost function value, given that reflectivities at different frequencies, polarizations, and observation geometries possess different magnitudes and sensitivities.The adaptive feature of ASA takes into account the different sensitivities of the response to each state parameter, making the optimization process more efficient and less time-consuming.

D. Retrieval Accuracy Assessment Methodology
We evaluated the retrieved SM profiles against the POMEbased SM profiles used to generate the simulated noisy observations.Performance is evaluated in terms of the root-meansquare error (RMSE) of SM retrievals at a given depth z, where N SM is the total number of retrieved SM values at a depth z (from all locations and times under consideration) The RMSE of the entire SM profile is then given by where l is the threshold depth, which is defined as the depth up to which the SM values have been included in the calculation of retrieval error and N z is the number of depths within between the surface and l.We adopted the concept of the threshold depth used in the AirMOSS mission [6] to account for the fact that the RMSE generally increases when deeper layers are included in the RMSE calculation.Similarly, the RMSE of VWC is also evaluated against the true VWC value computed from NDVI using Equation (4).For efficient comparison of retrieval accuracy determined by various SoOp-R parameters, we set a baseline RMSE such that the RMSE of SM profiles up to 50 cm deep is less than 0.06 m 3 /m 3 , which we have assumed to be the maximum acceptable error for useful science measurements.Although the baseline RMSE has been chosen arbitrarily by the author as a working requirement, it is more desirable that the accuracy requirements are determined by applications utilizing SM observations of a specific depth.Such requirements for RZSM observations can be derived by the observation system simulation experiment (OSSE), which is beyond the scope of this article.

E. Bivariate Model for Retrieval Error
In our sensitivity analysis, we linearly approximated the relationship between the retrieval error (RMSE) and each of the two error sources, i.e. measurement error and observation timing gap.Straight lines where x is either the measurement error parameter (η) or observation time window (T ), which will be introduced in Section III-B, were fitted to the RMSEs of the top 50 cm SM profiles and VWC.Note that b SM,η is theoretically equal to b SM,T because two RMSEs from each model should agree when time-coincident and error-free observations (η = 0 and T = 0) are used in the retrieval.Similarly, b VWC,η should equal to b VWC,T .They can be considered as a "representativeness" error containing contributions from other error sources such as the vegetation model and the inversion algorithm.These linear models were then combined into a bivariate model using the square root of the sum of the squares (RSS) to avoid double counting the representativeness error: (17) where b SM and b VWC were assumed to be b SM,η and b VWC,η , respectively.This model will be validated with the simulation with both errors later in Section IV-H.

III. SOOP-R SYSTEM PARAMETERS
We evaluated the sensitivity of the retrieval error (RMSE) to the following system parameters: frequency combination, polarization combination, measurement error, temporal discrepancy between observations, and inclination of the receiver orbit.This was done by varying three different parameters, (1) frequency combination, (2) observation error, and (3) observation time difference, independently.Each of these studies was performed for linear and circular polarization and two realistic orbit inclinations.

A. Frequency and Polarization
We utilized up to four different microwave frequencies and two types of receiver antenna polarizations in our simulations.The specific frequencies were selected because they are currently in use by operational satellite systems providing frequent coverage to most of the Earth: Orbcomm (137 MHz) [11], the US Navy MUOS (255 MHz, 370 MHz) [13], and GNSS (1575.42MHz) [14].All of these systems use right-hand circular polarization (RHCP).As such, the signals reflected at the specular point will have a large left-hand circular polarization (LHCP) component and the signal-tonoise ratio (SNR) measured by an RHCP antenna in space will be very small compared to the other polarizations.It has also been shown that the RHCP observations do not provide significant benefit in estimating RZSM [17].Therefore, we selected LHCP for our single-polarization observable and an orthogonal pair of linear polarizations (H&V) for our dualpolarization observable.

B. Observation Time Window
Multi-frequency SoOp-R observations from different transmitting satellites are highly unlikely to be obtained at the same time within the same region of interest.Temporal discrepancy between the observations will adversely affect the retrieval accuracy.To investigate this effect, we introduce a new SoOp-R system parameter, observation time window, T , illustrated in Figure 1.Only observations made within a period [t SMi − T /2, t SMi + T /2] (red dots) are used as input to the inversion algorithm in estimating the SM profile at time t SMi .
A retrieval period, P , is used as one of input parameters independent of the observation time window to set the number of SM profiles to be retrieved in our simulation.Due to the different spatiotemporal coverage of transmitting satellites, observations at all the four frequencies may not be available in a given observation time window.In selecting a given frequency combination, we will limit observations to only those times in which all observations needed for the multifrequency retrieval are visible within the observation time window, which could reduce the available coverage.This effect is not evaluated in the present study.

C. Incidence Angle
In SoOp-R system, the only design parameters that determine the observational geometry are the orbital elements of the receiver, given the non-cooperative nature of the sources.The time-varying geometry results in a multi-angular observation with each frequency, in general, observed at different incidence angles but would have some correlation since they are observed from the same receiver location.To find a realistic set of incidence angles for multi-frequency retrieval, we simulated the observation geometry of the three transmitting constellations and a constellation of 8 (chosen arbitrarily) receivers in two different representative orbits, the International Space Station (ISS)'s orbit (400 km, 51.6°inclination) and a polar orbit (400 km, 98°inclination), over a 7-day period.The receivers were equally spaced in true anomaly within the same orbital plane.The number of transmitters trackable by a receiver was assumed to be unlimited.We placed the 3-km resolution Equal-Area Scalable Earth (EASE) Grid 2.0 [57] over the globe and identified cells that include at least one specular point from every transmitting constellation.For each cell flagged as "full frequency", the set of incidence angles corresponding to each frequency was grouped into a vector, θ = [θ 137 , θ 370 , θ 1575 ], where each element corresponds to an incidence angle of a signal at 137, 370, and 1575.42MHz, observed simultaneously.Note that the incidence angle of 255 MHz signal is always the same as that of 370 MHz signal because they are transmitted from the same MUOS satellite.If a cell had multiple observations at the same frequency, one of the incidence angles was randomly chosen and any θ containing an incidence angle greater than 60°was discarded.The 7-day simulation produced an ensemble of 34,320 and 12,787 full-frequency cells for the ISS and the polar orbit, respectively.
Figure 2 presents the histograms of this ensemble for both orbits.The main difference between the two histograms is that the incidence angle of the MUOS signal is centered near the receiver's orbital inclination for the ISS case.As MUOS is in a geosynchronous orbit and its footprint extends up to ±65°latitude, a receiver can make more observations along its inclination as long as it orbits within the MUOS's footprint.The number of full frequency observations in the ISS orbit is notably higher than in the polar orbit.Most Orbcomm satellite orbits have an inclination of 45°and GPS has a global coverage, so there is little difference between these two corresponding distributions.
In order to incorporate this geometric statistics into SM retrieval, we randomly assign θ(k) to the hourly SM data as a known input for the truth model and the retrieval algorithm.Retrieval simulations for the two different orbits were conducted independently and compared to each other.
It is important to note that the purpose of these geometry simulations for this study was to generate a realistic distribution of correlated incidence angles to select from, in generating synthetic observations.Earth coverage, another important consideration in orbit selection, is beyond the scope of this article.

D. Ancillary Parameters
Other important ancillary parameters are soil texture, which determines the soil dielectric properties, and surface roughness, which affects the coherence of the reflected signal.The in-situ clay content in Table I used in the truth model is also given to the retrieval algorithm as input to the Mironov model.The surface roughness under vegetation is assumed to be smooth and follows Kirchhoff approximation for a Gaussian surface with the RMS height of 1 cm.

IV. SIMULATION STUDY
In this section, we describe our simulation procedure and present results of sensitivity studies performed using it.

A. Simulation Environment Setup
Figure 3 depicts the simulated environment consisting of a grass canopy over a flat ground.The boundary conditions and mean SM content are obtained from the in-situ measurements illustrated by the blue dots at five depths (5, 10, 20, 50, and 100 cm).POME is used to obtain a multilayered SM profile discretized into 1-cm interfaces from 5 cm to 100 cm, depicted by a black curve within the soil layer.SM layers from 0 cm Fig. 3. Representation of simulation environment for this article.The grass canopy is modeled as a layer of randomly distributed elliptical disks over a smooth surface.The ground environment is represented as a multilayered SM profile based on the POME model, depicted by the black curve within the soil layer.The SM profile is calculated from the in-situ measurements represented by the blue dots at discrete depths.

Parameters
Grass blade a i,0 Scale parameter s i to 4 cm are assumed to be identical to the SM layer at 5 cm depth.We assume a homogeneous soil texture composition and smooth subsurface layer roughness.The surface roughness represented by the RMS height of 1 cm is included in the surface reflectivity at the vegetation-soil interface, z = 0.
As all the vegetation cover of the selected USCRN sites are categorized as grasslands, a layer of randomly distributed grass blades is placed above the ground.The grass blade is modeled as a thin elliptical disk using vegetation parameters listed in Table III.The nominal values of the major axis, minor axis, thickness, and volumetric moisture content of a grass blade, number density of grass blades, and layer thickness are taken from [38].The scale parameters were arbitrarily selected such that the canopy height is less than ∼1 m and the volumetric moisture is less than ∼0.65 g/cm 3 .Applying the method explained in Section II-B2, the time-varying geometrical shape and distribution of the grass blades were computed via Equation (3) based upon the daily VWC computed from Equation (4) using the daily NDVI data interpolated from the 16-day NDVI product from MODIS [40].

B. Simulation Procedure
We devised three test cases that vary either the frequency combination, standard deviation of observation error, or observation time window, independently, for each of the two receiver polarizations (LHCP and H&V) and two orbital inclinations (51.6°and 98°).Monte Carlo simulations were conducted by computing large ensembles of simulated reflectivity from in-situ Θ v , under different variations of each of these three parameters and then directly computing the RMSE from Equation ( 12) and ( 13).We also performed retrieval simulations with specific parameter values selected based on the sensitivity analysis of each parameter to examine the complex interactions of different error sources.Linear models for RMSE prediction derived from the test cases will be applied and compared to these simulations.
Figure 4 illustrates the method and necessary inputs used to generate an ensemble of synthetic reflectivities at multiple frequencies.The POME-based SM profile Θ v from each of the selected USCRN stations, vegetation structure a, and a set of system parameters β i are used as inputs to the forward scattering model to compute reflectivity Γ veg i .Adding the modeled error ν i to this yields the calculated measurement of reflectivity Γ veg i .Time series of hourly synthetic reflectivities generated from the USCRN in-situ data from 2016 to 2020 were obtained through this approach.
The retrieval algorithm is shown in Figure 5.The state vector used in ASA consists of the 4-parameters defining the POME model and the VWC.Equation ( 9), in which the vegetation is characterized by VWC (though VOD), is used for the forward model.All simulations including the forward and inverse processes were performed on an end-to-end SoOp-R mission simulator [58].Details about each of the three test cases and their corresponding results are presented in the following subsections.

C. Effect of Frequency Combination
In this subsection, we evaluate the sensitivity of retrieval accuracy to different combinations of the observation frequencies, 137, 255, 370, and 1575.42MHz, denoted as a, b, c, and d, respectively.Simulated retrievals are generated using all four frequencies and compared to those using all possible dual, triple, and quadruple combinations of them.To minimize the computational time and maintain the statistical significance, we used one-year data for this case and for Section IV-D.4,025 SM profiles generated from the in-situ data of the 7 USCRN stations in 2016, sampled at 12-hour intervals (P = 12 hr), were used.We assumed a "perfect" instrument with no measurement errors (η = 0) and that all observations were timecoincident (T = 0).Only the frequency combinations were varied, all other parameters were fixed.
Figure 6 depicts the retrieval RMSEs for the SM profile down to 50 cm, SM values sampled at the depths of 5, 10, 20, 50, and 100 cm, and VWC, as a function of frequency combination.The top row of subplots show results using a single LHCP antenna and the bottom row of subplots shows those for a pair of linearly polarized (H&V) antennas.The left column of subplots shows the results for the set of incidence angles representative of the ISS orbit at 51.6°inclination and the right column of shows corresponding results for the incidence angles representative of the polar orbit at 98°i nclination.The red dotted line indicates the baseline RMSE of SM profiles for the top 50 cm soil.
As expected, the the lowest overall RMSE for both SM (all levels) and VWC is realized when observations from all four frequencies are combined.In this case, when the H&V-pol antenna is used, the RMSEs at 5 and 10 cm depths are less than 0.02 m 3 /m 3 and the RMSEs at 20, 50, and 100 cm are less than 0.04, 0.06, and 0.10 m 3 /m 3 , respectively.Frequency combinations abcd, abc, and abd can achieve the baseline RMSE of less than 0.06 m 3 /m 3 for the top 50 cm soil, when using the H&V-pol antenna.In general, all combinations of three frequencies provide lower RMSE than the combinations of two frequencies.These general results are expected, as more frequencies provide more information for the retrieval.The signal at 1575.42 MHz (d) is found to play an important role in estimating VWC.Combinations without d have a higher RMSE than those with d.A high RMSE in VWC results in a high RMSE in SM profiles.A lower frequency generally provides more information on deeper layers when it is combined with a high frequency, as the high frequency helps fix the VWC and the surface SM, while the low frequency provides sensitivity to lower depths.While the use of four frequencies would offer the least RMSE, when the 137 MHz (Orbcomm) signal is not available, 255 MHz and/or 370 MHz (MUOS) can be combined with GNSS to lower the RMSE.MUOS is in a geosynchronous orbit, whereas Orbcomm is in LEO, providing different coverage of the Earth, so the option to use either of them and still provide low RMSE suggests that a dual-linear pol antenna could improve the global coverage, using MUOS in regions where Orbcomm is not visible.
The RMSE of the full frequency case (RMSE abcd ) and the changes in RMSE for each of the frequency combinations relative to the abcd case are provided in Table IV.We use abcd in all of the remaining simulations to investigate the variation of RMSE due to measurement error and observation time delay.The variation will be linearly approximated and the coefficients calculated will be provided in the next subsections.Thus, in future design trade studies, one can roughly predict the retrieval error for a specific frequency combination with the linear fit scaled by a given RMSE ratio.

D. Effect of Measurement Error
Within this subsection, the relationship between the measurement error and the retrieval accuracy is examined.As introduced in Section II-B4, the measurement error of each SoOp-R observation is modeled as Gaussian random error whose standard deviation is set to be a fraction of the true reflectivity at the corresponding SoOp-R configuration.η was varied from 0% to 20%, while the other parameters are fixed.Each element of observation vector Γ veg includes a simulated independent zero-mean Gaussian noise.We used time-coincident (T = 0) observations from all four frequencies (abcd) generated from the same data set used in Section IV-C.
Figure 7 presents the RMSEs of the retrieved SM and VWC as a function of the standard deviation of measurement errors.RMSE appears to have an approximate linear relationship to the relative magnitude of the observation error, η.Straight lines introduced in Section II-E were fitted to the RMSEs of the top 50 cm SM profiles and VWC.The coefficients of the linear equation are summarized in Table V.When the measurement is "perfect", the retrieval error of SM profile and VWC are approximately 0.080 m 3 /m 3 and 0.184 kg/m 2 for the LHCP antenna and 0.040 m 3 /m 3 and 0.066 kg/m 2 for the H&V-pol antenna, respectively.The retrieval errors then increase steadily as the measurement error increases, regardless of polarization and inclination, with a slope of a SM ≈ 0.630 m 3 /m 3 and a VWC ≈ 1.447 kg/m 2 for the LHCP antenna and a slope of a SM ≈ 0.430 m 3 /m 3 and a VWC ≈ 1.040 kg/m 2 for the H&V-pol antenna (average values over receiver's inclinations).The increase in the retrieval error due to the measurement error can be roughly predicted through this approximate linear relationship.When the measurement error exceeds a certain level, the RMSE of SM profiles no longer increases significantly and flattens out.It is more obvious in the deeper layers.This implies that the retrieval accuracy loses its sensitivity to measurement errors as the error grows.Note that the maximum measurement error used in the linear regression is 10%.
The slope of SM RMSE is the steepest at 20 cm depth where the inflection point of the POME model is located, indicating the estimation of the inflection point is sensitive to measurement errors.Comparing the overall slopes of SM RMSEs at 5, 10, and 50 cm depths, the deeper the layer, the steeper the slope.This shows that the retrieval of RZSM is more susceptible to measurement errors than that of surface SM, making it more difficult to estimate.Due to the complexity of SoOp interactions within the multilayered SM structure, even a small increase in measurement errors creates large uncertainties in the retrieval of RZSM.We attribute the decrease in slope for the 100 cm depth to the insensitivity of observations at any of the four frequencies to the SM at 100 cm depth.
In order to achieve the baseline RMSE for the top 50 cm SM profile, the H&V-pol antenna should be used and the standard deviation of measurement errors should be less than 4%.

E. Effect of Observation Timing Gaps
For our third study, we examine the relationship between the observation time window and retrieval accuracy.All four Fig. 8. RMSEs of retrieved SM profile and VWC as a function of observation time window.The RMSEs were calculated for the top 50 cm SM profile, SM at the depths of 5, 10, 20, 50, and 100 cm, and VWC.The top row uses a single LHCP antenna while the bottom row uses a dual linearly polarized antenna.The left column uses the set of incidence angles generated from the ISS orbit at 51.6°inclination while the right column uses the set of incidence angles generated from the polar orbit at 98°inclination.The altitude for both orbits is 400 km.The red dotted line represents the baseline RMSE of SM profiles for the top 50 cm of soil.frequencies (abcd) are used to retrieve a single SM profile and VWC.The time of each observation is randomly selected with equal probability (uniform distribution) within an observation time window, T , and combined in a single retrieval referenced to the center of the time window (refer to Figure 1).For example, given an observation time window of 12 hours, the SM profile at 6 AM is estimated from four frequencies, each observed at a random time between 12 AM and 12 PM.Times are discrete with increments of 1 hour as we utilize the hourly in-situ data from the USCRN stations to generate the simulated observations.The observation time window is varied from 0 (timecoincident) to 7 days while the other parameters are fixed.No measurement error is added (η = 0).For a fair comparison of retrieval accuracy in different observation time windows, target SM profiles to be retrieved must be identical.Considering that the maximum observation time window is 7 days (168 hours), we selected the SM profiles with an interval of 7 days (P = 7 days) for all comparisons.To generate statistically valid RMSE, we used all the five-year SM profiles at the 7 USCRN stations from 2016 to 2020 in this case.The number of the profiles used to calculate the RMSE in each case ranges from 1,133 to 1,280.
Figure 8 depicts the RMSEs of retrieved SM and VWC as a function of the observation time window.A straight line is also fitted to the RMSEs of SM profiles and VWC and the observation time window in hours, with the coefficients listed in Table V.If the observation time window is large, it is more likely that reflectivities measured from other SM profiles and VWC that are temporally more distant from the target SM profile and VWC will be used in the retrieval.Therefore, the RMSEs of both SM profiles and VWC increase as the observation time window increases, regardless of polarization and inclination.The slopes of the linear equation are a SM ≈ 2.410 • 10 −4 m 3 /(m 3 hr) and a VWC ≈ 5.990 • 10 −4 kg/(m 2 hr) for the LHCP antenna and a SM ≈ 2.784 • 10 −4 m 3 /(m 3 hr) and a VWC ≈ 3.310 • 10 −4 kg/(m 2 hr) for the H&Vpol antenna (average values over receiver's inclinations).The approximate linear relationship between the retrieval error and the observation time window can be used to roughly estimate the increase in the retrieval error due to the increase in the observation time difference between subsequent observations at different frequencies.Note that the maximum observation time window used in the linear regression is 72 hours.
To meet the baseline RMSE for the top 50 cm SM profile, the H&V-pol antenna should be used and the observation time window should be less than 72 hours.

F. Effects of Polarization
Both options for receiver antenna polarizations (LHCP and dual H&V) were evaluated in each of the three comparisons in Figures 6 to 8. Comparing the results in Figure 6, the H&Vpolarized antenna generally provides a lower RMSE of SM profiles and VWC compared to the LHCP antenna.This is more evident when the number of frequencies is reduced.The use of two polarizations increases the amount of information available in the retrieval, improving the retrieval accuracy.It is noteworthy that, when comparing the VWC cases where 1575.42MHz is not included, the dual polarization compensates for loss of information and mitigates errors of both VWC and SM profiles.
Similar behavior is also observed in Figure 7 and 8.The overall RMSE of H&V is lower than LHCP.When the measurement error rises from 0% to 4%, the performance degradation of the LHCP antenna is more severe than that of the H&V-pol antenna.As shown in Table V, the slope of the LHCP antenna is greater than that of the H&V-pol antenna, except for a SM for time window at 98°inclination.The RMSE at 5 cm remains below 0.04 m 3 /m 3 until the standard deviation of measurement errors reaches 4% and 15%, respectively, for the LHCP antenna and H&V-pol antenna.In particular, the VWC estimates are significantly improved by approximately a factor of two in any case.This indicates that the additional polarization can compensate for an increase in measurement errors, thereby improving the accuracy of SM and VWC observations to some extent.Overall, the dual linearly polarized antenna can offer more accurate estimates of SM profiles and VWC than the single LHCP antenna.

G. Effects of Orbital Inclination
In the three comparisons, two constellation orbits at different inclinations were used only to generate a realistic ensemble of incidence angles, θ.No noticeable difference between the results for the different constellations can be found from Figure 6 to Figure 8.This suggests that the inclination of the receiver's orbit does not significantly affect the retrieval accuracy, in terms of the distribution of incidence angles.Keep in mind that Earth coverage and revisit time were not evaluated in this study but these still depend strongly on the receiver orbit.Our results could be interpreted as showing that a range of receiver orbits could be considered to meet coverage and revisit requirements without having an effect on retrieval accuracy.

H. Effects of Combined Errors
In all the three test cases, we used ideal SoOp-R parameters except for the variable of interest to investigate the effect of an individual SoOp-R parameter.In this subsection, we examine the retrieval simulation with realistic SoOp-R parameters to account for the complex interaction of different error sources in a realistic spaceborne SoOp-R mission. 1) Validation of bivariate model: First, we validate the bivariate model introduced in Section II-E with simulations encompassing the effects of both measurement error and observation time delay.As an example for validation, we picked the H&V-pol antenna and receiver's inclination of 98°t o run simulations with a range of measurement errors and observation time window.All the four frequencies were selected as our illuminating sources.The corresponding coefficients of the linear models, listed in Table V, were used in Equation ( 16) and ( 17) to calculate the model-based RMSEs of top 50 cm SM and VWC, respectively.Figure 9 shows the comparison of calculated RMSE (blue) and simulated RMSE (red) for the top 50 cm SM and VWC.The bivariate model agrees with the simulation with RMSE of 0.0043 m 3 /m 3 for SM and that of 0.0189 kg/m 2 for VWC, which is negligible for most applications.Interestingly, the slope (of both model and simulation) in each direction gradually decreases as the value of the other parameter increases.This implies the effect of measurement error becomes less as the observation time window increases and vice versa.
2) For different polarizations and inclinations: We now look at the effect of different polarizations and inclinations on the retrieval error and validate the bivariate model.The value of each parameter was selected based on the value that satisfies the baseline RMSE using the H&V-pol antenna with the four frequencies in the above sensitivity analysis.We assumed that the standard deviation of reflectivity errors at each frequency is 4% of its reflectivity value.The observation time window was chosen to be 12 hours (T = 12 hr).Thus, SM profiles with an interval of 12 hours (P = 12 hr) were selected as the target set of SM profiles to be retrieved.All SM profile data from the 7 USCRN stations between 2016 and 2020 were used in this simulation, resulting in 18,235 profiles on average retrieved in each case.
Figure 10(a) depicts the RMSE of retrieved SM on the xaxis as a function of soil depth on the y-axis.The RMSE of VWC retrieved with the SM profiles is shown in Figure 10(b).A total of four simulations were performed for combinations of two polarizations (LHCP and H&V-pol) and two receiver orbits (51.6°and 98°inclinations).The circle marker represents the prediction of RMSE based on the bivariate model.Each marker in Figure 10(a) has the same color as its corresponding case.It shows that the bivariate model predicts the retrieval error very close to the simulation.
The RMSE curve of SM profiles can be divided into three parts along the soil depth.The first part is from 0 cm to 5 cm, where the RMSE is constant due to the assumption that the SM layers are identical within this range.The second part is from 5 cm to 20 cm, which is above the inflection point of the POME model located at a depth of 20 cm.The third part is from 20 cm to 100 cm, which is below the inflection point.The RMSE above the inflection point grows faster than the RMSE below the inflection point as the soil depth increases.This implies that the retrieval performance may be improved by adjusting the depth and SM value of the inflection point as tuning parameters of the inversion algorithm.
In terms of polarization and inclination, the retrieval performance of the H&V-pol antenna is superior to that of the LHCP antenna in the entire column of SM profiles and in VWC.The receiver's inclination, however, does not affect the retrieval performance significantly.These confirm the conclusions of Section IV-F and IV-G.The RMSEs of SM profile for the top 50 cm depth are approximately 0.110 m 3 /m 3 and 0.058 m 3 /m 3 for the LHCP antenna and the H&V-pol antenna, respectively.The RMSEs of VWC are around 0.233 kg/m 2 and 0.092 kg/m 2 for the LHCP antenna and the H&V-pol antenna, respectively.Moreover, the RMSE of SM profile reaches 0.06 m 3 /m 3 at depths of approximately 13 cm and 56 cm for the LHCP antenna and H&V-pol antenna, respectively.The addition of polarization significantly improves the retrieval accuracy of deeper layers.We attribute the retrieval error primarily to the insensitivity of signals to deeper layers, attenuation due to the vegetation, uncertainties of the vegetation parameterization, and inaccuracy inherent in the inversion algorithm.
3) For different frequency combinations: For this analysis, we look at the effect of different combinations of frequencies on the retrieval error and validate the bivariate model.We picked the H&V-pol antenna and receiver's inclination of 98°a gain and ran retrieval simulations with combinations of four and three frequencies (abcd, abc, abd, acd, and bcd) and other SoOp-R parameters selected above (η = 4%, T = 12 hr).
The RMSE ratios given in Table IV were used to scale the coefficients of the bivariate model.Figure 11 shows the RMSE of retrieved SM profiles and VWC for different frequency combinations.The circle marker represents the prediction of RMSE based on the bivariate model.Each marker color in Figure 11(a) is the same as that of the corresponding case.While the prediction for abd is very close to the simulation, the other three combinations (abc, acd, and bcd) show some errors.Particularly, the lower simulated RMSE of acd and bcd indicates that the correlation between two errors (measurement error and observation time delay) is strong at the given observational configuration.In other words, this means that if the effect of one error is large, the effect of the other error is offset.Thus, the bivariate model overestimates the retrieval error.More accurate predictions can be obtained by fitting the linear model for each frequency combination directly, rather than scaling.
In terms of the effect of the frequency combinations, the combination of four frequencies provides the best performance of the retrieval of both SM profiles and VWC.Among the combinations of three frequencies, abc showed the lowest retrieval error of SM profiles, followed by abd, acd, and bcd.For the VWC retrieval, abc performs the worst among others.
We can see that the exclusion of the highest frequency d leads to the degradation of VWC retrieval.Overall, given the highest frequency d, a lower error can be obtained when a combination of lower frequencies is used together.The inclusion of the lowest frequency a has the greatest impact on improving retrieval errors in the deeper layers.Given two frequencies, bc, bd, and cd, adding a further reduces RMSE of the top 50 cm SM by roughly 0.028, 0.025, and 0.01 than adding d, c, and b, respectively.These results are consistent with those in Section IV-C.

V. DISCUSSION
In this section, we summarize significant findings and discuss several suggestions and limitations.The sensitivity of retrieval error of RSZM and VWC has been examined through simulations that vary each of the SoOp-R system parameters, including frequency, polarization, orbital inclination, observation time window, and measurement error.The observations regarding frequency and polarization in our simulation study are in agreement with the results from [17] where the Cramer-Rao lower bound was used to determine the lower bound of the estimation variance for RZSM modeled by two SM slabs.The results presented in this article, however, are useful in that the end-to-end retrieval algorithm is applied to the large set of the actual SM data, which is utilized to develop multilayered SM profiles.The multistatic configuration and temporal discrepancy of multi-SoOp-R observation were also taken into account in the simulations to mimic realistic spaceborne SoOp-R observation.
We observed that the additional frequencies generally reduced the retrieval error.Having a high frequency along with a low frequency is important for estimating RZSM by fixing the surface SM content and VWC in the retrieval algorithm.In addition, lower RMSE was achieved when lower frequencies were used in the frequency combination.This suggests that low frequencies play the most important role in detecting RZSM.
As for polarization, it is observed that the H&V-pol antenna provides more accurate estimates of SM profiles and VWC than the LHCP antenna by approximately a factor of two.Considering that circular polarization involves horizontal and vertical polarizations, the high retrieval accuracy achieved by the H&V-pol antenna is largely due to the addition of measurements used in the retrieval algorithm.
With respect to orbital inclination, the selection of receiver's orbital inclination does not appear to have a significant impact on the retrieval accuracy.This means that as long as signals of interest are visible from the receiver's orbit, the orbital inclination of a receiver can be chosen according to other things such as the required spatial coverage.Since the coverage of Orbcomm and MUOS extends to mid-latitudes [58], we can leverage them to improve the retrieval accuracy of SM profiles in low-and mid-latitudes, where most regions of interest for SM observations lie.Other illuminating sources would be required to cover high-latitude regions.
The retrieval error increases as the maximum latency between observations at different frequencies being combined increases.The temporal resolution requirement, usually justified by temporal variation of SM content, is an important factor in determining retrieval accuracy in multi-SoOp-R observation.The observation time window introduced in our study can be used to determine the temporal resolution of multi-frequency observations.For example, having the observation time window of 12 hours in Section IV-H will offer one SM observation every 12 hours, which is subdaily time scale.This temporal resolution can be achieved by a large number of small SoOp-R satellites in a constellation.
Measurement errors induced by instruments and scattering environments, indeed, will degrade retrieval accuracy.A highaccuracy instrument and a precise error model will be needed to mitigate this impact.
The bivariate model for predicting the retrieval error was introduced as a function of measurement error and observation time window and was validated with the simulations.We found a correlation between the two errors that as the effect of observation time delay increases, the effect of measurement error becomes less important and vice versa.
Several assumptions were made in this article.First, the SoOp-R measurement error was assumed to be additive zeromean Gaussian noise whose standard deviation is a fraction of the reflectivity.This may not be true for all frequencies in different scattering environments, which necessitates a different error model that adequately characterizes measurement noise.Second, considering the wavelength of the signals used in this study, we supposed that the reflected signals were dominantly coherent on the smooth surface (RMS height of 1 cm).However, this assumption may not be valid under rough surface and complex topography conditions.In this case, incoherent signals may need to be included in the retrieval algorithm.The effects of surface roughness and topography on retrieval accuracy also need to be dealt with in depth.Third, since the spatial resolution of coherent reflection is determined by the first Fresnel zone defined by the wavelength of the signal, they are not equal to each other at different frequencies.Thus, we assumed that the retrieval was performed over a region where the footprints of these observations overlap.Lastly, we assumed that the inflection point of the dynamic SM profile of the POME model was fixed at 20 cm depth in any case, which is not true for the actual SM profiles.This could deteriorate the performance of the retrieval algorithm.The depth of the inflection point should be adaptively searched in the inversion algorithm for each SM profile.
We neglected some effects other than the SoOp-R parameters in our analysis.Given an observation time window, the temporal variability of SM profiles in a local area affects the retrieval accuracy.High retrieval error can be caused by high temporal variability of SM profiles, which is affected by meteorological and physical factors such as precipitation and soil texture [59].Analyzing these effects depending on different climate types and soil hydraulic properties is beyond the scope of this article.
In addition, the RMSE presented here also includes errors induced by the inversion algorithm, which can be further improved by addressing the details described above.Defining a retrieval confidence flag can be useful for the quality control of retrieval products.Incorporating other hydrologic variables such as precipitation and evapotranspiration observed from other sensors may be helpful to estimate the mean SM content used in the POME model in the inversion algorithm.

VI. CONCLUSION
In this article, we presented the forward and inverse algorithm to retrieve a SM profile and VWC with a multi-SoOp-R observation and analyzed the sensitivity of retrieval accuracy to SoOp-R system parameters such as frequency combination, measurement error, time delays between measurements, receiver polarization, and receiver orbital inclination.The insitu SM data from the selected USCRN stations were used to construct multilayered SM profiles modeled by the POME model.A dynamic vegetation structure model was applied to represent seasonal variations of VWC for grasslands.The synthetic reflectivities were generated through the SCoBi-Veg model.Using the synthetic observations, the retrieval of SM profiles and VWC was performed by adaptive simulated annealing with the least squares cost function.A variety of aspects of multi-SoOp-R observation has been addressed in the sensitivity analysis, which is useful for developing requirements for a future SoOp-R mission.We recommend that both high and low frequencies should be used together to achieve high retrieval accuracy of SM profiles and VWC.Specifically, the GPS L1 signal (1575.42MHz) plays an important role in limiting uncertainties from vegetation and surface SM and the Orbcomm signal (137 MHz), the lowest frequency in the SoOp of our interest, provides sensitivity to the deepest soil layers.We found that the relationship between retrieval error and reflectivity measurement error is approximately linear.The retrieval error also linearly increases as the time delays between combined observations increase.A correlation between the measurement error and the observation time window is found.The dual linearly polarized antenna would provide more accurate estimates than the LHCP antenna.The receiver's orbital inclination does not significantly affect the retrieval accuracy.Given the reflectivity error of 4% and the observation time window of 12 hours, the four frequencies at 137, 255, 370, and 1575.42MHz provided the lowest retrieval error, followed by the combinations of three frequencies excluding 1575.42,370, 255, and 137 MHz, respectively, in that order.These SoOp-R system parameters should be carefully considered in the design of spaceborne SoOp-R missions.The bivariate model derived from the sensitivity analysis will provide a rough estimation of the retrieval error of a SM profile and VWC for future design studies.

Fig. 1 .
Fig. 1.Illustration of concept of observation time window T .An observation time window T is defined as a period of time during which SoOp-R observations are collected for retrieving a SM profile.SoOp-R observations within T (red dots) are used in the retrieval of SM profile at t SM i .

Fig. 4 .
Fig. 4. Method to generate an ensemble of synthetic reflectivity vectors, Γ veg , each consisting of a set of noisy observations for N obs combinations of frequencies and polarizations (i = 1, . . ., N obs ).

Fig. 5 .
Fig. 5. Inversion algorithm to retrieve a SM profile Θv and VWC W .The state vector, consisting of the 4-parameter state x for POME model and VWC, is iteratively estimated via adaptive simulated annealing.The tilde ( ) is used to indicate the calculated value while the hat ( ) is used to indicate the estimated value.

Fig. 6 .
Fig. 6.RMSEs of retrieved SM profile and VWC as a function of frequency combination.Frequency combination is denoted by the combination of a, b, c, and d which correspond to 137, 255, 370, and 1575.42MHz, respectively.The RMSEs were calculated for the top 50 cm SM profile, SM at the depths of 5, 10, 20, 50, and 100 cm, and VWC.The top row uses a single LHCP antenna while the bottom row uses a dual linearly polarized antenna.The left column uses the set of incidence angles generated from the ISS orbit at 51.6°inclination while the right column uses the set of incidence angles generated from the polar orbit at 98°inclination.The altitude for both orbits is 400 km.The red dotted line represents the baseline RMSE of SM profiles for the top 50 cm of soil.

Fig. 7 .
Fig.7.RMSEs of retrieved SM profile and VWC as a function of standard deviation of Gaussian random error added to the reflectivity in percent (%).The RMSEs were calculated for the top 50 cm SM profile, SM at the depths of 5, 10, 20, 50, and 100 cm, and VWC.The top row uses a single LHCP antenna while the bottom row uses a dual linearly polarized antenna.The left column uses the set of incidence angles generated from the ISS orbit at 51.6°inclination while the right column uses the set of incidence angles generated from the polar orbit at 98°inclination.The altitude for both orbits is 400 km.The red dotted line represents the baseline RMSE of SM profiles for the top 50 cm of soil.
Fig.10.RMSE of retrieved SM profiles and VWC for selected SoOp-R system parameters, which are a period of 12 hours, an observation time window of 12 hours, a standard deviation of measurement errors of 4%, and four frequencies at 137, 255, 370, and 1575.42MHz.Individual results are depicted for the combination of two polarizations (LHCP and H&V-pol) and two receiver orbital inclinations (51.6°and 98°).The circle marker represents the prediction of RMSE based on the bivariate model.Each marker in Figure10(a) has the same color as its corresponding case.

TABLE II VALUES
OF PARAMETER b f,q USED FOR RETRIEVAL.THE TRANSMITTER POLARIZATION STATE IS RIGHT-HAND CIRCULAR POLARIZATION (RHCP) WHILE THE RECEIVER POLARIZATION, INDICATED BY THE SUBSCRIPT q, ARE LEFT-HAND CIRCULAR POLARIZATION (LHCP), HORIZONTAL POLARIZATION (H-POL), AND VERTICAL POLARIZATION (V-POL),

TABLE IV RETRIEVAL
ERRORS OF THE FOUR FREQUENCY CASE (RMSE abcd ) AND RATIOS OF RMSE FOR THE FREQUENCY COMBINATIONS RELATIVE TO THE abcd CASE

TABLE V COEFFICIENTS
OF STRAIGHT LINES FITTED TO THE RETRIEVAL ERRORS.SEE SECTION IV-D AND SECTION IV-E.