High Spatial and Temporal Soil Moisture Retrieval in Agricultural Areas Using Multi-Orbit and Vegetation Adapted Sentinel-1 SAR Time Series

Abstract

The retrieval of soil moisture information with spatially and temporally high resolution from Synthetic Aperture Radar (SAR) observations is still a challenge. By using multi-orbit Sentinel-1 C-band time series, we present a novel approach for estimating volumetric soil moisture content for agricultural areas with a temporal resolution of one to two days, based on a short-term change detection method. By applying an incidence angle normalization and a Fourier Series transformation, the effect of varying incidence angles on the backscattering signal could be reduced. As the C-band co-polarized backscattering signal is prone to vegetational changes, it is used in this study for the vegetational correction of its related backscatter ratios. The retrieving algorithm was implemented in a cloud-processing environment, enabling a potential global and scalable application. Validated against eight in-situ cosmic ray neutron probe stations across the Rur catchment (Germany) as well as six capacitance stations at the Apulian Tavoliere (Italy) site for the years 2018 to 2020, the method achieves a correlation coefficient of R of 0.63 with an unbiased Root Mean Square Error of 0.063 m³/m³.

1. Introduction

By using 60% to 90% of the total available water, agricultural systems are being the main consumers of fresh water resources on a global scale [1]. In this regard, they are heavily affected by the increasing impacts of climate change on the available water resources, while adding further pressure on it with the rising global demand for agricultural products, and the accompanying intensification of their production [2,3,4]. As irrigated agricultural areas produce nearly twice as much yield compared to rain-fed areas, the demand for irrigation water will consequently grow in the future. It is estimated that until 2080, the demand will increase by over 50% for developing countries and by approximately 16% for developed countries compared to the year 2000, while the largest change will occur in Africa and Latin America with 300% and 119%, respectively [5]. Detailed knowledge about soil moisture, being a key parameter in the agricultural sector, is therefore crucial for mitigating these effects on both local and regional scales. This information is used to estimate the efficiency of irrigation, determine its impact on local water resources, improve regional and local water allocation and hydrological models, and estimate the soil carbon storage capacity [6,7]. Nevertheless, high-resolution surface soil moisture data for regional and local monitoring (down to managing precision farming level) are still challenging to obtain [8].

By using Synthetic Aperture Radar (SAR) satellite missions, this knowledge gap can be filled [9,10,11]. Providing a cloud- and weather-independent monitoring of the Earth’s surface, SAR observations are suitable for regional and local soil moisture estimations, however, with a global extent [12]. While currently available soil moisture products from satellite missions such as Advanced Scatterometer (ASCAT) [13], Advanced Microwave Scanning Radiometer—Earth Observing System/2 (AMSR-E/AMSR2) [14,15], Soil Moisture and Ocean Salinity (SMOS) [16], and Soil Moisture Active Passive (SMAP) [17] are limited due to their coarse resolution, various studies evaluate the potential of downscaling these products using additional remote sensing or ground reference data [18,19,20,21,22]. Choi and Hur [18] used Moderate Resolution Imaging Spectroradiometer (MODIS) multi-spectral data to enhance the AMSR-E soil moisture product from 25 km to 1 km spatial resolution. For the successive mission, AMSR2, Fang et al. [23] used as well MODIS data, to produce a 1 km soil moisture product, while Piles et al. [24] and Molero et al. [19] used MODIS data to downscale SMOS 40 km resolution soil moisture product to a finer resolution of 1 km. However, the use of additional optical satellite data limits these methods to cloud-free conditions, which can lead to large gaps within monitoring, especially in tropical and subtropical regions. Using only cloud- and weather-independent radar data for downscaling, Das et al. [22] used Sentinel-1A and Sentinel-1B SAR data to enhance the SMAP radiometer resolution, increasing the resolution from 36 km to 1 km, even though the temporal resolution is degraded from 3 to 12 days due to sparse spatial overlapping of both missions. By taking advantage of recent machine learning approaches, SAR-based soil moisture estimation can be adapted to local, non-linear relationships between backscattering signals and soil moisture [25]. Hachani et al. [26] trained an artificial neural network (ANN) using Sentinel-1 and auxiliary in-situ data (e.g., soil moisture, soil texture, topography, vegetation measurements) to estimate soil moisture within an arid region with a resolution of 20 m. Datta et al. [27] compared different machine learning approaches (Support vector machine (SVM), K-Nearest Neighbor (KNN), and Random Forest (RF)) to estimate soil moisture over bare soils with a spatial resolution of 20 m. Using only Sentinel-1 backscatter time series, Bauer-Marschallinger et al. [28] presented a 1-km resolution soil moisture product, using a well-established change detection algorithm operated in the Sentinel-1 Data Cube at TU Wien. Using representative backscattering values from the driest and wettest conditions at a grid point, the relative soil moisture is retrieved here by scaling the individual backscattering value between these minimum and maximum values. Most recently, Balenzano et al. [29] published a volumetric soil moisture data set based on Sentinel-1 time series using a change detection method, with a spatial resolution of 1 km and a temporal resolution of 6 days.

However, current methods have limited applicability in agriculture due to their coarse spatial or temporal resolution [30]. When using physical backscatter models to decompose the SAR signal, the lack of publicly available full-polarimetric SAR observations often limits their use; in addition, they require a lot of surface information. For empirical models, their site dependence is an obstacle to their representation in other areas due to the large variability within the agricultural landscape [31]. In this regard, short-term change detection approaches are promising because they are independent of surface parameters that do not change significantly relative to soil moisture (e.g., soil roughness, soil texture, vegetation water content) and can be applied on the freely available dual-polarized C-band Sentinel-1 time series. Nevertheless, especially for agricultural areas with distinct bare soil and vegetation conditions, the influence of changing vegetation still limits their application. In this regard, new methods for addressing the effect of changing vegetational conditions need to be developed.

We developed a new method to address this issue by minimizing the time period between consecutive SAR observations. This is achieved by extending the alpha approximation approach by Balenzano et al., 2011 to multi-orbit acquisition time series. Moreover, we developed a method to detrend the vegetational effect on the related backscatter ratios. The method is established within the cloud-processing platform Google Earth Engine (GEE), using temporal and spatial high-resolution Sentinel-1 time series. Its integration into cloud computing environment allows a fast and applicable way for on-demand computation of soil moisture for individual time periods and areas of interest on a global scale.

2. Study Area

Two different study sites, the Rur catchment in western Germany and the Apulian Tavoliere site in southern Italy, were selected for evaluation. The Rur catchment is located in the federal state of North-Rhine Westphalia in the West of Germany, while small parts of it are also situated in Belgium and the Netherlands. It has a total area of 2354 km², with an elevation of 680 m above sea level (asl) in the south and 30 m asl in the north (Figure 1) [32]. The northern, flat parts of the catchment are dominated by crop cultivation, with mainly sugar beet, maize, and winter cereals [33]. In the southern region, located in the low mountain range Eifel, the main landcovers are pastures and forest of deciduous and coniferous trees [32]. Looking at the soil types within the Rur catchment, at the most northern parts, Fimic Anthrosols and Dystric Cambisols soils with a loamy sand texture are present, while the other northern parts are characterized mostly by Cumulic Anthrosols and Haplic Luvisols with both silt loam textures. In the southern regions, Fluvisols, Gleysols, Eutric Cambisols, and Stagnic Gleysols with a silt loam texture are present [33]. While the mean field capacity in the northern part is greater than 200 mm, the soils in the southern part are less productive with field capacities ranging between 50 and 150 mm [32]. The mean annual precipitation ranges between 650 mm to 850 mm in the northern part. As the Rur catchment is located in the west wind zone, precipitation shadowing effects caused by the High Fen Mountains can be observed in the southeastern part of the catchment, causing lower annual precipitation of around 700 mm compared to the 1200 mm present in western luv site of the catchment [32,34]. In the southern Eifel region, most of the rainfall occurs in winter when evapotranspiration rates are at their lowest, which results in a large proportion of precipitation rapidly becoming runoff. Contrary to this, the northern part of the Rur catchment receives the highest rainfall during the summer, when potential evapotranspiration rates are high, which can effectively buffer precipitation events [33]. Previous studies revealed lower mean soil moisture values in the northern part of the Rur catchment, which could be attributed to the lower precipitation rates. The observed spatial variability in this region, being dominated by agricultural practices/farming and a very flat terrain, were found to be driven by vegetation, due to the temporal differences of evapotranspiration rates. Furthermore, varying management practices (e.g., planting and harvesting dates), and soil parameters (e.g., porosity and field capacity) were seen as the main reasons. In the southern region, higher mean soil moisture values were observed. Here, the highest variation occurs within sub-catchments dominated by forest stands, due to the heterogeneous topography and higher precipitation levels. For grassland areas, which are also characteristic for the southern part, also high mean soil moisture values, but with a much lower standard deviation could be observed [33]. During the observation period from 2018 to 2020, a long persisting drought period was present over central Europe, also leading to soil moisture deficits in the Rur catchment. For further information, we refer to the European Drought Observatory reports as well as the monthly drought monitoring from Helmholtz Centre for Environmental Research [35,36]. The Rur catchment also includes the Selhausen test site, which is due to the very good availability of in-situ data used as a validation site for several air- and space-borne microwave soil moisture products and listed as a super-site by the Committee on Earth Observation Satellites Land Product Validation Subgroup (CEOS LPV) [37,38,39,40,41,42].

The second study area is located at the Segezia experimental farm (41°22′16″N 15°29′30″E) with an area of approximately 4 km², located in the Cervaro basin (Northwest part of the Puglia region, Southern Italy) (Figure 2). The main crops are cereals (durum wheat, barley, oat) as well as pastures (natural or annual grass of mixed cereals and legume crops). With its semiarid Mediterranean climate, the summer is hot and dry while the winter is short and temperate. The annual precipitation is around 550 mm, being mostly concentrated in winter months. The topography is characterized by a gentle slope, ranging between 130 m and 200 m asl. The soil texture is sandy clay loam [43]. The study area is part of the greater Capitana plain with an area of approximately 4000 km², located around the city of Foggia and is characterized by intensive and irrigated agricultural management by farms with an average size of 20 ha.

3. Data

3.1. Sentinel-1 C-Band SAR

A total of 709 dual-polarized (VV + VH) Sentinel-1A and Sentinel-1B scenes for the Rur catchment (Figure 3) and 361 scenes for the Apulian Tavoliere site are used for the period between January 2018 to December 2020, covering all or parts of the study area. The polar-orbiting twin satellite constellation monitors the Earth’s surface, using 5400 MHz C-band SAR instruments, while following the same orbital plane 180° apart. In this regard, a combined revisit time of 6 days can be achieved [44,45]. Depending on the latitude, even denser time series are possible, as multiple orbits can cover a given area. In this study, four orbits (Orbit 139, Orbit 88, Orbit 15, Orbit 37) are covering the Rur catchment area, recording it in alternating incidence angles (Desc.: 43.8°, Asc.: 30.2°, Asc.: 35.4°, Desc.: 39.5°). Hence, a temporal resolution of one to two days is achieved here. For the Apulian site, two orbits are used (Orbit 124, Orbit 146), recording in alternating incidence angles (Desc.: 38.9°, Asc.: 32.3°), achieving a temporal resolution of five to two days. The scenes are in Interferometric Wide-Swath Mode (IW) and Ground Range Detected High Resolution (GRDH) format [46]. The IW mode captures a 250 km swath with a spatial resolution of 5 m by 20 m by combining three sub-swaths, using Terrain Observation with Progressive Scans SAR (TOPSAR). For GRDH product, the scenes are resampled into 10 m by 10 m pixel spacing. The combined time series from Sentinel-1A and Sentinel-1B are obtained and processed using the Google Earth Engine (GEE) web platform [47]. The platform provides Sentinel-1A and Sentinel-1B scenes, which are preprocessed by the following steps, using the Sentinel-1 Toolbox SNAP software: thermal noise removal, radiometric calibration, terrain correction using Shuttle Radar Topography Mission (SRTM) Version 3.0 Global 1 arc second dataset (SRTMGL1) [48,49]. The backscattering intensities are acquired as Digital Numbers (DN).

3.2. CORINE Land Cover Data

The Coordination of Information on the Environment Land Cover (CORINE) dataset provides a land cover classification for Europe with a resolution of 100 m, grouping the land types into five major groups: artificial surfaces (I), agricultural areas (II), forests and semi-natural areas (III), wetlands (IV), and water bodies (V), which are further characterized within 44 individual subcategories in total. The minimum mapping unit is 25 ha with a mapping accuracy of greater than 85%. The dataset is based on Sentinel-2 and Landsat-8 image classification using computer-assisted visual image interpretation tools (CAPI) [50,51]. Focusing on agricultural areas for estimating soil moisture, all areas within the group “agricultural areas” are considered, while all other categories were masked. Thereby, all subcategories within the category “agricultural land” are considered as one class. The CORINE land cover dataset only distinguishes between individual cultivation types when there is long-term land use (e.g., olive trees, fruit trees, rice fields, pastures). Areas with crop rotations are placed in the subcategory “complex cropping patterns”. In this respect, no information on the specific crop type is available for this study.

3.3. OpenLandMap Data

The OpenLandMap Dataset is a global gridded dataset, providing soil information with a spatial resolution of 250 m for six different soil depths (0, 10, 30, 60, 100 and 200 cm), based on 350,000 observations from United States Department of Agriculture (USDA) database in combination with ancillary remote sensing data using a machine learning approach [52]. In this study, information about the percentual sand and clay content is used for the inversion of soil moisture to dielectric constant and vice versa, using the formula developed by Hallikainen [53]. Furthermore, field capacity (volumetric %) at 33 kPa for the 0 cm soil depth is used, for serving as an initial soil moisture starting point. For a detailed overview, the soil properties for different soil depths are listed in

Table 1. Overview of CRNP stations.

			Soil Depth [m]	Clay [%]	Sand [%]	SOC [g/kg]	Bulk Density [kg/m³]
Rotation (mainly Sugar beet, Potato, Maize, Cereals)	Aachen	RU_C_006	0	24.5	18.7	23	1276.3
			0.10	24.6	18.4	20	1280.0
		6.0275, 50.7985	0.30	26.3	18.8	10	1425.9
			0.60	28.9	18.3	5	1490.7
	Gevenich	RU_BCK_002	0	22.8	21.1	15	1312.2
			0.10	22.9	21.2	13	1323.4
		6.3235, 50.9892	0.30	25.6	21.3	7	1420.5
			0.60	27.7	20.9	2	1482.2
	Merzenhausen	ME_BCK_001	0	15.9	23.6	20	1306.1
			0.10	16.2	23.6	15	1349.3
		6.2974, 50.9303	0.30	17.4	23.4	5	1453.2
			0.60	18.2	24.1	4	1494.3
	Selhausen	SE_C_001	0	17.1	20.2	11	1315.7
			0.10	17.1	20.2	12	1321.9
		6.4471, 50.8659	0.30	19.2	20.8	7	1458.2
			0.60	20.6	20.9	0	1497.2
Meadow	Kall	RU_C_005	0	25.0	22.9	36	1222.7
			0.10	25.0	22.8	38	1236.2
		6.5264, 50.5013	0.30	26.8	23.3	10	1365.6
			0.60	29.4	23.0	7	1414.9
	Kleinhau-Hürtgenwald	RU_C_007	0	18.5	36.9	43	1093.0
			0.10	18.4	36.6	42	1123.9
		6.3720, 50.7224	0.30	19.1	37.1	20	1227.4
			0.60	20.8	39.1	9	1346.4
	Rollesbroich	RO_C_001	0	19.2	35.2	43	1139.6
			0.10	19.4	35.2	40	1153.5
		6.3042, 50.6219	0.30	19.7	35.6	24	1271.7
			0.60	21.2	36.7	10	1393.9
	Schönes-eiffen	RU_BCDKR_001	0	22.6	33.6	59	1054.1
			0.10	22.7	33.3	60	1095.9
		6.3755, 50.5149	0.30	24.3	34.5	27	1179.1
			0.60	25.3	35.8	16	1357.2

and

Table 2. Overview of Apulian Tavoliere site.

			Soil Depth [m]	Clay [%]	Sand [%]	SOC [g/kg]	Bulk Density [kg/m³]
Wheat	Apulian Tavoliere	6.0275, 50.7985	0	17.8	41.2	87	958.5
			0.10	17.9	41.2	82	1038.9
			0.30	18.7	41.4	24	1117.4
			0.60	19.8	43.8	13	1307.9

for both Rur and Apulian test site. As the proposed soil moisture retrieval approach is implemented within the Google Earth engine, the OpenLandMap dataset already available within it offers a good compromise between high spatial resolution and global coverage. However, this choice does not mean that regional or local soil datasets with higher resolution or accuracy may not also be suitable.

3.4. Global Land Data Assimilation System (GLDAS)

The Global Land Data Assimilation System (GLDAS) is a near-real terrestrial modeling system, incorporating ground-based and spaceborne observations to provide information on surface states and fluxes on a global scale [54]. It was developed by the National Aeronautics and Space Administration (NASA) Goddard Space Flight Center and the National Oceanic and Atmospheric Administration National Centers for Environmental Prediction (NCEP). In this study, the GLDAS Noah Land Surface Model L4 3 hourly 0.25 × 0.25 degree V2.1 (GLDAS_NOAH025_3H_2.1) product was used, which is providing data from January 2000 to present with a spatial resolution of 0.25° [55]. The model uses National Oceanic and Atmospheric Administration (NOAA)/Global Data Assimilation System (GDAS) atmospheric analysis fields [56], the disaggregated Global Precipitation Climatology Project (GPCP) V1.3 Daily Analysis precipitation fields [57,58], and the Air Force Weather Agency’s Agricultural Meteorological modeling system (AGRMET) radiation fields as forcing input variables [59].

3.5. Cosmic-Ray Neutron Probe (CRNP) Stations

Eight cosmic-ray neutron probe (CRNP) stations are spatially distributed on agricultural areas within the Rur catchment to monitor soil moisture. Measuring hourly neutron intensity, the stations are using either CRS-1000 or CRS-2000/B cosmic-ray neutron probe (Hydroinnova LLC, Albuquerque, NM, USA), providing soil moisture information within a circular area of 150 m to 200 m diameter around the CRNP within a soil depth between 15 cm to 70 cm. The CRS-1000 probes are filled with ³He gas, while the CRS-2000/B probes are filled with ¹⁰BF₃ enriched gas [60]. Since the neutrons are partially absorbed by the gas as they pass through the probes, electrons are produced that are attracted to an anode. The resulting electric currents are then amplified, detected, and counted by a pulse module [61]. For soil moisture estimation, the CRNP probes were calibrated using the N₀ method [62], using a one-time representative in-situ soil moisture measurement within the CRNP footprint. Even though the method does not account for biomass correction, it gives reasonable results with an RMSE $\le$0.033 cm³/cm³ [63]. In this study, the daily average was used for comparison.

3.6. Capacitance Stations

In February 2014, eleven stations equipped with an EM50 datalogger were installed over the Segezia site, with an average spacing of approximately 500 m (Figure 2). They record volumetric soil moisture values every 15 min using two 5TM and two 10HS probes, which are installed horizontally at 0.025 m, 0.10 m, 0.20 m and 0.40 m depths. Furthermore, the 5TM device also measures soil temperature. The capacitance sensors were calibrated by thermogravimetric method collecting soil moisture samples throughout a year, and it is periodically controlled. Moreover, intensive campaigns collecting vegetation and soil parameters in coincidence with radar and optical acquisitions (e.g., PRISMA, COSMO-SkyMed, RADARSAT-2, Sentinel-1/2, 9LSO-2, SAOCOM) are frequently carried out. In this study, the capacitance stations at points 2, 3, 5, 7, 9, and 10 are used, having the most consistent and continuous data for the years 2018 to 2020. For validation, soil moisture values from 0.025 m and 0.10 m are used, matching with each orbit overpass time.

4. Methods

The soil moisture estimation workflow is completely implemented within the cloud-processing platform Google Earth Engine (GEE), allowing to process CPU-intense time series data sets on a catchment scale in considerably shorter time. The workflow can be divided into two major parts: (1) preprocessing of Sentinel-1 SAR data and (2) soil moisture estimation using the alpha approximation method (Figure 4). As the algorithm is constructed in a customizable way, both the observation period as well as the study region can be changed at a later stage. The output can be both visualized directly within the browser-based GEE platform and downloadable for local post-processing. In this study, the output was resampled to a resolution of 200 m, thus aggregating all introduced input grids (Sentinel-1, CORINE, OpenLandMap, GLDAS) into the same 200 m × 200 m pixel grid by spatial averaging. The Sentinel-1 data was down-sampled, while the coarser datasets (CORINE, GLDAS, OpenLandMap) were up-sampled. In this regard, the down-sampling of Sentinel-1 data to a resolution of 200 m further reduced the effect of speckle noise in the retrieval and averaged out the spatial heterogeneity, as a coarser scale reduces the effect of heterogeneous surface parameters (e.g., soil roughness, plant canopy structure and vegetation water content). On the other hand, the relationship with soil moisture is strengthened as it adjusts more uniformly across space [64]. Nevertheless, the different original resolutions of the used data sets can have varying influence on the resulting soil moisture. Regarding the CORINE data set, both an unnecessary loss of “valid” as well as not correctly masked Sentinel-1 pixels in the spatial domain are possible. In terms of the coarser GLDAS and OpenLandMap data set, multiple Sentinel-1 pixels are assigned to the same soil moisture starting and soil moisture range value, depreciating the correlation to the spatial soil moisture distribution that can possibly be achieved by using high-resolution Sentinel-1 time series.

4.1. Preprocessing

4.1.1. Masking

Two masking steps precede the workflow. In a first step, a land cover mask is applied to each Sentinel-1A and Sentinel-1B scene, as this study focuses on soil moisture estimation within agricultural areas. Using a binary mask created from the CORINE Land Cover Data, where non-agricultural areas are set to zero and agricultural areas are set to one, masked pixels are excluded from further processing, reducing the amount of computational resources in the first place. In this regard, the method is not applied to areas, where C-band is not able to penetrate through the vegetation cover and would give unreasonable results, e.g., forest areas [65]. In the next step, a soil freeze and snow cover filter were applied, as both events influence the retrieved dielectric constant and lead to wrong soil moisture estimations using a short-term change detection method. The snow filtering was done using the MODIS Terra Snow Cover Daily Global dataset [66], while the freeze filtering was done using the COM-C Land Surface Temperature (V2) dataset [67], filtering out scenes, where the air temperature at 2 m height is below 3 °C, as soil freezing is possible [68]. The filtered pixels are replaced by calculating the temporal mean between 15 days before and 15 days after the snow or freeze event, to ensure a continuous time series also during the winter period. Since heavy convective precipitation events could also potentially affect the backscattering signal, there could be theoretical possibility that a change in backscatter is related to such specific meteorological conditions. As these effects usually only occur under heavy convective storm events (with over 50 mm/h), these conditions only prevail in approximately 0.1% of areas on a global scale, with the majority of these areas located in the tropics [69]. Since the likelihood of such events occurring during SAR observations in both study areas is even lower, it is not considered in this study. However, if the method is applied on a global scale, it could become a crucial factor, particularly within the tropics.

4.1.2. Spatial Filtering

To reduce the influence of speckle on the change detection-based soil moisture estimation, a denoising filter was applied to the co-polarized band of each Sentinel-1 SAR scene. While improving the image quality, changes in backscattering signals between two consecutive images are becoming more related to changes in soil moisture, rather than due to the speckle effect. To preserve the edges of individual agricultural fields, a focal median filter with a 50 m radius circular kernel (according to the area of around 78 resolution cells) was applied to all scenes in this study. For each pixel, it calculates the median from its surrounding pixels values and rewrites it into the original pixel. While its fast performance is suitable for processing large time series, its better ability to preserve the edges of objects makes it more suitable for field scale soil moisture estimation compared to the focal mean filter [70].

4.1.3. Incidence Angle Normalization

As the images within the time series were recorded with multiple incidence angles, they must be normalized to reduce the effect of incidence angle on backscattering signal. In this regard, the SAR images are normalized to a reference incidence angle of 40°, using scene-based linear regression [71]:

$${\sigma }_0^{norm}={\sigma }_0^i-\beta \left({\theta }_i-{\theta }_{ref}\right)$$

(1)

with ${\sigma }_0^{norm/i}$ being the normalized and original backscattering signal, ${\theta }_{ref/i}$ being the reference and original incidence angle and $\beta$ as the slope derived from linear regression between ${\sigma }_0^i$ and ${\theta }_i$. An additional correction of azimuthal angle was not considered, as azimuthal anisotropy is mainly caused by the orientation of steep slopes, as present in higher and lower mountain ranges. Focusing on agricultural areas, predominantly located on flat surfaces, the influence of the azimuth angle is substantially lower [71]. Via the investigation of the influence of azimuthal angle over winter wheat fields, it was found negligible within a dense time series [72].

4.1.4. Fourier Transformation

Incidence angle normalization adjusts the backscattering signal within the individual scenes but does not address the different backscatter distributions of the four orbits, caused by the varying influence of soil and plant parameters as well as crop row orientation on the incidence angle effect [73,74,75]. Therefore, the time series still contains backscatter changes unrelated to soil moisture, which need to be removed. Knowing the frequency of the signal caused by the changing sequence of incidence angles in the backscatter time series and determined by the Sentinel-1 revisit time of 12 days, the time series of ${\sigma }_0^{norm}$ is transformed into a Fourier Series excluding these frequencies by:

$${\sigma }_0^F\left(\mathrm{t}\right) ~ \frac{a_0}{2}+{{\displaystyle \sum }}_{k=1}^n(a_k\cos \left(k\omega t\right)+b_k\sin \left(k\omega t\right))$$

(2)

with

$$\omega =\frac{2\pi }{T}$$

(3)

$$a_0=\frac{2}{T}{{\displaystyle \int }}_{-T/2}^{T/2}\mathrm{f}\left(t\right)dt$$

(4)

$$a_k=\frac{2}{T}{{\displaystyle \int }}_{-T/2}^{T/2}\mathrm{f}\left(t\right)\cos \left(k\omega t\right)dt$$

(5)

$$b_k=\frac{2}{T}{{\displaystyle \int }}_{-T/2}^{T/2}\mathrm{f}\left(t\right)\sin \left(k\omega t\right)dt$$

(6)

For each year, the combined backscattering time series was transformed into a Fourier Series, with $T$= 365 and $n$= 24, resulting in excluding signals with a frequency smaller than 15 days.

4.1.5. Vegetation Correction

In terms of the first-order radiative transfer model, the total backscattering signal can be written as the sum of the bare soil component and the vegetation affected scattering component:

$${\sigma }_0^F= \left(f_{b }\times {\sigma }_{soil}^F\right)+\left(1-f_{b }\right)\times \left({\sigma }_{veg}^F+{\gamma }^2{\sigma }_{soil}^F+{\sigma }_{int}^F\right)$$

(7)

with ${\sigma }_{soil}^F$, ${\sigma }_{veg}^F$ and ${\sigma }_{int}^F$ being the soil surface, vegetation and vegetation—surface interaction backscattering component, respectively, ${\gamma }^2$ being the vegetation attenuation factor and $f_{b }$ being the bare soil fraction [76]. Applying a short-term change detection approach, the ratio of two consecutive backscatter signals are as follows:

$$\frac{{\sigma }_0^{F2}}{{\sigma }_0^{F1}}=\frac{(f_{b }\times {\sigma }_{soil}^{F2})+\left(1-f_{b }\right)\times \left({\sigma }_{veg}^{F2}+{\gamma }^2{\sigma }_{soil}^{F2}+{\sigma }_{int}^{F2}\right)}{(f_{b }\times {\sigma }_{soil}^{F1})+\left(1-f_{b }\right)\times \left({\sigma }_{veg}^{F1}+{\gamma }^2{\sigma }_{soil}^{F1}+{\sigma }_{int}^{F1}\right)}$$

(8)

For areas dominated by agricultural crop production, the surface changes from bare soil conditions to fully vegetated conditions during one growing cycle. For agricultural areas, this leads to different equations during bare soil and vegetation period, assuming that the bare soil fraction changes from 0 to 1:

$$\frac{{\sigma }_{baresoil}^{F2}}{{\sigma }_{baresoil}^{F1}}=\frac{ {\sigma }_{soil}^{F2}}{ {\sigma }_{soil}^{F1}}$$

(9)

and

$$\frac{{\sigma }_{vegetation}^{F2}}{{\sigma }_{vegetation}^{F1}}=\frac{{\sigma }_{veg}^F+{\gamma }^2{\sigma }_{soil}^{F2}+{\sigma }_{int}^F}{{\sigma }_{veg}^F+{\gamma }^2{\sigma }_{soil}^{F1}+{\sigma }_{int}^F}$$

(10)

In this regard, a change in soil moisture leads to different changes in the total backscatter signal when comparing bare soil and vegetation periods:

$$\frac{ {\sigma }_{soil}^{F2}}{ {\sigma }_{soil}^{F1}} \ne \frac{{\sigma }_{veg}^F+{\gamma }^2{\sigma }_{soil}^{F2}+{\sigma }_{int}^F}{{\sigma }_{veg}^F+{\gamma }^2{\sigma }_{soil}^{F1}+{\sigma }_{int}^F}$$

(11)

Here, the related parameters of bare soil fraction $f_{b }$, vegetation attenuation ${\gamma }^2$, as well as vegetation and vegetation-surface backscatter component ${\sigma }_{veg}^{Fourier}$ and ${\sigma }_{int}^{Fourier}$, causing the difference in backscatter sensitivity of the total backscatter signal, are all related to vegetation. For simplification, we therefore assume that a systematic change in the backscatter ratios between bare soil and vegetation conditions can be observed and described by a linear relationship using an appropriate estimator for the transition between bare soil and vegetation conditions. Comparing the backscatter signals from both agricultural areas as well as meadow areas for the Rur catchment, having similar temporal soil moisture dynamics and soil moisture ranges the intensity of the co-polarized backscatter signal was chosen as an estimator for differentiating the influence of vegetation on the related backscatter ratios. As shown in Figure 5, the backscatter signal from agricultural areas ranges between 0.025 and 0.200, with significantly higher values within the bare soil period, compared to a range of 0.025 and 0.120 for the meadow stations. Especially during the shooting-phase in spring, the increasing vegetational effect on the backscatter signal leads to sharply declining total backscattering values, becoming partly uncorrelated to the temporal behaviour of soil moisture. Looking at the agricultural areas, this results in overall minima in backscattering values within the period between March and May, where the influence of changing vegetation is high on the total backscattering signal, while the overall minimum of in-situ measured soil moisture is occurring during August and September. As previous studies showed, the range of backscatter changes decreases with increased vegetation. Indicating a decreasing sensitivity of the backscattering signal to changes in soil moisture with growing vegetation, this relationship can be described by an empirical function [77,78,79]. As the vegetation influence on the backscattering ratios can vary quite significantly, depending on the crop type and the vegetation parameters (e.g., plant height, plant geometry/structure, biomass, vegetation water content) [42,80], we choose a linear regression on pixel-basis to remove this trend by:

$${\left(\frac{ {\sigma }_0^{F2}}{ {\sigma }_0^{F1}}\right)}^{detr.}= \left(\frac{ {\sigma }_0^{F2}}{ {\sigma }_0^{F1}}\right)-m\left( {\sigma }_0^{F1}-{\sigma }_{mean}^{F1}\right)$$

(12)

with $m$ and ${\sigma }_{mean}$ being the regression slope and mean backscattering value. Applying a pixel-based linear regression, no information about the specific crop type is used. Nonetheless, it is assumed that all crops affect the corresponding backscatter ratios during the transition from bare soil to fully vegetated conditions, and this is the dominating factor for a systematic change in the backscatter ratios during the season, as vegetation sensitive C-band is used. By resampling pixels to a coarser resolution of 200 m, each pixel furthermore contains different plant types and their individually mixed vegetation-induced attenuation behavior. The linear regression is only applied to stations with distinct bare soil and vegetation conditions during the growing season, with a yearly backscattering range greater than 0.10, which was empirically chosen for the study area.

4.2. Soil Moisture Estimation

4.2.1. Alpha Approximation

By using a change detection approach on a dense time series of SAR data, the changes in backscattering signal between two consecutive observations can be related to a change in soil moisture, assuming that the temporal variability of this parameter is highest compared to other surface parameters, such as soil roughness, vegetation height, vegetation water content, and vegetation geometry, and the recording parameters frequency and incidence angle are constant within the time series [81]. Thus, the ratio of two consecutive co-polarized backscattering signals $\left({\sigma }_0^{\left(2\right)}/{\sigma }_0^{\left(1\right)}\right)$ can be expressed as function of soil dielectric constants $\left({\epsilon }^{\left(1,2\right)}\right)$ and incidence angle $\left(\mathsf{\theta }\right)$ [82]:

$${\left(\frac{ {\sigma }_0^{F2}}{ {\sigma }_0^{F1}}\right)}^{detr.}\approx {\left|\frac{{\alpha }_{PP}^{\left(2\right)}\left({\epsilon }^{\left(2\right)},\theta \right)}{{\alpha }_{PP}^{\left(1\right)}\left({\epsilon }^{\left(1\right)},\theta \right)}\right|}^2$$

(13)

with

$$\left|{\alpha }_{VV}\left(\epsilon ,\theta \right)\right|=\left|\frac{\left(\epsilon -1\right)\left(si{n}^2\theta -\epsilon \left(1+si{n}^2\theta \right)\right)}{(\epsilon \cos \theta +\sqrt{\epsilon -si{n}^2\theta )²}}\right|$$

(14)

For solving the equation, it can be resolved to zero, giving [81]:

$$\left|{\alpha }_{PP}^{\left(1\right)}\left(\epsilon ,\theta \right)\right|-\sqrt{\frac{{\mathsf{\sigma }}_0^{\left(1\right)}}{{\mathsf{\sigma }}_0^{\left(2\right)}}}\times \left|{\alpha }_{PP}^{\left(2\right)}\left(\epsilon ,\theta \right)\right|=0$$

(15)

Applying this equation to all consecutive observation pairs of the SAR time series, an underdetermined linear equation system can be created, with N unknown variables and N−1 equations, depending on an initial soil moisture value to be analytically solvable [83]:

$$\left[\begin{array}{cccccc}1& -\sqrt{{\left(\frac{ {\sigma }_0^{F1}}{ {\sigma }_0^{F2}}\right)}^{detr.}}& 0& ⋮& 0& 0\\ 0& 1& -\sqrt{{\left(\frac{ {\sigma }_0^{F2}}{ {\sigma }_0^{F3}}\right)}^{detr.}}& ⋮& 0& 0\\ \cdots & \cdots & \cdots & ⋰& \cdots & \cdots \\ 0& 0& 0& ⋮& -\sqrt{{\left(\frac{ {\sigma }_0^{FN-2}}{ {\sigma }_0^{FN-1}}\right)}^{detr.}}& 0\\ 0& 0& 0& ⋮& 1& -\sqrt{{\left(\frac{ {\sigma }_0^{FN-1}}{ {\sigma }_0^{FN}}\right)}^{detr.}}\end{array}\right]\times \left[\begin{array}{c}\left|{\alpha }_{PP}^{\left(1\right)}\left({\epsilon }_1,\theta \right)\right|\\ \left|{\alpha }_{PP}^{\left(2\right)}\left({\epsilon }_2,\theta \right)\right|\\ \cdots \\ \left|{\alpha }_{PP}^{\left(N\right)}\left({\epsilon }_N,\theta \right)\right|\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ \cdots \\ 0\end{array}\right]$$

(16)

In addition, in this study, it is proposed to rescale the backscattering signal to dB, as this is expected to reduce possible biases in the final soil moisture estimates.

4.2.2. Soil Moisture to Dielectric Constant Inversion

As the alpha approximation method requires dielectric constant values ($\mathsf{ϵ})$, the initial in-situ volumetric soil moisture values $(m_v)$ need to be converted as well as the derived dielectric constant values from alpha approximation must be converted back into soil moisture. Using the soil texture class data set, this is done by applying the following equation by Hallikainen et al. [53]:

$$ϵ=\left(a_0+a_{1}S+a_{2}C\right)+\left(b_0+b_{1}S+b_{2}C\right)m_v+\left(c_0+c_{1}S+c_{2}C\right)m_v^2$$

(17)

with $a_{\mathrm{i}},b_i, {\mathrm{c}}_{\mathrm{i}}$ being frequency dependent constants, while $S$ and $C$ are the sand and clay percentages. In this study, we used constants for 5.3 GHz, matching Sentinel-1 C-band frequency of 5.4 GHz, provided by Chen et al. [84]. In the last step, the minimum and 95th percentile of retrieved soil moisture values are rescaled between the observed minimum and maximum of the region retrieved from GLDAS dataset.

5. Results and Discussion

The effect of the individual processing steps on both the backscattering time series as well as the resulting soil moisture estimations are presented and discussed in this chapter. Focusing on the main processing steps relevant for preparing the multi-orbit SAR time series for the short-term change detection approach under changing vegetational influence, the incidence angle normalization, Fourier Series transformation, and vegetation correction are discussed in detail. Finally, the soil moisture product is assessed in terms of its spatial and temporal accuracy over the Rur catchment and the Apulian Tavoliere test site. For the evaluation metrics, we applied both the squared Pearson correlation coefficient (R²) as well as unbiased Root Mean Square Error (uRMSE) as described by Gruber et al., 2020 [85].

5.1. Incidence Angle Normalization

Looking at the average backscattering signals for each orbit exemplarily over the Rur catchment, there is a clear trend visible of decreasing backscattering signals due to increasing incidence angles (

Table 3. Comparison between original median backscatter values and incidence angle normalized median backscatter values from individual orbits.

	Median Backscatter Value				Incidence Angle Normalized Median Backscatter Value
Orbit	88	15	37	239	88	15	37	239
VV	0.147	0.112	0.107	0.091	0.111	0.095	0.105	0.106

, Figure 6). The normalized signals do not show such a trend for the respective orbits. The incidence angle normalization reduces the range of observed median backscattering intensities related to changing incidence angles by approximately 70% for the co-polarized signal. Nevertheless, the probability distribution function shows that the different distributions of backscattering values, depending on the incidence angle, are still present. While the orbits 88, 37 and 139 show a stronger bimodal distribution in co-polarization, orbit 15 has a more down-skewed distribution. In this regard, the influence of both incidence angle as well as ascending and descending orbit is still present and is similar to previous studies [71,86]. As the influence depends on surface and vegetation parameters, e.g., surface roughness, vegetation biomass, canopy height, plant geometry and phenology, a crop specific normalization could improve the results. Nevertheless, the results imply further processing steps to level out the different backscattering signals for a short-term change detection approach.

5.2. Fourier Series Transformation

Although the time series was incidence angle normalized, co-polarized signals are still differently affected by the incidence angle over the season. Using the CRNP station Aachen (RU_C_006) as an example, the co-polarized signal is mostly affected during the non-vegetation period from October to March, as observed in previous studies [86], resulting in a soil moisture range of 0.75 m³/m³ for this period and 0.43 m³/m³ for the other months (Figure 7).

Looking at the in-situ soil moisture range of these periods, the lowest variability can be observed during the non-vegetation times from October to March with a mean soil moisture range of 0.120 m³/m³, while more variations occur during the remaining months with around 0.282 m³/m³. As the highest variability in estimated soil moisture is observable in the periods with the lowest in-situ measured soil moisture variability, it is likely caused by the stronger incidence angle effect on non-vegetated surfaces, as previously observed [87]. By calculating the Fourier Series from the normalized time series, the influence of alternating incidence angle on the estimated soil moisture is reduced, visible at the reduced variability looking at the related boxplots (Figure 7). To evaluate the sensitivity of the chosen minimum frequency, soil moisture is estimated from Fourier Series excluding frequencies higher than 10 days, 15 days, and 22 days. Here, both the squared Pearson correlation coefficient (R²) as well as unbiased Root Mean Square Error (uRMSE) [85] are improved due to the Fourier Series transformation. While the soil moisture estimation from incidence angle normalized backscattering time series achieves an R² of 0.10 with an uRMSE of 0.12 m³/m³, the estimated soil moisture from the Fourier Series transformed time series performs better with R² values between 0.20 to 0.23 and uRMSE between 0.07 m³/m³ and 0.08 m³/m³. Looking at the soil moisture range derived from Fourier Series transformed time series, for the non-vegetation period, the estimated soil moisture values have a range of 0.23 m³/m³, 0.24 m³/m³, and 0.22 m³/m³, while for the vegetation period they have a range of 0.22 m³/m³, 0.20 m³/m³, and 0.19 m³/m³. In this regard, all Fourier Series transformed backscattering time series are able to mitigate the effect of varying incidence angle. Especially during the non-vegetation period, the range could be reduced by up to 70%, even though it is still two times higher than observed from in-situ data. Regarding the sensitivity to the chosen maximum frequency, no significant influence on the retrieved soil moisture is observed.

By excluding frequencies higher than the revision time, some information about soil moisture is naturally lost in the backscatter time series. In principle, a diurnal cycle can be observed for soil moisture, with decreasing soil moisture one to two hours before sunrise and increasing soil moisture about two to four hours before sunset [88]. Since all available orbits are used, the original SAR time series includes soil moisture information from both morning (around 05:00 a.m. to 06:00 a.m.) and evening acquisitions (around 04:40 p.m. to 05:10 p.m.) for both study areas. Excluding high frequencies, information about the sub-daily dynamics of soil moisture is therefore filtered out. Nevertheless, the dynamics of changing vegetational water content and dew cycles during the day are also filtered out, which are also affecting the backscatter signal of C-band. Excluding high temporal frequencies also reduces the representation of short-term wetting and drying of topsoil due to precipitation events, while retaining weekly and seasonal patterns. Since the backscattering signal is not only affected by precipitation through soil moisture increase but also by intercepting water, these effects are also reduced. In doing so, a tradeoff can be made between losing temporal information about soil moisture and reducing backscatter changes unrelated to soil moisture. Thus, the exclusion of frequencies higher than the revisit time mainly links to the effect of the changing angle of incidence rather than the loss of soil moisture information. Looking at the power spectral density, calculated using Welch’s average periodogram [89], a small peak at 3 days frequency is observable, which is the period between ascending and descending scenes (Figure 7). The next increase in power compared to the in-situ measured soil moisture power spectral density curve is starting from 5 days frequencies onward, which is again the period between ascending and descending scenes, while it is also peaking at 6 days frequency, which is the period between matching orbits of Sentinel-1 A and B. As the power spectral density of the Fourier Transformed soil moisture time series is reduced for these frequencies, it is in line with filtering out backscattering signals related to alternating incidence angles.

5.3. Vegetational Correction

The fraction of bare soil and vegetation affected backscatter as well as the changes in the vegetation scatter contribution during the growing cycle lead to different correlations between the change of backscatter and change of soil moisture. Using the co-polarized backscatter signal as an estimator for bare soil and vegetation conditions during the growing season for agricultural areas, the consecutive backscatter ratios are adjusted through a linear regression (Section 4.1.4). In this regard, the effect of increasing vegetation scattering contribution within the total backscatter signal can be corrected. In Figure 8, the backscatter ratios from Fourier Series transformed and vegetational adapted time series of CRNP station SE_C_001 are shown. At the agricultural crop dominated SE_C_001 station, a small positive trend is present, indicating a decline in backscatter ratios caused by increasing vegetation scatter contribution. By removing this trend, the backscatter ratios related to lower backscatter value are increased, while the ratios related to higher backscatter values are decreased. By grouping the backscatter ratios into three vegetational periods of no-vegetation (January, February, October, November, and December in light grey), intermediate-vegetation (Mar, Apr, Aug, and Sep in dark grey), and full-vegetation cover (May, June, and July in black), most of the increase in backscatter ratios fall into the full-vegetation period. A similar correction can be observed within the non-vegetation period, leveling down the related backscatter ratios. Normalizing the backscatter ratios to the mean backscatter value, the ratios within the intermediate-vegetation period are not changed considerably.

Figure 9 shows the soil moisture from Fourier Series transformed and vegetation corrected backscatter time series, compared to the in-situ measured soil moisture for the CRNP station RU_C_006. In general, an increase in estimated soil moisture can be observed within the period from February to May, while a decrease is observed between June and October. From November onward, both soil moisture estimations are aligning. This leads to a better match between the in situ and estimated soil moisture and corrects for the relative underestimation of soil moisture in the first half of the year, dominated by growing vegetation, relative to the overestimation of soil moisture in the second half of the year, characterized by vegetation maturing and drying. In this respect, the temporal course of the measured and the estimated soil moisture are in better agreement. In particular, the estimated soil moisture minimum of the vegetation-corrected time series is in August and thus agrees with the values measured at the station and not in May as in the uncorrected soil moisture estimate. For this example, the soil moisture estimation without vegetation correction achieves an R² of 0.54 with an uRMSE of 0.044 m^3/m³, while the vegetation corrected ones lead to an increase in R² to 0.80 with an uRMSE of 0.029 m³/m³.

5.4. Effect of Individual Processing Steps on Soil Moisture Estimation

To evaluate the influence of the individual processing steps, the resulting soil moisture values are compared with the in-situ measured values from the eight CRNP stations within the Rur catchment as well as six capacitance station at 0.025 m and 0.1 m soil depth within the Apulian Tavoliere site. To relate the estimated soil moisture to the CRNP data and capacitance measurements, the mean estimated soil moisture around each location within a diameter of 200 m was calculated and was chosen to match the detection radius of the CRNP stations. We focus on the main processing steps that were used to create the temporally and spatially high-resolution time series from multi-orbit Sentinel-1 scenes. In this regard, the estimated soil moisture from incidence angle normalisation (IA), incidence angle normalisation and Fourier Series (IA + FS) and incidence angle normalisation, Fourier Series and vegetational adapted (IA + FS + VA) time series are compared (

Table 4. Comparison between overall soil moisture retrievals from incidence angle normalized (IA), IA and Fourier Series transformed (IA + FS), and IA + FS and vegetational adapted (IA + FS + VA) SAR time series for each catchment and soil depth, with better performance indicated by a darker background color.

		Mean R²			Mean uRMSE
	Test Site	IA	IA + FS	IA + FS + VA	IA	IA + FS	IA + FS + VA
2018	Rur	0.36	0.47	0.58	0.056	0.052	0.046
	Apulian Tavoliere 0.025 m	0.36	0.39	0.44	0.058	0.058	0.056
	Apulian Tavoliere 0.1 m	0.37	0.39	0.42	0.060	0.059	0.057
2019	Rur	0.27	0.48	0.55	0.058	0.045	0.042
	Apulian Tavoliere 0.025 m	0.15	0.16	0.49	0.081	0.081	0.059
	Apulian Tavoliere 0.1 m	0.36	0.36	0.39	0.065	0.067	0.066
2020	Rur	0.47	0.57	0.68	0.054	0.047	0.041
	Apulian Tavoliere 0.025 m	0.19	0.17	0.29	0.074	0.069	0.063
	Apulian Tavoliere 0.1 m	0.29	0.27	0.39	0.069	0.066	0.063

and

Table 5. Statistical comparison between the individual processing steps in regard to the correlation and error of the related estimated and in-situ measured soil moisture for each individual station within the Rur catchment. Darker colors indicate a better performance.

			R²			uRMSE [vol. %]
		CRNP	IA	IA + FS	IA + FS + VA	IA	IA + FS	IA + FS + VA	SM Range	VV Range
Crop dominated	2018	RU_C_006	0.12	0.23	0.51	6.54	6.14	4.67	26.12	0.11
		RU_BCK_002	0.42	0.23	0.52	5.76	6.43	5.03	26.21	0.20
		ME_BCK_001	0.28	0.38	0.56	5.74	5.36	4.47	23.40	0.14
		SE_C_001	0.46	0.45	0.60	5.55	5.54	4.79	26.10	0.15
	2019	RU_C_006	0.37	0.38	0.51	4.25	4.51	3.85	21.69	0.11
		RU_BCK_002	0.29	0.47	0.60	5.70	4.92	4.35	26.24	0.20
		ME_BCK_001	0.18	0.39	0.52	5.89	5.02	4.53	23.55	0.17
		SE_C_001	0.06	0.36	0.53	9.10	4.02	3.33	19.45	0.16
	2020	RU_C_006	0.54	0.54	0.80	4.21	4.43	2.89	23.14	0.16
		RU_BCK_002	0.66	0.44	0.82	4.50	5.41	3.02	24.65	0.20
		ME_BCK_001	0.62	0.42	0.87	4.38	5.44	2.66	24.99	0.16
		SE_C_001	0.05	0.67	0.39	10.64	6.26	8.02	37.13	0.17
Meadow dominated	2018	RU_C_005	0.61	0.73	0.73	3.66	3.12	3.10	22.23	0.07
		RU_C_007	0.28	0.65	0.65	6.55	5.34	5.34	30.86	0.08
		RO_C_001	0.36	0.53	0.52	5.48	4.86	4.87	30.88	0.08
		RU_BCDKR_001	0.31	0.58	0.58	5.61	4.46	4.44	28.73	0.09
	2019	RU_C_005	0.45	0.63	0.63	4.32	3.82	3.82	19.49	0.06
		RU_C_007	0.31	0.56	0.56	5.52	4.46	4.43	24.20	0.08
		RO_C_001	0.25	0.51	0.51	5.93	4.78	4.81	25.66	0.05
		RU_BCDKR_001	0.24	0.55	0.55	5.79	4.77	4.76	21.38	0.05
	2020	RU_C_005	0.55	0.76	0.76	4.11	3.22	3.21	22.35	0.08
		RU_C_007	0.36	0.66	0.73	6.19	4.47	4.12	27.73	0.10
		RO_C_001	0.48	0.51	0.51	4.88	4.76	4.79	24.52	0.07
		RU_BCDKR_001	0.49	0.55	0.55	4.27	3.91	3.91	24.26	0.08

). Since the agricultural crops at the respective CRNP sites change annually in some cases, the methods are evaluated for the individual years 2018, 2019 and 2020, using the coefficient of determination (R²) as well as the unbiased Root Mean Square Error (uRMSE). It needs to be mentioned that by comparing soil moisture estimations from C-band SAR observations to in-situ measurements from CRNP station, different soil depths are considered. While C-band only represents the first upper centimeters of soil, CRNP represents a soil depth between 15 cm to up to 70 cm, depending on the soil moisture condition. In general, a strong correlation between surface and root zone soil moisture was observed, when comparing SMOS near-surface soil moisture measurements to root-zone in-situ data [90]. Nonetheless, under extreme dry conditions, remotely sensed soil moisture shows a reduced variability, due to the limited storage capacity of the upper soil layer compared to the root zone soil moisture [91]. On the other hand, when comparing point data from capacitance sensors from 0.025 m and 0.1 soil depths, the spatial soil heterogeneity can result in soil moisture values, which are not representative for the area covered by the radar signal. Especially when losing contact to the surrounding soil, capacitance probes can lead to unreasonable values. In this regard, CRNP measurements are more suitable, providing an integrated soil moisture value over its footprint area. However, in this study, no significant difference between CRNP and capacitance measurements were found, regarding the performance of the introduced soil moisture estimation, whether due to different soil depth nor different spatial footprints of both in-situ instruments.

Using only the incidence angle normalisation, both low to high correlations between in-situ measured and estimated soil moisture are achieved for the Rur catchment as well as Apulian Tavoliere site. In general, the Rur catchment had mean R² values of 0.36–0.47 and mean uRMSE of 0.054–0.058 m³/m³ for the individual years (Table 4).

Almost no correlation is observed at the CRNP station SE_C_001 in 2019 and 2020 with R² of 0.06 and 0.07, while the highest correlations are observed in 2020 at the CRNP stations RU_BCK_002 and ME_BCK_001 with R² values of 0.66 and 0.62, respectively. The highest uRMSE is observed in 2019 and 2020 at station SE_C_001 with 0.091 m³/m³ and 0.106 m³/m³, respectively. The lowest uRMSE values are observed in 2020 at the station RU_C_006 with 0.042 m³/m³ and RU_BCDKR_001 with 0.043 m³/m³ (Table 5). At the Apulian Tavoliere site, mean R² values of 0.15–0.36 and mean uRMSE of 0.058–0.081 m³/m³ for a soil depth of 0.025 m, and mean R² values of 0.29–0.37 and mean uRMSE of 0.060–0.069 m³/m³ for a soil depth of 0.1 m are achieved (Table 4). The lowest correlation is found at 0.025 m soil depths at station 9 and 5 in 2020, with R² values of 0.01 and 0.09, respectively, while the highest correlations can be found at 0.1 m soil depth at station 10 with 0.68 in 2019 and station 5 with 0.51 in 2018. Here, the uRMSE is highest in 0.025 m soil depth at station 10 in 2019 with 0.112 m³/m³, while it is lowest at the same year and station but in 0.1 m soil depth with 0.044 m³/m³ (

Table 6. Statistical comparison between the individual processing steps in regard to the correlation and error of the related estimated and in-situ measured soil moisture for each individual station at the Apulian Tavoliere site. Darker colors indicate a better performance.

	Soil Depth	0.025 m						0.1 m
		R²			uRMSE [vol. %]			R²			uRMSE [vol. %]
	TDR	IA	IA + FS	IA + FS + VA	IA	IA + FS	IA + FS + VA	IA	IA + FS	IA + FS + VA	IA	IA + FS	IA + FS + VA
2018	Station_2	0.42	0.50	0.30	5.04	4.65	5.42	0.33	0.42	0.30	6.64	6.15	7.06
	Station_3	0.35	0.34	0.39	5.38	5.77	5.61	0.28	0.27	0.36	5.71	6.09	5.74
	Station_5	0.43	0.44	0.60	6.73	6.69	5.82	0.51	0.47	0.54	5.46	5.70	5.32
	Station_7	0.41	0.42	0.53	6.00	6.10	5.25	0.45	0.43	0.51	5.94	6.16	4.48
	Station_9	0.23	0.30	0.47	5.52	5.11	5.28	0.32	0.44	0.41	5.92	5.46	5.97
	Station_10	0.36	0.36	0.39	6.33	6.39	6.26	0.33	0.33	0.37	5.96	6.02	5.90
2019	Station_2	0.14	0.15	0.39	7.50	7.49	5.92	0.17	0.20	0.46	7.30	7.05	5.96
	Station_3	0.36	0.35	0.55	5.91	6.16	4.79	0.33	0.34	0.60	6.04	6.18	4.49
	Station_5	0.10	0.11	0.21	9.31	9.36	8.46	0.23	0.23	0.10	7.23	7.38	8.58
	Station_7	0.14	0.15	0.79	8.15	8.06	3.64	0.55	0.56	0.37	6.35	6.33	7.61
	Station_9	0.10	0.11	0.34	6.37	6.45	5.96	0.21	0.19	0.42	8.07	8.26	6.78
	Station_10	0.08	0.09	0.64	11.24	11.32	6.49	0.68	0.64	0.40	4.34	4.89	6.26
2020	Station_2	0.13	0.17	0.02	7.12	6.28	8.52	0.12	0.16	0.49	7.24	6.38	5.30
	Station_3	0.13	0.05	0.30	7.97	7.90	6.03	0.43	0.34	0.43	6.03	5.92	5.32
	Station_5	0.04	0.03	0.19	8.17	7.67	6.23	0.20	0.20	0.22	6.47	5.91	5.66
	Station_7	0.37	0.41	0.46	6.75	5.79	4.93	0.50	0.44	0.43	6.26	6.22	5.93
	Station_9	0.01	0.01	0.50	8.55	8.04	5.52	0.34	0.28	0.52	7.26	7.46	5.67
	Station_10	0.45	0.38	0.26	6.01	5.84	6.58	0.16	0.15	0.25	8.40	7.89	7.18

). Even though this could be caused by differing soil moisture within the related soil depths, most likely it is due to corrupted measurements, as other stations do not show such a significant difference.

As described in Section 5.1, the influence of the alternating incidence angles on the backscattering signal is still present in the incidence angle normalised time series, affecting the resulting soil moisture estimation. Grouped in the three different vegetation periods, the different range of in-situ measured and estimated soil moisture is displayed (Figure 10). For the Apulian Tavoliere site at 0.025 m soil depth, the estimated soil moisture range matches the observed one. However, for the Rur catchment, there was an overestimation of the minimum soil moisture, while for the Apulian Tavoliere site at 0.1 m, there is an underestimation. This suggests that the scaling from GLDAS may not be optimal but can be further improved. Furthermore, the estimated soil moisture still showed extreme low and high values due to the incidence angle still effect present in the backscatter time series. Crop-dominated stations showed significant variation in correlation and error between individual years due to the influence of incidence angle, which depends on plant geometry, plant phase, and vegetation water content, changing both between the year as well as within a year [72,84]. In contrast, the soil moisture estimation at the meadow station (RU_C_005) is stable over the years, with only slight changes in correlation and uRMSE. Dominated by meadows with short grass vegetation, it results in a more stationary incidence angle effect during the growing period. As the incidence angle effect is in general not as prominent due to the short vegetation cover, for these types of landcovers, an incidence angle normalized backscattering time series from multi-orbit Sentinel-1 scenes can already estimate the soil moisture sufficiently during the whole year.

The soil moisture estimation from incidence angle normalized and Fourier Series transformed time series (IA + FS) generally improves in terms of both correlation coefficient and uRMSE. The mean R² from individual CRNP stations is between 0.47 and 0.57, while the corresponding uRMSE values range between 0.045 m³/m³ and 0.052 m³/m³ (Table 4). The strongest improvement is observed in 2019, with a 79% increase in R² and 22% decrease in uRMSE. However, Fourier transformation led to reduced correlations in four cases, all being observed at crop dominated stations (Table 5). At the Apulian Tavoliere site, the mean R² from the individual capacitance stations at 0.025 m range between 0.16 and 0.39 with corresponding uRMSE between 0.058 m³/m and 0.081 m³/m³. At 0.1 m soil depth, R² values between 0.27 and 0.39 with corresponding uRMSE values between 0.059 m³/m³ and 0.067 m³/m³ are achieved.

Focusing on the alternating incidence angles, the Fourier Series transformation further reduces its effect on the resulting soil moisture. The minimum and maximum estimated soil moisture values at the CRNP stations in the three vegetation periods range from 0.186 m³/m³ to 0.563 m³/m³, from 0.162 m³/m³ to 0.481 m³/m³, and from 0.145 m³/m³ to 0.412 m³/m³ (Figure 11). While still exceeding the observed range of soil moisture values in the non-vegetation period, there is a clear reduction compared to the range of IA soil moisture estimation. The transformation into the Fourier Series is especially addressing the maximum soil moisture values. Looking at the varying correlation between in-situ measured soil moisture and estimated soil moisture for the individual years at each CRNP station, the additional Fourier Series transformation could reduce the observed range of R² during the observation period at each station in the mean by around 16%, compared to the IA soil moisture. At the Apulian Tavoliere site, a similar behavior can be observed, with a decrease in the upper extreme soil moisture values.

The addition of vegetation correction improves soil moisture estimation at the CRNP stations with mean R² values ranging between 0.55 and 0.68 and corresponding uRMSE values between 0.041 m³/m³ and 0.046 m³/m³. The improvements in R² range from 15% to 105% compared to the previous soil moisture estimations without vegetation correction, and uRMSE is reduced by up to 27%. The soil moisture estimations for the soil depth of 0.025 m and 0.1 m at the Apulian Tavoliere site improved in general as well. For the 0.025 m soil depth, R² values range between 0.29 and 0.49, with the corresponding uRMSE between 0.056 m³/m³ and 0.063 m³/m³. For the 0.1 m soil depth, R² values range between 0.39 and 0.42, with corresponding uRMSE between 0.057 m³/m³ and 0.066 m³/m³. In this regard, the improvements ranged from 14% to 217% in R² and from 1% to 28% in uRMSE compared to the backscatter time series without vegetational correction.

Regarding the correlation coefficient at the crop dominated stations, the vegetation correction improved the observed R² values in all but one case (Table 5). At station SE_C_001 in 2020, the IA + FS soil moisture estimations are performing better with an R² of 0.67 compared to 0.39 from the IA + FS + VA soil moisture estimation. For uRMSE, a similar behavior is present at that station. At the Tavoliere site in Puglia, the increase in performance due to vegetation adaptation is most evident when compared to soil moisture values at 0.025 m depth, while it is not as pronounced at 0.1 m soil depth. Looking at the meadow dominated CRNP stations, the correlations between IA + FS and IA + FS + VA are mostly the same, as the vegetational correction is not applied but for one station. Having a backscatter intensity range greater than 0.10, the vegetation correction was applied and increased the correlation from R² = 0.66 to R² = 0.73 at station RU_C_007 in 2020. As the footprint of the CRNP stations as well as the 200 m pixels incorporates multiple agricultural fields, a change in the field diversity could be responsible for the change in backscatter range compared to the previous years.

With an overall R of 0.68, R² of 0.46 and uRMSE of 0.056 m³/m³ at the Rur catchment, as well as an overall R of 0.56 (0.63), R² of 0.31 (0.40) and uRMSE of 0.065 m³/m³ (0.064 m³/m³) at the Apulian Tavoliere site at 0.025 m (0.1 m) soil depth, the proposed method for estimating soil moisture achieves consistent results (Figure 12). Here, the highest correlation can be observed at the CRNP stations, followed by the capacitance station at 0.1 m soil depth. The performance for the shallow 0.025 cm soil depth is lower, likely being caused by the higher spatial variability of surface soil moisture. As topsoil measurements tend to be more heterogeneous and may lack spatial representativeness for a 40.000 m² pixel, while measurements from deeper soil layers, being less variable, represent more the average dynamics and achieve better correlations with spatial SAR measurements. On the other hand, this could also be caused by the exclusion of high frequencies from the Fourier transformed backscatter time series. As the topsoil shows higher temporal soil moisture dynamics, the low-pass filtering effect could have resulted in lower correlation for this soil moisture depth, even though it is the most comparable in terms of penetration depth of the SAR signal. Reflecting more weekly and seasonal changes, soil moisture measurements from deeper depth could thus have a higher correlation to the low-passed SAR time series. Since at the CRNP station, both the areal measurement of soil moisture and the low-frequency dynamics of soil depth from 0.15–0.70 cm are combined, this could possibly be the reason for the highest correlation and the lowest uRMSE between estimated and in-situ measured soil moisture observed here.

For comparison, we also calculated the correlation and error between used GLDAS and in-situ measured soil moisture derived from CRNP stations in the Rur catchment. Evaluating GLDAS for all stations and all years in the Rur catchment, an R value of 0.69 with an uRMSE of 0.061 m³/m³ was achieved, compared to R = 0.68 and uRMSE = 0.056 m³/m³, achieved with the newly proposed approach. In this sense, similar results were obtained, with a significantly higher spatial resolution (GLDAS = 27 km, our approach 0.2 km). The overall performance at the Apulian Tavoliere site is also comparable to the published SMOSAR soil moisture product validation study over seven soil moisture test sites located across Europe, North America, and Australia for the period between 2015 to 2018 [29]. They were using only a single orbit time series with 1068 individual scenes to estimate soil moisture with a spatial resolution of 1 km and a temporal resolution of 6 to 12 days. Validated against the spatial average of the Apulian site on the network scale, an R value of 0.71 with an RMSE of 0.064 m³/m³ was achieved. As the spatial averaging increases the methods sensitivity to soil moisture, averaging our final soil moisture product over the six capacitance stations further increases the correlation, R, to a value of 0.68 (0.69) with an RMSE = 0.050 m³/m³ (0.053 m³/m³) for the soil depth of 0.025 m (0.1 m). Spatially averaging across the eight CRNP stations in Rur catchment on a daily basis, R can be further increased to a correlation of R = 0.79 with a RMSE = 0.039 m³/m³, while having 709 averaged soil moisture estimations. In this respect, the proposed method can further improve the short-term change detection approach, by using smaller temporal intervals between each observation. Furthermore, by applying vegetational correction, the influence of the changing scattering contribution of vegetation can also be managed.

As shown, major improvements were achieved in the intermediate and full vegetational period, being addressed by the vegetational correction. In this regard, the presented method is able to estimate soil moisture also during the vegetation period, by reducing the effect of varying incidence angles and vegetation scattering using the temporally dense Sentinel-1 co-polarized C-band time series. It enables a continuous soil moisture monitoring of agricultural areas during the whole crop cycle, representing soil moisture patterns on a regional scale, and thus can fill the gap of recent SAR based soil moisture products, regarding their lacking spatial or temporal resolution for a regional soil moisture monitoring. By using the co-polarized backscatter signal as an indicator for the vegetational influence on the related backscatter ratios for agricultural areas, which have a distinct bare soil and vegetation cover period, the limitations of C-band vegetation sensitivity, leading to biased soil moisture information within the backscattering signal, can be addressed. The observation gaps, which occur when masking out pixels dominated by volume backscattering as proposed by the SMOSAR soil moisture method, can be avoided, especially as these gaps occur during the most important phenological stages (shooting, flowering, and ripening of crops) in the growing period. While being able to improve soil moisture estimation under changing vegetation conditions in general, the retrieval of soil surface information using C-band is physically limited. In this regard, depending on the crop type and the vegetation density, for some periods during the growing cycle, C-band can become completely insensitive to soil moisture [42]. Even though the resampling to a coarser resolution of 200 m averages out the influence of vegetation, there might be pixels, where no correlation between a change in backscatter and a change in soil moisture is present. Hence, further research needs to be carried out, developing a (multi-sensor) method more adapted to the individual crop types and having deeper penetration capabilities through the canopy.

While the method performs well in both the overall mean and the growing season, the lowest correlations and highest uRMSE are observed in the non-growing season, where wet soil moisture conditions are prevailing. As previous studies show, the alpha approximation method has the highest bias within the soil moisture range from 0.30–0.55 m³/m³, requiring further improvements especially for wet soil moisture conditions [29]. In addition, the assumption, that the initial soil moisture condition is on the field capacity level is not fully matching for the years 2019 and 2020 at the Rur catchment, as the prevailing drought conditions in these two years are causing soil moisture values lower than field capacity in both January months. By using adapted initial soil moisture conditions, the problem of setting appropriate initial soil conditions should be addressed. In this regard, an iterative approach, as proposed by Balenzano et al., 2011 [81], could mitigate the effect of the influence of pixels not suitable to be reference for the initial soil moisture. Using a moving window across the time series, soil moisture is estimated for a shorter period, while the mean of the resulting soil moisture estimation serves as initial soil moisture values for the further iteration. By limiting the calculated mean values of soil moisture to the maximum value derived from the field capacity when used as initial soil moisture in the next iteration, the high number of outliers can also be addressed. Nevertheless, this comes with an increase in computational costs, reducing the applicability of the method for regional monitoring. Here, further research is needed, optimizing the iteration period regarding its improvement and computational time.

6. Conclusions

The study presents a methodology for estimating volumetric soil moisture content [m³/m³] under changing vegetation using multi-orbit co-polarized Sentinel-1 C-band time series. The method was tested using TERENO sites in the Rur catchment, Germany (Figure 13) as well as the Apulian Tavoliere site, Italy. The high revisit time of the Sentinel-1 A and B satellites enables a temporal resolution of the resulting soil moisture product of one to two days when using all available orbits. The method was evaluated using a spatial resolution of 200 m in accordance with the footprint radius of the eight in-situ soil moisture sites collecting cosmic ray neutron probe (CRNP) measurements for the years 2018 to 2020. In addition to the SAR recordings, auxiliary data of soil texture and field capacity were used as initial starting value and soil moisture to dielectric constant conversion. The presented method is based on a change in detection approach, using the short-term change between two consecutive SAR recordings caused solely by changes in soil moisture. The influence of other surface parameters such as soil roughness, vegetation biomass, plant height, or phenological stage can be disentangled from the change in backscattering signal to some extent, as they are associated with a slower rate of change. In this regard, establishing a temporal high resolution time series from multiple orbits increases the validity of this algorithm pre-requisite. To reduce the effect of alternating incidence angles on the backscattering signal introduced by using the multi-orbit time series, an incidence angle normalization was applied, and an additional Fourier Series transformation was developed. Furthermore, the effect of changing vegetation scattering contribution on the backscattering signal was addressed. The backscatter ratios were corrected for vegetation influence based on the related co-polarized backscatter time series as an indicator for the vegetation backscatter contribution at pixel scale. In the study, the effect on the individual processing steps on the backscattering time series as well as on the resulting soil moisture estimation are analyzed and discussed. Furthermore, the final soil moisture product was validated with in-situ measurements for two different sites with two different in-situ measuring techniques as well as three different soil depths.

Each of the applied preprocessing steps significantly increased the overall fitness of the estimated soil moisture, both in terms of correlation coefficient as well as regarding uRMSE. While the incidence angle normalization led to an alignment of the mean backscattering values between individual orbits, the effect of incidence angle is still present in the individual backscattering signal distribution. By using the Fourier Series transformation, dismissing frequencies higher than the revisit time of the Sentinel-1A and Sentinel-1B, the change in backscattering signal caused from varying incident angles, can be excluded, resulting in higher correlation as well as reduced uRMSE at most individual CRNP and capacitance stations. Adding the newly developed vegetation correction further increase the correlation between estimated and in-situ measured soil moisture, especially at the crop dominated stations.

To provide an outlook, further research should be carried out adapting the initial soil moisture condition for drought periods as well as evaluating other globally available datasets for this purpose. Especially, the assumption that soil moisture conditions in the Rur catchment in January are near the field capacity level is not true throughout the drought years of 2019 and 2020. Furthermore, the criteria for applying the vegetational adaption was chosen empirically based on the characteristics of the study areas. For a global application, further research should be carried out on establishing a more generalized criterion.

By establishing the soil moisture retrieval algorithm within a cloud-processing environment, the SAR-based soil moisture estimation takes a step towards a globally applicable. The scaling of the obtained soil moisture information in terms of adaptable areas, time periods as well as spatial and temporal resolution leads to improvements in its usability within the different operational sectors, e.g., for management decisions in agriculture on a field basis or for drought monitoring on a regional level. Moreover, the temporal resolution of Sentinel-1 satellite recordings will continuously improve with the upcoming launches of Sentinel-1C and Sentinel-1D, despite the failure (23 December 2021) of Sentinel-1B. In this regard, high temporal resolution time series of 1 to 2 days using only single orbit acquisitions as well as sub-daily time series using all available orbits for estimating soil moisture will be possible in the near future. A combination with upcoming launches of L-band SAR missions like NISAR, ALOS-4 and ROSE-L and Tandem-L as well as P-band missions like BIOMASS will offer great potential for improving soil moisture estimation under dense vegetation cover as well as for different soil depths.