Spatial and Temporal Distribution of Soil Moisture at the Catchment Scale Using Remotely-Sensed Energy Fluxes

Abstract

Despite playing a critical role in the division of precipitation between runoff and infiltration, soil moisture (SM) is difficult to estimate at the catchment scale and at frequent time steps, as is required by many hydrological, erosion and flood simulation models. In this work, an integrated methodology is described to estimate SM at the root zone, based on the remotely-sensed evaporative fraction (Λ) and ancillary information on soil and meteorology. A time series of Terra MODIS satellite images was used to estimate SM maps with an eight-day time step at a 250-m spatial resolution for three diverse catchments in Europe. The study of the resulting SM maps shows that their spatial variability follows the pattern of land cover types and the main geomorphological features of the catchment, and their temporal pattern follows the distribution of rain events, with the exception of irrigated land. Field surveys provided in situ measurements to validate the SM maps’ accuracy, which proved to be variable according to site and season. In addition, several factors were analyzed in order to explain the variation in the accuracy, and it was shown that the land cover type, the soil texture class, the temporal difference between the datasets’ acquisition and the presence of rain events during the measurements played a significant role, rather than the often referred to scale difference between in situ and satellite observations. Therefore, the proposed methodology can be used operationally to estimate SM maps at the catchment scale, with a 250-m spatial resolution and an eight-day time step.

1. Introduction

Soil moisture (SM) plays a critical role in the division of precipitation between runoff and infiltration, as well as the evapotranspiration and development of vegetation. It therefore plays an important role in the catchment water budget. It is a required input for various models and studies related to the water and carbon cycle, soil erosion, crops development, flood forecasting and drought monitoring. SM can be easily measured in situ by collecting gravimetric samples of undisturbed soil or using time domain reflectometry. These methods can provide accurate estimates over the entire root zone, but they are point measurements and cannot account for the high spatial soil moisture variability at large scales [1]. At the catchment level, soil moisture can be estimated from water balance models, which, however, require rigorous parameterization and are data intensive [2,3,4].

Remote sensing methods offer the advantage of providing spatially-distributed SM estimates. Satellite images of various wavelengths have been used to estimate SM, capitalizing on the advantages of remote sensing, i.e., low cost data for large areas, wall-to-wall coverage and easy production of time series, including access to historical data. Several SM products are produced at a large scale (25–50 km). The NOAA SSM/I (Special Sensor Microwave/Imager) product is the daily SM product for land regions produced by the NOAA team at a 25-km spatial resolution [5]. The ESA SMOS (Soil Moisture and Ocean Salinity) mission was launched in 2009 in order to provide world maps of SM with a frequency of 1–3 days at a 35–50-km resolution [6]. The ASCAT (Advanced Scatterometer) on board MetOp satellites has provided surface soil moisture at a coarse resolution (25–50 km) with a daily global coverage since 2007 [7]. The recently-launched (January 2015) Soil Moisture Active Passive satellite (SMAP) provided measurements of land surface SM and freeze-thaw state with near-global revisit coverage every 2–3 days. However, the above-mentioned products have low resolution, and thus, methods for their downscaling have been proposed to improve their usability [8,9].

Several physical, semi-empirical and empirical models have been proposed for estimating SM using microwave satellite images at high resolution; but, their application is usually limited to small homogeneous areas, and they are heavily influenced by the presence of vegetation [10,11,12]. So far, little attention has been paid to the use of the surface temperature of vegetated surfaces for the determination of subsurface SM [13]. A first approximation of volumetric SM using SEBAL (surface energy balance algorithms for land) was made by Bastiaanssen et al. [14]. They introduced an empirical relationship between the evaporative fraction and the volumetric SM content. The equation was used to estimate the moisture content of the entire root zone in the irrigated soils of the Lerma-Chapala basin in Mexico. Scott, Bastiaanssen and Ahmad [13] modified this relationship using the saturated soil water content in order to derive the SM content in the root zone. When testing this method on irrigated plains in central Pakistan and Mexico, they obtained an average error of 0.035 cm³ cm⁻³, and 90% of the cases had an error of less than 0.07 cm³ cm⁻³ and an RMS error of 0.05 cm³ cm⁻³. The method was further developed to estimate phreatic water storage in the unsaturated zone with the combined use of biannual piezometric observations [15].

Root zone SM has largely been the main focus of scientific research rather than the topsoil surface [16]. Manfreda et al. [17] concluded that the relationship between the volumetric SM at the surface to that in deeper layers of soil is not always clear, as SM is not uniformly distributed with depth. This implies that prediction of SM in the deep layer that is based on the superficial SM carries a certain level of uncertainty. Thus, direct estimation of root zone SM rather than the surface SM, usually offered by remote sensing, has increased benefit. On the other hand, Bastiaanssen, Molden and Makin [14] point out that methods for surface SM estimation using airborne microwave radiometry, though complex, are more accurate than the root zone moisture estimation ones at low spatial resolution.

The methodologies and products for mapping SM with remote sensing either use scatterometer data at sparse resolutions, or high resolution microwave images for sub-catchment areas, or are fusions of remote sensing and hydrological models [2,8,9,18]. The aim of this work is to describe an integrated methodology to use remotely-sensed energy balance fluxes and ancillary information on soil and meteorology in order to provide a time series of soil moisture maps at the catchment scale operationally. Specific objectives include the analysis of the effect of land cover type, soil texture and other environmental parameters. The methodology was applied in several diverse European catchments for a hydrological year, and in situ measurements were used to validate its accuracy.

2. Materials and Methods

2.1. Study Areas

The following study areas were selected to demonstrate the applicability of the methodology in diverse hydrological conditions (Figure 1).

Figure 1. Location map of study areas, main land cover types and measurement locations. SM, soil moisture.

Figure 1. Location map of study areas, main land cover types and measurement locations. SM, soil moisture.

2.1.1. The Nestos River

The Nestos is a transboundary river between Greece and Bulgaria. Its watershed covers an area of 6100 km². Water is used for energy production and irrigation. The land cover is mainly natural vegetation (60%), and the rest is agricultural land. The climate is Mediterranean-continental, with a mean annual precipitation of 750 mm. Soil texture is mostly medium and fine along the river valley.

2.1.2. Rijnland

Rijnland is a low-lying area of land surrounded by dikes with an artificially-controlled water system using canals, pumps and weirs (polder) in the western part of the Netherlands, bordering the North Sea. The total area is about 1000 km², of which the majority (72%) is occupied by low-lying land reclamation areas, while 15% is covered by free-draining areas and 8% by dunes. The area consists of urban and rural sectors committed to horticulture, agriculture and grasslands. The climate is temperate-maritime, and the annual precipitation is 960 mm. Soil texture is coarse on the northwest side along the coast, medium in the central area and medium-fine in the southern section of the study area.

2.1.3. The Tamega River

The Tamega is a transboundary river between Portugal and Spain and is a tributary of River Douro. Its watershed has an area of approximately 3270 km² and drains to the Torrao reservoir. The main land cover types are forest (65%) and agriculture (31%). The climate is warm temperate-Mediterranean, and the mean annual precipitation is 850 mm. Soil texture is medium throughout the study area.

2.2. Data and Pre-Processing

Several products from the MODIS (Moderate Resolution Imaging Spectroradiometer) instrument aboard the Terra (formerly known as Earth Observing System—EOS AM) satellite were used:

• MODIS surface reflectance (MOD09Q1, MOD09A1) provides an estimate of the surface spectral reflectance at a 250- and a 500-m spatial resolution, as would have been measured at ground level in the absence of atmospheric scattering or absorption.
• MODIS global land surface temperature (LST) and emissivity (MOD11A2) are estimated using the split-window technique on the thermal bands of MODIS at a 1-km spatial resolution.

The above-mentioned products are 8-day composites, meaning that each pixel contains the best possible observation during an 8-day period as selected on the basis of high observation coverage, low view angle, the absence of clouds or cloud shadow and aerosol loading. Additional data provided with these products include quality assessment, the day of the year for the pixel, along with solar, view and zenith angles. The quality assessment information of each MODIS Land HDF-EOS (Hierarchical Data Format-EOS) data product was utilized in order to remove low quality pixels using the LDOPE software tools provided by the MODIS Land Quality Assessment [19]. The products mentioned above are validated Stage 2, meaning that accuracy has been assessed over a widely-distributed set of locations and time periods via several ground truth and validation efforts [20]. For instance, the accuracy of MODIS LST is better than 1 °C, when validated with daily in situ data [21,22] and better than 2 °C using nighttime ground measurements [23].

Saturated water content (θs) was used as an input dataset for the estimation of SM. θs was estimated by Joint Research Centre (JRC)/Soil Bureau using pedotransfer functions as described in Tóth et al. [24]. The prediction method was used depending on the type of available input information, which were the Database of Hydraulic Properties of European Soils (HYPRES) [25], the Hungarian Detailed Soil Hydrophysical Database (MARTHA) [26] and the Soil Map of the Region of Eastern Macedonia and Thrace [27] datasets. The accuracy of estimations reached an RMSE (vol%) of 3.437 and a Pearson correlation coefficient of 0.857 (p = 0.01) with the measured data.

An updated map of land cover was created per study area. Multispectral SPOT XS and Landsat 5 TM satellite images were used, together with Web Map Service data (Google Earth) and in situ observations, in an automated supervised classification to produce the map of the main land cover types [28].

Raster meteorological data were utilized instead of point measurements from ground meteorological stations, in order to account for the large spatial variability within each catchment and the sparse coverage of ground station observations. The raster meteorological data were provided by the Center for Weather Forecast and Climate Analysis/Brazilian National Institute for Space Research (CPTEC/INPE). The CPTEC/INPE configured its regional Eta model to run for the study areas with a 5-km horizontal resolution and boundary conditions from CPTEC/INPE Global Model, T213L42 version. The initial conditions were obtained from the GFS (Global Forecast System) numerical model of the NCEP/NOAA (National Centers for Environmental Prediction/National Oceanic and Atmospheric Administration). The ETA model was executed daily at 00 and 12 Zulu time (GMT) and provided forecasts with a one-hour time step up to 3 days ahead. The 24-hour forecasts were used to avoid the transient errors from the initial stabilization of the model results before getting optimum accuracy. Table 1 lists the meteorological parameters that were used in further analyses.

Rainfall data were gathered from at least one rain gauge lying within each study area. Precipitation heights were recorded daily, with the exception of the Chrysoupoli meteorological station in Nestos that recorded the presence of precipitation events daily. Thus, the information from this station was used qualitatively.

Table 1. List of meteorological parameters that were used.

Parameter Name	Unit
Surface downward short wave flux	(W/m2)
2 m relative humidity	(%)
2 m temperature	(K)
10 m u wind (zonal component)	(m/s)
10 m v wind (meridional component)	(m/s)

Table 1. List of meteorological parameters that were used.

Parameter Name	Unit
Surface downward short wave flux	(W/m2)
2 m relative humidity	(%)
2 m temperature	(K)
10 m u wind (zonal component)	(m/s)
10 m v wind (meridional component)	(m/s)

Field surveys were conducted before and during the study period to collect in situ measurements of soil moisture in order to assess the accuracy of the results. A stratified random sampling was adopted for reducing the number of samples, while at the same time maintaining the desired level of accuracy. The strata were defined from overlay analysis of broad vegetation and soil texture categories from recent satellite images and soil parameter maps. The number of sample sites per stratum was determined as a function of SM variability and the desired accuracy in the estimations [29]. Their location was determined randomly using a GIS function. In each sample location, soil moisture measurements were conducted on at least 15 sub-locations, spread around the sample location, aiming to cover the variability within the satellite pixel [30]. The soil moisture was measured with a FieldScout TDR probe. The instrument is based on time domain reflectometry, where the apparent dielectric constant of the soil matrix (Ka) is determined using an electrical field and is empirically related to the volumetric soil moisture content. The TDR was previously calibrated using 50 soil core samples covering a wide range of soil types and moisture contents, where soil moisture was measured with the gravimetric method (R² = 0.857, p < 0.01). Soil moisture was recorded using the 7.5–20-cm probes, and basic soil observations were noted.

2.3. Methods for Estimating Soil Moisture with Thermal Infrared Satellite Images

Soil moisture in the root zone (θ_rz) was estimated using the following equation [13]:

$${\theta }_{rz}={\theta }_{sat}\text{ exp}\left(\frac{\Lambda -1}{0.42}\right)$$

(1)

where θ_sat is the saturated water content and Λ is the evaporative fraction. θ_sat is a normalization term, so that when the relative soil moisture ($\frac{{\theta }_{rz}}{{\theta }_{sat}}$) reaches 1 (soils fully saturated), the soil moisture θ_rz equals the soil porosity.

The evaporative fraction (Λ) is the ratio between the latent heat flux (λE) that is the energy consumed for evapotranspiration and the net available energy that is the net radiation (R_n) minus soil heat flux (G₀).

$$\Lambda =\frac{\lambda E}{R_n-G_0}$$

(2)

Under non-advective conditions, the Λ values ranges usually between 0 and 1, which represents zero to maximum evapotranspiration [31].

Τhe latent heat flux (λE) is derived using the energy balance equation:

$$\lambda E=R_n-G_0-H$$

(3)

The net radiation represents the actual energy available at the surface resulting from the total incoming short and longwave radiation and its loss through various processes (evapotranspiration, heat production, metabolic activities) [32].

$$R_n=\left(1-a\right)K_{in}+L_{in}-L_{out}-\left(1-{\epsilon }_0\right)L_{in}$$

(4)

where α is the albedo, ε₀ is the emissivity, K_in is the shortwave input radiation and L_in is the incoming longwave radiation defined as:

$$L_{in}=(1.08{(-\ln {\tau }_{sw})}^{0.265})5.67{10}^{-8}{T}_a^4$$

(5)

L_out corresponds to the outgoing longwave radiation and is calculated as a function of the surface temperature (T_s) and ε₀:

$$L_{out}=\left({\epsilon }_{0}5.67{10}^{-8}{T}_s^4\right)$$

(6)

The soil heat flux G₀ refers to the heat stored in the soil and the vegetation. Having estimated R_n, G₀ can be easily deduced through the following expression:

$$G_0=R_n\left(T_s-273\right)\left(0.0032\frac{a}{0.9}+0.0062{\left(\frac{a}{0.9}\right)}^2\right)\frac{1-0.978NDV{I}^4}{a}$$

(7)

In the above equations, T_s, ε₀, α and the Normalized Difference Vegetation Index (NDVI) are intermediate variables computed using the MODIS data using conventional equations. The shortwave transmissivity (τ_sw), T_a and K_in are derived from meteorological data output from the weather forecast model.

The sensitive heat flux (H) corresponds to the heat loss to the air due to the temperature difference, as shown by the equation below:

$$H=\frac{C_p\rho \mathrm{d}T}{r_{ah}}$$

(8)

where C_p is the air specific heat, ρ is the air density, r_ah is the aerodynamic resistance to heat transport and dT is the difference between soil and air temperatures. The accurate computation of H using raster meteorological data is detailed in [33].

The latent heat flux, thus the evaporative fraction, is relative to the whole root zone, so root depth and soil depth maps were used to estimate the θ_sat information, using a linear relation.

In order to solve the energy balance equation, a variation of the SEBAL algorithm was used. Alexandridis et al. [34] developed the “ITA-Water” tool (Integrated Thermodynamic Algorithms for Estimating Agricultural Water Use) for the implementation of the SEBAL approach using Terra MODIS composite images. Based on this tool and considering the requirements for using spatially-distributed meteorological data, the tool evolved into the “ITA-MyWater”, and its functionalities were extended to estimate the soil moisture at the root zone [33].

3. Results and Discussion

3.1. Spatial Patterns of Soil Moisture

Four SM maps per site were selected from the time series to show the spatial variability of each study area at various seasons (Figure 2). One-way analysis of variance (ANOVA) on a per-pixel basis was used to examine the effects of land cover and soil texture on the selected SM maps for each basin and to interpret its spatial patterns [35,36]. The results are presented as SM mean and standard deviation for each class, as well as the level of significance of the differences between classes (Table 2 and Table 3).

Table 2. SM means and standard deviations per land cover type (cm³·cm⁻³).

Land cover	Nestos	Rijnland	Tamega
Rainfed crops	0.11/0.05 a	0.22/0.09 a	0.10/0.04 a
Irrigated crops	0.16/0.08 b	n/a	0.16/0.07 b
Broadleaved forests	0.20/0.06 c	0.18/0.07 b	0.15/0.07 b
Coniferous forests	0.22/0.07 c	0.18/0.06 b	0.12/0.03 c
Shrubs	0.15/0.05 b	n/a	0.11/0.05 a
Pastures	0.11/0.06 a	0.20/0.07 a	0.12/0.06 c

Note: Different letters give significant differences within each column (p < 0.001).

Table 2. SM means and standard deviations per land cover type (cm³·cm⁻³).

Land cover	Nestos	Rijnland	Tamega
Rainfed crops	0.11/0.05 a	0.22/0.09 a	0.10/0.04 a
Irrigated crops	0.16/0.08 b	n/a	0.16/0.07 b
Broadleaved forests	0.20/0.06 c	0.18/0.07 b	0.15/0.07 b
Coniferous forests	0.22/0.07 c	0.18/0.06 b	0.12/0.03 c
Shrubs	0.15/0.05 b	n/a	0.11/0.05 a
Pastures	0.11/0.06 a	0.20/0.07 a	0.12/0.06 c

Note: Different letters give significant differences within each column (p < 0.001).

Table 3. SM means and standard deviations per soil texture class (cm³·cm⁻³).

Soil Texture	Nestos	Rijnland	Tamega
Fine	0.19/0.06 a	0.23/0.11 a	n/a
Medium	0.18/0.07 b	0.22/0.09 b	0.13/0.06
Coarse	0.16/0.06 c	0.18/0.07 c	n/a

Note: Different letters give significant differences within each column (p < 0.001).

Table 3. SM means and standard deviations per soil texture class (cm³·cm⁻³).

Soil Texture	Nestos	Rijnland	Tamega
Fine	0.19/0.06 a	0.23/0.11 a	n/a
Medium	0.18/0.07 b	0.22/0.09 b	0.13/0.06
Coarse	0.16/0.06 c	0.18/0.07 c	n/a

Note: Different letters give significant differences within each column (p < 0.001).

Figure 2. Spatial variability of SM in the study areas: (a) Nestos 24 November 2012, 1 January 2013, 9 May 2013, 18 June 2013; (b) Rijnland 7 October 2012, 25 January 2013, 18 February 2013, 5 August 2013; (c) Tamega 7 October 2012, 24 November 2012, 1 January 2013, 4 August 2013. No-data pixels appear in white.

Figure 2. Spatial variability of SM in the study areas: (a) Nestos 24 November 2012, 1 January 2013, 9 May 2013, 18 June 2013; (b) Rijnland 7 October 2012, 25 January 2013, 18 February 2013, 5 August 2013; (c) Tamega 7 October 2012, 24 November 2012, 1 January 2013, 4 August 2013. No-data pixels appear in white.

The spatial variability of SM in the Nestos basin (Figure 2a) follows the meteorological conditions, as well as the main geomorphological features. Lower values of SM were observed in the areas covered by rainfed crops, shrubs and pastures (Table 2). Significantly higher values were estimated in the forested areas, which appear towards the elevated middle parts of the study area. For the selected dates of Figure 2a, SM was high in most parts on 25 November 2012 following a few rain events and higher on 9 May 2013 after four consecutive rainy days. SM was lower than 0.15 cm³ cm⁻³ in most parts of the study site, with the exception of the river’s narrow valley, which is characterized by higher values of SM even during dry periods, being at the lowest end of the natural drainage system and dominated by medium to fine clays (Table 3). SM values were low in the northern part of the study site on 18 June 2013, but very high in the eastern and southern ends. This map captures the peak of the irrigation season and shows high values of SM at the irrigation network that has been developed in the river’s floodplain at the southern end of the study area.

The spatial variability of SM in Rijnland is guided by the meteorological conditions, as well as the main soil types and shallow water table regime, closely connected to a dense surface water network of lakes and canals (Figure 1). SM is higher throughout the catchment, as compared to the other study areas. The natural vegetation of the north and western end of the study area have significantly lower SM than the pastures at the southeast (Table 2). The areas surrounding the water bodies at the central and western side display the highest SM on all occasions, being dominated by medium to fine texture soils (Table 3). SM on 7 October 2012 (Figure 2b) was high in the study area (0.20–0.35 cm³·cm⁻³) following a few days of rainfall. Similarly high SM appeared on 25 January 2013 in the central and northern parts of the study area (0.18–0.30 cm³·cm⁻³), following a smaller rain event. SM values were much lower (0.15–0.25 cm³·cm⁻³) on 28 February 2013, which followed several dry days. On 5 August 2013, SM was low in the north (0.09–0.22 cm³·cm⁻³) and higher at the southern half of the study area (0.26–0.38 cm³·cm⁻³), following a local rain event, as reported by the meteorological station.

The spatial variability of SM in the Tamega basin follows the meteorological conditions, as well as the main geomorphological features and the high land cover variability (Figure 1). Lower values appear at the mountain tops, where sparse vegetation, barren land and rocks dominate, and significantly higher values at the narrow river valley towards the centerline of the study area (Table 2). The soil texture is uniform, thus appearing to have no impact on the spatial pattern of SM (Table 3). SMs on 7 October 2012 and 1 January 2013 (Figure 2c) were low overall (0.08–0.18 cm³·cm⁻³ and 0.13–0.2 cm³·cm⁻³, respectively), following a dry period, with the exception of the land along the river, where SM was very high (0.23–0.30 cm³·cm⁻³ and 0.31–0.38 cm³·cm⁻³, respectively). On 24 November 2012, SM was very high in the northwest end of the study area, where a recent rainfall had occurred (evident around the cloud-covered pixels), while it was lower at the southern parts. On 4 July 2013, SM was very low throughout the study area (0.08–0.15 cm³·cm⁻³), with the exception of the agricultural land on the southern end, where SM was kept close to 0.20 cm³·cm⁻³ via irrigation.

3.2. Temporal Dynamics of Soil Moisture

Figure 3 shows the time series of SM at selected locations, one per land cover type (Figure 1). SM of rainfed crops in Nestos (Figure 3a) was governed by the temporal distribution of rain events: it showed an undulating pattern with higher values (>0.2 cm³·cm⁻³) during the rainy seasons of autumn (October–November) and spring (February–March). During the dry summer season, SM fell to lower values (0.07–0.17 cm³·cm⁻³), which is, nevertheless, after the harvest of the winter cereals. SM of irrigated crops was similar to the SM of rainfed crops during the wet months (October–May) and then remained high (0.15–0.3 cm³·cm⁻³) during the irrigation season (June–September). SM of broadleaved forests was high during the wet winter season (0.17–0.35 cm³·cm⁻³) and dropped to lower values (0.13–0.25 cm³·cm⁻³) during the dry summer months. The high values during the wet season were due to the north-facing slopes and the higher altitude at which the broadleaved forests appear, where evaporation was lower. SM of coniferous forests was stable throughout the year (0.1–0.3 cm³·cm⁻³), with slightly higher values during the wet spring and late summer, probably due to the frequent storms at the high altitudes in which this vegetation type appears. SM of shrubs and pastures fluctuated in accordance with the temporal distribution of rain events, mostly during the wet season, similar to rainfed crops. SM of pastures was lower than the others with the highest temporal fluctuation, possibly because they are located in shallow soils.

SM in Rijnland was higher, compared to the other study areas, with higher minima (Figure 3b). SM temporal patterns of broadleaved forest, coniferous forest and pastures were similar with a higher temporal variation (standard deviations ranging from 0.06 to 0.09 cm³·cm⁻³). SM of rainfed crops was the highest of all examined land cover types (annual average 0.22 cm³·cm⁻³). SM of all land cover types related directly to the pattern of rainfall during most of the year, with the exception of summer, when the phreatic water level was artificially risen (by either letting in additional water or applying less drainage) and fed the soil moisture, even at times of lower rainfall. Similar temporal patterns had been experienced in a location to the east of Rijnland; however, the recorded SM had reached higher values in the first quarter of 2002 [37]. The exceptionally high soil moisture of mid-September 2013 was due to the heavy rainfall of these days, which had influenced the SM of all examined land cover types.

In Tamega (Figure 3c), there was higher SM and variability (0.05–0.35 cm³·cm⁻³) during the winter–spring wet period (December–March), which dropped during the dryer summer months (0.08–0.15 cm³·cm⁻³). SM rose again in September after the small, but frequent rain events. SM of rainfed crops, coniferous forests, shrubs and pastures followed the pattern of rainfall, with peaks during the wet season. Irrigated crops followed a similar pattern during the winter; however, they displayed constantly higher values (0.14–0.26 cm³·cm⁻³) during the summer irrigation season (June–September). Broadleaved forests also displayed exceptionally high values during the summer months, due to the high altitude and sporadic storms. Despite the frequent precipitation events during November, all land cover types had low SM, due to its low intensity. Similar temporal patterns had been noted for the forested area adjacent to the east catchment for the previous hydrological year, when some extreme SM values had been noted during winter 2011 [8].

Figure 3. Time-series of SM at the selected locations per land cover and rainfall in the study areas: (a) Nestos; (b) Rijnland; and (c) Tamega.

Figure 3. Time-series of SM at the selected locations per land cover and rainfall in the study areas: (a) Nestos; (b) Rijnland; and (c) Tamega.

3.3. Validation and Analysis of Environmental Parameters

Comparison of the remotely-sensed estimations of SM was performed with in situ measurements of SM for the top-soil acquired during the field surveys. The correlation coefficient, the Student’s t-test and the RMSE were calculated. The results are presented in Table 4.

For Nestos, the relationship between remotely-sensed SM estimates and in situ measurements was high and statistically significant during the summer period and low during the winter season (r = 0.62 and 0.36, respectively). It was not possible to establish any relation during the field survey of April 2013, as the satellite images were influenced by consistently cloudy weather. The RMSE was low throughout the seasons, ranging from 0.037 to 0.055 cm³·cm⁻³. The mean differences were mostly positive, revealing an underestimation of remotely-sensed SM as compared to in situ measurements.

Table 4. Correlation coefficients, RMSE and mean differences of field-surveyed vs. remote sensing data.

Study Area	Date of Field Survey	Count	r	RMSE (cm3·cm−3)	Mean difference (cm3·cm−3)
Nestos	13–15 July 2011	20	0.62 *	0.037	0.004
Nestos	3–5 November 2012	22	0.36	0.041	0.057 **
Nestos	9–10 April 2013	27	-	0.055	0.043 *
Nestos	2–3 July 2013	24	0.42 *	0.039	−0.042 **
Rijnland	17–19 July 2012	17	0.37	0.076	0.061 **
Rijnland	26–29 June 2013	25	0.31	0.087	0.091 **
Rijnland	23–26 September 2013	22	0.52 *	0.086	0.039
Tamega	26–29 September 2011	28	0.43 *	0.015	−0.062 **
Tamega	20–23 May 2013	38	0.31	0.023	0.044 **
Tamega	2–5 September 2013	38	-	0.024	−0.070 **

Note: * Statistically significant at the 0.05 level; ** statistically significant at the 0.001 level.

Table 4. Correlation coefficients, RMSE and mean differences of field-surveyed vs. remote sensing data.

Study Area	Date of Field Survey	Count	r	RMSE (cm3·cm−3)	Mean difference (cm3·cm−3)
Nestos	13–15 July 2011	20	0.62 *	0.037	0.004
Nestos	3–5 November 2012	22	0.36	0.041	0.057 **
Nestos	9–10 April 2013	27	-	0.055	0.043 *
Nestos	2–3 July 2013	24	0.42 *	0.039	−0.042 **
Rijnland	17–19 July 2012	17	0.37	0.076	0.061 **
Rijnland	26–29 June 2013	25	0.31	0.087	0.091 **
Rijnland	23–26 September 2013	22	0.52 *	0.086	0.039
Tamega	26–29 September 2011	28	0.43 *	0.015	−0.062 **
Tamega	20–23 May 2013	38	0.31	0.023	0.044 **
Tamega	2–5 September 2013	38	-	0.024	−0.070 **

Note: * Statistically significant at the 0.05 level; ** statistically significant at the 0.001 level.

In Rijnland, the relationship between SM estimates and measurements was high in June and July (r = 0.31 and 0.37, respectively) and even higher and statistically significant in September (r = 0.52). The RMSE and mean differences were amongst the highest across the study sites (<0.087 cm³·cm⁻³ and <0.091 cm³·cm⁻³, respectively). The mean differences were consistently positive, revealing an underestimation of remotely-sensed SM as compared to in situ measurements.

In Tamega, the relationship between SM estimates and measurements was high and statistically significant in September 2011 and medium in May 2013, but it was not possible to establish any relation during September 2013, when soil moisture was too low, almost reaching the low measurable range of the TDR. The RMSE was low (<0.024 cm³·cm⁻³), and the mean differences were medium to low, but of a variable sign. Similar accuracy was reported for the catchment adjacent to Tamega and the previous hydrological year, where errors were lower during the spring and summer seasons [8].

The relationships between SM estimates and in situ measurements were variable. This can be explained by several factors, such as the temporal difference between the datasets, rain events during the measurements and the scale difference. The effect of these factors was evaluated using one-way analysis of variance (ANOVA) of the residuals between measured and estimated SM, for a selected field survey, per study area. The results are presented as p-values and the level of significance of the parameter’s effect, where low p-values denote a high confidence that there is an actual difference caused by the parameter (Table 5).

MODIS data were acquired daily for each location and then put together in eight-day composites, which means that each pixel of the SM map may refer to a different day within the respective eight-day period. Therefore, the estimated SM may be as much as eight days apart from the measured SM. At the same time, soil moisture has been reported to vary greatly, even daily, depending on topography, soil properties and meteorological conditions [38].

The temporal difference between the measured and estimated SM ranged from −2 to 3 days in Nestos, from 3 to 5 in Rijnland and from −3 to 1 days in Tamega (negative temporal difference means that the estimation preceded the measurement). This effect has significant influence at the Rijnland and Tamega sites (p = 0.003 and 0.026, respectively), but not at Nestos (p = 0.389), possibly because the meteorological conditions have increased the influence of temporal variation of SM in the first two sites.

Field surveys required from 2 to 4 days to cover each case study catchment. During this time, there was a possibility of rain on the whole or part of the site. This was expected to influence the quality of the measurements, especially if the rain event were between the times of in situ measurement and satellite overpass.

During the selected field surveys, there was a single rain event in Tamega and multiple ones in Rijnland, while there was dry weather in Nestos. The rain events seem to significantly influence the quality of the results in Rijnland (p = 0.023), but not in Tamega (p = 0.114), probably because the rain event in Tamega was very light.

The field surveys provided point measurements, spread around the sampling location to describe at an adequate level the local variability of SM. However, the comparison with the large MODIS pixels (250 m) was expected to show a deviation, especially at the most heterogeneous locations. The level of heterogeneity of land cover at each location was evaluated with the landscape metric Shannon’s Diversity Index (SHDI) [39]. The SHDI of each location was compared to the residuals to identify the level of influence of the scale difference between the point measurements and the MODIS observations, providing no significant influence on the quality of the results (p > 0.501 at all sites). This verified that the measurement locations were correctly selected in large homogeneous areas.

Table 5. Analysis of variance and the effect of various parameters on the residuals at all sites (p-values).

Study Area	Date of Field Survey	Temporal Difference	Rain Event	Location Heterogeneity
Nestos	13–15 July 2011	0.389	-	0.618
Rijnland	17–19 July 2012	0.003 *	0.023 *	0.712
Tamega	26–29 September 2011	0.026 *	0.114	0.501

Note: * Statistically significant at the 0.05 level.

Table 5. Analysis of variance and the effect of various parameters on the residuals at all sites (p-values).

Study Area	Date of Field Survey	Temporal Difference	Rain Event	Location Heterogeneity
Nestos	13–15 July 2011	0.389	-	0.618
Rijnland	17–19 July 2012	0.003 *	0.023 *	0.712
Tamega	26–29 September 2011	0.026 *	0.114	0.501

Note: * Statistically significant at the 0.05 level.

4. Discussion

The temporal variability of SM is routinely measured by automated soil moisture probes, which is useful for a field or other homogeneous area. A time series of spatially-distributed measurements for the entire catchment is, however, not feasible. The importance of this work consists in providing continuous time series of SM maps at the catchment scale. The main advantages of the methodology are that it avoids some basic problems reported in SM mapping with remote sensing, such as the bias of dense vegetation and surface roughness of microwave images, and the influence of surface reflectance variability of the visible and near-infrared images [40,41].

Input MODIS products surface reflectance and LST were eight-day composites, which offer the advantage of improved data quality at the pixel level. However, some cloud effects cannot be successfully removed, such as the cloud shadow [42] and artificial heterogeneity [43]. Persistent cloud cover for more than eight days could be dealt with by using longer compositing periods (e.g., 16 days or one month); however, there is a risk of losing important phenological changes of crops with short growing seasons [44,45].

The partitioning of surface fluxes, in this case Λ, was used as an indicator of SM based on the argument that the observed difference between the actual and potential evapotranspiration is caused by the available water in the root zone [13]. The importance of SM on the evapotranspiration process of plants has been reported positively in dry and negatively in wet climates [37]. Nevertheless, similar methods have provided successful results in various climates (Pakistan, Mexico, Brazil, Ghana and others), at extreme conditions (near the permanent wilting point of crops or at field capacity), using Landsat and MODIS satellite images [13,46,47].

Deviations of remote sensing estimates from in situ measurements were comparable to those obtained in the Iberian peninsula by the disaggregation of SMOS data to 500 m [8] and with those reported for the SMOS products [48]. A similar range of accuracy was reported by using surface energy fluxes in Pakistan and Mexico [13]. The almost consistent underestimation of SM by the satellite estimations as compared to in situ measurements was also noted by Sánchez–Ruiz, Piles, Sánchez, Martínez–Fernández, Vall–llossera and Camps [8], while Fang and Lakshmi [9] produced mixed results.

The observed higher temporal variability during the rainy season could be due to errors introduced by the frequent cloud cover and variable atmospheric conditions, the low performance of the MODIS compositing algorithm, as well as the uncertainty of the meteorological model used for input. Thus, there are occasions when the compositing algorithm underperforms and errors are introduced in the remote sensing estimations [49]. Nevertheless, using composited data is a routine method to overcome frequent cloud coverage and to minimize the gaps in time series, and techniques have been proposed to estimate the introduced error [50]. Akuraju et al. [51] point out some limitations in the use of an exponential model to fit soil moisture in the root zone to Λ. Indeed, the relation is sensitive to vegetation biomass, with a correlation significantly higher when there is enough biomass. In addition, Λ values depend on the ET limiting conditions (water or energy conditions); thus, Akuraju et al. (2013) advise to adapt the model to the NDVI and net radiation values, whether they are above or below certain thresholds. A time lag between rainfall events and SM increase can be noted in some cases, similar to related literature [8,41]. This could be due to the compositing within eight days of the satellite observations and the meteorological data. The in situ measurements were performed with a TDR that measures topsoil moisture (0–7.5 and 0–20 cm) and were compared to the estimated root zone SM in the study areas, which can vary from 0–0.3 m for grassland to 0–2 m for trees. Although they can vary differently in time, due to their non-linear relation, these two parameters can be compared, since each test site is located within a certain geoclimatic environment [16], under homogeneous bioclimatic conditions [52]. The comparison is more efficient in wet conditions [35]. Ford et al. [53] demonstrated that there is generally a strong relationship between near-surface (5–10 cm) and root zone (25–60 cm) soil moisture in a study conducted in an extensive area (Great Planes, USA), while Ragab [18] used an empirical relationship to estimate root-zone soil moisture from surface soil moisture. Using a limited number of in situ observations to validate a remote sensing product is acceptable, considering the associated cost of wide-scale field surveys and the spatial representativeness of the sampling design [54,55].

The presented methodology is a semi-automated procedure. Several tools have been prepared to facilitate the operation by the analyst [33,34,56]. It uses freely available data at no cost from operational satellites, meteorological models and soil data. Alternative sources of satellite data can be used, such as Aqua MODIS, NOAA AVHRR or other satellites with visible, near and thermal infrared bands.

The presented eight-day SM maps are essential for water and land management studies. They are useful as inputs in meteorological, hydrological, gas and energy exchange models. Yet, SM is very dynamic and influenced by multiple factors, thus difficult to model [57]. They can be directly inserted into models or they can provide a means for calibration and operational adjustment (through data assimilation) of the spatial distribution of hydrological parameters [2,3,4,18].

5. Summary and Conclusions

An integrated methodology was presented that connects the evaporative fraction, as an indicator of SM, with MODIS satellite images, which offer frequent, global, validated and free of cost data and, therefore, allow for operational monitoring applications.

The methodology was tested in three diverse European environments. The resulting maps were successfully validated using in situ measurements, slightly underestimating SM, but the RMSE was comparable to similar attempts found in the literature.

Several factors were analyzed to explain the variation in the accuracy, and it was shown that the land cover type, the soil texture class, the temporal difference between the datasets acquisition and the presence of rain events during the measurements played a significant role, unlike the scale difference between in situ and satellite observations.

Therefore, the proposed methodology can be used operationally to estimate SM maps at the catchment scale, with a 250-m spatial resolution and eight-day time step.