Sentinel-1 soil moisture content and its uncertainty over sparsely vegetated fields

Soil moisture content (SMC) retrievals from synthetic aperture radar (SAR) observations do not exactly match with in situ references due to imperfect retrieval algorithms, and uncertainties in the model parameters, SAR observations and in situ references. Information on the uncertainty of SMC retrievals would contribute to their applicability. This paper presents a methodology for deriving the SMC retrieval uncertainty and decomposing this in its constituents. A Bayesian calibration framework was used for deriving the total uncertainty and the model parameter uncertainty. The methodology was demonstrated with the integral equation method (IEM) surface scattering model, which was employed for reproducing Sentinel-1 backscatter (σ0) observations and the retrieval of SMC over four sparsely vegetated fields in the Netherlands. For two meadows the calibrated surface roughness parameter distributions are remarkably similar between the ascending and the descending Sentinel-1 orbits as well as between the two meadows, and yield consistent SMC retrievals for the calibration and validation periods (RMSDs of 0.076 m3 m 3 to 0.11 m3 m 3). These results are promising for operational retrieval of SMC over meadows. In contrast, the surface roughness parameter distributions of two fallow maize fields differ significantly and the surface roughness conditions changing over time result in less consistent SMC retrievals (calibration RMSDs of 0.096 m3 m 3 and 0.13 m3 m 3 versus validation RMSDs of 0.26 m3 m 3). The SMC retrieval uncertainty derived with the Bayesian calibration successfully reproduces the uncertainty estimated empirically using in situ references. The main uncertainty originates from the in situ references and the Sentinel1 observations, whereas the contribution from the surface roughness parameters is relatively small. The presented research yields further insights into the surface roughness of agricultural fields and SMC retrieval uncertainties, and these insights can be used to guide SAR-based SMC product developments.


Introduction
The soil moisture content (SMC) is a key state variable in climatological, meteorological, hydrological and ecological processes. Its control on the exchanges of water and energy at the land surface plays an important role in the development of climate and weather systems (Global Climate Observing System, 2016;Massari et al., 2014;Seneviratne et al., 2010). In addition, it is important for the partitioning of rainfall in infiltration and runoff (Beck et al., 2009;Massari et al., 2014;Wanders et al., 2014), regarded as an indicator for the onset of droughts (Miralles et al., 2016;Seneviratne et al., 2010;Vautard et al., 2007), and essential for vegetation growth (Feddes et al., 1976;Ines et al., 2013). Hence, information about the SMC would benefit a number of applications.
Backscatter (σ 0 ) observations by synthetic aperture radar (SAR) instruments can be used to estimate the SMC at much finer scale, even up to agricultural field scale (e.g. Amazirh et al., 2018;El Hajj et al., 2017;Lievens and Verhoest, 2012;Su et al., 1997). Bauer-Marschallinger et al. (2019) developed an operational 1 km resolution SMC product from Sentinel-1 SAR σ 0 observations, based on a change detection algorithm that assumes static surface roughness and vegetation conditions. However, at the field scale this assumption is unlikely to hold because spatial surface roughness and vegetation effects on the σ 0 are not averaged out over a large area (Bauer-Marschallinger et al., 2019). In those situations, the relation between the σ 0 signal and SMC must be separated from the effects of surface roughness and vegetation before the SMC can be estimated reliably (Kornelsen and Coulibaly, 2013;Paloscia et al., 2013;Verhoest et al., 2008). Physically based scattering models, such as the integral equation method (IEM) for surfaces (Fung et al., 1992) and the Tor Vergata model for vegetation (Bracaglia et al., 1995), simulate the scattering contributions from soil-vegetation systems based on prescribed electromagnetic characteristics. This supports the application of these models to various site conditions and sensor configurations (Paloscia et al., 2013;Petropoulos et al., 2015), the understanding of backscattering processes (Baghdadi et al., 2002;Balenzano et al., 2012;Wang et al., 2018) and the propagation of uncertainty sources (Satalino et al., 2002;Van der Velde et al., 2012a).
Surface scattering models, including the frequently-used IEM model, simulate the scattering of electromagnetic waves from a surface and are used to estimate the σ 0 from soils (Ulaby and Long, 2014). The surface roughness essentially governs the σ 0 response and, thus, the sensitivity to SMC. The parameterisation of the surface roughness is, therefore, an important input. Measuring the surface roughness was part of many field campaigns, such as EMAC'94 (Su et al., 1997), FLOODGEN in 1994, 1998(Baghdadi et al., 2004, Orgeval'94 (Zribi et al., 1997), OPE3 in 2002 (Joseph et al., 2010), SMAPVEX12 (McNairn et al., 2015), SMAPVEX16-IA (Hornbuckle et al., 2017) and SMAPVEX16-MB (McNairn et al., 2016). However, Baghdadi et al. (2002), Baghdadi et al. (2004) and Su et al. (1997) have shown that the IEM model does not accurately reproduce σ 0 observations using measured surface roughness parameters. Lievens et al. (2011) and Verhoest et al. (2008) attributed this to both uncertainties in the surface roughness measurements and simplifications in the representation of surfaces.
A pragmatic approach for applying surface scattering models to land surfaces is considering the surface roughness parameters as 'effective parameters', obtained by model calibration instead of field measurements (Baghdadi et al., 2002;Lievens et al., 2011;Lievens and Verhoest, 2012;Rahman et al., 2008;Su et al., 1997;Verhoest et al., 2008;Verhoest et al., 2007). The calibration of the surface roughness parameters is accomplished by searching for a parameter set that results in a match between σ 0 observations and model simulations. Subsequently, the calibrated surface roughness parameters can be used to retrieve SMC from other σ 0 observations and/or on other fields (Su et al., 1997).
In addition to the surface roughness parameterisation, the SMC estimates from σ 0 observations will contain uncertainties specific for the selected retrieval algorithm (De Lannoy et al., 2014;Pathe et al., 2009) and due to uncertainty in the σ 0 observations Pathe et al., 2009). The σ 0 observations contain uncertainty from calibration uncertainties, sensor instabilities and speckle effects, which are together referred to as radiometric uncertainty . Furthermore, SMC references are required for the calibration of scattering models and the validation of the SMC retrievals. The SMC references are typically obtained from in situ measurements. This introduces uncertainties due to a SMC probe's measurement uncertainty (Cosh et al., 2005) and spatial scale mismatches with satellite-observed SMC (Western and Blöschl, 1999;Cosh et al., 2006). A horizontal spatial scale mismatch between the SMC at an in situ monitoring station and fieldaveraged SMC originates from differences in land cover, soil texture and structure, and local features such as nearby ditches and subsurface drainage pipes. A vertical spatial scale mismatch originates from the (Sentinel-1 C-band) σ 0 observations having a sampling depth that varies from the surface to a depth of 1 cm to 10 cm (Nolan and Fatland, 2003;Ulaby et al., 1996), whereas in practice SMC measurements at 5 cm or 10 cm depth, with an influence zone of e.g. 4 cm above and below the probe (Benninga et al., 2018), often have to be adopted for calibration and validation purposes (e.g. Bauer-Marschallinger et al., 2019;Chan et al., 2018;Kornelsen and Coulibaly, 2013;Pathe et al., 2009;Van der Velde et al., 2015).
Information on the uncertainty of SMC retrievals is essential to assess their reliability and for their applicability, for example for the assimilation of SMC retrievals into land surface models and for combining SMC products (Pierdicca et al., 2014;Verhoest et al., 2007;De Lannoy et al., 2014). Verhoest et al. (2007) estimated the uncertainty of SMC retrievals with the IEM model by defining uncertainty distributions for the surface roughness parameters. As a result of assumed surface roughness parameter uncertainties of ±7.5%, ±15% and ±25%, Verhoest et al. (2007) reported SMC retrieval uncertainties (standard deviations) of 0.023 m 3 m − 3 , 0.041 m 3 m − 3 and 0.060 m 3 m − 3 , respectively. Vernieuwe et al. (2011) continued on the study by Verhoest et al. (2007) by considering the correlation between the parameters based on a synthetically generated surface roughness data set. Doubková et al. (2012) and Pathe et al. (2009) estimated the uncertainty of SMC retrievals from parameter uncertainty assumptions and the radiometric uncertainty. Pulvirenti et al. (2018) defined fuzzy logic rules in order to assign a degree of uncertainty (low, medium, high) to each SMC retrieval. These previous studies, however, relied on assumptions regarding the uncertainty of model parameters for the estimation of the SMC retrieval uncertainty. This is reflected in the applied calibration methods in general, which ignore uncertainties and aim for one optimal parameter set that results in a match between observations and simulations (e.g. Joseph et al., 2008;Lievens et al., 2011;Verhoest et al., 2007).
Bayesian calibration approaches allow for the derivation of parameter distributions and the separation of parameter uncertainty from the total simulation uncertainty, based on statistical assumptions of which the validity can be verified (Barber et al., 2012;De Lannoy et al., 2014;Haddad et al., 1996;Notarnicola et al., 2006;Notarnicola and Posa, 2004;Pierdicca et al., 2014). For example, using semi-empirical Oh surface scattering models (Oh et al., 1992;Oh et al., 2002), Haddad et al. (1996), Pierdicca et al. (2014) and Pierdicca et al. (2010) formulated Bayesian frameworks for the retrieval of surface roughness parameters and SMC along with estimates of their retrieval uncertainty.
Bayesian frameworks cannot be solved analytically for highly nonlinear models (Vrugt, 2016), such as physically based scattering models. To provide an efficient solution for such models, the DiffeRential Evolution Adaptive Metropolis package (DREAM; Vrugt (2016)) implements a multi-chain Markov chain Monte Carlo simulation algorithm for generating samples from the posterior distributions that describe the parameter uncertainty and the total simulation uncertainty. De Lannoy et al. (2014) used DREAM to calibrate a radiative transfer model for simulating SMOS L-band brightness temperatures and to estimate the uncertainty of the parameters and the total simulation uncertainty.
In this study, the uncertainties involved in surface scattering model simulations and SMC retrievals were investigated. We focused on the calibration of the IEM surface roughness parameters, and, therefore, used Sentinel-1 σ 0 observations and SMC measurements from sparsely vegetated fields, namely two meadows and two fallow cultivated parcels. The Bayesian calibration was performed with DREAM. The paper extends on previous research on SMC retrieval from SAR σ 0 observations by (1) adopting a Bayesian calibration framework for deriving the uncertainty of the IEM surface roughness parameters and the total uncertainty, (2) assessing the derived SMC retrieval uncertainty against the uncertainty estimated empirically using in situ references, and (3) decomposing the total uncertainty in its four constituents.

Definitions of uncertainties
The standard deviation is selected as uncertainty measure. The standard deviation of the differences between two data sets, such as SMC retrievals and references, is often referred to as the unbiased root mean square deviation (uRMSD; Kerr et al. (2016)): where N stands for the number of match-ups between estimates (Y e ) and references (Y r ), t stands for the observation number and the bars denote the means of Y e and Y r . A SMC retrieval, its total uncertainty and constituents are illustrated in Fig. 1. The surface roughness parameters for retrieving the SMC, as well as the parameter uncertainty (U p ) and total uncertainty were derived with Bayesian calibrations, using DREAM as described in Section 4. We refer to the total uncertainty that is derived with the Bayesian calibration as U total− B . The surface roughness parameter set with the highest posterior probability, also referred to as the 'maximum a posteriori' (MAP; (Vrugt, 2016;De Lannoy et al., 2014;Lu et al., 2017)), was used for the optimal SMC retrieval. U total− B should be of similar magnitude as the empirical uncertainty of SMC retrievals for cases that the Bayesian calibration was statistically valid (De Lannoy et al., 2014). The empirical SMC retrieval uncertainty can be calculated with Eq. (1) using in situ references. The U total− B and U p are visualized by two histograms in Fig. 1, which partly overlap and show that the distribution of the U total− B is wider than the distribution of the U p . This is expected as U p is one of the constituents of the total uncertainty.
The other constituents are inherent to the in situ references and the satellite observations, namely the measurement uncertainty of the station probes providing the in situ references (U sp ), the in situ references' uncertainty attributable to a spatial scale mismatch with Sentinel-1 observed SMC (U s,S1 ), and Sentinel-1's radiometric uncertainty (U S1 ). Fig. 1 illustrates that U sp and U s,S1 apply to the in situ references, and U S1 applies to the SMC retrievals. In Section 3, the U sp , U s,S1 and U S1 are quantified.
Next to its estimation with the Bayesian calibration, the total uncertainty can be found by combining its constituents. This is referred to as U total− C and can be calculated by following the addition rule for variances (Moore et al., 2017): where Cov stands for the covariance terms between the uncertainty constituents. The constituents are assumed to be uncorrelated, whereby Cov reduces to 0. The relative contributions of U sp , U s,S1 , U S1 and U p can then be calculated in a similar fashion as was done in Van der Velde et al. (2012b): RC s,S1 = U s,S1 and With Eq.
(2) we can evaluate to what extent U total− C explains U total− B , and with Eqs. (3)-(6) we can assess their individual relative contributions.

Study region, fields and periods
The SMC measurements used as references were collected by monitoring stations in the Twente region, located in the eastern part of the Netherlands (Fig. 2). The Twente region is flat with some elevated glacial ridges and it has a temperate oceanic climate with a Cfb Köppen-Geiger climate classification (Beck et al., 2018). The SMC monitoring stations in this region are collectively known as the Twente network (Dente et al., 2011(Dente et al., , 2012Van der Velde et al., in review, 2019).
Stations are installed at the border of fields for safety and continuity reasons. Adjacent to monitoring stations, we selected two meadows (hereafter field I and II) and two fallow cultivated fields (field III and IV) as study fields for which we collected additional field measurements in total on 87 occasions. The study fields and the locations of the field measurements are shown in Fig. 3. Field I and III are adjacent to the same monitoring station. The study fields have loamy sandy surface layers. Supplement 1 details the study fields' surface layer soil textures and bulk densities from the soil physical map of the Netherlands (BOFEK2012; Wösten et al. (2013)). Table 1 lists the study periods. We used the winter season of October 2016 -March 2017 for the calibration of the surface roughness parameters and the winter season of October 2017 -March 2018 for the validation of σ 0 simulations and SMC retrievals. The study periods are taken outside the growing season, from October (after harvesting and other agricultural practices) to March (before ploughing and sowing), so that the fields were fallow or covered with non-growing sparse In situ reference Reference's measurement uncertainty (U sp ) Reference's spatial mismatch uncertainty (U s,S1 ) Optimal or 'MAP' soil moisture content retrieval Sentinel-1's radiometric uncertainty (U S1 )

Model parameter uncertainty (U p )
Total uncertainty (U total ) Fig. 1. The SMC total retrieval uncertainty and its constituents, which are quantified in this study. The arrows represent one standard deviation. The two histograms partly overlap.
vegetation and no agricultural practices were applied during the study periods. In between the winter seasons, several agricultural practices are applied on cultivated fields (field III and IV), such as sowing, harvesting, manuring and ploughing. On meadows (field I and II) typically no ploughing is applied and the surface roughness is expected to change little. Table 2 lists the land covers at the locations of the SMC monitoring stations and on the study fields during the calibration and validation period. Field I and II are covered with grass which is virtually static and sparse during winters: field measurements with the LI-COR LAI-2000 (LI-COR, 1992) indicated leaf area indices (LAI) of 1.1 m 2 m − 2 and 1.3 m 2 m − 2 outside the growing season versus maximums of 8.0 m 2 m − 2 and 7.7 m 2 m − 2 in the growing season for field I and II, respectively (Benninga et al., 2020a). Field III and IV were fallow with remaining maize stubble in the winter of 2016/2017. In the winter of 2017/2018 field III was again fallow with maize stubble, whereas field IV was used to grow winter wheat. Table 2 indicates that also the land cover at station IV changed between the calibration and the validation period. This station was installed on 20 May 2016. In the first period after installation the land cover at the station's location was similar to the study field. The second year it was covered with grassy vegetation because the area with the station probes was no longer directly subjected to agricultural practices.

Soil moisture content references
The SMC monitoring stations are equipped with 5TM probes (METER Group, 2019) installed at nominal depths of 5 cm, 10 cm, 20 cm, 40 cm and 80 cm, of which the readings are stored every 15 min. We used the 5 cm SMC measurements collected at Sentinel-1 overpass times as the in situ references. The probes at 5 cm depth provide an integrated measurement over a soil depth of 1 cm to 9 cm (Benninga et al., 2018).
The SMC measurements inside the study fields revealed an inconsistency in the station measurements of field IV. During the period May 2016 -November 2016 the station measurements had a bias of − 0.024 m 3 m − 3 with respect to the field measurements, whereas for the period April 2017 -September 2017 the bias increased to − 0.12 m 3 m − 3 (see Table S3). This likely is a consequence of the change in the land cover at the station's location from fallow with maize stubble to grassy vegetation (see Table 2 and Section 3.1).

Measurement uncertainty
A soil-specific calibration function was developed for the station probes of the Twente network (Van der Velde et al., in review, 2019).
The calibration accuracy, quantified by Eq. (1) between the calibrated probe measurements and gravimetrically determined volumetric SMC (GVSMC) references, is 0.027 m 3 m − 3 . We adopted this value as general measure for the U sp .

Spatial scale mismatch uncertainty
The horizontal and vertical spatial scale mismatches have a systematic and a variable impact on differences between the SMC references and Sentinel-1 observed SMC. The systematic component is a bias which will implicitly be accounted for via the calibration of the surface roughness parameters. The variable component is U s,S1 , which contributes to the uncertainty of Sentinel-1 SMC retrievals.
The U s,S1 was quantified by Eq. (1) between the station measurements and the spatial mean of the 0 cm -5 cm layer SMC measurements inside field I -IV. Supplement 2 provides further information on the estimation of U s,S1 . The values for U s,S1 in Table S3 demonstrate that adopting the station measurements as reference for the Sentinel-1 SMC retrievals introduces a significant amount of uncertainty, varying from 0.036 m 3 m − 3 to 0.068 m 3 m − 3 . We adopted the mean of 0.051 m 3 m − 3 over field I -IV as the common measure for U s,S1 .

Data processing
The Sentinel-1 constellation provides images in C-band (5.405 GHz), over land in Interferometric Wide Swath (IW) mode at VV and VH polarization. We only used the Sentinel-1 σ 0 observations in VV polarization because of the higher expected sensitivity to SMC than the VH polarization (e.g. Amazirh et al., 2018;El Hajj et al., 2017) and because the definitions of the surface roughness parameters in the IEM model are different for VV and VH due to underlying assumptions. The radiometric accuracy is specified at 1 dB (three standard deviations). After multilooking, the Ground Range Detected (GRD) High Resolution (HR) product has a resolution of 20 m × 22 m (4.4 equivalent number of looks) (Torres et al., 2012;Bourbigot et al., 2016).
To obtain Sentinel-1 backscatter (σ 0 ), we downloaded Level-1 GRD HR IW Sentinel-1 images from the Copernicus Open Access Hub (Copernicus, 2019) and processed them using the following operations in ESA's Sentinel Application Platform (SNAP) V6.0 (European Space Agency, 2019): (1) Apply Orbit File, (2) Thermal Noise Removal, and (3) Range Doppler Terrain Correction, including radiometric normalization to σ 0 (in m 2 m − 2 ) with projected local incidence angles on a geographic grid (WGS84) with a pixel spacing of 9.0E-5 • (equivalent to 10 m × 6.1 m at the study region's latitude). Subsequently, the Sentinel-1 σ 0 observations were averaged over the study fields, excluding a 20 m distance from the borders of the fields and 40 m from trees and buildings to avoid possible influences of features outside the fields (see the net area in Table 2). The last step was to convert the σ 0 values to decibel (dB). Table 3 specifies the orbits that cover the study region. Sentinel-1A provides images since 3 October 2014 and Sentinel-1B since 28 September 2016. The combination of Sentinel-1A and Sentinel-1B gives a temporal resolution of 1.5 days over the study region. However, frozen conditions, wet snow and intercepted rain can disturb σ 0 observations and we masked the Sentinel-1 observations for these weather-related surface conditions with the masking rules presented in Benninga et al. (2019), which are summarized in Supplement 3. Furthermore, in situ references that decreased during frozen soil periods (Van der Velde et al., in review, 2019) were removed, and from 18 January 2018 to 16 March 2018 the SMC monitoring station adjacent to field II was malfunctioning and no references are available for this period. Table 4 lists the number of Sentinel-1 observations with a matching in situ reference, before and after masking for the above-mentioned weather-related surface conditions.

Radiometric uncertainty
Sentinel-1's radiometric uncertainty (s S1 , in dB) was estimated by the standard deviation of Sentinel-1 σ 0 observations from a target which is assumed time-invariant . This resulted in a second-order power function between s S1 and the surface area over which the Sentinel-1 σ 0 observations are averaged. The SMC retrieval uncertainty due to the s S1 (being U S1 ) is then derived through combination with the σ 0 to SMC sensitivity, which follows from simulations with the IEM model.

Surface scattering model application
IEM is a physically based surface scattering model (Fung et al., 1992)  that has widely been used to simulate the σ 0 from bare and sparsely vegetated land surfaces. Readers are referred to Ulaby and Long (2014) for more background on the IEM model and to Kornelsen and Coulibaly (2013) for a discussion of previous studies in which IEM was used. Vegetation effects are not accounted for by the IEM model, and accordingly, we limited the calibration and validation periods to the fallow or non-growing sparse vegetation conditions outside the growing season (see Section 3.1). The applicability of the IEM model to sparse grass covers is justified by the results of Van der Velde and Su (2009) (Knyazikhin et al., 1999;Tesemma et al., 2014), which is comparable to our LAI measurements outside the growing season. Moreover, the SMC retrieval uncertainties attributable to vegetation effects were found to be fairly small compared to uncertainties caused by the surface roughness parameterisation (Van der Velde et al., 2012a).
The IEM model requires parameterisations on the dielectric and geometric properties of the land surface. The dielectric properties were estimated with the Mironov soil dielectric mixing model (Mironov et al., 2009). SMC and the soil textures from Supplement 1 served as input to the Mironov model. The geometry of the land, also known as the surface roughness, is parameterised by the root mean square surface height (s), the autocorrelation length (c l ) and an autocorrelation function. The exponential autocorrelation function was selected, because it is viewed as most appropriate for agricultural fields (Ulaby and Long, 2014;Verhoest et al., 2008). Callens et al. (2006) have shown that changes in surface roughness due to heavy rainfall are limited when no agricultural practices were applied recently. The study periods were taken such that no agricultural practices were applied within them (see Section 3.1) and, therefore, surface roughness can be assumed time-invariant. This assumption is discussed in Section 5.3.
Agricultural surfaces are generally anisotropic (Verhoest et al., 2008) and the fields are viewed from a different direction in Sentinel-1's ascending and descending orbits, causing that the surface roughness would be different for these orbits. This is especially expected for the fallow maize fields, which have tillage rows (see Fig. 3). The meadows do not have a clear row structure. The anisotropy of the study fields was considered by separating the calibration of s and c l for the Sentinel-1 σ 0 observations made in the ascending and in the descending orbits. By combining the two ascending orbits (15 and 88) and the two descending orbits (37 and 139) respectively, the surface roughness parameters were calibrated on two incidence angles (see Table 3) and the varying SMC conditions encountered during the calibration period.

Bayesian model calibration
Bayesian model calibration derives posterior parameter distributions conditioned on prior parameter distributions (prior) and the match between model simulations and reference data (likelihood), by solving Bayes' rule (Vrugt, 2016): where p(θ|z) is the resulting posterior probability density function (PDF) of the parameters (θ) given the reference data (z). The likelihood function evaluates how well the model reproduces z given θ, by describing the PDF of the residuals between simulations and references. The generalized likelihood function, derived by Schoups and Vrugt (2010), offers a wide flexibility in heteroscedasticity, distribution and autocorrelation of the residuals. The likelihood model parameters have to be inferred jointly with the model parameters or must be given a fixed value. The validity of the residual model can be verified with a residual analysis. For more background on residual analysis, readers are referred to Lu et al. (2017), Scharnagl et al. (2011), Schoups and Vrugt (2010) and Thyer et al. (2009).

Application of DREAM
We adopted a simple implementation of the generalized likelihood function, and assumed homoscedastic, Gaussian and uncorrelated residuals. These assumptions are often made and convenient to use (e.g. Lu et al., 2017;Raj et al., 2018;Scharnagl et al., 2011), and lead to the common standard least squares approach (Schoups and Vrugt, 2010). In Section 5.1, the validity of the residual model is verified with a residual analysis. The standard deviation of the residuals (σ 0 ) has to be inferred with the calibration. The two surface roughness parameters in combination with σ 0 bring the dimensionality (number of unknowns) at three.
The validity of the IEM model is limited to medium surface roughness conditions with ks⩽3, where k is the free-space wavenumber (Baghdadi et al., 2004;Su et al., 1997). For the wavelength of Sentinel-1,    Table 4 The number of Sentinel-1 observations for which a matching in situ reference is available. Total  After masking  Total  After masking   I  87  50  96  54  II  74  42  63  40  III  92  54  96  54  IV  79  55  98  56 this corresponds to a maximum s of 2.68 cm. For c l no IEM validity domain has been formulated. Calibration ranges of 0.2 cm to 400 cm have been used, and the resulting calibrated c l values ranged from 1.4 cm to 13 cm for maize and bare agricultural fields (Joseph et al., 2010;Joseph et al., 2008;Lievens et al., 2011;Satalino et al., 2002;Verhoest et al., 2007) and from 0.2 cm to 7 cm for a mosaic of grasslands and wetlands (Van der Velde et al., 2012a). Non-informative (uniform) priors are preferred for scientific objectivity (Lunn et al., 2013;Notarnicola and Posa, 2004;Notarnicola et al., 2006). We defined the prior distributions as uniform distributions with the ranges (0.1 cm, 2.68 cm)

Field Calibration Validation
for s and (0.1 cm, 100 cm) for c l . The prior distribution of σ 0 is defined as a uniform distribution with ranges (0 dB, 2 dB). We used the standard DREAM settings (Vrugt, 2016), with ten Markov chains. A burn-in of 50% of the realizations is recommended to allow initialization to the posterior parameter distributions (Vrugt, 2016). Convergence of the chains was assessed by the multivariate Gelman-Rubin convergence diagnostic R d , where R d below 1.2 indicates convergence (Brooks and Gelman, 1998;Vrugt, 2016), and by visual inspection of the mixing of the Markov chains (Raj et al., 2018;Vrugt, 2016). 7000 realizations per chain appeared sufficient to reach convergence after 50% of the realizations, which results in 35000 samples describing the posterior parameter distributions.

Soil moisture content retrieval
The MAP surface roughness parameter set was used for the optimal σ 0 simulations and SMC retrievals. These were evaluated against the Sentinel-1 σ 0 observations and in situ SMC references, respectively, with the root mean square deviation (RMSD), the unbiased RMSD (uRMSD) and the Pearson correlation coefficient (r P ), defined in Supplement 4.
For the retrieval of SMC from the Sentinel-1 σ 0 observations, we generated look-up tables of σ 0 simulations for SMC values ranging from 0.01 m 3 m − 3 to 0.75 m 3 m − 3 , with an increment of 0.001 m 3 m − 3 , and combinations of soil textures, incidence angles and surface roughness parameter sets. A SMC retrieval is then taken equal to the SMC value for which the minimum difference between σ 0 simulations and a Sentinel-1 σ 0 observation is found.
For deriving U total− B , we generated 1000 σ 0 residual samples from the skew exponential power distribution that underlies the likelihood function (Schoups and Vrugt, 2010), using the σ 0 that was found for the MAP surface roughness parameter set. The resulting 1000 SMC retrievals with the MAP surface roughness parameters, after superimposing the σ 0 residual samples on a Sentinel-1 σ 0 observation, describe U total− B . For the computation of U p , we randomly sampled 1000 surface roughness parameter sets from their posterior distributions and derived 1000 SMC retrievals.

Residual analysis
The residual analysis plots for the Bayesian calibrations are included in Supplement 5.1. The figures (a) in Figs. S1-S4 show for fields I to IV that the residual variances are generally independent of the simulated σ 0 , which justifies the use of a homoscedastic residual model. The figures (b) show that the deviations from the theoretical quantiles for a Gaussian distribution are only substantial for a few σ 0 simulations at the tails and not systematic among the calibration cases, so we accepted the validity of the Gaussian residual distribution. Only for field IV a number of outliers can be observed, which is further discussed with regard to the σ 0 simulations in Section 5.3.
The calibration cases show some autocorrelation, with mean values of 0.40 at a lag of one time step and 0.28 at a lag of two time steps (figures (c) in Figs. S1-S4). In the Bayesian calibration of process models, such as rainfall-runoff models (Schoups and Vrugt, 2010) and terrestrial ecosystem models (Lu et al., 2017), autocorrelation in the residuals can be accounted for with autoregressive residual models. However, the IEM model does not contain state variables. Using autoregressive residual models, therefore, does not change the posterior parameter distributions nor the residual analysis plots of our calibration results. In Supplement 6 this is demonstrated by showing for field I the calibration results obtained with a first-order and a second-order autoregressive residual model. Supplement 5.2 contains the residual analysis plots for the validation period. Figs. S5-S8 show that the homoscedastic Gaussian residual model is generally also valid for the validation period. The quantile-quantile plots already give an outlook on the performances of the σ 0 simulations and SMC retrievals in the validation period. Regarding field I, the quantile-quantile plot (Fig. S5b) is steeper (larger dispersion) than the quantile-quantile plot for the calibration period and it reveals a bias (compare to the plot's origin, (0, 0)). Hence, a 'slightly' degraded performance and a bias are expected for the validation period. For field II the quantile-quantile plots for the calibration and the validation period (Figs. S2b and S6b) are comparable, so we expect similar performances. For field III (Fig. S7b) and field IV (Fig. S8b) steep lines and high biases are observed, suggesting worse performances for the validation period.

Posterior parameter distributions
Figs. 4-7 show the posterior parameter sets and the MAP s and c l for fields I to IV. Joseph et al. (2008), Lievens et al. (2011), Rahman et al. (2008 and Verhoest et al. (2007) already reported that multiple optimal combinations of s and c l are possible. The scatter plots in Figs. 4-7 illustrate that the posterior distributions of the surface roughness parameters actually cover a large part of the solution space, and that the s and c l are highly correlated (Spearman's rank correlation coefficient, r S , is 0.97 to 1.0). Individual values of s or c l , therefore, do not contain much information about the surface roughness. For example, s values of 0.5 cm to 1.5 cm are found in the posterior parameter distributions of all the fields. Hence, both s and c l or a ratio between them should be used to characterize the roughness of a surface. From Figs. 4-7 it is clear that the relation between s and c l is non-linear and that the simple s/c l ratio will not suffice. Instead, it is approximately a square root relation and the parameter Z s = s 2 /c l (Zribi and Dechambre, 2002) is suitable for characterizing the roughness of the surfaces.
For the meadows (Figs. 4 and 5), the ascending and the descending orbits' posterior distributions coincide. In other words, the surface roughness is similar for the Sentinel-1 σ 0 observations made in the ascending and in the descending orbits. This is an indication that the meadows have an isotropic surface roughness, at least in Sentinel-1's ascending and descending orbit viewing directions. Therefore, we also calibrated the surface roughness parameters with the Sentinel-1 σ 0 observations from both passes combined, of which the results are also shown in Figs. 4 and 5. The parameter sets obtained from the combined calibration were used in the remainder of this paper.
Furthermore, the posterior parameter distributions of the two meadows are quite similar. The MAP s values are 0.16 cm and 0.18 cm, and the c l values are 1.31 cm and 1.49 cm for field I and field II, respectively. In Section 5.3 we discuss the cross-validation results of the MAP surface roughness parameters of field I applied to retrieve the SMC for field II, and vice versa.
For the fallow maize fields (Figs. 6 and 7), the ascending and the descending orbits' posterior distributions are different. This was expected, as these fields do have an anisotropic surface due to tillage rows and these are viewed from different angles in the ascending and descending orbits.

Retrievals
The MAP SMC retrievals, U p and U total− B are plotted as time series in Fig. 8, and Table 5 lists the performance metrics of the MAP SMC retrievals for the calibration and the validation period. Time series and performance metrics of the forward σ 0 simulations, using the SMC references and calibrated surface roughness parameters as input to the IEM model, are shown in Fig. 9 and Table 6 respectively.

Meadows
The performance of the meadows' MAP σ 0 simulations is comparable for the calibration and the validation period. This indicates that the surface roughness remained similar, which can be explained by the fact that no ploughing was applied on the meadows. The increase in the empirical uncertainty (uRMSD, Eq. (1)) of the SMC retrievals can be explained by the wetter conditions during the validation period. IEM model simulations show that the σ 0 to SMC sensitivity diminishes with   increasing SMC, see for example Fig. 3 in Altese et al. (1996) and the results in Benninga et al. (2019), which results in larger SMC deviations for equal σ 0 deviations under wetter conditions. Because of this, the SMC retrieval uncertainty distributions in Fig. 8 are also wider at higher SMC and they are skewed towards the higher SMC levels.
The posterior surface roughness parameter distributions and the MAP values are quite similar for the two meadows. To further verify this, we performed a cross-validation by retrieving the SMC for field II using the MAP surface roughness parameters of field I, and vice versa. Table 5 lists the SMC retrieval performances, and Supplement 7 includes the SMC and σ 0 time series figures. The calibration has aimed to optimize the RMSD of the σ 0 simulations, so it could be expected that the RMSD for the calibration period is higher using the MAP parameter set of the other meadow. In general, the SMC retrieval performances are comparable using the MAP surface roughness parameter sets obtained for the other meadow.

Fallow fields
The validation performances for the fallow fields are considerably worse than the calibration performances. Field III was fallow with maize stubble during both the calibration and the validation period, but the surface roughness is likely to be different due to agricultural practices in between. Furthermore, Figs. 8c (and 9c) show in the validation period three distinctive periods: from 14 October 2017 to 10 November 2017 with low SMC retrievals (high σ 0 simulations), from 15 November 2017 to 14 January 2018 with high SMC retrievals (low σ 0 simulations), and from 15 January 2018 to 21 March 2018 with low SMC retrievals (high σ 0 simulations). However, the field was harvested in September/October 2017 and not ploughed after May 2017, so changes in surface roughness due to agricultural practices and heavy rainfall were not expected. This result demonstrates that even under these circumstances, the fallow field cannot be simulated with a single set of surface roughness parameters.
Besides, Sentinel-1 σ 0 observations of up to -3 dB in the calibration period and -2 dB in the validation period were acquired for field III, with a maximum of -1.62 dB. These Sentinel-1 observations all originate from relative orbit number 15. The IEM model with an exponential autocorrelation function for the surface roughness cannot reproduce such high σ 0 observations with any set of surface roughness parameters. Therefore, we omitted the Sentinel-1 σ 0 observations of field III acquired in relative orbit 15 for further analysis in this study. An additional calibration was performed exclusively on the Sentinel-1 observations acquired in relative orbit 88, of which the results are presented in Supplement 8. The posterior parameter distributions (Fig. S14) and the σ 0 simulation performances (Fig. S17) for this calibration are, however, similar to the calibration on both ascending orbits. Using the original surface roughness parameter sets and omitting the orbit 15 Sentinel-1 σ 0 observations does improve the σ 0 simulation and SMC retrieval performances (Tables 5 and 6).
For field IV, the validation performances are more degraded than for field III. This can be explained by the different land covers in the calibration and the validation period (Table 2) and by the bias in the in situ references during the validation period (Section 3.2). Fig. 10a shows the residuals of the MAP SMC retrievals with the original references and with the references corrected for the bias of -0.12 m 3 m − 3 . Part of the residuals can indeed be explained by this bias. However, still three periods can be distinguished in the residuals: between 17 October 2017 and the sowing of the winter wheat on 10 November 2017 the RMSD against the bias-corrected references is 0.10 m 3 m − 3 (13 observations), between 13 November 2017 and 6 December 2017 the RMSD is smallest with a value of 0.050 m 3 m − 3 (10 observations), and after 15 December 2017 the RMSD is 0.19 m 3 m − 3 (32 observations). The development of the winter wheat vegetation on this field during the validation period does not have a large effect, as this should otherwise be visible as a gradual trend in the residuals extending to April 2018. Moreover, at the end of the validation period the wheat cover is still sparse, as is shown in Fig. 10b. A number of heavy rainfall events occurred between 6 December 2017 and 15 December 2017 (in total 64 mm). Callens et al. (2006) demonstrated that rainfall smoothens the surface and reduces the surface roughness on recently tilled fields. Indeed, the Sentinel-1 σ 0 observations being lower after 15 December 2017 is in accordance with a reduced surface roughness.
A number of outliers were observed in the residual analysis plots of the calibrations on field IV. As visualized in Fig. S4, the σ 0 simulations in the first part of the calibration period, between 14 October 2016 and 6 November 2016, hold the largest residuals and all these residuals are on one side of the quantile-quantile plots. This indicates that the surface roughness conditions has changed within the calibration period. As discussed in the previous paragraphs, for both fallow fields the same holds between the calibration and the validation period as well as within the validation periods.

Note on the soil moisture content references
It should be noted that the SMC references extend to higher levels than saturated SMCs generally observed. BOFEK2012 (Wösten et al., 2013) lists saturated SMC values of 0.44 m 3 m − 3 -0.45 m 3 m − 3 for the surface layers (0 cm to 23 cm depth) of field I to IV. These values are exceeded by the station SMC measurements. This can be partly attributed to higher organic matter content and root density near the soil surface. Organic matter increases SMC values especially in sandy and loamy sandy soils (Minasny and McBratney, 2018). In addition, local soil variability is not captured by BOFEK2012, local soil variability is not    . 9. Sentinel-1 σ 0 observations and simulations. The parameter and total simulation uncertainty are visualized by the 95% confidence interval.

Table 6
Performance metrics of the MAP σ 0 simulations against the Sentinel-1 σ 0 observations. considered in the probes' calibration function, and roots and macropores in the probes' influence zone can increase measured SMC (Benninga et al., 2018). However, even with consideration of these factors, the very high SMC measurements, especially for field II, seem unrealistic in absolute sense. Nevertheless, the correlations between the station and field measurements, listed in Table S3, are high. It can, therefore, be expected that the station measurements capture the temporal variability of the adjacent field's SMC.
The absolute SMC measurement values may still deviate from realistic values. This will affect the surface roughness parameters obtained by the calibration, and for the SMC retrieval over independent periods or fields it may be necessary to apply an unbiasing procedure. This is reflected in the meadows' cross-validation results, see Tables 5 and 6: the RMSDs, which include the bias in the mean, are generally higher than the original calibration metrics, whereas the r P and uRMSDs, which exclude this bias, are comparable. Fig. 11 shows U total− B in comparison to the uncertainty of the MAP SMC retrievals estimated empirically using the SMC references and Fig. 12 shows U total− B relative to U total− C , for bins of SMC references. The empirical uncertainty is quantified with Eq. (1), but without removing the bias for each bin separately to preserve the integrity of the time series' PDFs. For field I and II both the calibration and the validation period are included (Figs. 11a and b). Since it was found that the parameters calibrated for the cultivated fields (field III and IV) are invalid for the validation period, the latter period is not included in Figs. 11c-f. As a consequence of that and because the ascending and descending orbits are separated for field III and IV, the total number of pairs is larger for field I and II and for visualization purposes the number of pairs per bin in Figs. 11 and 12 is ten for field I and II and five for field III and IV.

Retrieval uncertainty
The increasing empirical uncertainty and U total− B with increasing SMC in Fig. 11 are explained by the diminishing σ 0 to SMC sensitivity with increasing SMC, as was discussed in Section 5.3.1. Both the increasing trend and the magnitude of the empirical uncertainty are rather closely approximated by U total− B . In other words, the SMC retrieval uncertainty derived with the Bayesian calibration does successfully reproduce the uncertainty estimated empirically. This does, however, not hold for field IV. As explained in Section 5.3.2, the IEM model does not correctly reproduce the σ 0 of field IV within the calibration period with a single set of surface roughness parameters. As a consequence, the likelihood function implementation with a homoscedastic residual standard deviation is not valid over the complete calibration period. Fig. 12 shows that the combination of U sp , U s,S1 , U S1 and U p , i.e. U total− C , approximately explains U total− B , except for field IV again. Fig. 12 also shows the relative squared contributions of U sp , U s,S1 , U S1 and U p , namely RC sp , RC s,S1 , RC S1 and RC p . The RC p is relatively small and constant across the investigated SMC domain, with an average of 13% over the SMC domain and fields I to III. This means that the U p increases with SMC because the total uncertainty increases with increasing SMC. From the assumption that U sp and U s,S1 are equal to 0.027 m 3 m − 3 and 0.051 m 3 m − 3 across the entire SMC domain follows that their relative contributions (RC sp and RC s,S1 ) decrease with increasing SMC because the total uncertainty increases with SMC. The average RC sp and RC s,S1 decrease, respectively, from 13% and 46% at a SMC of 0.26 m 3 m − 3 to 4% and 15% at a SMC of 0.53 m 3 m − 3 . The average RC S1 increases from 31% at a SMC of 0.26 m 3 m − 3 to 67% at a SMC of 0.53 m 3 m − 3 , which is explained by the increasing U S1 with increasing SMC . The U S1 is found to be the dominant driver for the increasing SMC retrieval uncertainty with increasing SMC. For field III, the RC S1 is even larger than for the other fields (at an equal SMC level) because of field III's smaller surface area.

Conclusions
The total uncertainty and its constituents were investigated for SMC retrievals from Sentinel-1 σ 0 observations over four sparsely vegetated fields (two meadows and two fallow cultivated fields). A Bayesian framework was used for calibrating the surface roughness parameters that are input to the IEM surface scattering model, and for deriving the parameter and total uncertainty distributions. Subsequently, these distributions were used to retrieve the SMC and its uncertainty, and the relative contributions of four uncertainty sources were evaluated. This resulted in the following conclusions: 1. The simplest implementation of the likelihood function, using a homoscedastic Gaussian residual model, describes the simulation residuals. An exception is when the IEM model is not capable of reproducing the Sentinel-1 σ 0 observations in a calibration or validation period with a single set of surface roughness parameters. 2. The surface roughness parameters (s and c l ) are highly correlated, with Spearman's rank correlation coefficients (r S ) of 0.97 to 1.0. The s and c l have approximately a square root relation and the parameter Z s = s 2 /c l , which was already introduced in Zribi and Dechambre (2002), is shown to be suitable for characterizing the roughness of the surfaces. This result also implies that it is valid to fix one of the parameters s or c l for simplifying the calibration while still acquiring the same posterior Z s distribution. 3. For the two meadows the surface roughness parameter distributions coincide for Sentinel-1's ascending and descending orbits, despite the different directions from which Sentinel-1 views the fields in these passes. Furthermore, the surface roughness parameter distributions of the two meadows are quite similar. In contrast, for the two fallow fields the surface roughness parameter distributions depend on the pass direction and the distributions differ between the two fields. This is attributed to the anisotropic nature of these surfaces caused by tillage rows. Fig. 11. The empirical uncertainty of SMC retrievals and U total− B , for bins of SMC references. The number of pairs per bin is ten for field I (a) and II (b) and five for field III (c-d) and IV (e-f).
4. The performance of the SMC retrievals for the calibration period, expressed by the RMSD, is between 0.076 m 3 m − 3 and 0.13 m 3 m − 3 . The validation results for an independent period confirm that, for the meadows, the surface roughness parameters can be used across years. For the fallow fields, however, the surface roughness conditions change; not only between the calibration and the validation period, but even within single winter periods.
5. The total SMC retrieval uncertainty derived with the Bayesian calibration successfully reproduces the uncertainty estimated empirically using in situ references, including the trend of increasing uncertainty with increasing SMC. 6. The in situ references' measurement uncertainty (U sp ) and spatial scale mismatch uncertainty (U s,S1 ), the SMC retrieval uncertainty due to Sentinel-1's radiometric uncertainty (U S1 ) and the parameter Fig. 12. U total− B and its four constituents relative to U total− C . The bins are the same as in Fig. 11. uncertainty (U p ) constitute the total uncertainty. The main uncertainty originates from the in situ references and the Sentinel-1 σ 0 observations, whereas the contribution from the surface roughness parameters is relatively small.
The two meadows' coinciding surface roughness parameter distributions for the ascending and descending orbits, their similar surface roughness and consistent SMC retrievals for the calibration and validation period are promising results for operational retrieval of SMC over meadows. The value of such a SMC product would be substantial as meadows cover a major portion of the land in use for agriculture, e.g. 71% in the study region Twente and 55% in the Netherlands in 2017 (Ministerie van Economische Zaken, 2017). Therefore, further research to the selection of a common surface roughness parameter set for meadows and the associated retrieval uncertainty would be interesting.
To improve the performance of the Sentinel-1 SMC retrievals it will be essential to reduce the in situ references' uncertainties and the radiometric uncertainty. The references' uncertainties can be reduced by averaging multiple spatially distributed measurements. Reducing the impact of radiometric uncertainty can be achieved by accepting a coarser spatial resolution or a further improvement of the SAR image processing.
By the Bayesian calibration of the IEM model, further insights into the surface roughness of agricultural fields and SMC retrieval uncertainties have been derived. These insights can be used to guide SARbased SMC product developments. Moreover, the study shows the utility of Bayesian calibration approaches for deriving such new insights and the presented methodology may serve as an example for the Bayesian calibration of other scattering model applications.

CRediT authorship contribution statement
Harm-Jan F. Benninga

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.