# Sentinel-1 soil moisture at 1 km resolution: a validation study

Anna Balenzano, Francesco Mattia, Giuseppe Satalino, Francesco Lovergine, Davide Palmisano, Jian Peng, Philip Marzahn, Urs Wegmüller, Oliver Cartus, Katarzyna Dąbrowska-Zielińska, et al.

► To cite this version:

Anna Balenzano, Francesco Mattia, Giuseppe Satalino, Francesco Lovergine, Davide Palmisano, et al.. Sentinel-1 soil moisture at 1 km resolution: a validation study. Remote Sensing of Environment, 2021, 263, pp.112554. 10.1016/j.rse.2021.112554. hal-04829319

HAL Id: hal-04829319

https://hal.science/hal-04829319v1

Submitted on 11 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sentinel-1 soil moisture at 1km resolution: a validation study

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

Anna Balenzano<sup>a, *,</sup> Francesco Mattia<sup>a</sup>, Giuseppe Satalino<sup>a</sup>, Francesco P. Lovergine<sup>a</sup>, Davide Palmisano<sup>a</sup>, Jian Peng<sup>b, #</sup>, Philip Marzahn<sup>b</sup>, Urs Wegmuller<sup>c</sup>, Oliver Cartus<sup>c</sup>, Katarzyna Dabrowska-Zielinska<sup>d</sup>, Jan P. Musial<sup>d</sup>, Malcolm W. J. Davidson<sup>e</sup>, Valentijn R.N. Pauwels<sup>f</sup>, Michael H. Cosh<sup>g</sup>, Heather McNairn<sup>h</sup>, Joel T. Johnson<sup>i</sup>, Jeffrey P. Walker<sup>f</sup>, Simon H. Yueh<sup>j</sup>, Dara Entekhabi<sup>k</sup>, Yann H. Kerr<sup>l</sup> and Thomas J. Jackson<sup>g</sup>

<sup>a</sup> National Research Council of Italy (CNR), Institute for Electromagnetic Sensing of the Environment (IREA), UOS Bari, Italy

<sup>b</sup> Ludwig-Maximilians Universität München (LMU), Department of Geography, Munich, Germany

<sup>c</sup> Gamma Remote Sensing Research and Consulting AG (GAMMA), Gümligen, Switzerland

<sup>d</sup> Institute of Geodesy and Cartography (IGiK), Remote Sensing Centre, Warsaw, Poland

<sup>e</sup> European Space Agency, Mission Science Division, Noordwijk, The Netherlands

<sup>f</sup> Monash University, Department of Civil Engineering, Clayton, Victoria, Australia

<sup>g</sup> USDA-ARS Hydrology and Remote Sensing Laboratory, Beltsville, Maryland, USA

<sup>h</sup> Agriculture and Agri-Food Canada (AAFC), Ottawa, Ontario, Canada

<sup>i</sup> Ohio State University, Dep. of Elect. and Computer Engineering, Columbus, Ohio, USA

<sup>j</sup> Jet Propulsion Laboratory (JPL), California Institute of Technology, Pasadena, California, USA

<sup>k</sup> Massachusetts Institute of Technology (MIT), Civil and Environmental Engineering, Cambridge, Massachusetts, USA

<sup>l</sup> Centre d'Etudes Spatiales de la BIOsphere (CESBIO), Toulouse, France

* Corresponding author at: National Research Council of Italy, Institute for Electromagnetic Sensing of the Environment
E-mail address: anna.balenzano@cnr.it. (A. Balenzano)

# Now at Helmholtz Centre for Environmental Research (UFZ), Department of Remote Sensing and Leipzig University, Remote Sensing Centre for Earth System Research, Leipzig, Germany.

## Abstract

26

27

28

29

30

31

32

33

This study presents an assessment of a pre-operational soil moisture product at 1 km resolution derived from satellite data acquired by the European Radar Observatory Sentinel-1 (S-1), representing the first space component of the Copernicus program. The product consists of an estimate of surface soil volumetric water content Θ [m³/m³] and its uncertainty [m³/m³], both at 1 km. The retrieval algorithm relies on a time series based Short Term Change Detection (STCD) approach, taking advantage of the frequent revisit of the S-1 constellation that performs C-band Synthetic Aperture Radar (SAR) imaging. The performance of the S-1 Θ product is estimated through a direct comparison between 1068 S-1 Θ images against in situ Θ measurements acquired
by 167 ground stations located in Europe, America and Australia, over 4 years between January 34
2015 and December 2020, depending on the site. The paper develops a method to estimate the 35
spatial representativeness error (SRE) that arises from the mismatch between the S-1 Θ retrieved at 36
1 km resolution and the in situ point-scale Θ observations. The impact of SRE on standard 37
validation metrics, i.e., root mean square error (RMSE), Pearson correlation (R) and linear 38
regression, is quantified and experimentally assessed using S-1 and ground Θ data collected over a 39
dense hydrologic network (4 – 5 stations/km²) located in the Apulian Tavoliere (Southern 40
Italy). Results show that for the dense hydrological network the RMSE and correlation are ~0.06 41
m³/m³ and 0.71, respectively, whereas for the sparse hydrological networks, i.e., 1 station/km², 42
the SRE increases the RMSE by ~0.02 m³/m³ (70% Confidence Level). Globally, the S-1 Θ 43
product is characterized by an intrinsic (i.e., with SRE removed) RMSE of ~0.07 m³/m³ over the 44
Θ range [0.03, 0.60] m³/m³ and R of 0.54. A breakdown of the RMSE per dry, medium and wet Θ 45
ranges is also derived and its implications for setting realistic requirements for SAR-based Θ 46
retrieval are discussed together with recommendations for the density of in situ Θ observations. 47

48 **Keywords**

49 Soil Moisture, High Resolution, Sentinel-1, Synthetic Aperture Radar (SAR), Spatial Representativeness
50 Error (SRE), Validation.
51

52 **1 Introduction**

Measurements of Earth's surface soil moisture (Θ) at global scales and at spatial resolutions
of 20 × 20 km² or coarser are currently provided as products of the Soil Moisture and Ocean
Salinity (SMOS) mission of the European Space Agency (ESA) (Kerr et al., 2010), the Soil
Moisture Active Passive (SMAP) mission of the National Aeronautics and Space Administration
(NASA) (Entekhabi et al., 2010), and the Advanced SCATterometer (ASCAT) system aboard the
Meteorological Operational (MetOp) platform of the European Organisation for the Exploitation of
Meteorological Satellites (EUMETSAT) (Wagner et al., 2013). Satellite-measured Θ has proved
useful for improving understanding of the global water and energy cycles (McColl et al., 2017; Seneviratne et al., 2010) and strengthening land applications such as large scale hydrological modelling (Heimhuber et al., 2017; Lievens et al., 2016), numerical weather prediction (NWP) (Dharssi et al., 2011; Rodríguez-Fernández et al., 2019), flood forecasting and drought monitoring and prediction (Mishra et al., 2017; Nicolai-Shaw et al., 2017; Wanders et al., 2014). Despite the usefulness of existing products, significant interest remains in improving the spatial resolution of Θ products to extend and facilitate applications such as mapping the impact of irrigation on local water budgets, assessing the impact of local Θ variability on atmospheric instability and improving NWP and hydrological modelling at regional scales (Dorigo et al., 2017; Peng et al., 2021). In response to these science and application needs, a number of recent studies have proposed techniques to downscale microwave Θ products at low resolution using optical and thermal data (see Peng et al., 2017; Sabaghy et al., 2018 for review), although these approaches are subject to corruption by cloud cover that impacts the optical imagery. An additional approach has become possible following the launch of the European Radar Observatory Sentinel-1 (S-1), developed in the framework of the Copernicus programme (www.copernicus.eu). S-1 systematically provides C-band Synthetic Aperture Radar (SAR) imagery from two identical spacecraft, (S-1 A & S-1 B), at high spatial and moderate temporal (6-day exact repeat cycle) resolutions with a sustained observation strategy for the next decades which foresees first the S-1 C & S-1 D satellites from 2022 onwards and then the S-1 Next Generation satellites from 2028 onwards (Torres et al., 2020, 2012). Spaceborne SAR sensors are currently the most suitable systems to retrieve Θ at high spatial resolution at spatial scales ranging from local to regional and continental. In the past, the use of SAR systems - such as the Advanced Synthetic Aperture Radar (ASAR) aboard the ENVIRONMENTAL SATellite (ENVISAT) - for Θ retrieval has been hampered by their inappropriate observational characteristics, particularly their long revisit times. Also, the validation of past SAR retrieval algorithms has been limited to relatively short campaigns often adopting different strategies for the
in situ sampling of Θ, although other validation approaches have been also used (e.g. Das et al., 2014; Tomer et al., 2015).

This paper presents a pre-operational Θ product, derived from VV&VH S-1 observations at 1 km resolution and its validation status. The retrieval uses a time series method introduced in Balenzano et al., (2011), further developed in Balenzano et al., (2013) and applied in consecutive papers (e.g., Al-khaldi et al., 2019; Iacobellis et al., 2013; Ouellette et al., 2017). In this study, the algorithm has been consolidated and extensively validated, and technical challenges for optimized processing of a product at regional/continental scales, 1 km resolution and 6-12 day revisit, have been addressed. In this respect, the algorithm is considered mature for uptake and use in operational environments such as Copernicus Services operated by the European Union, or other services developed by space agencies or downstream users. The long-term continuity of S-1 data plays a critical and positive role supporting this development. The technique exploits the frequent revisit of S-1 to realize a time series based Short Term Change Detection (STCD) algorithm applicable for bare and vegetated areas dominated by soil attenuated scattering. The code implementing the algorithm is referred to as SMOSAR ("Soil Moisture retrieval from multi-temporal SAR data"). There are two main differences compared to previous papers addressing Θ retrieval from S-1 data (e.g., Bauer-Marschallinger et al., 2019, 2018; El Hajj et al., 2017; Paloscia et al., 2013; Pulvirenti et al., 2018). The first is that the developed S-1 Θ product includes uncertainty information in terms of the Θ standard deviation provided as a coregistered layer at the same resolution and unit. It is noted that, unlike previous studies that used the propagation of uncertainties from SAR observations to Θ retrieved values - see Gruber et al., (2020) for a critical review - SMOSAR provides the observed standard deviation of Θ at 1 km as measure of the uncertainty. Such information allows the imaged areas to be discriminated into different levels of uncertainty (Merchant et al., 2017) that responds to the needs of data assimilation (Pan and Wood, 2010). Second, an extensive validation study of the product was also conducted. The study adopts procedures and metrics recommended by the
Committee on Earth Observation Satellites (CEOS) Working Group on Calibration and Validation (WGCV) (Montzka et al., 2020). In particular, the implemented multi-scale validation activity consists of a comparison of S-1 Θ estimates with in situ observations collected over seven cal/val sites located in the USA, Canada, Australia and Europe. Per each site, the objective has been to analyse S-1 time series in a timeframe of 4 years between January 2015 and December 2020.

A crucial aspect in the time series comparison of satellite estimated Θ against in situ Θ observations is the spatial mismatch between the point-scale (~0.1 m) in situ measurements and the satellite estimates retrieved at resolutions ranging from tens of kilometres (e.g., SMAP, SMOS, ASCAT) to hundreds of meters (e.g., S-1 Θ), which generates the so-called spatial representativeness error (SRE). This issue has previously not been considered consistently in the validation of satellite Θ products at high resolution, e.g., ≤ 1 km. Little effort has been dedicated both to quantify the corresponding SRE and to set up cal/val sites dedicated to high-resolution Θ retrieval. In this study, an emphasis is placed on addressing the SRE for S-1 Θ retrievals at 1 km. Two measures were developed. The first consists of setting up a core validation site, located in the Apulian Tavoliere (Italy), that is characterized by a dense network of ground stations, i.e., ~4 stations/km² (Balenzano et al., 2014). Such a dense network allows upscaling of the in situ Θ observations at ~1 km with a relatively low SRE (e.g., SRE < 0.03 m³/m³). The second measure develops a method for modelling the SRE across scales and the evaluation of its impact on the standard metrics at ~1 km resolution.

The paper is organized as follows. In Section 2, the test sites and ground data are described. Section 3 then summarizes the S-1 data and the low-resolution satellite Θ products analysed in the study. In Section 4, the Θ retrieval algorithm and the S-1 Θ product at 1 km are presented. Sections 5 and 6 illustrate the validation approach and results. Finally, conclusions are drawn in Section 7, including recommendations on SAR-derived Θ product validation requirements.
134

## 2 Test sites and ground data

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

The SMOSAR performance assessment was carried out over the following seven sites: Apulian Tavoliere (Italy), Red de Estaciones de Medición de Humedad del Suelo - REMEDHUS (Spain), Hydrological OBservatory and Exploratorium - HOBE (Denmark), Yanco (New South Wales, Australia), Little Washita (Oklahoma, USA), Texas Soil Observation Network - TxSON (Texas, USA) and Elm Creek (Manitoba, Canada). The sites cover a broad range of agronomic, hydrologic and climatic conditions, are instrumented with a network of calibrated ground stations continuously measuring surface soil moisture (0.05 m depth), and are routinely observed with S-1 data. The REMEDHUS, Elm Creek and HOBE data sets were collected from the International Soil Moisture Network (ISMN) (Dorigo et al., 2011), while Little Washita and Yanco data were downloaded from the Agricultural Research Service Micronet (Starks et al., 2014) (ars.mesonet.org), and the OzNet hydrological monitoring network (www.oznet.org.au), respectively; TxSON (www.beg.utexas.edu/research/programs/txson) data were shared through scientific collaborations. TxSON and Yanco use a nested design, replicating the soil moisture measurements at 3 km and 9 km inside their extent in support of the SMAP's Cal/Val Program (Caldwell et al., 2019; Yee et al., 2016). The Apulian Tavoliere network is hosted by Consiglio per la Ricerca in Agricoltura (CREA) in the experimental farm of Segezia and operated and maintained by the Italian National Research Council (CNR).

152

153

154

155

156

157

158

Table 1 summarizes the main features of the hydrological networks in terms of spatial extent, number of stations (S), spatial density (S/km²) and analyzed period. At the Apulian Tavoliere site, the Segezia experimental farm is mostly cropped with cereals (i.e., wheat, barley and oat). The area is quite homogeneous apart from one station deployed in a pasture field and one station located in an olive grove with cover crops, usually wheat. The 11 stations are situated in 10 S-1 Θ pixels of 520 m (i.e., the pixel spacing is approximately half the resolution) and cover a total area of ~1.6 × 1.6 km² (Fig. 1). This site features the only high-density network available in this study,
i.e., 4.3 stations/km². For this reason, it is considered as a core validation site as compared to the "low density" validation sites characterized by 1 station/km² (i.e. maximum 1 station is available for the comparison against the S-1 Θ at 1 km). The climate zone, the land cover classification, and the mean soil texture are also reported in Table 1. The agronomic, hydrologic and climatic conditions of the various sites range from semi-arid croplands/grasslands areas in Australia, Spain and Italy, to humid-subtropical grassland/shrubland areas in Texas, rangelands/croplands in Oklahoma and wet-all-seasons croplands in Denmark and Canada. Regarding the soil texture, there is an important sand component over the HOBE site as well as over REMEDHUS and Yanco. The soil texture with a 250 m pixel spacing over the experimental areas was obtained from the International Soil Reference Information Centre (ISRIC) (Hengl et al., 2017). The sites are mostly flat or gently undulating, with REMEDHUS and TxSON having the highest topography variability (mean and standard deviation of topography are shown in Table 1).

The extent of the validation dataset differs from site to site based on the number of the stations per network and the availability of the S-1 time series coverage. A minimum number of one hundred S-1 images per site was required. For the European and Australian sites, it was achieved in the time frame between 2015 and 2018; whereas for the remaining non-European sites the time series extend up to 2020. Quality control on the ground data consisted of excluding Θ values that showed anomalously low variations over the entire study period: four stations over REMEDHUS that recorded mean Θ values of ~0.026 m³/m³ and standard deviation 0.024 m³/m³; two stations over HOBE that recorded mean Θ values of ~0.524 m³/m³ and standard deviation ~0.048 m³/m³. Furthermore, Θ measures lower than 0.03 m³/m³ were removed from the data sets because this is the typical level of the calibration error of ground stations (Rowlandson et al., 2013). The number of measurements excluded is 191 out of 2389 for the REMEDHUS site, 198 out of 3515 for the Yanco site and 61 out of 2045 for the Little Washita. The maximum threshold for Θ was established at 0.60 m³/m³ as proposed in (Dorigo et al., 2013), and Θ values
184 > 0.60 $m^3/m^3$ (16 in total over the entire dataset) were excluded. Finally, $\Theta$ values measured during frozen soil conditions were excluded. The selection was carried out by the quality flag provided in the ISMN dataset, which identified the frozen soils using the soil temperature information. In particular, it resulted that the Canadian site was severely affected, and therefore the $\Theta$ values between October/November and March were discarded.

<table><caption>Table 1. Spatial extent (L), number of stations (S) available and spatial density (i.e., S/km²), Koeppen and Geiger climate classification (Rubel et al., 2017), land cover (LC), mean soil texture and topography variability (Digital Elevation Model (DEM) from Shuttle Radar Topography Mission, 30m) of the test sites. The climate labels are: BSk (arid, steppe, cold arid), Cfa (warm temperature, fully humid, hot summer), Dfb (Snow, fully humid, warm summer). For the soil texture and DEM both the mean and standard deviation at a site scale are reported.</caption><thead><tr><th>Site</th><th>Apulian Tavoliere (Italy)<br>(Balenzano et al., 2014)</th><th>Elm Creek (Canada)<br>(McNairn et al., 2015)</th><th>Little Washita (Oklahoma)<br>(Cosh et al., 2006)</th><th>HOBE (Denmark)<br>(Bircher et al., 2012)</th><th>REMEDHUS (Spain)<br>(Martinez F. and Ceballos, 2005)</th><th>TxSON (Texas)<br>(Caldwell et al., 2019)</th><th>Yanco (Australia)<br>(Smith et al., 2012)</th></tr></thead><tbody><tr><td>L [km²]</td><td>1.6 * 1.6</td><td>17 * 17</td><td>25 * 25</td><td>30 * 30</td><td>35 * 35</td><td>36 * 36</td><td>60 * 60</td></tr><tr><td>S</td><td>11</td><td>9</td><td>20</td><td>30</td><td>20</td><td>40</td><td>37</td></tr><tr><td>Spatial density [S/km²]</td><td>4.30</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.02</td><td>0.03</td><td>0.01</td></tr><tr><td>Climate zone</td><td>BSk</td><td>Dfb</td><td>Cfa</td><td>Dfb</td><td>BSk</td><td>Cfa</td><td>BSk</td></tr><tr><td>LC</td><td>Croplands</td><td>Croplands</td><td>Rangelands<br>Croplands</td><td>Croplands</td><td>Croplands</td><td>Grasslands<br>Shrublands</td><td>Croplands<br>Grasslands</td></tr><tr><td>Clay</td><td>10.32±2.84</td><td>29.59±7.59</td><td>17.97±1.65</td><td>8.14±1.85</td><td>19.65±2.06</td><td>25.76±2.25</td><td>35.32±4.70</td></tr><tr><td>Sand [%]</td><td>44.73±8.59</td><td>41.41±13.52</td><td>41.18±4.35</td><td>75.35±4.12</td><td>50.50±3.78</td><td>41.12±4.38</td><td>53.37±4.62</td></tr><tr><td>DEM [m]</td><td>152.3±12.6</td><td>229.1±6.5</td><td>366.9±27.9</td><td>99.9±20.0</td><td>821.2±61.2</td><td>503.0±58.4</td><td>140.1±10.2</td></tr><tr><td>period</td><td>Jan15-Dec18</td><td>Jun16-Dec20</td><td>Apr16-Dec20</td><td>Jan15-Dec18</td><td>Jan15-Dec18</td><td>Apr16-May20</td><td>Jan15-Dec18</td></tr></tbody></table>

Fig. 2 shows $\mu_{obs} = E(\Theta_{obs})$ (solid line bars), where $E(\cdot)$ is the spatio-temporal sample mean operator, and standard deviation, $\delta_{\Theta_{obs}}$, (dashed line bars) of the $\Theta$ distribution ($0.03 m^3/m^3 \le \Theta \le 0.60 m^3/m^3$) for each site (hereafter identified by AT=Apulian Tavoliere, EC=Elm Creek, T=TxSON, H= HOBE; Y=Yanco, LW=Little Washita, R=REMEDHUS). The yellow and black lines display the 75th (Q75) and 25th (Q25) percentiles of $\Theta$ values. The total number of the $\Theta$
measurements (M) over each test site available for the comparison after the quality selection is also reported in the legend. Although the sites were selected with different climatic conditions to cover as much as possible the Θ variability, the validation Θ data set is not equally-distributed over the Θ range, i.e., 75% of Θ measurements are approximately within 0.30 m³/m³. In particular, the lowest Q75 was observed over REMEDHUS and Yanco, i.e., approximately 0.20 m³/m³. This may be due to a combination of the semi-arid climate and the high sand component of soil texture, which reduces the water retention capacity of the soil (Montzka et al., 2018). Also for the fully-humid sites, such as HOBE, the high percentage of sandy soil likely explains the limited Θ values observed, as compared, for example, to Elm Creek. Additionally, HOBE shows the lowest δ<sub>Θobs</sub>, which is conversely the highest over Yanco.

<figure><img src="image_2.png" /></figure>

210

211

212

213

<figure><img src="image_3.png" /><figcaption>Fig. 1. Left panel: Location of the hydrological network (in red) at the Apulian Tavoliere site (Apulia region, Southern Italy). Right panel: distribution of the stations (red points) at the Segezia experimental farm (black line) over the S-1 soil moisture product grid with 520 m grid spacing (yellow lines).</figcaption></figure>

214

215

216

217

218

<caption>Fig. 2. Mean (solid line bars), μ<sub>obs</sub>, and standard deviation, δ<sub>Θobs</sub>, (dashed line bars) of the Θ distribution per site (AT=Apulian Tavoliere, EC=Elm Creek, T=TxSON, H= HOBE; Y=Yanco, LW=Little Washita, R=REMEDHUS). The yellow and black diamonds display the 75<sup>th</sup> (Q75) and 25<sup>th</sup> (Q25) percentiles. M is the total number of Θ measurements per site.</caption>
## 3 Sentinel-1 data collection

Time series of S-1 A & S-1 B Interferometric Wide (IW) Ground Range Detected (GRD) High Resolution (HR) data (Bourbigot et al., 2016) were collected. These data were pre-processed to obtain time series of calibrated, co-registered, geocoded and temporally filtered (Quegan and Yu, 2001) stacks of VV and VH backscatter coefficients at 40 m pixel size (roughly corresponding to ~100 m resolution) and with an equivalent number of looks (which is the ratio of the square of the backscatter spatial mean to the corresponding variance) ≈100.

Table 2 summarizes the S-1 acquisition pass and time, Relative Orbit Number (RON), number (N) of S-1 images and mean incidence angle (ϑ). The TxSON and Yanco sites were imaged at lower ϑ than Little Washita and Elm Creek. The Apulian Tavoliere site was covered by both descending and ascending S-1 tracks at different ϑ. To derive a robust calibration curve of the STCD algorithm over the entire S-1 swath (see Section 4.3.1), the S-1 tracks were selected to cover as much as possible the S-1 ϑ range, i.e., ~29° – 46°. It is noted that the S-1 time series over the non-European sites are characterized by 12-day revisit time, while for the European sites, S-1 A & S-1 B time series with 6-day revisit are available from the end-September 2016. However, a limited number of gaps in the S-1 time series occurred. Finally, it is worth noting that in case of Elm Creek, despite the availability of 84 S-1 images, only 64 were considered, excluding dates with frozen soils.

<table><caption>Table 2. Number (N) of the S-1 Interferometric Wide Swath (IW) Ground Range Detected (GRD) acquired over the experimental sites and for which also the ground data are available (AT=Apulian Tavoliere, EC=Elm Creek, T=TxSON, H= HOBE; Y=Yanco, LW=Little Washita, R=REMEDHUS). A, D, RON and UTC indicate the ascending or descending S-1 acquisition pass, Relative Number Orbit and the Coordinated Universal Time, respectively.</caption><thead><tr><th>Site</th><th>AT</th><th>T</th><th>Y</th><th>AT</th><th>H</th><th>R</th><th>LW</th><th>EC</th></tr></thead><tbody><tr><td>PASS<br>RON<br>(UTC)</td><td>A<br>146<br>(16:48:45)</td><td>A<br>107<br>(00:34:43)</td><td>D<br>118<br>(19:31:51)</td><td>D<br>124<br>(05:02:53)</td><td>D<br>139<br>(05:40:08)</td><td>D<br>154<br>(06:24:48)</td><td>A<br>107<br>(00:36:46)</td><td>A<br>136<br>(00:23:25)</td></tr><tr><td>N</td><td>183</td><td>97</td><td>112</td><td>173</td><td>165</td><td>161</td><td>113</td><td>64 out of<br>85</td></tr><tr><td>ϑ</td><td>33.4°</td><td>35.5°</td><td>37.2°</td><td>38.4°</td><td>38.8°</td><td>39.6°</td><td>41.4°</td><td>44.0°</td></tr></tbody></table>
242

# 4 Sentinel-1 Θ product

243

## 4.1 The Concept of the Algorithm

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

The implemented S-1 Θ retrieval algorithm transforms a dense or quasi-dense time series (i.e., 6- or 12-day revisit) of N dual-polarized S-1 IW images at 40 m pixel size (~100 m resolution) into N-Θ maps (Balenzano et al., 2011; Balenzano et al., 2013; Ouellette et al., 2017) at 520 m pixel size (~1 km resolution). The premise for the algorithm is that Θ changes take place at relatively short temporal scales (i.e., a few days or less), whereas changes associated with other surface parameters affecting the radar backscatter, such as soil roughness, canopy structure and vegetation biomass, are typically characterized by significantly longer temporal scales (e.g., a few weeks). Consequently, a SAR change detection approach (e.g., Rignot and Van Zyl, 1993) is expected to track changes in Θ only, since other parameters affecting the radar backscatter can be considered constant. Of course, the shorter the SAR revisit, the better the assumptions of the algorithm. This is a key difference over, for instance, the approach developed at the Vienna University of Technology (TU WIEN) (i.e., Bauer-Marschallinger et al., 2018, 2019), which in contrast requires a very long time series to estimate extreme values of Θ. The name of the algorithm - short term change detection (STCD) – reflects the importance of working on a dense time series of SAR data. To this regard, it is noted that Θ retrieval by the STCD algorithm will benefit from the launch of S-1 C, which will jointly operate with S-1 A & S-1 B (Torres et al., 2020) at least initially. Moreover, the development of new concepts of geostationary SAR platforms with a hyper-temporal resolution (e.g., Hobbs et al., 2019) will allow in the near future the full exploitation of the potential of time series retrieval approaches like STCD.

263

264

265

266

A second pillar of the algorithm is acknowledgment that the SAR signal at C-band does not always penetrate the vegetation layer, which is necessary to sense Θ. The lack of sensitivity to Θ is characteristic for mature dense forests, where the C-band SAR signal interacts primarily with the tree crown (i.e., volume scattering), whereas the signal scattered from the soil (i.e., either attenuated
surface scattering or double bounce) is not significant (Quegan et al., 2000). A distinctive radar feature of these targets is a high level of cross-polarized backscattering arising from the multiple reflections characteristic of volume scattering. In the case of agricultural or short vegetated areas (e.g. grassland and herbaceous cover), the interaction between the C-band radar signal and crops can significantly vary with the crop canopy structure and with the plant water content (i.e., fresh biomass). These canopy characteristics are strongly related to the phenological stage and ultimately to the plant development stage (Cookmartin et al., 2000; Khabbazan et al., 2019; Macelloni et al., 2001; McNairn and Brisco, 2004; Moran et al., 2012; Palmisano et al., 2020; Picard et al., 2003; Saich et al., 2000; Le Toan et al., 1997). Under these circumstances, it is clear that before performing a quantitative retrieval of Θ at SAR C-band, masking is required for those surfaces characterized by volume scattering.

The mathematical framework for the STCD algorithm is provided by a first-order approximation of the radiative transfer (RT) theory, which expresses the total backscatter of a vegetated surface as a superposition of three terms: the attenuated soil backscatter, the volume contribution and the soil-vegetation interaction (Tsang et al., 2001). In this context, the STCD algorithm adopts two main approximations:

• the first is that STCD applies only to bare or vegetated soils dominated by attenuated surface backscattering ($\sigma_0$), which at VV polarization can be expressed as in (1), i.e., volume scattering and soil-vegetation interaction are neglected

$$ \sigma_0 \approx \sigma_0^s \cdot \tau^2 = |\alpha_{VV}(\varepsilon, \vartheta)|^2 \cdot \Omega(\nu, \vartheta, \chi) \cdot \tau^2, \quad (1) $$

where $\tau^2$ is the two-way vegetation attenuation and $\sigma_0^s$ is the ground backscatter. The latter is written as the product of a term, $\alpha_{VV}(\varepsilon, \vartheta)$, which represents the influence of the surface permittivity ($\varepsilon$) and incidence angle ($\vartheta$) and a term, $\Omega(\nu, \vartheta, \chi)$, which represents the influence of the soil roughness, depending on the surface roughness power spectrum, $\chi(\cdot)$, the SAR frequency, $\nu$, and $\vartheta$. This
factorization is consistent with surface scattering models like the Small Perturbation Model and Small Slope Approximation (Voronovich, 1994) for which the expression of the reflection coefficient in VV polarization is

$$ |\alpha_{VV}(\varepsilon, \vartheta)| = \left| \frac{(\varepsilon-1)(\sin^2 \vartheta - \varepsilon(1+\sin^2 \vartheta))}{(\varepsilon \cos \vartheta + \sqrt{\varepsilon - \sin^2 \vartheta})^2} \right| \quad (2) $$

The second is that the backscatter ratio between two subsequent SAR acquisitions, at $DoY(i)$ and $DoY(i + 1)$, depends only on the ratio between the surface reflection coefficients of the two correspondent dates, such that

$$ \frac{(\sigma_0)_{DoY(i+1)}}{(\sigma_0)_{DoY(i)}} \approx \frac{(\sigma_0^s \cdot \tau^2)_{DoY(i+1)}}{(\sigma_0^s \cdot \tau^2)_{DoY(i)}} \approx \frac{|\alpha_{VV}(\varepsilon, \vartheta)|^2_{DoY(i+1)}}{|\alpha_{VV}(\varepsilon, \vartheta)|^2_{DoY(i)}}, \quad (3) $$

which requires that the roughness and vegetation parameters in (1) do not change between the two acquisition dates. The approximation in (3) was first proposed in Balenzano et al. (2011), under the name of the “alpha approximation”. The code implementing the STCD algorithm is referred to as SMOSAR (Soil Moisture retrieval from multi-temporal SAR data) (Balenzano et al., 2013). Fig. 3 is a schematic of the logic implemented in SMOSAR. The input is a time series of N (N=4) S-1 IW images at 40m pixel, which is firstly masked (using static land cover and dynamic S-1 VH observations). Then, it is transformed into Θ maps at 40m pixel, through the retrieval module (the soil texture maps are used to convert the dielectric constant into Θ). The final step is a low pass filter and resampling module that deliver Θ mean and standard deviation at 520m pixel (corresponding to half the spatial resolution of the Θ maps, which is ~1km). The retrieval module requires a calibration parameter, which can be updated as external information. The aforementioned processing chain is described in detail in Sections 4.2-4.4.

Finally, it is noted that the S-1 data acquired from different orbits are not handled by this version of the code, which means that in (3) the incidence angle is assumed not to change between two subsequent S-1 images. In this respect, an advanced version of SMOSAR accepting as input S-1 time series acquired from different orbits is under development.
<figure><img src="image_4.png" /><figcaption>Fig. 3. SMOSAR soil moisture (Θ) retrieval algorithm schema. Input and output data and main modules, i.e., masking and retrieval blocks, are drawn.</figcaption></figure>

**4.2 Location masking**

The need for limiting the S-1 Θ retrieval to surfaces over which the C-band SAR signal shows good sensitivity to Θ requires a masking process before applying the retrieval. This task is implemented as a two-step process in SMOSAR. The first uses the quasi-static ESA CCI land cover at 300 m spatial resolution (v2.0.7; Product User Guide, 2017) to mask forests, urban areas, water bodies and permanent snow and ice. The second step consists of a dynamic masking of the vegetation, which exploits a classification algorithm developed by Satalino et al. (2014). The method uses the S-1 VH observations to separate the radar response of seasonal crops into volume and soil attenuated scattering. The areas dominated by volume scattering are obscured, whereas those dominated by soil attenuated scattering are left unmasked. The rationale of the approach is that the higher the VH level, the higher the volume contribution. The detailed procedure implemented in SMOSAR is reported in Appendix B, as Supplementary material.
## 4.3 The retrieval of Θ and its variance

In SMOSAR, the variable that is initially retrieved is the absolute value of the alpha coefficient, $|α_{VV}(ε, ϑ)|$, which is subsequently inverted into the relative dielectric constant, $ε$, and then into Θ using the soil texture information. Mathematically at pixel scale and at a given date, the retrieved Θ content, $Θ_{retr}$, and its variance, $δ_{F}^2$, can be expressed by the Taylor uncertainties propagation

$$ Θ_{retr} = F(ε(α_{VV})) $$

(4)

$$ δ_F^2 = \left( \left( \frac{\partial F}{\partial ε} \right) (ε) \cdot \frac{\partial ε}{\partial α_{VV}} (α_{VV}) \right)^2 \cdot δ_{α_{VV}}^2 $$

(5)

where $F(·)$ is the function relating Θ to $ε$ (e.g., Hallikainen, 1985), $ε(α_{VV})$ is the analytical inversion of the reflection coefficient in (2) and $δ_{α_{VV}}^2$ is the variance of the random variable $|α_{VV}(ε, ϑ)|$.

### 4.3.1 The maximum likelihood estimator of $|α_{VV}(ε, ϑ)|$

The retrieval method described in Balenzano et al. (2013) is here summarized. It is a time series approach that applies the approximation in (3) to $N$ subsequent dates. For the sake of simplicity, the square roots of the quantities in (3) are considered and the following ratio defined as $\hat{S}_{ij} = \sqrt{(\sigma_0)_{DoY(i)}/(\sigma_0)_{DoY(j)}}$. It is therefore possible to write a linear underdetermined stochastic system of $N-1$ equations with $N$ unknowns, $|α_{VV}(ε, ϑ)|$:

$$ \bar{\bar{A}} \vec{\alpha}_{VV} = \begin{bmatrix} -\hat{S}_{21} & 1 & 0 & 0 & \dots & 0 & 0 \\ 0 & -\hat{S}_{32} & 1 & 0 & \dots & 0 & 0 \\ 0 & 0 & -\hat{S}_{43} & 1 & \dots & 0 & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & 0 & 0 & \dots & -\hat{S}_{N(N-1)} & 1 \end{bmatrix} \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \dots \\ \alpha_N \end{bmatrix} = \vec{0} \quad (6) $$
where $\bar{\bar{A}}$ is a full row rank ($N-1) \cdot N$ matrix, $\vec{\alpha}_{VV} = [\alpha_1, \alpha_2, ..., \alpha_N]$ is a $N$-dimensional vector and $\vec{0}$ is the $(N-1)$-dimensional null vector. The solution of the system can be expressed as a function of a free parameter ($0 < \lambda < \infty$)

$$ \vec{\alpha}_{VV} = \hat{\lambda} \cdot [\hat{S}_{1N}, \hat{S}_{2N}, ..., \hat{S}_{NN}]. \qquad (7) $$

A set of linear constraints can be applied to the solution of the underdetermined system

$$ 0 < \alpha_{min} \leq |\alpha_{VV}(\varepsilon, \vartheta)|_{DoY(i)} \leq \alpha_{max} < \infty \quad i = 1 \dots N, \qquad (8) $$

and noting that for a number of looks of the S-1 images much larger than 10, as is always the case in this study, the probability density function (pdf) of each $\hat{S}_{iN}$ element in (7) can be approximated well by a normal distribution (Lee et al., 1994; Oliver and Quegan, 1998), the maximum likelihood (ML) criterion provides an optimal estimate for the $\lambda$ parameter as

$$ \hat{\lambda} = \max \left( \frac{\alpha_{min}}{\hat{S}_{iN}} \right) \quad i = 1, \dots, N. \qquad (9) $$

It is noted that the ML solution in this case is the minimum norm solution which is equal to the least square solution (LS) (e.g., Ouellette et al., 2017). In summary, the algorithm transforms time series of S-1 observations into time series of $\Theta$ absolute values, under the assumption that additional information is available in terms of an estimate of the minimum value of the alpha coefficient, $|\alpha_{VV}(\varepsilon, \vartheta)|$ during the $N$ S-1 acquisitions.

There are various options to ascertain the value of $\alpha_{min}$ required in (9). For instance, an estimate of $\alpha_{min}$ at low resolution (e.g., ~40 km) can be obtained from $\Theta$ operational products, e.g., SMOS, SMAP, ASCAT, etc. Such an option was implemented in Ouellette et al. (2017) and Al-Khaldi et al. (2019). Another option is to use in situ data to drive the retrieval (e.g., Palmisano et al., 2020). In this study, a calibration curve expressing S-1 VV observations versus $|\alpha_{VV}(\varepsilon, \vartheta)|^2$ values at low resolution was adopted. The rationale is that the spatial average at coarse scale reduces the influence of surface parameters characterized by a high spatial frequency (~0.1 km), such as roughness, crop
canopy structure and vegetation water content, while strengthening the relationship with Θ, which adjusts steadily in space (Macelloni et al., 1999). The calibration curve was first implemented in an ESA feasibility study (Mattia et al., 2011), and improved by using a subset of 1/3 of the total couples available in Table 2. A data set of S-1 VV observations and |α<sub>VV</sub>(ε,ϑ)|<sup>2</sup> observed at site scales was built. The Apulian Tavoliere data set was not used to identify the calibration curve, because of its limited extent. S-1 observations were expressed in terms of the γ coefficient (i.e., γ = σ<sub>0</sub>/cos(ϑ)) rather than σ<sub>0</sub> in order to mitigate the effect of difference in ϑ from site to site (Table 2). The S-1 observations were first masked, as described in Section 4.2, in order to select the area dominated by the soil attenuated scattering (1), then γ was estimated over the unmasked areas. The Θ measurements of the hydrological networks were temporally collocated with the S-1 acquisitions and averaged at the site scale. The |α<sub>VV</sub>(ε,ϑ)|<sup>2</sup> coefficients were derived from the mean Θ values considering the mean ϑ (Table 2) and the soil texture (Table 1) specific for each test site. Finally, the linear relationship |α<sub>VV</sub>(ε,ϑ)|<sup>2</sup> versus γ was identified at coarse scale and used to derive α<sub>min</sub> during the Θ retrieval process. The calibration curve was applied for the Θ retrieval over all validation sites in Table 1, and it is expected to improve with time as new S-1 observations are integrated. For this reason, in the SMOSAR algorithm there is an option to update the parameters of the calibration.

It is also noted that even though $\alpha_{min}$ in (9) is derived at a coarse scale, the time series approach enables resolving at high resolution the Θ fields undergoing a different temporal evolution. This is because the maximum condition in (9) is enforced at the pixel scale and, therefore, the indexing of $max(\frac{1}{\hat{S}_{iN}})$ in (7) changes with the local temporal evolution of the backscatter.

### 4.3.2 The Θ retrieval error

The expression of the retrieved variance $\delta_F^2$ in (5) depends on two main terms: i) the partial derivatives $\frac{\partial F(\epsilon)}{\partial \epsilon}$ and $\frac{\partial \epsilon}{\partial \alpha_{VV}}(\alpha_{VV})$ and ii) the variance of the reflection coefficient, $\delta_{\alpha_{VV}}^2$. The
computation of the partial derivative can be carried out analytically (the computation is straightforward, though quite lengthy and tedious and therefore not reported here). The estimate of $\delta_{\alpha_{VV}}^2$ is detailed in the Appendix A. The result (see eq. (A12)) indicates an upper bound that can be expressed as

$$ \delta_{\alpha_{VV}}^2 \le [ \delta_{stat}^2 + \delta_{cal}^2 ], \qquad (10) $$

where $\delta_{stat}^2$ and $\delta_{cal}^2$ arise from the propagation of the measurement error affecting the terms $\hat{S}_{iN}$ ($i = 1, \dots, N$) in (7) and the error affecting the estimate of the parameter $\alpha_{min}$ in (9). In Appendix A, the analytical expressions of these two terms are derived and discussed. Here, it is noted that on top of the two aforementioned error sources, a third contribution, $\delta_{mod}^2$, accounting for possible failures of the approximations reported in (1) and (3), needs to be included in the error budget. For instance, (3) assumes that between two subsequent S-1 acquisitions the only surface parameter that changes is $\Theta$. In reality, there is always a certain probability that the roughness and/or vegetation parameters also change and this probability is expected to increase with the revisit time. Therefore, the total variance can be expressed as

$$ \delta_{retr}^2 = \delta_{F}^2 + \delta_{mod}^2 = \delta_{stat}^2 + \delta_{cal}^2 + \delta_{mod}^2. \qquad (11) $$

In principle, the term $\delta_{mod}^2$, called the model error, can be characterized experimentally in those cases where $\delta_{stat}^2 \approx \delta_{cal}^2 \approx 0$. In many circumstances, the $\delta_{mod}^2$ term can be the dominant contribution in (11). In this respect, a thorough investigation of the impact of the S-1 time revisit on the RMSE and R will be conducted in a future study, in which constant observation conditions but the time revisit are kept.

## 4.4 The output product

Once the coefficient $|\alpha_{VV}(\varepsilon, \vartheta)|$ on each date is retrieved, $\varepsilon$ can be analytically derived and then $\Theta$ estimated by inverting the Hallikainen et al. (1985) empirical model. To this regard, SMOSAR includes the global gridded ISRIC soil texture at 250 m spatial resolution (Hengl et al., 2017).
The last step in SMOSAR is a low pass filter, with a kernel of $W \times W$ pixels (W=13), applied to the $\Theta$ maps at 40 m pixel size. The advantage is twofold. First, the uncertainty on the $\Theta$ retrieved is reduced and, second, the impact of errors due to abrupt changes of vegetation and/or soil roughness, which normally take place at field scales and can be wrongly interpreted as $\Theta$ changes, are mitigated. The mitigation would probably increase when averaging over larger areas. However, 1 km resolution is a tradeoff between the need of reducing the presence of biases and preserving a high resolution in the final $\Theta$ product. A side effect of the masking process is that a number of pixels in each 1 km$^2$ resolution cell could be masked and therefore null. To handle this aspect, the adopted rule is that if the ratio between the remaining pixels over the total pixels is less than 33%, than the $\Theta$ estimate for that kernel window is set to null. Finally, a resampling of the $\Theta$ map is performed (using the same resampling parameter W), changing the pixel size from 40 m to 520 m, which corresponds to a spatial resolution of approximately 1 km. The standard deviation associated with the mean $\Theta$ value at 1 km resolution is also estimated and delivered as a companion layer.

As an example, Fig. 4 (upper panel) shows a prototype of a 1 km S-1 $\Theta$ 6-day composite (24 overlapping S-1 descending tracks composed of 165 frames) over the Mediterranean basin as well as the 1 km S-1 $\Theta$ standard deviation 6-day composite (lower panel). The main river basins are delineated by black lines. The $\Theta$ patterns related to precipitation (light to dark blue) and drying of the soil (yellow to orange and then red) are visible. The $\Theta$ uncertainty is generally high over areas with a contrasting topography and very wet conditions. The S-1 $\Theta$ product prototype at the Mediterranean scale has been generated for one year (December 2017 – 2018) and its extension, systematic production and distribution (WebGIS Service) is currently under study.

### 4.4.1 The sliding window processing

SMOSAR processes the S-1 time series in a continuous chain using a sliding-window of four S-1 images each time (N = 4). As soon as a new S-1 image is acquired (e.g. image $k^{th}$ on date(k)), it is processed together with the previous (N - 1) images. As a result, a time series of N-Fast Delivered
(FD) $\Theta$ maps referring to [date($k$ - (N - 1)),..., date($k$)] is obtained (horizontal black box at time $T$ in Fig. 5). When subsequent S-1 images are processed (from time $T+1$ up to $T+N-1$), multiple estimates of $\Theta$ maps are obtained on the same date($k$) (i.e., those on the same column in Fig. 5). The FD $\Theta$ maps on the same date are partly correlated and can be averaged to reduce their total variability. The averaged output is called the Precision $\Theta$ product. It is noted that the temporal standard deviation can be an indicator of the extent to which the vegetation and/or surface roughness were stationary during the 4 subsequent S-1 acquisitions. This indicator is, however, not analyzed in the present study.

In the event that a gap in S-1 acquisitions occurs, a new processing chain is started. The first and last S-1 images of the processed time series are prone to larger errors because there are no multiple $\Theta$ estimates on the same date to be averaged. The Precision $\Theta$ product is the product validated in the following Sections.

1km S-1 Θ 6-day composite

<figure><img src="image_5.png" /></figure>

1km S-1 Θ standard deviation 6-day composite

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 5 Multiple Fast Delivered (FD) Θ maps produced by the processing of N S-1 data applied continuously to the pipeline of S-1 images (N=4) and Precision Θ product derived by averaging the FD images for the same date.</figcaption></figure>

# 5 Methods for local validation

The analysis focused on the comparison of time series of retrieved (Θ<sub>retr</sub>) versus observed Θ values (Θ<sub>obs</sub>) collected over the validation sites. The comparison was performed at two spatial scales: 1 km and the network scale. To match the spatial and temporal time series of retrieved and in situ measured Θ, the ground stations closest to the centroids of the S-1 Θ retrieved grid were selected and mean Θ measurements within 1 hour (from 1 to 4 samples according to the temporal sampling of the hydrologic networks) before the S-1 acquisition time were considered.

The metrics selected for Θ validation are standard statistical scores, such as bias, root mean square error (RMSE), unbiased RMSE (ubRMSE), Pearson correlation coefficient (R) and its significance (p), standard deviation (δ), and the parameters of the linear regression, i.e., slope (β) and intercept, of retrieved versus observed Θ (Entekhabi et al., 2010). However, the interpretation of these metrics should be approached with caution when independent (Θ<sub>obs</sub>) and dependent (Θ<sub>retr</sub>)
variables are affected by similar error levels, as it is the case in this study. Indeed, the measurement errors for $\Theta_{obs}$ reduce the magnitude of the observed correlation between the independent and dependent variables, bias the estimate of the slope towards zero (Kelly, 2007) and increase the RMSE (Dorigo et al., 2015). Under these circumstances, the Weighted Least Square (WLS) method (York et al., 2004) is most suited to evaluate the linear regression parameters and $R$, instead of the Ordinary Least Square (OLS). The WLS code implemented in (Thirumalai et al., 2011) was used in the analysis. It is noted that WLS also accounts for the presence of heteroscedastic errors (unequal variability across the $\Theta$ range) both in the dependent and independent variables (Cantrell, 2008; Thirumalai et al., 2011). The rationale of WLS is that data with the least errors have the greatest influence on the slope and intercept of the fitted line (i.e., weights are proportional to the inverse of the variance of the data values). Therefore, the analysis of the error sources and the quantification of the various contributions is crucial for the implementation of the WLS. For the independent variable, the most important source of error is the SRE. In preparation for the validation activity of the SMAP mission various experimental sites have been identified, or established or upgraded to deploy an appropriate number of ground stations to ensure an accurate estimate of the average $\Theta$ (Colliander et al., 2017). This analysis, however, depends on the resolution of the EO system. For this reason, a similar effort should be undertaken at higher resolution, e.g. ~1 km. In this study, an experimental analysis was carried out over the Apulian Tavoliere (Italy) core site, as well as a modelling analysis to provide to SRE the appropriate weights across a number of spatial scales.

510

## 5.1 Spatial representativeness error

Each technique measuring $\Theta$ is characterized by its own “support”, which is the effective area that each measurement represents (Western and Blöschl, 1999). In this respect, the support of in situ observations is, in general, much smaller than that of satellite $\Theta$ retrieved products. To allow meaningful comparisons it is necessary to sample and then average a number ($S$) of independent $\Theta_{obs}$ over an area comparable to the resolution cell of the satellite product. The spatial

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510
representativeness error (SRE) is the margin of error in estimating the mean Θ value of that area, at a specific confidence level (CL) and using S independent point-scale Θ<sub>obs</sub> observations. Its mathematical expression is:

$$
\begin{equation}
\begin{array}{@{}l@{}}
  \tilde{\delta}_{SRE} = z_{\alpha/2} \cdot \tilde{\delta}_{\Theta_{obs}} / \sqrt{S} \\
  \quad (12)
\end{array}
\end{equation}
$$

It is noted that Θ is a multiscale, heteroscedastic process (Das et al., 2010; McColl et al., 2017; Western and Blöschl, 1999), so that $\tilde{\delta}_{\Theta_{obs}}$ is a function both of the sampling scale (L) and $\tilde{\mu}_{obs}$, i.e., $\tilde{\delta}_{\Theta_{obs}} = L(L, \tilde{\mu}_{obs})$. To estimate it, Gilbert, (1987) recommends characterizing first its coefficient of variation ($CV_L$), then $\tilde{\delta}_{\Theta_{obs}}$ can be obtained as a product, i.e., $CV_L \cdot \tilde{\mu}_{obs}$. This is beneficial because $CV_L$ usually shows less variability than $\tilde{\delta}_{\Theta_{obs}}$, and its experimental relationship with $\tilde{\mu}_{obs}$ can be fitted using an exponential law depending on two parameters, $k_1$ and $k_2$ such that

$$
\begin{equation}
\begin{array}{@{}l@{}}
  CV_L = \frac{\tilde{\delta}_{\Theta_{obs}}}{\tilde{\mu}_{obs}} = k_1 \cdot e^{-k_2 \tilde{\mu}_{obs}}, \\
  \quad (13)
\end{array}
\end{equation}
$$
$\tilde{\mu}_{obs} \cdot e^{-k_2 \tilde{\mu}_{obs}}$ is that the coordinates of their maximum, i.e., $(\tilde{\mu}_{obs}^{max}, \tilde{\delta}_{\Theta_{obs}}^{max})$, are simply expressed in terms of $k_1$ and $k_2$ parameters:

$$ \left\{ \begin{array}{l} \tilde{\mu}_{obs}^{max} = 1/k_2 \\ \tilde{\delta}_{\Theta_{obs}}^{max} = k_1/(ek_2) \end{array} \right. , \qquad (14) $$

where $e$ is Euler's number. This feature can be further exploited to predict the dependence of $\tilde{\delta}_{\Theta_{obs}}$ on $\tilde{\mu}_{obs}$ at any extent scale $L$, as will be shown in the next Section.

### 5.1.1 Scaling of $\Theta$ variability at 1km and SRE quantification

The multi-scale nature of $\Theta_{obs}$ suggests that the relation between its sample mean and standard deviation, i.e., $\tilde{\mu}_{obs}$ and $\tilde{\delta}_{\Theta_{obs}}$, and the extent scale $L$ is approximated well by a power-law. In (Famiglietti et al., 2008), this property was assessed in a Log - Log plot between $\tilde{\delta}_{\Theta_{obs}}$ and the extent scale, $L$. Since for each scale, $\tilde{\delta}_{\Theta_{obs}}$ depends on $\tilde{\mu}_{obs}$, as shown in Fig. C1 (in the Supplementary material), Famiglietti et al. (2008) plotted the mean standard deviation at each scale and suggested that a fit of the plot could be used to estimate the mean variance at a particular scale. Following this line of reasoning, the maximum value of $\tilde{\delta}_{\Theta_{obs}}$ at each scale is considered in this study, rather than its mean value. In other words, a power law relation between $(\tilde{\mu}_{obs}^{max}, \tilde{\delta}_{\Theta_{obs}}^{max})$ and $L$ was sought such that

$$ \left\{ \begin{array}{l} \tilde{\mu}_{obs}^{max} = \alpha_1 \cdot L^{\beta_1} \\ \tilde{\delta}_{\Theta_{obs}}^{max} = \alpha_2 \cdot L^{\beta_2} \end{array} \right. , \qquad (15) $$

where $\alpha_1$ ($\alpha_2$) and $\beta_1$ ($\beta_2$) are fitting parameters, the latter related to the fractal dimension (i.e., Mandelbrot and Van Ness, 1968). For each scale (site) reported in Table C1 (in the Supplementary material), the couples $(\tilde{\mu}_{obs}^{max}, \tilde{\delta}_{\Theta_{obs}}^{max})$ were estimated. To increase the number of the fitting points the curves reporting $\tilde{\delta}_{\Theta_{obs}}$ vs $\tilde{\mu}_{obs}$ in Famiglietti et al., (2008) and Jacobs et al., (2004) were also included. Additionally, the $\Theta$ variability at the 3 km and 9 km nested grids of the Yanco and TxSON networks was also considered. The total number of the fitting points is 24. Fig. 6 (left panel) shows the $\text{Log}_{10}(\tilde{\mu}_{obs}^{max})$ (orange points) and the $\text{Log}_{10}(\tilde{\delta}_{\Theta_{obs}}^{max})$ (green points) versus the
Log₁₀(L) respectively. The fitting parameters α₁, β₁ and α₂, β₂ are reported in the caption. The vertical dotted line identifies the Log₁₀(δ̃<sub>Θ<sub>obs</sub></sub><sup>max</sup>) and Log₁₀(μ̃<sub>obs</sub><sup>max</sup>) at 1 km. The coefficient of determination for the fit of μ̃<sub>obs</sub><sup>max</sup> and δ̃<sub>Θ<sub>obs</sub></sub><sup>max</sup> are ~0.6 and ~0.8, respectively. In both cases, the correlation is highly significant confirming not only that most of the variability of Θ<sub>obs</sub> over the various sites can be explained with the extent of the site, but also that it is possible to predict the variability of Θ<sub>obs</sub> at a particular scale. For instance, using (15) (μ̃<sub>obs</sub><sup>max</sup>, δ̃<sub>Θ<sub>obs</sub></sub><sup>max</sup>) can be estimated at L = 1 km and then using (14), it is possible to derive the associated δ̃<sub>Θ<sub>obs</sub></sub> = CV<sub>L</sub> · μ̃<sub>obs</sub>. In particular, for L = 1 km: k₁ = 0.686 and k₂ = 4.328(m³/m³)⁻¹.

Moreover, using (12), the behaviour of δ̃<sub>SRE</sub> as a function of μ̃<sub>obs</sub> and for any S can be estimated. Fig. 6 (right panel) shows δ̃<sub>SRE</sub> at 1 km scale, 70% CL, for S=1 and 4. For S=4, δ̃<sub>SRE</sub> is within the typical calibration error of ground station probes, i.e., ~0.03 m³/m³. For S=1, δ̃<sub>SRE</sub> is significantly higher, hence in comparing retrieved Θ values at 1 km with that observed by a single station, the δ̃<sub>SRE</sub> cannot be disregarded. The outcome of this Section is that in the absence of an adequate number of stations, a viable alternative is to predict δ̃<sub>SRE</sub> and take this into account in the validation metrics, as shown in the next Section.

<figure><img src="image_8.png" /></figure>

583
Fig. 6. Left panel: Log$_{10}$(δ̃<sub>obs</sub><sup>max</sup>) vs Log$_{10}$(L) (green squares). The fitting parameters are α₂ = 0.023 m³/m³; β₂ = 0.132, R² = 0.81, p < 0.01. Log$_{10}$(μ̃<sub>obs</sub><sup>max</sup>) vs Log$_{10}$(L) (orange circles). The fitting parameters are α₁ = 0.085 m³/m³; β₁ = 0.144, R² = 0.60, p < 0.01 (N=24). The vertical line indicates the Log$_{10}$(μ̃<sub>obs</sub><sup>max</sup>) and Log$_{10}$(δ̃<sub>obs</sub><sup>max</sup>) at 1 km. Right panel: spatial representativeness error (δ̃<sub>SRE</sub>) as a function of μ̃<sub>obs</sub> at 70% CL, at 1 km scale and S=1 station (blue line) and S=4 stations (red line).

## 5.2 Error budget

According to (11), the total Θ retrieval error (δ<sub>retr</sub><sup>2</sup>) can be split into three main terms, namely:
the statistical (δ<sub>stat</sub><sup>2</sup>), the calibration (δ<sub>cal</sub><sup>2</sup>) and the model (δ<sub>mod</sub><sup>2</sup>) error. Conversely, the error
affecting the Θ observations (δ<sub>obs</sub><sup>2</sup>) consists of two independent contributions: the first due to the
SRE (δ<sub>SRE</sub><sup>2</sup>) and the second due to the sensor calibration error δ<sub>sensor</sub><sup>2</sup>. Therefore, it can be written
as

$$ \delta_{\text{obs}}^2 = \delta_{\text{SRE}}^2 + \delta_{\text{sensor}}^2, \qquad (16) $$

where in most cases δ<sub>sensor</sub><sup>2</sup> can be disregarded compared to δ<sub>SRE</sub><sup>2</sup> because its RMSE level is
~0.03 m³/m³ (Rowlandson et al., 2013) or even lower (Coopersmith et al., 2016). As a result, the
RMSE is:

$$ \begin{align*} 600 \quad rmse^2 &= E[(\theta_{retr} - \theta_{obs})^2] = E[(\theta'_{retr} + \epsilon_{retr}) - (\theta'_{obs} + \epsilon_{obs})]^2 = \\ 601 \quad &= E[(\theta'_{retr} - \theta'_{obs})^2] + \delta_{retr}^2 + \delta_{obs}^2 = rmse_{intr}^2 + \delta_{obs}^2 \approx rmse_{retr}^2 + \delta_{SRE}^2, \quad (17) \\ 602 \quad &\phantom{=} E[(\theta'_{retr} - \theta'_{obs})^2] + \delta_{retr}^2 + \delta_{obs}^2 = rmse_{retr}^2 + \delta_{SRE}^2 \quad (17) \\ 603 \quad &\phantom{=} E[(\theta'_{retr} - \theta'_{obs})^2] + \delta_{retr}^2 + \delta_{obs}^2 = rmse_{retr}^2 + \delta_{SRE}^2 \quad (17) \\ 604 \quad &\phantom{=} E[(\theta'_{retr} - \theta'_{obs})^2] + \delta_{retr}^2 + \delta_{obs}^2 = rmse_{retr}^2 + \delta_{SRE}^2 \quad (17) \\ 605 \quad &\phantom{=} E[(\theta'_{retr} - \theta'_{obs})^2] + \delta_{retr}^2 + \delta_{obs}^2 = rmse_{retr}^2 + \delta_{SRE}^2 \quad (17) \\ 606 \quad &\phantom{=} E[(\theta'_{retr} - \theta'_{obs})^2] + \delta_{retr}^2 + \delta_{obs}^2 = rmse_{retr}^2 + \delta_{SRE}^2 \quad (17) \\ 607 \quad &\phantom{=} E[(\theta'_{retr} - \theta'_{obs})^2] + \delta_{retr}^2 + \delta_{obs}^2 = rmse_{retr}^2 + \delta_{SRE}^2 \quad (17) \\ 608 \quad &\phantom{=} E[(\theta'_{retr} - \theta'_{obs})^2] + \delta_{retr}^2 + \delta_{obs}^2 = rmse_{retr}^2 + \delta_{SRE}^2 \quad (17) \\ 609 \quad &\phantom{=} E[(\theta'_{retr} - \theta'_{obs})^2] + \delta_{retr}^2 + \delta_{obs}^2 = rmse_{retr}^2 + \delta_{SRE}^2 \quad (17) \\ 610 \quad &\phantom{=} E[(\theta'_{retr} - \theta'_{obs})^2] + \delta_{retr}^2 + \delta_{obs}^2 = rmse_{retr}^2 + \delta_{SRE}^2 \quad (17) \\ 611 \quad &\phantom{=} E[(\theta'_{retr} - \theta'_{obs})^2] + \delta_{retr}^2 + \delta_{obs}^2 = rmse_{retr}^2 + \delta_{SRE}^2 \quad (17) \end{align*} $$

where ε<sub>retr</sub> and ε<sub>obs</sub> are zero-mean random errors, with variance δ<sub>retr</sub><sup>2</sup> = E(ε<sub>retr</sub><sup>2</sup>) and δ<sub>SRE</sub><sup>2</sup> = E(ε<sub>SRE</sub><sup>2</sup>), while θ'<sub>retr</sub> and θ'<sub>obs</sub> are two random variables with their intrinsic variability. Eq. (17) assumes that error cross-correlation and error orthogonality (i.e., homoscedasticity) components are zero or perfectly compensating (Gruber et al., 2020). Disregarding δ<sub>SRE</sub><sup>2</sup>, the RMSE between θ<sub>retr</sub> and θ<sub>obs</sub> is the intrinsic mean square error, RMSE<sub>intr</sub><sup>2</sup>.

For sites with a low density of ground stations at ~1 km resolution, a workable approach is to use the procedure outlined in Section 5.1 to predict δ̃<sub>SRE</sub> in (12) on each date for an arbitrary number of stations, S. Then δ<sub>SRE</sub><sup>2</sup> in (17) can be approximated by the temporal mean of the retrieved
612 $\tilde{\delta}_{SRE}$. Finally, the RMSE$_{intr}$ can be retrieved using (17). In the case of high-density sites,
613 RMSE$^2 \approx$ RMSE$^2_{intr}$, as $\delta^2_{SRE}$ is expected to be negligible. For biased estimates of $\Theta_{retr}$, i.e.,
614 $\mu_{obs} - \mu_{retr} = bias$, the relation between $ubRMSE^2$ and RMSE$^2$ is (Entekhabi et al., 2010):
615 $ubRMSE^2 = RMSE^2 - bias^2$.

(18)

## 6 Validation results

### 6.1 Validation at 1km

619 The performance of the STCD algorithm was investigated under two different experimental
620 conditions:

621 • over the core site, i.e., the Apulian Tavoliere, where the error on the independent variable
622 ($\Theta_{obs}$) is much lower than the retrieval error, i.e., $\delta^2_{obs} \ll \delta^2_{retr}$

623 • over the low-density sites, where $\delta^2_{obs} \sim \delta^2_{retr}$ and therefore, it is necessary to estimate the
624 $\delta_{obs}$ and evaluate RMSE$_{intr}$ through (17).

#### 6.1.1 Validation over the Apulian Tavoliere core site

625 The relatively high density of ground stations deployed on this site (i.e., ~4 stations/km$^2$)
626 enables the use of the standard metrics to characterize the performance of the STCD algorithm. This
627 site provides the opportunity to investigate experimentally the dependence of the RMSE on the
628 number of stations (from 1 to 11) used to estimate the spatial mean $\Theta$, i.e., $\tilde{\mu}_{obs}$. To investigate this,
629 first the time series of $\Theta_{retr}$ values averaged at site scale (i.e., 1.6 km × 1.6 km; see Fig. 1),
630 i.e. $\tilde{\mu}_{retr}$, was derived using the S-1 time series of ascending A146 orbit (Table 2). Next, the ground
631 stations were randomly aggregated in 11 non-overlapping groups (G), including S = 1, 2, 3, 4 ... 11
632 stations as illustrated in Fig. 7 (left panel). The aggregated stations are shown in the same colour.
633 For each group and on each date, the correspondent $\Theta_{obs}$ values were averaged and compared
634
against $\tilde{\mu}_{retr}$. For those groups with the same number of stations $S$, the related RMSEs were evaluated and averaged and the standard deviation was calculated. In Fig. 7 (right panel), the dotted line shows the mean RMSE and standard deviation (up to $S = 5$) versus the number of the stations.

The higher the value of $S$, the lower the mean RMSE as well as the associated standard deviation.

The mean RMSE decreases from $0.078 \text{ m}^3/\text{m}^3$ for 1 station to the minimum RMSE, $0.064 \text{ m}^3/\text{m}^3$, for $S = 11$. However, below ~0.06 $\text{m}^3/\text{m}^3$, the RMSE is minimally reduced, likely because the $\delta_{SRE}$ reaches a value lower than $\delta_{sensor}$ (see Fig. 6, right panel). Likewise, the observed threshold level of RMSE, i.e., ~0.06 $\text{m}^3/\text{m}^3$ can be considered as a good estimate of the STCD retrieval error.

An overview of the STCD algorithm performance is displayed in Fig. 8 (left panel), which reports the data of Fig. 7 in the form of a scatter plot, for $S = 11$. Three outliers (residual errors > 3 standard deviations (Komorowski et al., 2016)) are reported in black circles. A further aspect that was investigated is the variability of the RMSE as a function of $\Theta$ ranges, when $\delta_{SRE}$ is minimized.

Fig. 8 (right panel) illustrates the distribution of RMSE (green bars) and ubRMSE (blue bars) per $\Theta$ interval. To balance the number of points per $\Theta$ interval, $\Theta_{obs}$ were aggregated in uneven intervals ($\text{m}^3/\text{m}^3$), i.e., [0.05,0.20], [0.20,0.30] and [0.30,0.55]. The (ub)RMSE increases as a function of $\Theta$ and its minimum value in the first interval is $\approx 0.05 \text{ m}^3/\text{m}^3$. The relative error is ~20% of the mean of the two subsequent intervals. Under these circumstances, the adoption of a relative error for the requirements on SAR $\Theta$ retrieval products seems more realistic than an absolute threshold.

The Appendix A further elaborates on why the (ub)RMSE increases with $\Theta$.
<figure><img src="image_9.png" /><figcaption>Fig. 7. Left panel: Number of stations (S) per each group (G). Stations S with the same colour belong to the same group G. Right panel: RMSE between ascending 1.6 km*1.6 km S-1 Θ and Θ measured by 1 station or averaged from 2 up to 11 stations as a function of the number of the stations within the Apulian Tavoliere core test site.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Left panel: Scatter plot (Dates=183) between Θ derived from the S-1 ascending track (A146) and the Θ values averaged over the 11 stations at the Apulian Tavoliere site (1.6 km *1.6 km). The Ordinary Least Square (OLS) fit (in black), as well as the statistical scores, are reported. Three outliers are in black circles. Right panel: Distribution of RMSE (green bars) and ubRMSE (blue histogram bars), as defined in (18), per Θ interval (without the three outliers).</figcaption></figure>

Fig. 9 compares the temporal behavior of $\tilde{\mu}_{retr}$ (red line) and $\tilde{\mu}_{obs}$ (blue line) over the Apulian Tavoliere site. The dashed areas represent the daily spatial standard deviations. The daily precipitation available from a meteo station 10 km distant from the site is also reported. Few comments are in order:
• In general, the level and the temporal evolution of the continuous red and blue lines are in good agreement (R=0.71), in particular $\Theta_{retr}$ captures fairly well the wetting and dry-downs observed in situ. Nevertheless, the retrieved standard deviation is important for medium-high values of $\Theta$, while it is fairly small in dry conditions. This is in agreement with the prediction of the statistical error of STCD reported in Appendix A (see Fig. A1). As a consequence, the RMSE in the spatial domain can be highly affected by biases for wet rather than for dry surfaces.

• Some nuisance in the $\Theta$ retrieval performance can be observed when the STCD physical approximations (e.g., constant soil roughness and vegetation between two S-1 acquisitions) are not fulfilled. This is the case of the three outliers reported in Fig. 9, as red points. For instance, between July and August 2015 various fields of the farm were ploughed and then arrowed (i.e. smoothed). The change in soil roughness and consequent drastic change of backscatter between August, 11 and 23 was interpreted by SMOSAR as due to a change of $\Theta_{retr}$, thus leading to the $\Theta$ overestimation. Similarly, the $\Theta$ underestimation in April is likely related to the rapid growth of wheat canopy ("stem elongation" phase) which characterized several fields of the farm. Conversely, the two outliers occurred on September 29 and October, 23 2017 are probably due to a high spatial/temporal variability of the precipitation fields in the area. Indeed, although precipitation events were recorded by the meteo station and a consequent increase of $\Theta$ retrieved values is observed, the ground stations in average did not measure a significant increase of $\Theta$.
<figure><img src="image_11.png" /></figure>

Fig 9 Time-series comparing S-1 soil moisture product with respect to the site observations averaged at the network scale. The in situ average, $\tilde{\mu}_{obs}$, is the blue continuous line and the S-1 average, $\tilde{\mu}_{retr}$, is the red line. The shaded areas represent the daily soil moisture standard deviation. Daily precipitation from a meteo station 10 km far from the site is indicated by the black line. The three outliers over the Apulian Tavoliere are reported as red points.

### 6.1.2 Validation over low-density test sites

Apart from the Apulian Tavoliere, the density of ground stations of all the sites can be considered to be 1 station/km² (Table 1). As a consequence, $\delta_{SRE}$ is expected to contribute significantly to the RMSE. Therefore in this section, $\delta_{SRE}$ is estimated for each site and accounted for to compute the intrinsic RMSE, according to (17) as described in Section 5.2.

Fig. 10 (upper panel) reports per site, the (ub)RMSE, the intrinsic RMSE, $\delta_{SRE}$ at 70% CL and bias. The total number of points and outliers (in brackets) are also reported. $\delta_{SRE}$ ranges between 0.050 $m^3/m^3$ and 0.056 $m^3/m^3$ for the low-density experimental sites and decreases with $\mu_{obs}$. In general, the $\delta_{SRE}$ increases the RMSE by approximately 0.02 $m^3/m^3$ and the intrinsic RMSE ranges between ~0.06$m^3/m^3$ and ~0.08 $m^3/m^3$. Positive biases (i.e., $\Theta$ underestimation) are mostly observed; the highest value ~0.035 $m^3/m^3$ was detected over the Yanco and Txson sites. This effect should be assessed over longer time series and eventually reduced by improving the calibration of the retrieval algorithm. Indeed, although the impact of biases is in average modest in specific areas/times can be fairly high. An anomalous behaviour is observed over Elm Creek, which
shows a RMSE and bias significantly higher than the remaining results. For a further insight, the Elm Creek dataset was split into two parts: spring and summer-autumn, excluding the frozen soil periods. While the summer-autumn retrieval performance is in line with that observed over the other sites (see EC (SuAu) in Fig. 10), the statistical scores estimated in spring are the worst of the entire data set. Fig. 10 (lower panel), shows the correlation R (orange points) observed over each site, including EC (SuAu), versus the standard deviation of Θ<sub>obs</sub>, i.e., δ<sub>Θ<sub>obs</sub></sub> (see Fig. 2). It is noted that 61% of the variability of R is explained by δ<sub>Θ<sub>obs</sub></sub>. The correlation for the whole Elm Creek dataset is reported as a separated black point. It can be observed that R is significantly reduced for Elm Creek if spring Θ time series is added in the comparison. An interpretation of this phenomenon is provided in section 6.2.

724

<figure><img src="image_12.png" /></figure>

725
<figure><img src="image_13.png" /><figcaption>Fig. 10. Performance metrics over the low-density hydrological networks. Upper panel: ubRMSE (blue bars), RMSE (green bars), intrinsic RMSE (red bars), $\delta_{SRE}$ (yellow bars) according to (17) and bias (white bars). The total number of point and removed outliers (in brackets) are reported (AT=Apulian Tavoliere, EC=Elm Creek, EC (SuAu)=Elm Creek for the season summer and autumn, T=TxSON, H= HOBE; Y=Yanco, LW=Little Washita, R=REMEDHUS). Lower panel: ordinary least square Pearson correlation vs $\delta_{\Theta_{obs}}$. Linear fit is also reported R$^2_{fit}$=0.61.</figcaption></figure>

### 6.1.3 Overall results at 1km

An evaluation of the performance of $\Theta_{retr}$ over all the experimental test sites was carried out and the impact of the errors on the dependent and independent variables on the fitting parameters was assessed. Furthermore, in case of the Apulian Tavoliere sites only the S-1 ascending time series was considered, as well as only the summer-autumn periods for the Elm Creek dataset. The total number of $\Theta_{obs}$ and $\Theta_{retr}$ pairs is $M = 15057$.

First, Table 3 summarizes the statistical scores. Accounting for the $\delta_{SRE}$ implies that the intrinsic RMSE is ~0.02 $m^3/m^3$ lower than the RMSE. Then, the WLS was applied to mitigate the effect of the errors on the dependent and independent variables. In this regard, the S-1 $\Theta$ standard deviation, provided with the retrieved S-1 mean $\Theta$ at 1 km resolution, was used as an estimate of the uncertainty associated with the S-1 $\Theta$. Conversely, (12) was used as $\Theta_{obs}$ uncertainty in the WLS for each station on each day. Table 3 reports the WLS R and linear fitting parameters and shows a clear improvement with respect to the correspondent OLS (Réjou-Méchain et al., 2014). Finally,
Fig. 11 illustrates the distribution of intrinsic RMSE (red bars) and RMSE (green bars) per $\Theta_{obs}$ ranges. The RMSE at 1 km is generally higher than the one reported in Fig. 8, because the $\Theta_{obs}$ are not averaged at network scale and therefore include the RSE. Indeed, the intrinsic RMSE is comparable with that observed over the core validation site, per $\Theta$ intervals up to 0.30 $m^3/m^3$. Conversely, for very wet surfaces both RMSEs are higher than those observed in Fig. 8. The reason is that the moderate positive bias, observed in average in Fig. 10 (top panel), is significantly higher in the third interval [0.30 $m^3/m^3$– 0.55 $m^3/m^3$]. Likely, this is the result of various effects. First, the spatial variability of $\Theta_{retr}$ that is larger for wet and very wet soils (see Fig. A1) directly amplifies the bias. Besides, there is a spurious effect because the high tail of the $\Theta_{retr}$ distribution in the second interval [0.20 $m^3/m^3$– 0.30 $m^3/m^3$] largely falls in the third interval further skewing the distribution toward lower values. Finally, the calibration of the algorithm needs to be improved particularly for very wet surfaces.

<table><caption>Table 3. Overall performance metrics over the experimental sites. R refers to the Pearson correlation ($p$ < 0.01). The parameters of the ordinary least squares (OLS) and weighted least squares (WLS) are shown.</caption><thead><tr><th>M<br>(out)</th><th>(ub)RMSE<br>m<sup>3</sup>/m<sup>3</sup></th><th>μ<sub>obs</sub>±<br>δ<sub>θ<sub>obs</sub></sub><br>m<sup>3</sup>/m<sup>3</sup></th><th>μ<sub>retr</sub>±<br>δ<sub>θ<sub>retr</sub></sub><br>m<sup>3</sup>/m<sup>3</sup></th><th>Bias<br>m<sup>3</sup>/m<sup>3</sup></th><th>δ<sub>SRE</sub><br>(70%CL)<br>m<sup>3</sup>/m<sup>3</sup></th><th>intrinsic<br>(ub)RMSE<br>m<sup>3</sup>/m<sup>3</sup></th><th>OLS<br>linear fit<br>R</th><th>WLS linear<br>fit<br>R</th></tr></thead><tbody><tr><td>15057<br>(82)</td><td>0.088<br>(0.085)</td><td>0.188±<br>0.090</td><td>0.168±<br>0.073</td><td>0.021</td><td>0.053</td><td><strong>0.070</strong><br>(0.067)</td><td>y=0.377*x<br>+ 0.097<br>R=0.46</td><td>y=0.810*x +<br>0.017 and<br><strong>R=0.54</strong></td></tr></tbody></table>

762
763
764
<figure><img src="image_14.png" /><figcaption>Fig. 11. Distribution of RMSE (green bars) and intrinsic RMSE (red bars) according to (17) per $\Theta_{obs}$ interval.</figcaption></figure>

## 6.2 Validation at the network scale

Fig. 12 reports the validation results at the network scale, i.e., both $\Theta_{retr}$ and $\Theta_{obs}$ were spatially averaged at the network scale on each S-1 acquisition date, $\tilde{\mu}_{obs}$ and $\tilde{\mu}_{retr}$. Results of the whole (EC) and the summer-autumn Elm Creek (EC (SuAu)) dataset are reported separately both for the RMSE and the Pearson correlation. The estimated RMSE (upper panel) ranges between approximately $0.050m^3/m^3$ and $0.066m^3/m^3$, excluding EC. The Pearson correlation (lower panel) is reported as a function of the standard deviation of the temporal series of the daily spatial mean $\tilde{\mu}_{obs}$, i.e. $\delta_{\tilde{\mu}_{obs}}$. R can still be ordered by $\delta_{\tilde{\mu}_{obs}}$ as for Fig. 10 (lower panel), but R at the network scale, ranging between 0.35 and 0.86, generally increases due to the spatial average process which enhances the sensitivity to $\Theta$. Besides, higher correlation (> 0.60) are observed for the sites imaged at incidence angle lower than 35° (yellow points). In particular, the lower correlation coefficient for the Apulian Tavoliere descending (D) time series relative to the ascending (A) passes is likely due to the higher S-1 incidence angle (Palmisano et al., 2020).
<figure><img src="image_15.png" /><figcaption>Fig. 12. Statistical scores at the site scale. Upper panel: ubRMSE (blue bars), RMSE (green bars), and bias (white bars) per site. The total number of compared points is reported. Lower panel: Pearson correlation vs the observed $\delta_{\tilde{\mu}_{obs}}$. Linear fits are also reported, $R^2_{fit}=0.57$. For the Apulian Tavoliere site, the metrics for both the S-1 ascending (A) and descending (D) tracks are shown (AT=Apulian Tavoliere, EC=Elm Creek, EC (SuAu) =Elm Creek in summer-autumn, T=TxSON, H= HOBE Y=Yanco, LW=Little Washita, R=REMEDHUS). Sites imaged at incidence angle lower than 35deg are shown identified by the yellow points.</figcaption></figure>

Fig. 13 and Fig. 14 show the temporal behaviour of retrieved and observed $\Theta$ averages, together with their standard deviations, at the Yanco and Elm Creek network scales. Over the Yanco site, the daily precipitation is available for several sub-areas, therefore, they were averaged on each date to obtain an estimate representative of the network. In particular, Fig. 13 confirms that STCD (red
line) reproduces fairly well the changes of Θ due to precipitation events and dry down cycles, over the entire period of approximately 4 years (R=0.86). It is noted that Yanco, unlike Segezia, cannot be considered a homogeneous site in terms of soil properties and vegetation cover. As a consequence, the standard deviation of Θ<sub>retr</sub> observed in Fig. 13 (red shaded area) is likely driven by land cover, topography and soil conditions.

The Elm Creek site was selected because it is a challenging case for the retrieval algorithm, due to a number of site factors. First, Fig. 14 shows the presence of large gaps of data, mainly in winter, which are due to frozen soil conditions. In addition, missing S-1 data in July and August 2017 inhibited the retrieval during the active growing season, a period of favourable conditions for estimating soil moisture. Moreover, a systematic underestimation of Θ in springtime is observed. This underestimation is related to a biased estimation of the α<sub>min</sub> parameter at a coarse scale, which propagated the bias at high resolution. In May an important drop in backscatter (i.e., 4 – 5 dB with respect to April) was observed at large scale in the Manitoba region. Spring wheat, canola and corn, which are the main crops of the area, are typically seeded in the first 10 days of May (Powers, 2021). As a result, the status of the soil roughness changes in large areas of Manitoba. During this active spring period, roughness diminishes as farmers prepare a smooth soil seedbed. As crops emerge, canopies in fill with varying canopy structures and scattering mechanisms. The S-1 response increases on average and is again more distinct among fields. Currently, this dynamic is not handled well by the SMOSAR calibration (at least over the spring period). Under these conditions, a further understanding and adaptation of the code is required.
<figure><img src="image_16.png" /><figcaption>Fig 13 Time-series comparing S-1 soil moisture data product with respect to Yanco site observations averaged at the network scale. The in situ average is the blue continuous line and the S-1 average is the red line. The shaded areas represent the daily standard deviation. Daily precipitation averaged at network scale is indicated by the black line.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig 14 The same as Fig. 13 but for the Elm creek site. Periods with frozen soils show no data.</figcaption></figure>

Finally, Fig. 15 compares $\Theta_{retr}$ against $\Theta_{obs}$ at the network scale for all the validation sites (EC (SuAu) included) and the S-1 acquisitions. The statistical scores are also reported. The overall correlation and RMSE are 0.67 and 0.058 $m^3/m^3$, respectively. Three outliers were removed, as in Fig. 8. Overall, the bias and the $\delta_{SRE}$ are quite small (e.g., $\delta_{SRE} \ll 0.03m^3/m^3$) for the seven
experimental sites. Therefore, $\delta_{SRE}$ has a marginal impact on the RMSE and the OLS method was adopted for estimating the linear fit in Fig. 16, i.e., $y = 0.735x + 0.045$.

<figure><img src="image_18.png" /><figcaption>Fig. 15. Site scale comparison between $\Theta$ retrieved from S-1 and observed over AT=Apulian Tavoliere (A146), AT=Apulian Tavoliere (D124), EC (SuAu) =Elm Creek in summer-autumn, T=TxSON, H= HOBE, Y=Yanco, LW=Little Washita, R=REMEDHUS sites. The comparison includes 1068 dates. Three outliers (>3 standard deviations), i.e. black points, are also reported.</figcaption></figure>

837

## 7 Conclusions

This study presents an extended validation of a pre-operational surface soil volumetric water content $\Theta$ product [$m^3/m^3$] at ~1 km resolution derived from VV&VH S-1 observations. The VH S-1 channel is used for the dynamic masking of vegetation, while the $\Theta$ retrieval is based on the VV S-1 observations. Only static information about land cover and soil texture is needed for $\Theta$ retrieval in addition to the S-1 backscatter. The S-1 $\Theta$ retrieval algorithm consists of a time series based short term change detection approach. This approximation decouples the effect of $\Theta$ from that of the other surface parameters (e.g., surface roughness and vegetation) on the SAR signal, hence significantly simplifying the $\Theta$-retrieval and expediting the processing. The strength of the retrieval algorithm is its conceptual simplicity and its robustness as the $\Theta$ estimates depend on a single free parameter, i.e., $\alpha_{min}$. Conversely, the accuracy of the S-1 $\Theta$ product depends on: i) the statistical

830

831

832

833

834

835

836

838

839

840

841

842

843

844

845

846

847

848
variability of the backscatter temporal ratio between subsequent S-1 acquisitions, and ii) the calibration of the retrieval algorithm needed to identify $\alpha_{min}$. The latter is a continuous process improving with the length of the time series of S-1 and in situ data analysed (in particular, wet and very wet soil conditions are still underrepresented). An important aspect of the developed $\Theta$ product is that a co-registered standard deviation layer is also provided with each mean $\Theta$ product at 1 km. A quantitative analysis on the dependence of $\Theta$ uncertainty on land cover, topography, soil texture and mean $\Theta$, at increasing spatial scales, remains to be assessed.

In the study, the performance of the S-1 mean $\Theta$ product was estimated through direct comparison against in situ $\Theta$-observations recorded over 7 hydrologic networks with 167 ground stations, located in Italy, Spain, Denmark, Canada, USA and Australia covering 4 years between January 2015 and December 2020. An emphasis was also placed on addressing the SRE that arises from the mismatch between the S-1 $\Theta$ retrieved at ~1 km and the in situ point-scale $\Theta$-observations and its impact on the metrics. At the core validation site (4.3 stations/km²) in the Apulian Tavoliere (Italy), the in situ $\Theta$ observations were upscaled at 1 km to a relatively low SRE (i.e., SRE < 0.03 m³/m³, 70% C.L) and therefore it was considered as a benchmark for the S-1 $\Theta$ performance. The observed correlation between S-1 $\Theta$ and in situ observations is 0.71 and the RMSE is ~0.06 m³/m³ with a bias of ~0.01 m³/m³ over the $\Theta$ range [0.03, 0.60] m³/m³. The S-1 $\Theta$ time series over Southern Italy is available through the connected Data in Brief article (Balenzano et al, 2021).

A first recommendation steaming from this study is to establish cal/cal sites (either new or as an update of existing hydrologic networks) dedicated to SAR soil moisture retrieval having a minimum density of 4-5 stations/km². It is noted that such an estimate is independent on the characteristics of the SAR sensor as it is solely based on consideration of the statistics of the observed $\Theta$ fields. A second suggestion is to reconsider the requirement for SAR $\Theta$ product accuracy, usually expressed in terms of the RMSE. The motivation is that the SAR signal
uncertainty increases with its mean, and this statistical property is also observed on the derived Θ estimates. A more realistic requirement could therefore be to adopt a relative error (either unbiased or not, i.e., (ub)RMSE/mean). For instance, the results of this study would indicate a level of 20% for Θ equal or higher than 0.20 m³/m³. For lower Θ values, a constant threshold for (ub)RMSE of, e.g., 0.05 m³/m³ could be adopted.

In terms of outlook, despite the important progress that the launch of the S-1 constellation has brought for high-resolution Θ retrieval, technical and programmatic issues persist. The former mainly concern the temporal resolution of presently available satellite SAR data that is still non-optimal for Θ applications, as user requirements point to a temporal resolution of 1-2 days or less [e.g., www.wmo-sat.info]. To tackle this issue a programmatic effort for a coordinated acquisition plan of multi-mission SAR data, such as those acquired by S-1 and the recently launched RADARSAT Constellation Mission, or those provided by S-1 and the forthcoming EU L-band Radar Observation System for Europe (ROSE-L) system (Davidson et al., 2019), is needed. Success in this approach will require the development of retrieval algorithms combining C- and L-band SAR data (Zhu et al., 2019), which can also lead to simultaneous and consistent retrieval of Θ and vegetation water content at high temporal and spatial resolution with beneficial impact on coupled hydrology-crop growth models (Pauwels and Verhoest, 2007) and on a better understanding of the land-atmosphere interaction (Vereecken et al., 2010).

892

893 Acknowledgment

This research was supported by the Scientific Exploitation of Operational Missions (SEOM) program of the European Space Agency, through the project "Exploitation of S-1 for Surface Soil Moisture Retrieval at High Resolution (Exploit-S-1)" (contract 4000118762/16/I-NB). The authors express sincere thanks to the data providers and the International Soil Moisture Network for the REMEDHUS, HOBE and Elm Creek data and USDA Agricultural Research Service for the Little

898
Washita data. The authors are grateful to Dr Todd Caldwell for providing the TxSON data and Prof Jeffrey Walker for supplying the Yanco data, respectively. Dr Frederik Uldall is acknowledged for his support in preprocessing the Hobe data. The authors are indebted to Dr. Jarrett Powers for his support in the interpretation of the Elm Creek results and to Dr. Michele Rinaldi for his continuous assistance on the management of the Apulian Tavoliere hydrologic network. Finally, the authors are grateful to the Reviewers for their valuable comments which helped to improve the quality of this paper.

## APPENDIX A: Error Budget

The symbolic expression for the variance of the retrieved $Θ$, $δ_F^2$, at pixel scale and for a given date, is reported in (5), which shows that $δ_F^2$ depends on $δ_{|α_{VV}|}^2$. Here, the procedure for estimating $δ_{|α_{VV}|}^2$ is briefly sketched. The quantitative retrieval of $|α_{VV}(ε, ϑ)|$ is based on (7), which for a specific date $i$ can be written as

$$|\hat{\alpha}_{VV}| = \hat{\lambda} \cdot \hat{S}_{iN} \quad i = 1, \dots, N \qquad (A1)$$

where $\hat{\lambda}$ and $\hat{S}_{iN}$ are two random variables. $\hat{S}_{iN}$ is defined in (6); for $i = N - 1$ and $L \gg 10$ ($L$ is the number of looks) the probability density function of $\hat{S}_{iN}$ is normal and its moments' expression can be found in (Lee et al., 1994; Oliver and Quegan, 1998). In particular, using the Stirling's approximation of a gamma function ratio (Tricomi and Erdelyi, 1951), the estimate of its sample coefficient of variation (CV) decreases with $L$ as

$$\frac{\tilde{\delta}_{\hat{S}_{iN}}}{\tilde{\mu}_{\hat{S}_{iN}}} \approx \frac{1}{\sqrt{2L}} \qquad (A2)$$

where $\tilde{\delta}_{\hat{S}_{iN}}$ and $\tilde{\mu}_{\hat{S}_{iN}}$ are the spatial sample standard deviation and mean, computed over $L \gg 10$ independent samples.

For $i = N$, $\hat{S}_{NN}$ assumes constant values equal to 1 and, therefore, its variance is zero. The maximum likelihood estimate of $\hat{\lambda}$ is given in (9) and can be cast in the following form
$$922 \qquad \hat{\lambda} = \max\left(\frac{\alpha_{min}}{\hat{S}_{iN}}\right) = \hat{\alpha}_{min} \cdot \max(\hat{S}_{Ni}) = \hat{\alpha}_{min} \cdot \hat{S}_{Nj^*} \quad i = 1, \dots, N \qquad (A3)$$

where $\hat{S}_{Nj^*}$ has the same statistics of $\hat{S}_{Ni}$ in (A1). For the sake of notation simplicity, from now on $\hat{S}_{Nj^*}$ will be indicated as $\hat{S}^*$. In (A3), $\hat{\alpha}_{min}$ is obtained from the calibration curve, discussed in Section 4.1.2.1, and its variance is derived later on in this Appendix. The two random variables $\hat{S}^*$ and $\hat{\alpha}_{min}$ are independent, therefore the sample mean of $\hat{\lambda}$ is

$$927 \qquad \tilde{\mu}_{\hat{\lambda}} = \tilde{\mu}_{\hat{\alpha}_{min}} \cdot \tilde{\mu}_{\hat{S}^*}$$

(A4)

and the variance $\tilde{\delta}_{\hat{\lambda}}^2$ can be simply derived as

$$930 \qquad \tilde{\delta}_{\hat{\lambda}}^2 = \tilde{\mu}_{\hat{\alpha}_{min}}^2 \cdot \tilde{\delta}_{\hat{S}^*}^2 + \tilde{\mu}_{\hat{S}^*}^2 \cdot \tilde{\delta}_{\hat{\alpha}_{min}}^2$$

(A5)

From (A4) and (A5) and for a sufficiently large $L$, the sample squared CV is

$$933 \qquad \frac{\tilde{\delta}_{\hat{\lambda}}^2}{\tilde{\mu}_{\hat{\lambda}}^2} = \frac{\tilde{\delta}_{\hat{\alpha}_{min}}^2}{\tilde{\mu}_{\hat{\alpha}_{min}}^2} + \frac{\tilde{\delta}_{\hat{S}^*}^2}{\tilde{\mu}_{\hat{S}^*}^2} \approx \frac{\tilde{\delta}_{\hat{\alpha}_{min}}^2}{\tilde{\mu}_{\hat{\alpha}_{min}}^2} + \frac{1}{2L} \qquad (A6)$$

In the same vein as for (A5), we can estimate the variance $\delta_{|\alpha_{VV}|}^2$ of $|\alpha_{VV}|$ in (A1). It is, however, worth noting that in (A1), $\hat{\lambda}$ and $\hat{S}_{iN}$ are somehow correlated because both in $\hat{S}^*$ and $\hat{S}_{iN}$ appears the backscatter acquired on the last date $N$ in the time series. Under these circumstances, the mean of $|\alpha_{VV}|$ is

$$938 \qquad \tilde{\mu}_{|\alpha_{VV}|} = \tilde{\mu}_{\hat{\lambda}} \cdot \tilde{\mu}_{\hat{S}_{iN}} \cdot \left[ 1 + \tilde{\rho}_{\hat{\lambda}, \hat{S}_{iN}} \cdot \frac{\tilde{\delta}_{\hat{\lambda}}}{\tilde{\mu}_{\hat{\lambda}}} \cdot \frac{\tilde{\delta}_{\hat{S}_{iN}}}{\tilde{\mu}_{\hat{S}_{iN}}} \right] \le \tilde{\mu}_{\hat{\lambda}} \cdot \tilde{\mu}_{\hat{S}_{iN}} \cdot \left[ 1 + \frac{\tilde{\delta}_{\hat{\lambda}}}{\tilde{\mu}_{\hat{\lambda}}} \cdot \frac{1}{\sqrt{2L}} \right] \approx \tilde{\mu}_{\hat{\lambda}} \cdot \tilde{\mu}_{\hat{S}_{iN}} \qquad (A7)$$

where $\tilde{\rho}_{\hat{\lambda}, \hat{S}_{iN}}$ is the sample Pearson correlation coefficient between $\hat{\lambda}$ and $\hat{S}_{iN}$. Therefore, for a sufficiently large $L$, the variance of $|\alpha_{VV}|$ is bounded as

$$941 \qquad 0 \le \tilde{\delta}_{|\alpha_{VV}|}^2 \le \tilde{\mu}_{\hat{\lambda}}^2 \cdot \tilde{\delta}_{\hat{S}_{iN}}^2 + \tilde{\mu}_{\hat{S}_{iN}}^2 \cdot \tilde{\delta}_{\hat{\lambda}}^2 = \tilde{\mu}_{\hat{\lambda}}^2 \cdot \tilde{\mu}_{\hat{S}_{iN}}^2 \cdot \left[ \frac{\tilde{\delta}_{\hat{\lambda}}^2}{\tilde{\mu}_{\hat{\lambda}}^2} + \frac{\tilde{\delta}_{\hat{S}_{iN}}^2}{\tilde{\mu}_{\hat{S}_{iN}}^2} \right] \quad i = 1, \dots, N \qquad (A8)$$
942 Substituting (A6) and (A7) into (A8), we obtain

$$
943 \qquad \delta_{|\alpha_{VV}|}^2 \le \tilde{\mu}_{|\alpha_{VV}|}^2 \cdot \left[ \frac{\tilde{\delta}_{\tilde{\alpha}_{min}}^2}{\tilde{\mu}_{\tilde{\alpha}_{min}}^2} + \frac{1}{L} \right] \quad i = 1, \dots, N \tag{A9}
$$

944 To estimate $\frac{\tilde{\delta}_{\tilde{\alpha}_{min}}^2}{\tilde{\mu}_{\tilde{\alpha}_{min}}^2}$ it is necessary to consider the calibration curve of $|\hat{\alpha}_{VV}|^2$ versus $\gamma$, obtained at
945 the low resolution ($\gamma$ is the S-1 backscatter divided by the cosine of the incidence angle). The curve
946 is a first order polynomial $|\hat{\alpha}_{VV}|^2 = \hat{A} \cdot \hat{\gamma} + \hat{B}$, then $\hat{\alpha}_{min}$ is

$$
947 \qquad \hat{\alpha}_{min} = |\hat{\alpha}_{VV}|_{min} = \sqrt{\hat{A} \cdot \hat{\gamma}_{min} + \hat{B}} \tag{A10}
$$

948 The related CV is derived by using the propagation of uncertainty on $\hat{\alpha}_{min}$, considering $\hat{\gamma}_{min}$, $\hat{A}$ and
949 $\hat{B}$ all affected by errors.

$$
950 \qquad \frac{\delta_{\hat{\alpha}_{min}}^2}{\mu_{\hat{\alpha}_{min}}^2} \approx \frac{1}{4} \cdot \frac{(\mu_{\hat{\gamma}_{min}}^2 \cdot \delta_A^2 + \delta_B^2) + \mu_A^2 \cdot \delta_{\hat{\gamma}_{min}}^2}{[\mu_A \cdot \mu_{\hat{\gamma}_{min}} + \mu_B]^2} \qquad (A11)
$$

951 Substituting (A11) into (A9), it results

$$
952 \qquad \tilde{\delta}_{|\alpha_{VV}|}^2 \le \tilde{\mu}_{|\alpha_{VV}|}^2 \cdot \left\{ \frac{1}{4} \cdot \frac{(\mu_{\tilde{\gamma}_{min}}^2 \cdot \delta_A^2 + \delta_B^2) + \mu_A^2 \cdot \delta_{\tilde{\gamma}_{min}}^2}{[\mu_A \cdot \mu_{\tilde{\gamma}_{min}} + \mu_B]^2} + \frac{1}{L} \right\} = [\tilde{\delta}_{cal}^2 + \tilde{\delta}_{stat}^2] \tag{A12}
$$

954 In Fig. A1, it is reported the retrieved Θ ± δ<sub>F</sub> as a function of the true |α<sub>VV</sub>| (and Θ on the
955 secondary x-axis) for VV polarization and 30° incidence, where δ<sub>F</sub> is evaluated by (5), using (A12).
956 The predicted retrieval error increases from dry to wet surfaces. This is expected to be a general
957 property that can be easily understood considering that i) |α<sub>VV</sub>|<sup>2</sup> is almost linear with Θ (e.g., (Kim
958 and van Zyl, 2009)) and ii) the radar backscatter is proportional to |α<sub>VV</sub>|<sup>2</sup>, see (1). Therefore, the
959 standard error on Θ increases with the standard error on the radar backscatter, which in turn is
960 proportional to the backscatter itself. As a result, the higher the retrieved Θ, the higher its standard
961 error. In terms of the weight of δ̃<sub>stat</sub><sup>2</sup> versus δ̃<sub>cal</sub><sup>2</sup>, for Θ retrieved at a resolution of 1 km,
corresponding to an equivalent number of looks L~10⁴, the estimated statistical error is always
below 0.02 m³/m³. Therefore, the major contribution to the error comes from the calibration error,
as for instance δF is approximately 13% of Θ in average over the range of Θ variability.

<figure><img src="image_19.png" /><figcaption>Fig. A1. Example of the error budget for 1 km Θ retrieved at VV polarization and 30° incidence angle.</figcaption></figure>

968

969 References

970 Al-khaldi, M.M., Johnson, J.T., Brien, A.J.O., Balenzano, A., Mattia, F., Member, S., 2019. Time-
971 Series Retrieval of Soil Moisture Using CYGNSS. IEEE Trans. Geosci. Remote Sens. 57,
972 4322-4331. https://doi.org/10.1109/TGRS.2018.2890646

973 Balenzano, A., Mattia, F., Satalino, G., Lovergine, P.L., Palmisano D., Davidson, M.W.J., 2021.
974 Dataset of Sentinel-1 surface soil moisture time series at 1 km resolution over Southern Italy.
975 Data in brief.

976 Balenzano, A., Satalino, G., Iacobellis, V., Gioia, A., Manfreda, S., Rinaldi, M., De Vita, P.,
977 Miglietta, F., Toscano, P., Annicchiarico, G., Mattia, F., 2014. A ground network for SAR-
978 derived soil moisture product calibration, validation and exploitation in Southern Italy, in:
979 International Geoscience and Remote Sensing Symposium (IGARSS).
980 https://doi.org/10.1109/IGARSS.2014.6947206

981 Balenzano, A., Satalino, G., Lovergine, F., Rinaldi, M., Iacobellis, V., Mastronardi, N., Mattia, F.,
982 2013. On the use of temporal series of L- and X-band SAR data for soil moisture retrieval.
983 Capitanata plain case study. Eur. J. Remote Sens. 46. https://doi.org/10.5721/EuJRS20134643

984 Balenzano, A., Mattia, F., Satalino, G., Davidson, M.W.J., 2011. Dense temporal series of C- and
985 L-band SAR data for soil moisture retrieval over agricultural crops, in: IEEE Journal of
986 Selected Topics in Applied Earth Observations and Remote Sensing. pp. 439–450.

987 Bauer-Marschallinger, B., Freeman, V., Cao, S., Paulik, C., Schaufler, S., Stachl, T., Modanesi, S., 988 Massari, C., Ciabatta, L., Brocca, L., Wagner, W., 2019. Toward Global Soil Moisture 989 Monitoring with Sentinel-1: Harnessing Assets and Overcoming Obstacles. IEEE Trans. 990 Geosci. Remote Sens. 57, 520–539. https://doi.org/10.1109/TGRS.2018.2858004

991 Bauer-Marschallinger, B., Paulik, C., Hochstöger, S., Mistelbauer, T., Modanesi, S., Ciabatta, L., 992 Massari, C., Brocca, L., Wagner, W., 2018. Soil moisture from fusion of scatterometer and 993 SAR: Closing the scale gap with temporal filtering. Remote Sens. 10, 1–26. 994 https://doi.org/10.3390/rs10071030

995 Bircher, S., Skou, N., Jensen, K.H., Walker, J.P., Rasmussen, L., 2012. A soil moisture and 996 temperature network for SMOS validation in Western Denmark. Hydrol. Earth Syst. Sci. 16, 997 1445–1463. https://doi.org/10.5194/hess-16-1445-2012

998 Bourbigot, M., Johnsen, H., Piantanida, R., Hajduch, G., 2016. Sentinel-1 Product Definition. ESA 999 Doc. no. S1-RS-MDA-52-7440.

1000 Brocca, L., Melone, F., Moramarco, T., Morbidelli, R., 2010. Spatial-temporal variability of soil 1001 moisture and its estimation across scales. Water Resour. Res. 46, 1–14. 1002 https://doi.org/10.1029/2009WR008016

1003 Caldwell, T.G., Bongiovanni, T., Cosh, M.H., Jackson, T.J., Colliander, A., Abolt, C.J., Casteel, R., 1004 Larson, T., Scanlon, B.R., Young, M.H., 2019. The Texas Soil Observation Network: A 1005 Comprehensive Soil Moisture Dataset for Remote Sensing and Land Surface Model 1006 Validation. Vadose Zo. J. https://doi.org/10.2136/vzj2019.04.0034

1007 Cantrell, C.A., 2008. Technical Note: Review of methods for linear least-squares fitting of data and 1008 application to atmospheric chemistry problems. Atmos. Chem. Phys. 8, 5477–5487. 1009 https://doi.org/10.5194/acp-8-5477-2008

1010 Caves, R., “Final report: Technical assistance for the implementation of the AgriSAR 2009 1011 campaign,” ESA, Paris, France, Contract 22689/09, 2009.

1012 Colliander, A., Jackson, T.J., Bindlish, R., Chan, S., Das, N., Kim, S.B., Cosh, M.H., Dunbar, R.S., 1013 Dang, L., Pashaian, L., Asanuma, J., Aida, K., Berg, A., Rowlandson, T., Bosch, D., Caldwell, 1014 T., Caylor, K., Goodrich, D., al Jassar, H., Lopez-Baeza, E., Martínez-Fernández, J., González- 1015 Zamora, A., Livingston, S., McNairn, H., Pacheco, A., Moghaddam, M., Montzka, C., 1016 Notarnicola, C., Niedrist, G., Pellarin, T., Prueger, J., Pulliainen, J., Rautiainen, K., Ramos, J., 1017 Seyfried, M., Starks, P., Su, Z., Zeng, Y., van der Velde, R., Thibeault, M., Dorigo, W., 1018 Vreugdenhil, M., Walker, J.P., Wu, X., Monerris, A., O’Neill, P.E., Entekhabi, D., Njoku, 1019 E.G., Yueh, S., 2017. Validation of SMAP surface soil moisture products with core validation 1020 sites. Remote Sens. Environ. 191, 215–231. https://doi.org/10.1016/j.rse.2017.01.021

1021 Cookmartin, G., Saich, P., Quegan, S., Cordey, R., Burgess-allen, P., Sowter, A., 2000. Modeling 1022 microwave interactions with crops and comparison with ERS-2 SAR observations. IEEE 1023 Trans. Geosci. Remote Sens. vol. 38pp658-670.

1024 Coopersmith, E.J., Cosh, M.H., Bell, J.E., Crow, W.T., 2016. Multi-Profile Analysis of Soil 1025 Moisture within the US Climate Reference Network. Vadose Zo. J. 15, 1–8.

1026 Cosh, M.H., Jackson, T.J., Starks, P., Heathman, G., 2006. Temporal stability of surface soil 1027 moisture in the Little Washita River watershed and its applications in satellite soil moisture 1028 product validation. J. Hydrol. 323, 168–177. https://doi.org/10.1016/j.jhydrol.2005.08.020

1029 Crow, W. T., A.A. Berg, M.H. Cosh, A. Loew, B.P. Mohanty, R. Panciera, P. Rosnay, D. Ryu, J.P. 1030 Walker, 2012. Upscaling sparse ground-based soil moisture observations for the validation of 1031 coarse-resolution satellite soil moisture products. Reviews of Geophysics, vol. 50, no 3.

1032 Das, N.N., Entekhabi, D., Njoku, E.G., Shi, J.J.C., Johnson, J.T., Colliander, A., 2014. Tests of the 1033 SMAP combined radar and radiometer algorithm using airborne field campaign observations 1034 and simulated data. IEEE Trans. Geosci. Remote Sens. 52, 2018–2028. 1035 https://doi.org/10.1109/TGRS.2013.2257605
1036 Das, N.N., Mohanty, B.P., Njoku, E.G., 2010. Profile Soil Moisture Across Spatial Scales Under
1037 Different Hydroclimatic Conditions. Soil Sci. 175, 315-320.

1038 Davidson, M.W.J., Chini, M., Dierking, W., Djavidnia, S., Haarpaintner, J., Hajduch, G., Vaglio
1039 Laurin, G., Lavalle, M., Lopez Martinez, C., Nagler, T., Pierdicca, N., Su, B., 2019.
1040 Copernicus L-band SAR Mission Requirements Document.

1041 Dharssi, I., Bovis, K.J., Macpherson, B., Jones, C.P., 2011. Operational assimilation of ASCAT
1042 surface soil wetness at the Met Office. Hydrol. Earth Syst. Sci. 15, 2729-2746.
1043 https://doi.org/10.5194/hess-15-2729-2011

1044 Dorigo, W., Wagner, W., Albergel, C., Albrecht, F., Balsamo, G., Brocca, L., Chung, D., Ertl, M.,
1045 Forkel, M., Gruber, A., Haas, E., Hamer, P.D., Hirschi, M., Ikonen, J., Jeu, R. De, Kidd, R.,
1046 Lahoz, W., Liu, Y.Y., Miralles, D., Mistelbauer, T., Nicolai-shaw, N., Parinussa, R., Pratola,
1047 C., Reimer, C., Schalie, R. Van Der, Seneviratne, S.I., Smolander, T., Lecomte, P., 2017. ESA
1048 CCI Soil Moisture for improved Earth system understanding: State-of-the art and future
directions. Remote Sens. Environ. 203, 185-215. https://doi.org/10.1016/j.rse.2017.07.001
1049

1050 Dorigo, W., Xaver, A., Vreugdenhil, M., Gruber, A., Hegyiová, A., Sanchis-Dufau, A.D., Zamojski,
1051 D., Cordes, C., Wagner, W., Drusch, M., 2013. Global Automated Quality Control of In Situ
1052 Soil Moisture Data from the International Soil Moisture Network. Vadose Zo. J. 12,
vzj2012.0097. https://doi.org/10.2136/vzj2012.0097

1053
1054 Dorigo, W.A., Gruber, A., De Jeu, R.A.M., Wagner, W., Stacke, T., Loew, A., Albergel, C.,
1055 Brocca, L., Chung, D., Parinussa, R.M., Kidd, R., 2015. Evaluation of the ESA CCI soil
moisture product using ground-based observations. Remote Sens. Environ. 162, 380-395.
1056
1057 https://doi.org/10.1016/j.rse.2014.07.023

1058 Dorigo, W.A., Wagner, W., Hohensinn, R., Hahn, S., Paulik, C., Xaver, A., Gruber, A., Drusch, M.,
1059 Mecklenburg, S., Van Oevelen, P., Robock, A., Jackson, T., 2011. The International Soil
1060 Moisture Network: A data hosting facility for global in situ soil moisture measurements.
1061 Hydrol. Earth Syst. Sci. 15, 1675-1698. https://doi.org/10.5194/hess-15-1675-2011

1062 El Hajj, M., Baghdadi, N., Zribi, M., Bazzi, H., 2017. Synergic use of Sentinel-1 and Sentinel-2
1063 images for operational soil moisture mapping at high spatial resolution over agricultural areas.
1064 Remote Sens. 9, 1-28. https://doi.org/10.3390/rs9121292

1065 Entekhabi, B.D., Njoku, E.G., Neill, P.E.O., Kellogg, K.H., Crow, W.T., Edelstein, W.N., Entin,
1066 J.K., Goodman, S.D., Jackson, T.J., Johnson, J., Kimball, J., Piepmeier, J.R., Koster, R.D.,
1067 Martin, N., Mcdonald, K.C., Moghaddam, M., Moran, S., Reichle, R., Shi, J.C., Spencer,
1068 M.W., Thurman, S.W., Tsang, L., Zyl, J. Van, 2010. The Soil Moisture Active Passive
1069 (SMAP) Mission. Proc. IEEE 98.

1070 Entekhabi, D., Reichle, R.H., Koster, R.D., Crow, W.T., 2010. Performance Metrics for Soil
1071 Moisture Retrievals and Application Requirements. J. Hydrometeorol. 11, 832-840.
1072 https://doi.org/10.1175/2010JHM1223.1

1073 Entin, J., Robock, A., Vinnikov, K.Y., Hollinger, S.E., Liu, S., Namkhai, A., 2000. Temporal and
1074 Spatial Scales of Observed Soil Moisture Variations in the Extratropics. J. Geophys. Res. 105,
1075 11,865-11,877. https://doi.org/10.1029/2000JD900051

1076 Famiglietti, J.S., Devereaux, J.A., Laymon, C.A., Tsegaye, T., Houser, P.R., Jackson, T.J., Graham,
1077 S.T., Rodell, M., Oevelen, P.J. Van, 1999. Ground-based investigation of soil moisture
1078 variability within remote sensing footprints during the Southern Great Plains 1997 (SGP97)
1079 Hydrology Experiment. Water Resour. 35, 1839-1851.

1080 Famiglietti, J.S., Ryu, D., Berg, A.A., Rodell, M., Jackson, T.J., 2008. Field observations of soil
1081 moisture variability across scales. Water Resour. Res. 44, 1-16.
1082 https://doi.org/10.1029/2006WR005804

1083 Gilbert, R.O., 1987. Statistical Methods for Environmental Pollution Monitoring. Van Nostrand
1084 Reinhold Company Inc., New York.

1085 Gruber, A., Lannoy, G. De, Albergel, C., Al-yaari, A., Brocca, L., Calvet, J., Colliander, A., Cosh,
1086 M., Crow, W., Dorigo, W., Draper, C., Hirschi, M., Kerr, Y., Konings, A., Lahoz, W., Mccoll, K., Montzka, C., Muñoz-sabater, J., Peng, J., Reichle, R., Richaume, P., Rüdiger, C., 2020. Remote Sensing of Environment Validation practices for satellite soil moisture retrievals: What are ( the ) errors? Remote Sens. Environ. 244, 111806.
https://doi.org/10.1016/j.rse.2020.111806

1091 Hajnsek, I., R. Bianchi, M. Davidson, G. D'Urso, J. A. Gomez-Sanchez, A. Hausold, R. Horn, J. Howse, A. Loew, J. M. Lopez-Sanchez, R. Ludwig, J. A. Martinez-Lozano, F. Mattia, E. Miguel, J. Moreno, V. R. N. Pauwels, T. Ruhtz, C. Schmullius, H. Skriver, J. A. Sobrino, W. Timmermans, C. Wloczyk, and M. Wooding, "AGRISAR optical and radar campaign," ESA, Paris, France, Tech. Rep. Contract 19974/06/I/LG, 2007, Final Report

1096 Hallikainen, M.T., Ulaby, F.T., Dobson, M.C., El-rayes, M.A., Wu, L., 1985. Microwave dielectric behavior of wet soil-Part I: Empirical models and experimental observations. IEEE Trans. Geosci. Remote Sensing, vol. GE-23 ge-23, 25-34.

1099 Heimhuber, V., Tulbure, M.G., Broich, M., 2017. Modeling multidecadal surface water inundation dynamics and key drivers on large river basin scale using multiple time series of Earth-observation and river flow data. Water Resour. Res. 53, 1251-1269.
https://doi.org/10.1002/2016WR019858.Received

1103 Hengl, T., De Jesus, J.M., Heuvelink, G.B.M., Gonzalez, M.R., Kilibarda, M., Blagotić, A., Shangguan, W., Wright, M.N., Geng, X., Bauer-Marschallinger, B., Guevara, M.A., Vargas, R., MacMillan, R.A., Batjes, N.H., Leenaars, J.G.B., Ribeiro, E., Wheeler, I., Mantel, S., Kempen, B., 2017. SoilGrids250m: Global gridded soil information based on machine learning, PLoS ONE. https://doi.org/10.1371/journal.pone.0169748

1108 Hobbs, S.E., Guarnieri, A.M., Broquetas, A., Calvet, J., Chini, M., Ferretti, R., Nagler, T., Pierdicca, N., Wadge, G., 2019. G-CLASS: geosynchronous radar for water cycle science - orbit selection and system design. IET Int. Radar Conf. 2019, 7534-7537.
https://doi.org/10.1049/joe.2019.0601

1112 Iacobellis, V., Gioia, A., Milella, P., Satalino, G., Balenzano, A., Mattia, F., 2013. Inter-comparison of hydrological model simulations with time series of SAR-derived soil moisture maps. Eur. J. Remote Sens. 46, 739-757.

1115 Jacobs, J.M., Mohanty, B.P., Hsu, E.C., Miller, D., 2004. SMEX02: Field scale variability, time stability and similarity of soil moisture. Remote Sens. Environ. 92, 436-446.
https://doi.org/10.1016/j.rse.2004.02.017

1118 Kelly, B.C., 2007. Some Aspects of Measurement Error in Linear Regression of Astronomical Data. Astrophys. J. 665, 1489-1506. https://doi.org/10.1086/519947

1120 Kerr, Y.H., Waldteufel, P., Wigneron, J.-P., Delwart, S., Cabot, F., Boutin, J., Escorihuela, M.-J., Font, J., Reul, N., Gruhier, C., Juglea, S.E., Drinkwater, M.R., Hahne, A., Martín-Neira, M., Mecklenburg, S., 2010. The SMOS Mission: New Tool for Monitoring Key Elements of the Global Water Cycle. Proc. IEEE 98, 666-687. https://doi.org/10.1109/JPROC.2010.2043032

1124 Khabbazan, S., Vermunt, P., Steele-dunne, S., Arntz, L.R., Marinetti, C., Valk, D. Van Der, Iannini, L., Molijn, R., 2019. Crop Monitoring Using Sentinel-1 Data: A Case Study from The Netherlands. Remote Sens. 11, 1-24. https://doi.org/10.3390/rs11161887

1127 Kim, Y., van Zyl, J.J., 2009. A Time-Series Approach to Estimate Soil Moisture Using Polarimetric Radar Data. IEEE Trans. Geosci. Remote Sens. 47, 2519-2527.
https://doi.org/10.1109/TGRS.2009.2014944

1130 Kittler, J., Illingworth, J., 1986. Minimum Error Thresholding. Pattern Recognit. 19, 41-47.

1131 Komorowski, M., Marshall, D.C., Salciccioli, J.D., Crutain, Y., 2016. Secondary Analysis of Electronic Health Records. Second. Anal. Electron. Heal. Rec. 1-427.
https://doi.org/10.1007/978-3-319-43742-2

1134 Lee, J.S., Hoppel, K., Mango, S.A., Miller, A.R., 1994. Intensity and phase statistics of multi-look polarimetric SAR imagery. IEEE Trans. Geosci. Remote Sens. 32, 1017-1028.
1136 https://doi.org/10.1109/igarss.1993.322213.

1137 Le Toan, T., Ribbes, F., Wang, L., Floury, N., Ding, K., 1997. Rice Crop Mapping and Monitoring Using ERS-1 Data Based on Experiment and Modeling Results. IEEE Trans. Geosci. Remote Sens. 35, 41–56.

1140 Lievens, H., Lannoy, G.J.M. De, Bitar, A. Al, Drusch, M., Dumedah, G., Franssen, H.H., Kerr, Y.H., Tomer, S.K., Martens, B., Merlin, O., Pan, M., Roundy, J.K., Vereecken, H., Walker, J.P., Wood, E.F., Verhoest, N.E.C., Pauwels, V.R.N., 2016. Assimilation of SMOS soil moisture and brightness temperature products into a land surface model. Remote Sens. Environ. 180, 292–304. https://doi.org/10.1016/j.rse.2015.10.033

1145 Macelloni, G., Paloscia, S., Pampaloni, P., Marliani, F., Gai, M., 2001. The relationship between the backscattering coefficient and the biomass of narrow and broad leaf crops. IEEE Trans. Geosci. Remote Sens. 39, 873–884.

1148 Macelloni, G., Paloscia, S., Pampaloni, P., Sigismondi, S., De Matthaeis, P., Ferrazzoli, P., Schiavon, G., Solimini, D., 1999. The SIR-C/X-SAR experiment on Montespertoli: Sensitivity to hydrological parameters. Int. J. Remote Sens. 20, 2597–2612.

1151 Mandelbrot, B. and van Ness, J.W. (1968) Fractional Brownian Motions, Fractional Noises and Applications. SIAM Review, 10, 422-437. http://dx.doi.org/10.1137/1010093

1153 Martinez-Fernandez, J., Ceballos, A., 2005. Mean soil moisture estimation using temporal stability analysis. J. Hydrol. 312, 28–38. https://doi.org/10.1016/j.jhydrol.2005.02.007

1155 Mattia, F., Satalino, G., Pauwels, V.R.N., De Lathauwer, E., 2011. GMES Sentinel-1 Soil Moisture Algorithm Development, Final Report, ESA ESTEC Contract no 4000101352/10/NL/MP/ef.

1157 McColl, K.A., Alemohammad, S.H., Akbar, R., Konings, A.G., Yueh, S., Entekhabi, D., 2017. The global distribution and dynamics of surface soil moisture. Nat. Geosci. 10. https://doi.org/10.1038/NGEO2868

1160 McNairn, H., Brisco, B., 2004. The application of C-band polarimetric SAR for agriculture: A review. Can. J. Remote Sens 30, 525–542.

1162 McNairn, H., Jackson, T.J., Wiseman, G., Bélair, S., Berg, A., Bullock, P., Colliander, A., Cosh, M.H., Kim, S.B., Magagi, R., Moghaddam, M., Njoku, E.G., Adams, J.R., Homayouni, S., Ojo, E.R., Rowlandson, T.L., Shang, J., Goita, K., Hosseini, M., 2015. The soil moisture active passive validation experiment 2012 (SMAPVEX12): Prelaunch calibration and validation of the SMAP soil moisture algorithms. IEEE Trans. Geosci. Remote Sens. 53, 2784–2801. https://doi.org/10.1109/TGRS.2014.2364913

1168 Mengen, D., Montzka, C., Jagdhuber, T., Fluhrer, A., Brogi, C., Baum, S., ... & Vereecken, H. (2021). The SARSense Campaign: Air-and Space-Borne C-and L-Band SAR for the Analysis of Soil and Plant Parameters in Agriculture. Remote Sensing, 13(4), 825.

1171 Merchant, C. J., Paul, F., Popp, T., Ablain, M., Bontemps, S., Defourny, P., ... & Mittaz, J., 2017. Uncertainty information in climate data records from Earth observation. Earth System Science Data. 9(2), 511–527.

1174 Mishra, A., Vu, T., Valiya, A., Entekhabi, D., 2017. Drought monitoring with soil moisture active passive ( SMAP ) measurements. J. Hydrol. 552, 620–632. https://doi.org/10.1016/j.jhydrol.2017.07.033

1177 Molero, B., Leroux, D.J., Richaume, P., Kerr, Y.H., Merlin, O., Cosh, M.H., Bindlish, R., 2018. Multi-Timescale Analysis of the Spatial Representativeness of In Situ Soil Moisture Data within Satellite Footprints. J. Geophys. Res. Atmos. 123, 3–21. https://doi.org/10.1002/2017JD027478

1181 Montzka, C., Cosh, M., Bayat, B., Bitar, A. Al, Berg, A., Bogena, H.R., Bolten, J.D., Cabot, F., Caldwell, T., Chan, S., Colliander, A., Crow, W., Das, N., Lannoy, G. De, Dorigo, W., Steven, R., Gruber, A., Hahn, S., Jagdhuber, T., Jones, S., Kerr, Y., Kim, S., Koyama, C., Kurum, M., Lopez-baeza, E., Mattia, F., McColl, K.A., Mecklenburg, S., Mohanty, B., Neill, P.O., Or, D., Petropoulos, G.P., Piles, M., Reichle, R.H., Rodriguez-fernandez, N., Rüdiger, C., Scanlon, T.,
1186 Schwartz, R.C., Spengler, D., Srivastava, P., Suman, S., Schalie, R. Van Der, Wagner, W., 1187 Wegmüller, U., Wigneron, J., 2020. Soil Moisture Product Validation Good Practices Protocol 1188 Version 1.0. Good Practices for Satellite Derived Land Product Validation, 123, 1189 https://lpvs.gsfc.nasa.gov/PDF/CEOS_SM_LPV_Protocol_V1_20201027_final.pdf.

1190 Montzka, C., Rötzer, K., Bogena, H.R., Nilda Sanchez, Vereecken, H., 2018. A New Soil Moisture 1191 Downscaling Approach for SMAP, SMOS, and ASCAT by Predicting Sub-Grid Variability. 1192 Remote Sens. 10, 1–25. https://doi.org/10.3390/rs10030427

1193 Moran, M.S., Alonso, L., Moreno, J.F., Cendrero Mateo, M.P., Fernando De La Cruz, D., Montoro, 1194 A., 2012. A RADARSAT-2 quad-polarized time series for monitoring crop and soil conditions 1195 in Barrax, Spain. IEEE Trans. Geosci. Remote Sens. 50, 1057–1070. 1196 https://doi.org/10.1109/TGRS.2011.2166080

1197 Nicolai-Shaw, N., Zscheischler, J., Hirschi, M., Gudmundsson, L., Seneviratne, S.I., 2017. A 1198 drought event composite analysis using satellite remote-sensing based soil moisture. Remote 1199 Sens. Environ. 203, 216–225. https://doi.org/10.1016/j.rse.2017.06.014

1200 Oliver, C., Quegan, S., 1998. Understanding Synthetic Aperture Radar Images. Artech House Bost. 1201 42, 479.

1202 Ouellette, J.D., Johnson, J.T., Balenzano, A., Mattia, F., Satalino, G., Kim, S.-B., Dunbar, R.S., 1203 Colliander, A., Cosh, M.H., Caldwell, T.G., Walker, J.P., Berg, A.A., 2017. A Time-Series 1204 Approach to Estimating Soil Moisture from Vegetated Surfaces Using L-Band Radar 1205 Backscatter. IEEE Trans. Geosci. Remote Sens. 55. 1206 https://doi.org/10.1109/TGRS.2017.2663768

1207 Palmisano, D., Mattia, F., Balenzano, A., Satalino, G., Pierdicca, N., Monti-Guarnieri, A.V., 2020. 1208 Sentinel-1 Sensitivity to Soil Moisture at High Incidence Angle and the Impact on Retrieval 1209 over Seasonal Crops. IEEE Trans. Geosci. Remote Sens. DOI:10.1109/TGRS.2020.3033887.

1210 Pan, M., and Wood, E. F., 2010. Impact of accuracy, spatial availability, and revisit time of satellite-1211 derived surface soil moisture in a multiscale ensemble data assimilation system. IEEE Journal 1212 of Selected Topics in Applied Earth Observations and Remote Sensing, 3(1), 49-56.

1213 Pan, F., and Peters-Lidard, C.D., 2008. On the Relationship Between Mean and Variance of Soil 1214 Moisture Fields. Journal of the American Water Resources Association (JAWRA) 44(1):235- 1215 DOI: 10.1111/j.1752-1688.2007.00150.x

1216 Paloscia, S., Pettinato, S., Santi, E., Notarnicola, C., Pasolli, L., Reppucci, A., 2013. Soil moisture 1217 mapping using Sentinel-1 images: Algorithm and preliminary validation. Remote Sens. 1218 Environ. 134, 234–248. https://doi.org/10.1016/j.rse.2013.02.027

1219 Pauwels, V.R.N., Verhoest, N.E.C., 2007. Optimization of a coupled hydrology – crop growth 1220 model through the assimilation of observed soil moisture and leaf area index values using an 1221 ensemble Kalman filter. Water Resour. 43, 1–17. https://doi.org/10.1029/2006WR004942

1222 Peng, J., Albergel, C., Balenzano, A., Brocca, L., Cartus, O., Cosh, M.H., Crow, W.T., Dabrowska-1223 zielinska, K., Dadson, S., Davidson, M.W.J., Rosnay, P. De, Dorigo, W., Gruber, A., 1224 Hagemann, S., Hirschi, M., Kerr, Y.H., Lovergine, F., Mahecha, M.D., Marzahn, P., Mattia, 1225 F., Pawel, J., Preuschmann, S., Reichle, R.H., Satalino, G., Silgram, M., Bodegom, P.M. Van, 1226 Verhoest, N.E.C., Wagner, W., Walker, J.P., Wegmüller, U., Loew, A., 2021. A roadmap for 1227 high-resolution satellite soil moisture applications – confronting product characteristics with 1228 user requirements. Remote Sens. Environ. 252, 112162. 1229 https://doi.org/10.1016/j.rse.2020.112162

1230 Peng, J., Loew, A., Merlin, O., Verhoest, N.E.C., 2017. A review of spatial downscaling of satellite 1231 remotely sensed soil moisture. Rev. Geophys. 55, 341–366. 1232 https://doi.org/10.1002/2016RG000543

1233 Picard, G., Toan, T. Le, Mattia, F., 2003. Understanding C-band radar backscatter from wheat 1234 canopy using a multiple-scattering coherent model. IEEE Trans. Geosci. Remote Sens. 41, 1235 1583–1591.
1236 Powers, J., 2021, private communication.

1237 Product User Guide, 2017. ESA Land Cover CCI v2.0. Tech. Rep.

1238 Pulvirenti, L., Squicciarino, G., Cenci, L., Boni, G., Pierdicca, N., Chini, M., Versace, C., Campanella, P., 2018. A surface soil moisture mapping service at national (Italian) scale based on Sentinel-1 data. Environ. Model. Softw. 102, 13–28.
https://doi.org/10.1016/j.envsoft.2017.12.022

1242 Quegan, S., Toan, T. Le, Yu, J.J., Ribbes, F., Floury, N., 2000. Multitemporal ERS SAR Analysis Applied to Forest Mapping. IEEE Trans. Geosci. Remote Sens. 38, 741–753.

1243 Quegan, S., Yu, J.J., 2001. Filtering of Multichannel SAR Images. IEEE Trans. Geosci. Remote Sens. 39, 2373–2379.

1246 Réjou-Méchain, M., Muller-landau, H.C., Detto, M., Thomas, S.C., Toan, T. Le, Saatchi, S.S., Barreto-Silva, J.S., Bourg, N.A., Bunyavejchewin, S., Butt, N., Brockelman, W.Y., Cao, M., Cárdenas, D., Chiang, J.-M., Chuyong, G.B., Clay, K., Condit, R., Dattaraja, H.S., Davies, S.J., Duque, A., Esufali, S., Ewango, C., Fernando, R.H.S., Fletcher, C.D., Gunatilleke, I.A.U.N., Hao, Z., Harms, K.E., Hart, T.B., Hérault, B., Howe, R.W., Hubbell, S.P., Johnson, D.J., Kenfack, D., Larson, A.J., Lin, L., Lin, Y., Lutz, J.A., Makana, J.-R., Malhi, Y., Marthews, T.R., McEwan, R.W., McMahon, S.M., McShea, W.J., Muscarella, R., Nathalang, A., Noor, N.S.M., Nytch7, C.J., Oliveira, A.A., Phillips, R.P., Pongpattananurak, N., Punchi-Manage, R., Salim, R., Schurman, J., Sukumar, R., Suresh, H.S., Suwanvecho, U., Thomas, D.W., Thompson, J., Uríarte, M., Valencia, R., Vicentini, A., Wolf, A.T., Yap, S., Yuan, Z., Zartman, C.E., Zimmerman, J.K., Chave, J., 1Laboratoire, 2014. Local spatial structure of forest biomass and its consequences for remote sensing of carbon stocks. Biogeosciences 11, 6827–6840. https://doi.org/10.5194/bg-11-6827-2014

1259 Rignot, E.J.M., Van Zyl, J.J., 1993. Change detection techniques for ERS-1 SAR data. IEEE Trans. Geosci. Remote Sens. 31, 896–906.

1261 Rodríguez-Fernández, N., Rosnay, P. De, Albergel, C., Aires, F., Prigent, C., Kerr, Y., 2019. SMOS Neural Network Soil Moisture Data Assimilation in a Land Surface Model and Atmospheric Impact. Remote Sens. 11, 1–23. https://doi.org/10.3390/rs11111334

1264 Rowlandson, T.L., Berg, A.A., Bullock, P.R., Ojo, E.R.T., McNairn, H., Wiseman, G., Cosh, M.H., 2013. Evaluation of several calibration procedures for a portable soil moisture sensor. J. Hydrol. 498, 335–344. https://doi.org/10.1016/j.jhydrol.2013.05.021

1267 Rubel, F., Brugger, K., Haslinger, K., Auer, I., 2017. The climate of the European Alps: Shift of very high resolution Köppen-Geiger climate zones 1800–2100. Meteorol. Zeitschrift 26, 115–125. https://doi.org/10.1127/metz/2016/0816

1270 Sabaghy, S., Walker, P., Renzullo, L.J., Jackson, T.J., 2018. Spatially enhanced passive microwave derived soil moisture: Capabilities and opportunities. Remote Sens. Environ. 209, 551–580. https://doi.org/10.1016/j.rse.2018.02.065

1273 Saich, P., Borgeaud, M., Member, S., 2000. Interpreting ERS SAR Signatures of Agricultural Crops in Flevoland, 1993–1996. IEEE Trans. Geosci. Remote Sens. 38, 1993–1996.

1275 Satalino, G., Balenzano, A., Mattia, F., Davidson, M.W.J., 2014. C-band SAR data for mapping crops dominated by surface or volume scattering. IEEE Geosci. Remote Sens. Lett. 11, 384–388.

1278 Seneviratne, S.I., Corti, T., Davin, E.L., Hirschi, M., Jaeger, E.B., Lehner, I., Orlowsky, B., Teuling, A.J., 2010. Investigating soil moisture-climate interactions in a changing climate: A review. Earth-Science Rev. 99, 125–161. https://doi.org/10.1016/j.earscirev.2010.02.004

1281 Starks, P.J., C.A. Fiebrich, D.L. Grimsley, J.D. Garbrecht, J.L. Steiner, J.A. Guzman, and D.N. Moriasi, 2014. Upper Washita River experimental watersheds: meteorologic and soil climate measurement networks. Journal of Environmental Quality 43:1239-1249.

1284 Smith, A.B., Walker, J.P., Western, A.W., Young, R.I., Ellett, K.M., Pipunic, R.C., Grayson, R.B., Siriwardena, L., Chiew, F.H.S., Richter, H., 2012. The Murrumbidgee Soil Moisture
1286 Monitoring Network data set. Water Resour. Res. 48, 1-6.
1287 https://doi.org/10.1029/2012WR011976.
1288 Thirumalai, K., Singh, A., Ramesh, R., 2011. A MATLAB™ code to perform weighted linear
1289 regression with (correlated or uncorrelated) errors in bivariate data. J. Geol. Soc. India 77,
1290 377-380. https://doi.org/10.1007/s12594-011-0044-1
1291 Tomer, S.K., Bitar, A. Al, Sekhar, M., Zribi, M., Bandyopadhyay, S., Sreelash, K., Sharma, A.K.,
1292 Corgne, S., Kerr, Y., 2015. Retrieval and Multi-scale Validation of Soil Moisture from.
1293 Remote Sens. 7, 8128-8153. https://doi.org/10.3390/rs70608128
1294 Torres, R., Davidson, M.W.J., Geudtner, D., 2020. Copernicus Sentinel Mission at C- and L-band:
1295 Current Status and Future Perspectives. 2020 IEEE Int. Geosci. Remote Sens. Symp. 4055-
1296 4058.
1297 Torres, R., Snoeij, P., Geudtner, D., Bibby, D., Davidson, M., Attema, E., Potin, P., Rommen, B.,
1298 Floury, N., Brown, M., Traver, I.N., Deghaye, P., Duesmann, B., Rosich, B., Miranda, N.,
1299 Bruno, C., L'Abbate, M., Croci, R., Pietropaolo, A., Huchler, M., Rostan, F., 2012. GMES
1300 Sentinel-1 mission. Remote Sens. Environ. 120, 9-24.
1301 https://doi.org/10.1016/j.rse.2011.05.028
1302 Tricomi, F.G., Erdelyi, A., 1951. The asymptotic expansion of a ratio of gamma functions. Pacific
1303 J. Math. 1, 133-142.
1304 Tsang, L., Kong, J.A., Ding, K.-H., 2001. Scattering of Electromagnetic Waves: Theories and
1305 Applications. John Wiley & Sons, Inc.
1306 Vereecken, H., Kollet, S., Simmer, C., 2010. Patterns in Soil - Vegetation- Atmosphere Systems :
1307 Monitoring, Modeling and Data Assimilation. Vadose Zo. J. 9, 821-827.
1308 https://doi.org/10.2136/vzj2010.0122
1309 Voronovich, A.G., 1994. Wave Scattering from Rough Surfaces. Part of the Springer Series on
1310 Wave Phenomena book series (SSWAV, volume 17). Springer-Verlag, Berlin, Heidelberg.
1311 https://doi.org/10.1007/978-3-642-97544-8
1312 Wagner, W., Hahn, S., Kidd, R., Melzer, T., Bartalis, Z., Hasenauer, S., Figa-Saldaña, J., De
1313 Rosnay, P., Jann, A., Schneider, S., Komma, J., Kubu, G., Brugger, K., Aubrecht, C., Züger, J.,
1314 Gangkofner, U., Kienberger, S., Brocca, L., Wang, Y., Blöschl, G., Eitzinger, J., Steinnocher,
1315 K., Zeil, P., Rubel, F., 2013. The ASCAT soil moisture product: A review of its specifications,
1316 validation results, and emerging applications. Meteorol. Zeitschrift 22, 5-33.
1317 https://doi.org/10.1127/0941-2948/2013/0399
1318 Wanders, N., Karssenberg, D., Roo, A. De, Jong, S.M. De, Bierkens, M.F.P., 2014. The suitability
1319 of remotely sensed soil moisture for improving operational flood forecasting. Hydrol. Earth
1320 Syst. Sci. 18, 2343-2357. https://doi.org/10.5194/hess-18-2343-2014
1321 Wang, C., Zuo, Q., Zhang, R., 2008. Estimating the necessary sampling size of surface soil
1322 moisture at different scales using a random combination method. J. Hydrol. 352, 309-321.
1323 https://doi.org/10.1016/j.jhydrol.2008.01.011
1324 Western, A.W., Blöschl, G., 1999. On the spatial scaling of soil moisture. J. Hydrol. 217, 203-224.
1325 https://doi.org/10.1016/S0022-1694(98)00232-7
1326 Yee, M.S., Walker, J.P., Monerris, A., Rüdiger, C., Jackson, T.J., 2016. On the identification of
1327 representative in situ soil moisture monitoring stations for the validation of SMAP soil
1328 moisture products in Australia. J. Hydrol. 537, 367-381.
1329 https://doi.org/10.1016/j.jhydrol.2016.03.060
1330 York, D., Evensen, N.M., Martínez, M.L., De Basabe Delgado, J., 2004. Unified equations for the
1331 slope, intercept, and standard errors of the best straight line. Am. J. Phys. 72, 367-375.
1332 https://doi.org/10.1119/1.1632486
1333 Zhu, L., Walker, J.P., Ye, N., Rüdiger, C., 2019. Roughness and vegetation change detection: A
1334 pre-processing for soil moisture retrieval from multi-temporal SAR imagery. Remote Sens.
1335 Environ. 225, 93-106. https://doi.org/10.1016/j.rse.2019.02.027
1336

1337

1338

1339

1340

1341

1342

## 1343 APPENDIX B: Supplementary material

### 1344 Dynamic masking of vegetation

1345 The method for dynamic masking of vegetation adopted in SMOSAR was suggested by Satalino 1346 et al. (2014), who proposed to use the C-band VH backscatter temporal signature to classify the 1347 agricultural areas into crops dominated by volume and soil attenuated scattering. The underlying 1348 hypothesis is that the higher the volume contribution, the higher the level of VH backscatter. For 1349 example, it has been shown that there exist crops for which the volume scattering mechanism is 1350 dominant since early phenological stages (e.g. winter rape), whereas for other crops (e.g. winter 1351 wheat) the attenuated soil scattering remains the dominant mechanism throughout the growing 1352 season (Cookmartin et al., 2000; Picard et al. 2003). Satalino et al. (2014) carried out a literature 1353 review about the C-band scattering mechanisms of various crops and, in parallel, an analysis of long 1354 and dense time series of C-band backscatter collected over agricultural sites in Europe and North 1355 America during the ESA AgriSAR’06 and ‘09 campaigns (Hajnsek et al., 2007; Caves et al., 2009). 1356 A first outcome of the analysis has been to confirm that crop canopies characterized by large leaves 1357 or random branching structure, predicted to be dominated by volume scattering at C-band (e.g. 1358 sugar beet, potato, maize, onion, garlic etc), show a VH backscatter temporal behavior similar to 1359 that of winter rape. Conversely, those crops with canopies characterized by small stems and then
dominated by attenuated surface scattering mechanisms at C-band (e.g., barley, oat, alfalfa, bean,
grass etc.) present a VH backscatter temporal behavior similar to that of winter wheat. Of course,
there exists a large intra-class variability for both classes of crops. Nevertheless, both Khabbazan et
al., (2019) and Palmisano et al., (2020) recently confirmed that the temporal behaviour of the S-1
VH backscatter of the two classes of crops remains substantially different.

As an example, Fig. B1 shows a time-series of S-1 VH γ coefficients (i.e. backscatter normalized
for the cosine of the incidence angle) collected over one wheat (blue) and one winter rape (orange)
fields over the Selhausen site (Germany) in 2018 (Mengen et al., 2021). The shaded areas
correspond to the gamma standard deviation estimated over the field. The cross-pol radar response
of winter rape increases during the growing season, while towards the end of the season, when the
crop canopy dries and then tends to be transparent to the SAR signal, the radar response of winter
rape drops again to the soil backscatter level. Conversely, the VH level of winter wheat backscatter
remains generally low (usually not exceeding the level of -14 dB) throughout the growing season
and the behaviour below -14 dB is erratic and likely related to changes in the level of soil moisture.

<figure><img src="image_20.png" /><figcaption>Fig. B1 Temporal behaviour of the mean value of S-1 VH γ coefficient for winter rape (red) and winter wheat (blue) fields over Selhausen (Germany) in 2018 (Mengen et al., 2021). The constant threshold at -14dB is also shown (black dashed line). The shaded areas represent the intra-field standard deviation.</figcaption></figure>
The segmentation proposed by Satalino et al., (2014) is an adaptive thresholding approach applied to S-1 VH gamma coefficient (i.e. $\gamma_{VH}$). The algorithm implements an iterative solution of the Kittler-Illingworth (KI) method (Kittler and Illingworth, 1986), which is an adaptive scheme seeking for the optimal separation of two classes statistically described as a mixture of two Gaussian pdfs. The iterative solution of the KI method is implemented because it can use a quasi-optimal guess value for the threshold (approximately -14 dB) that was identified in the experimental analysis.

The procedure starts from the initial guess threshold $T_0$ and splits the areas remaining after the ESA CCI masking in two sets (see Fig. B2), i.e. pixels with $\gamma_{VH}$ values higher (h) and lower (l) than $T_0$. The average ($\tilde{\mu}_i$), the standard deviation ($\tilde{\delta}_i$) and the number of pixels ($P_i$) ($i = h, l$) of each set are computed. They completely define the Gaussian pdf of the two classes, which approximates the histogram of data. Then, the estimate of the optimal threshold at iteration k=1 ($T_k$), represented by the crossover of the two Gaussian pdfs, is obtained by solving the following equation:

$$ \left[ \frac{T_k - \tilde{\mu}_h}{\tilde{\delta}_h} \right] + 2 \ln(\tilde{\delta}_h) - 2 \ln(P_h) = \left[ \frac{T_k - \tilde{\mu}_l}{\tilde{\delta}_l} \right] + 2 \ln(\tilde{\delta}_l) - 2 \ln(P_l), \quad k > 0 \quad (B1) $$

The procedure can be iterated until the threshold value becomes stable (i.e. $T_n$). Finally, pixels having backscatter greater or lower than $T_n$ are separated in two classes. The number of iterations required to reach the nearly optimal threshold depends on the goodness of the initial threshold $T_0$ and on the actual data pdf. As a result, only those land surfaces dominated by soil attenuated scattering are left unmasked, and it is for these surfaces that the $\Theta$ retrieval algorithm is applied.
<figure><img src="image_21.png" /><figcaption>Fig. B2. Scheme of thresholding method for volume and surface attenuated crop classes.</figcaption></figure>

# APPENDIX C: Supplementary material

## Fitting the coefficient of variation over the experimental sites

To characterize the coefficient of variation (CV_L) at each test site a time series of continuous Θ measurements recorded by the hydrological networks over various periods from one to four years (Table C1) was considered. For each test site daily averages and standard deviations of Θ, i.e., μ̃<sub>obs</sub> and δ̃<sub>Θ<sub>obs</sub></sub>, recorded by the ground stations were calculated. Additionally at the Apulian Tavoliere network, ten intensive ground campaigns were also carried out to extend the estimate of the Θ statistics over small areas, with a size of ~0.1 km × 0.1 km. The CV_L was then fitted with the curve in (13) and the k₁ and k₂ parameters identified at the various extents. Results are shown in Table C1. As an example, Fig. C1 shows the fitted curve δ̃<sub>Θ<sub>obs</sub></sub> = CV_L · μ̃<sub>obs</sub> versus μ̃<sub>obs</sub> for Yanco, TxSON, and Apulian Tavoliere in continuous lines. As well, δ̃<sub>Θ<sub>obs</sub></sub> measured only during the S-1 acquisitions in the time frame reported in Table C1 are shown as squared points. The plotted curves of δ̃<sub>Θ<sub>obs</sub></sub> have a convex upward shape and with levels increasing with the extent, as expected (Crow
et al., 2012; Famiglietti et al, 2008). For instance over the Apulian Tavoliere, which is the smallest site, $\delta_{\Theta_{obs}}$ is within 0.04 $m^3/m^3$ whereas over Yanco, the largest site, it may reach 0.12 $m^3/m^3$. An additional modulation of the relation between $\delta_{\Theta_{obs}}$ and $\tilde{\mu}_{obs}$ is due to the different dry out processes that depend on the local soil texture and land cover (Crow et al., 2012; Pan and Peters-Lidard, 2008).

<table><caption>Table C1. Fitting parameters of (13) between the coefficient of variation, $CV_L$, and $\tilde{\mu}_{obs}$ over the experimental sites. L is the extent of the site.</caption><thead><tr><th>Site</th><th>Station</th><th>AT</th><th>EC</th><th>LW</th><th>H</th><th>R</th><th>T</th><th>Y</th></tr></thead><tbody><tr><td>L (km)</td><td>0.1</td><td>1.6</td><td>17</td><td>25</td><td>30</td><td>35</td><td>36</td><td>60</td></tr><tr><td>k<sub>1</sub></td><td>0.563</td><td>0.324</td><td>0.713</td><td>0.933</td><td>0.905</td><td>0.905</td><td>0.534</td><td>0.601</td></tr><tr><td>k<sub>2</sub></td><td>-5.477</td><td>-2.623</td><td>-2.700</td><td>-4.258</td><td>-3.262</td><td>-3.421</td><td>-2.010</td><td>-2.018</td></tr><tr><td>period</td><td>Intensive campaigns</td><td>Feb14-Jan18</td><td>Jan14-Dec16</td><td>Jan15-Sep17</td><td>Jan15-Feb17</td><td>Jan15-Dec17</td><td>Jan15-Mar18</td><td>Jan15-Sep17</td></tr></tbody></table>

1423

<figure><img src="image_22.png" /></figure>

1424

<caption>Fig. C1. Examples of standard deviation, $\delta_{\Theta_{obs}}$, vs mean $\Theta$, $\tilde{\mu}_{obs}$, in the range [0.03-0.60] $m^3/m^3$ over Yanco (Y), TxSON (T) and Apulian Tavoliere (AT) derived from fitting parameters in Table C1. Only the $\delta_{\Theta_{ob}}$ and $\tilde{\mu}_{obs}$ values during the S-1 acquisitions in the time frame reported in Table C1 are shown in squared points.</caption>

1425

1426

1427