Multi-Temporal Evaluation of Soil Moisture and Land Surface Temperature Dynamics Using in Situ and Satellite Observations

Abstract

Soil moisture (SM) is an important component of the Earth’s surface water balance and by extension the energy balance, regulating the land surface temperature (LST) and evapotranspiration (ET). Nowadays, there are two missions dedicated to monitoring the Earth’s surface SM using L-band radiometers: ESA’s Soil Moisture and Ocean Salinity (SMOS) and NASA’s Soil Moisture Active Passive (SMAP). LST is remotely sensed using thermal infrared (TIR) sensors on-board satellites, such as NASA’s Terra/Aqua MODIS or ESA & EUMETSAT’s MSG SEVIRI. This study provides an assessment of SM and LST dynamics at daily and seasonal scales, using 4 years (2011–2014) of in situ and satellite observations over the central part of the river Duero basin in Spain. Specifically, the agreement of instantaneous SM with a variety of LST-derived parameters is analyzed to better understand the fundamental link of the SM–LST relationship through ET and thermal inertia. Ground-based SM and LST measurements from the REMEDHUS network are compared to SMOS SM and MODIS LST spaceborne observations. ET is obtained from the HidroMORE regional hydrological model. At the daily scale, a strong anticorrelation is observed between in situ SM and maximum LST (R $\approx -$0.6 to −0.8), and between SMOS SM and MODIS LST Terra/Aqua day (R $\approx -$ 0.7). At the seasonal scale, results show a stronger anticorrelation in autumn, spring and summer (in situ R $\approx -$0.5 to −0.7; satellite R $\approx -$0.4 to −0.7) indicating SM–LST coupling, than in winter (in situ R ≈ +0.3; satellite R $\approx -$0.3) indicating SM–LST decoupling. These different behaviors evidence changes from water-limited to energy-limited moisture flux across seasons, which are confirmed by the observed ET evolution. In water-limited periods, SM is extracted from the soil through ET until critical SM is reached. A method to estimate the soil critical SM is proposed. For REMEDHUS, the critical SM is estimated to be ∼0.12 m${}^3$/m${}^3$, stable over the study period and consistent between in situ and satellite observations. A better understanding of the SM–LST link could not only help improving the representation of LST in current hydrological and climate prediction models, but also refining SM retrieval or microwave-optical disaggregation algorithms, related to ET and vegetation status.

1. Introduction

Soil moisture (SM) is an essential climate variable that plays an important role in both the Earth’s water and carbon cycles [1] . SM regulates the net ecosystem energy exchange at the land-atmosphere boundary layer over the continents through the heat fluxes (latent and sensible) and its feedback on precipitation. Also, it influences hydrological processes, runoff generation, flood or drought development, and agricultural productivity [2]. A better understanding of these processes and their couplings is needed for the improvement of meteorological and climate forecasts, and for a more sustainable use of water resources.

The scientific community considers L-band (1–2 GHz) the optimal bandwidth for measuring the Earth’s global surface SM, penetrating approximately the top 5 cm of the soil. At these frequencies, the amount of thermal radiation naturally emitted by the Earth’s surface depends mainly on its physical temperature and the soil dielectric constant (ϵ), which for the near surface soil layer is primarily a function of the soil water content [3]. Active sensors (radars) consist of a transmitter and a receiver that measures the signal that is reflected, refracted or scattered; while passive sensors (radiometers) consist only of a receiver that measures the radiation naturally emitted from the scene under observation. L-band radiometers are based on the large contrast between the dielectric properties of dry soil ($ϵ\sim$ 4) and water ($ϵ\sim$ 80); the Earth’s surface emissivity decreases with increasing SM [4,5]. Other advantages of L-band are: (i) it is significantly less affected by rain and atmospheric effects than higher microwave frequencies, with the atmosphere being nearly transparent [6]; (ii) the radiation emitted by the Earth’s surface can pass through sparse and up to moderate canopies (i.e., those with Vegetation Water Content (VWC) ≤ 5 kg/m${}^2$, corresponding to 70% of non-frozen land areas on Earth, excluding dense forest), but even vegetation with lower VWC attenuates the signal and, consequently, it needs to be accounted for in SM retrieval algorithms [7]; and (iii) since L-band is part of the microwave frequency range, the measurements are also independent of solar illumination.

Currently, there are two L-band satellite missions in orbit specifically devoted to monitoring the Earth’s global surface SM: (i) the Soil Moisture and Ocean Salinity (SMOS, 2009–2017) mission from the European Space Agency (ESA); and (ii) the Soil Moisture Active Passive (SMAP, 2015–2018) mission from the U.S. National Aeronautics and Space Administration (NASA).

SMOS, launched on 2 November 2009, is the first satellite to globally measure surface SM [8] and sea surface salinity [9]. Its single payload is the two-dimensional Microwave Imaging Radiometer with Aperture Synthesis (MIRAS), an interferometric radiometer which provides multi-angular (0°–65°) and full-polarimetric observations with a spatial resolution of ∼35–50 km [10].

SMAP was launched on 31 January 2015, to globally measure surface SM and the freeze/thaw state of the soil [11]. It includes a real aperture radiometer and a Synthetic Aperture Radar (SAR), providing full-polarimetric observations at a single look angle (40°) with spatial resolutions of ∼36 km and ∼3 km, respectively [12]. SMAP planned to provide a 9-km SM product using a disaggregation algorithm that combines active-passive observations [13]. However, due to a failure on the data transmission of the SMAP SAR on 7 July 2015, the duration of active observations is limited. Still, there are ten weeks of combined measurements to allow new insights in active-passive SM remote sensing.

Land surface temperature (LST) is also a key climate variable. Remotely sensed LST, acquired with thermal infrared (TIR) sensors, are routinely used in many operational applications, including weather forecasting. Polar orbiting satellites, such as NASA’s Terra (1999–2020) and Aqua (2002–2020) missions, provide global LST maps at 1 km every 1–2 day using the Moderate-resolution Imaging Spectroradiometer (MODIS) sensor [14]. Geostationary satellites (i.e., those with the same orbital period as the Earth’s rotation) are able to provide continental coverage at higher temporal resolution. The Spinning Enhanced Visible and InfraRed Imager (SEVIRI) sensor from the Meteosat Second Generation (MSG) mission (2002–2021) of ESA & European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) provides continental LST maps at 3 km every 15 min [15]. However, the effective temporal resolution of both polar and geostationary LST observations can be much lower than the theoretical value since optical observations, and particularly TIR observations, are masked in presence of cloud cover.

SM and LST are closely inter-related and their link is central to land-surface variables and processes such as soil emission, evapotranspiration (ET) and thermal inertia. Soil emission is influenced by the SM and soil physical temperature. Owing to this relationship, LST is a necessary input of any microwave SM retrieval algorithm [16,17,18]. LST and SM are also related through ET, which is a fundamental variable in hydrological and water/energy balance studies [19,20,21]. According to the latter, there are two different ET regimes: the energy-limited and the water-limited. In the energy-limited regime, SM is above a certain SM value (critical SM) and ET is mainly controlled by the available energy (not by the soil water content). In the water-limited regime, in contrast, SM is below the critical SM and acts as the main driver of the ET process. The critical SM is known to be bounded by the soil wilting point (WP), which is the lower limit of the plant available water, and the soil field capacity (FC), defined as the upper limit of SM available to plants. It typically lies between the 50% to 80% of FC [22]. The thermal inertia describes the resistance of a material to temperature variations and is proportional to the material’s thermal conductivity, density and specific heat capacity. In general, a decrease in SM produces an increase in LST diurnal range, and the maximum LST under moist conditions occurs later than under dry conditions (with a time lag within the day) [23]. This is because the specific heat capacity of water is typically higher than that of bare soil and, consequently, dry soils respond faster than wet soils to temperature variations.

Several studies have been focused on the indirect estimation of SM using TIR remote sensing [24]. Because the surface brightness temperature depends on its physical temperature and also on the surface water content and vegetation, most of these approaches are based on the so-called triangle or trapezoidal method of the LST/vegetation index space [25,26,27,28]. Usually, it is used with an ET model, e.g., the Atmosphere-Land Exchange Inverse (ALEXI) model [29], or the thermal inertia [30]. An intercomparison study showed that SM estimates from microwave and TIR provide complementary information [31]. Microwave retrievals are optimal for measuring accurate SM under all-weather conditions, but the spatial resolution obtained is limited (around tens of km). TIR retrievals, in turn, provide high resolutions (from tens of meters to several km), but are sensitive to atmospheric effects. To benefit from the advantages of both SM retrieval techniques, there are several SM downscaling or disaggregation algorithms based on the synergy of coarse resolution microwave and high resolution TIR observations [32,33,34]. A recent review describes these relationships and the synergistic use of passive microwave and optical Visible/InfraRed (VIS/IR) for multi-scale SM estimation [35].

Most radiometer-based and TIR-based SM retrieval algorithms employ the instantaneous LST in the inversion. However, some studies have revealed that a composite of maximum apparent temperature provides a closer representation of land surface conditions [36,37]. Also, the strongest correlation was observed between SMOS Level 2 (L2) SM and the MODIS LST acquisition closer to the daily maximum (e.g., Aqua day) in a previous study [38]. Consistently, better agreement to in situ was also observed in SMOS/MODIS SM disaggregated products when using LST at day times (closer to daily maximum) than at night times [39]. The above-mentioned studies suggest that daily maximum LST is more closely linked to instantaneous SM than collocated LST. In this regard, gaining understanding of the SM–LST link from observations could help improving the representation of LST in current hydrological and climate prediction models, and also refine SM retrieval or microwave-optical disaggregation algorithms, related to ET and vegetation status.

This study presents an evaluation of SM and LST dynamics using in situ and satellite observations over the central part of the river Duero basin in Spain. The study period comprises four years, from 2011 to 2014, and provides an assessment of the role of ET and thermal inertia physical processes into the SM–LST link. The in situ data include SM and LST observations from the Soil Moisture Measurement Stations Network of the University of Salamanca (REMEDHUS) [40]. The satellite data used are the SMOS Barcelona Expert Centre (BEC) Level 3 (L3) SM [41] and the MODIS Terra/Aqua LST. The correlation of instantaneous SM and a variety of LST-derived parameters is analyzed at the daily scale to better understand the fundamental link of the SM–LST relationship. The seasonal dependence of this relationship is evaluated by means of SM and LST temporal correlation using a modeled ET as supporting information. The Hydrological Model for Operational estimates of Recharge and actual Evapotransporation (HidroMORE) is used to generate the modeled ET [42,43]. Additionally, a method to estimate the critical SM from the SM–LST coupling/decoupling transitions has been developed and effectively applied to in situ and spaceborne measurements.

Section 2 provides a description of the in situ and satellite data used as well as a brief explanation of the ET model. Section 3 describes the applied methodology. Results and discussion are shown in Section 4. Main conclusions and contributions of this work are outlined in Section 5.

2. Data

2.1. In Situ Data

The REMEDHUS network [40] is located at the central part of the river Duero basin in Spain. It covers a semi-arid continental-Mediterranean agricultural region of 35 km × 35 km (41.1–41.5° N, 5.1–5.7° W), with an area of approximately 1300 km${}^2$. REMEDHUS is a member of the International Soil Moisture Network (ISMN) and includes 21 permanent stations providing hourly SM and LST measurements in the top 5 cm of the soil [44,45]. Both are acquired with Hydra Probes that measure with an accuracy of 0.003 m${}^3$/m${}^3$ for SM and 0.6 °C for LST [46]. The main land uses in this area are vineyard (stations E10, F6, H7, K6, J3, and L3), rainfed cereals or fallow (stations F11, H13, J12, J14, K9, K10, K13, L7, M5, M9, N5, and O7), and forest-pasture (stations H9 and M13). During the study period, station K4 changes its land use from vineyard (2011) to rainfed cereals or fallow (2012, 2013 and 2014). The SM and LST measurements of REMEDHUS from a particular year corresponding to stations over irrigated crops (K9 during 2012, K13 during 2013 and 2014) have been removed from the dataset. The reason for filtering irrigated crops is the quick SM fluctuations related with the high frequency of water supply (less than 12 h, specially during the summer), which are not captured by the SMOS 3-day repeat. Stations with corrupted LST (J12 during 2011) or corrupted SM (M13 during 2014) have also been discarded. In REMEDHUS there are also four automatic weather stations that measure air temperature, relative humidity, precipitation, solar radiation, reference or potential evapotranspiration ($E{T}_0$) and wind speed. The monthly accumulated $E{T}_0$ and the daily mean precipitation have also been used in this study. Figure 1 (left) shows a land use map of the REMEDHUS region, including the 21 SM and LST stations and the 4 weather stations, depicted with black triangles and white dartboards, respectively.

2.2. Satellite Data

The SMOS SM retrieval algorithm is based on the minimization of a cost function whose main component is the sum of the squared weighted differences between measured and modeled brightness temperatures ($T_B$). The target accuracy of the global SMOS L2 SM maps is 0.04 m${}^3$/m${}^3$ [17]. The SMOS BEC Level 3 (L3) SM v.1 products acquired over the Iberian Peninsula are freely distributed at BEC [47]. These products were obtained by quality-filtering and re-gridding from the operational ESA SMOS L2 v.5.51 products. Grid points affected by radio frequency interference (RFI) and/or SM with Data Quality Index (DQX, defined as the error standard deviation) greater than 0.07 m${}^3$/m${}^3$ are discarded. A DQX-inverse weighted average is applied to bin the data from the Icosahedral Snyder Equal Area (ISEA 4H9) grid to the global cylindrical 25-km Equal-Area Scalable Earth (EASE ML) grid. More details are available in [41]. A comprehensive validation of SMOS L2 and L3 SM products using two complementary small-scale and large-scale in situ networks (REMEDHUS and Inforiego, respectively) and a surface water balance model can be found in [48].

The MODIS LST v.5 products from Terra (MOD11A1) and Aqua (MYD11A1) satellites are provided by the U.S. Land Processes Distributed Active Archive Center (LP DAAC) [49]. They have 1 km of spatial resolution and a nominal accuracy of 1 °C under clear sky conditions. Higher error (from 4 °C to 10 °C) could be expected in presence of clouds and heavy aerosols.

Figure 1 also displays maps of 1-km MODIS LST Aqua day (top right) and 25-km SMOS BEC L3 SM (bottom right) over the Iberian Peninsula, both time-averaged during the study period. Note the different size of MODIS and SMOS pixels.

2.3. HidroMORE Model

The real or actual evapotranspiration ($E{T}_a$) during 2014 is obtained from the HidroMORE hydrological model [42,43]. It is an operational distributed model based on the water balance equation [50]. HidroMORE provides daily estimates of deep percolation (DP), water storage (WS) and ET after considering effective precipitation (P), surface runoff (RO) and irrigation (I):

$$P+I-RO-ET-DP=\pm \mathrm{\Delta }WS$$

(1)

where $\pm \mathrm{\Delta }WS$ is a residual of the balance and expresses the rate of change in soil moisture content or water storage [mm]. The calculation in Equation (1) accounts for the effective root zone.

The rationale of the $E{T}_a$ estimation is based on the dual crop coefficient-reference ET methodology [50]:

$$E{T}_a=E{T}_0(K_s{K}_{cb}+K_e)$$

(2)

where $E{T}_0$ stands for the reference or potential evapotranspiration, $K_s$ is the water stress coefficient, $K_{cb}$ is the basal crop coefficient, and $K_e$ is the soil evaporation coefficient. The $K_{cb}$ is calculated on the basis of the daily Normalized Vegetation Difference Index (NDVI), retrieved from temporal series of optical imagery (NDVI–$K_{cb}$ approach), $K_e$ is computed from the Food and Agriculture Organization (FAO) water balance procedure, and water stress conditions are applied to estimate $K_s$.

The $K_{cb}$ approach has been previously applied with different formulations (linear and exponential) and also using different indirect growing parameters, such as the Leaf Area Index (LAI), the Fractional Vegetation Cover (FVC) or the NDVI. Results show that the linear approach is adequate for SM estimation at the root zone, which is denoted by $\pm \mathrm{\Delta }WS$ in Equation (1). The $E{T}_a$ differences produced by the indirect growing parameter are lower than 3 mm of the total annual $E{T}_a$) [51]. Ground-based SM measurements from REMEDHUS have also been used for validating these results [52]. The higher correlations obtained between in situ SM and HidroMORE SM at the root-zone (R ≈ 0.78 to 0.95) indicate that the model was well calibrated and working correctly [51].

3. Methodology

3.1. Analysis of SM–LST Relationship

Ground-based SM measurements from the REMEDHUS network have been averaged to mimic the instantaneous SM acquired during SMOS morning/afternoon passes over the region. Ground-based LST measurements at a variety of acquisition times have also been selected, leading to six LST-derived parameters: (i) instantaneous $T_i$; (ii) daily mean $\overline{T}$; (iii) daily median $T_{med}$; (iv) daily maximum $T_{max}$; (v) daily minimum $T_{min}$; and (vi) diurnal range $\mathrm{\Delta }T=T_{max}-T_{min}$.

Regarding satellite data, the four 25-km SMOS L3 SM pixels covering the REMEDHUS area have been extracted in both morning and afternoon SMOS passes. The SMOS equatorial crossing times are: 6:00/18:00 UTC for morning (ascending)/afternoon (descending) orbits. In the case of MODIS LST, all 1-km MODIS LST pixels (≈552 to 576) within the corresponding SMOS pixel have been averaged in order to obtain four aggregated 25-km LST pixels per each MODIS platform (Terra and Aqua) and overpass (day and night). The MODIS equatorial crossing times are: 10:30/22:30 for Terra day (descending)/night (ascending), and 13:30/01:30 for Aqua day (ascending)/night (descending). The aggregated 25-km MODIS pixels with a LST lower than −5 °C have been filtered to screen out measurements affected by frozen soil. In addition, a 3-day averaging window has been applied to SMOS and MODIS data to reduce the noise in the time-series. This averaging does not affect results and neither conclusions obtained from satellite data analysis.

The SM–LST relationship has been assessed for: (i) in situ SM and LST from the REMEDHUS stations and (ii) satellite SMOS SM and MODIS LST over the REMEDHUS area. To do so, the temporal correlation of the two variables has been computed at daily and seasonal scales. As SM and LST annual cycles are expected to be highly related, their anomaly values have been used to detrend the time-series and perform a separate analysis. The anomaly has been obtained by subtracting the mean seasonal cycle of each variable, which represents the annual cycle. This mean seasonal cycle (or seasonal climatology) has been computed as the monthly mean value of the variable during the entire study period. In all cases, non-significant correlations (at the 95% of significance level, $p_{value}>0.05$) have been discarded. Normalized occurrence frequency density diagrams of the SM–LST space of both in situ and satellite data have also been computed.

3.2. Analysis of SM–LST Coupling/Decoupling and Critical SM

The temporal evolution of SM and LST has been analyzed using time-series of in situ and satellite observations together with the daily mean precipitation. The latter is computed from the average of the four weather stations located within the REMEDHUS area. A method to get an estimate of the critical SM (i.e., the SM value at the SM–LST coupling/decoupling transitions) is proposed. As a first step, both SM and LST are normalized as:

$$X_N=\frac{X-X_{min}}{X_{max}-X_{min}}$$

(3)

where X corresponds to the variable (SM or LST) and $X_{min}$ and $X_{max}$ are its minimum and maximum values along the entire time-series, respectively.

As a second step, the SM–LST transition points are identified as those in which the normalized SM and LST have a difference of less than 0.01. Then, the SM and LST values of these transition points are denormalized to compute their maximum, minimum, mean and standard deviation. The critical SM and its uncertainty can be inferred from these statistics.

3.3. Analysis of Water/Energy-Limited ET Regimes

In this study, 11 images from Landsat 8 at Level 1T during 2014 have been used in the HidroMORE simulation. Firstly, we have obtained reflectances at the top of the atmosphere (TOA). They have been converted to bottom of the atmosphere (BOA) reflectances through a relative atmospheric correction and the NDVI has been obtained and cloud-masked at a spatial resolution of 30 m. Later, the NDVI maps from Landsat together with other datasets, such as a land cover map at the same spatial resolution, climatic data and soil data, have been used as input for HidroMORE. The climatic data has been obtained from the REMEDHUS weather stations. The soil database came from the Sistema de Información Agroclimática para el Regadío (SiAR) service of the Spanish Ministry of Agriculture, Food and Environment [53].

A linear relationship between the basal crop coefficient $K_{cb}$ and NDVI has been computed. This NDVI–$K_{cb}$ approach has been applied in the HidroMORE simulation as follows:

$$K_{cb}=1.36·NDVI-0.03$$

(4)

The water balance has been calculated for each day and cell, and results at each REMEDHUS station location have been extracted. Runoff has not been considered in this simulation, owing that in this area, the topographic features and the soil properties make the infiltration the predominant process. HidroMORE considers the root-zone soil layer and calculates SM and $E{T}_a$, but our interest lies in the comparison between ET and the SM at the surface layer. Thus, the SM at root-zone from HidroMORE has not been used.

Since the SM–LST coupling/decoupling periods are an indicator of a different ET regime (water/energy-limited, respectively), the behavior of the ET has been studied over the year 2014, comparing the mean potential evapotranspiration ($E{T}_0$) of the four weather stations from the REMEDHUS network with the actual evapotranspiration ($E{T}_a$) obtained from the HidroMORE hydrological model. The correlation between the $E{T}_a$ from HidroMORE and SM from REMEDHUS has also been analyzed for the different ET regime periods during 2014.

The flowchart in Figure 2 graphically displays the data and methodology approach followed in this study.

4. Results and Discussion

4.1. Characterization of SM and LST over the Study Area

Figure 3 shows the diurnal and the annual cycles of in situ SM (blue line) and LST (red line) observations from the REMEDHUS network during the study period (2011–2014), averaged over all stations. The error bar indicates the standard deviation (std) among the 21 stations. In the diurnal cycle (left), the hourly mean LST exhibits a considerable diurnal range (from ∼13 °C to ∼19 °C), and the daily minimum and maximum LST ($T_{min}$ and $T_{max}$) occur at 7 h and at 16 h, respectively. These times approximately correspond to the SMOS morning and afternoon passes over the REMEDHUS region. Additionally, the LST variability among the stations is higher when the LST is high (std LST ∼2 °C) than when the LST is low (std LST ∼ 1 °C). By contrast, the hourly mean SM remains almost constant (∼0.13 m${}^3$/m${}^3$); its variation along the day is within the accuracy of the sensor, and the SM variability among the stations is low (std SM ∼ 0.09 m${}^3$/m${}^3$). The stability of the daily SM indicates that a unique SM value per day could be representative of the entire day, except when rainfall events occur. In this line, previous research studied the temporal evolution of SM in this region and verified its persistence for a few hours per day, and even during a few days [54]. Also, in a SMOS L2 SM validation study at the same area, no differences where detected between using the time-overpass or the daily average of the in situ SM values [40]. In the annual cycle (right), the monthly mean of SM and LST present opposite behaviors, with the lowest SM and the highest LST during summer months (June, July and August), and the highest SM and the lowest LST during winter months (December, January and February). The SM annual cycle ranges from ∼0.07 m${}^3$/m${}^3$ to ∼0.2 m${}^3$/m${}^3$ and the LST annual cycle ranges from ∼6 °C to ∼27.7 °C. The variability in SM is higher in winter (std SM ∼ 0.12 m${}^3$/m${}^3$) and lower in summer (std SM ∼ 0.05 m${}^3$/m${}^3$), whereas variability in LST is lower in winter (std LST ∼ 0.8 °C) and higher in summer (std LST ∼ 2 °C). These values are consistent with the temporal evolution of SM from the ground-based measurements from the last 15 years of data recorded in REMEDHUS [55] and the climate type of this region.

The possible relation of SM and LST dynamic ranges with soil texture is explored in Figure 4. For the case of SM (left), a very strong anticorrelation with the soil sand content (red asterisks, R = −0.93) and a strong correlation with the soil clay content (blue crosses, R = 0.86) is observed. Two stations (H9 and M13) present a higher SM dynamic range compared with the other stations. This can be explained by the fact that they are forest-pasture stations. Also, they are located at the bottom of valleys where the water table is shallow in winter, and where occasional flooding occurs [40]. In the case of LST (right), a non-significant correlation of LST dynamic range with the soil sand or clay content is found.

4.2. SM–LST Relationship at the Daily Scale

The correlation of SM acquired at SMOS morning/afternoon passes and LST-derived parameters is shown in Figure 5 (left). In the boxplot, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers, and the outliers are plotted individually using red crosses. Considering the entire study period, SM and LST variables are in general anticorrelated: if SM decreases, LST increases, and vice versa. This result agrees with those obtained in Figure 3 (right). There is a single outlier with a very low correlation (R ≈ −0.2 to 0.1) corresponding to H7, a vineyard station that recorded extremely dry conditions (with SM ranging from 0 to 0.08 m${}^3$/m${}^3$ along the four years). In the morning passes, the anticorrelation values (in median) sorted from strong to weak are obtained for: daily maximum LST and daily mean LST (R${}_{T_{max},\overline{T}}\approx$−0.71), daily median LST (R${}_{T_{med}}\approx$ −0.70), instantaneous LST and daily minimum LST (R${}_{T_i,T_{min}}\approx$ −0.68), and LST diurnal range (R${}_{\mathrm{\Delta }T}\approx$ −0.59). In the afternoon passes, the sorted R values are obtained for: daily median LST (R${}_{T_{med}}\approx$ −0.70), daily mean LST (R${}_{\overline{T}}\approx$−0.69), daily maximum LST and instantaneous LST (R${}_{T_{max},T_i}\approx$−0.68), daily minimum LST (R${}_{T_{min}}\approx$−0.66), and LST diurnal range (R ${}_{\mathrm{\Delta }T}\approx$−0.55). Note that the correlation obtained for instantaneous LST becomes equal to the correlation obtained for daily minimum LST in the morning, and equal to the correlation obtained for daily maximum LST in the afternoon. This could be due to $T_i$ being closer in time and similar to $T_{min}$ values at the morning, and to $T_{max}$ values at the afternoon passes, respectively. Considering the percentiles and whiskers, the strongest anticorrelation is obtained for the daily maximum LST ($T_{max}$) both in morning and afternoon passes. This can be explained by the time of $T_{max}$ being also the time of maximum potential evapotranspiration, i.e, when there is a higher atmospheric demand of water and the air turbulence (uplift) is maximum.

When analyzing the detrended in situ time-series, Figure 5 (right) shows that, as expected, correlations are weaker than when using the variables directly. In the morning passes, the sorted anticorrelation values (in median) from strong to weak are: LST diurnal range (R${}_{\mathrm{\Delta }T}\approx$−0.33), daily maximum LST (R${}_{T_{max}}\approx$−0.31), daily mean LST (R${}_{\overline{T}}\approx$−0.26), daily median LST (R${}_{T_{med}}\approx$−0.25), and instantaneous LST and daily minimum LST (R${}_{T_i,T_{min}}\approx$−0.20). In the afternoon passes, the sorted anticorrelation values are obtained for: LST diurnal range (R${}_{\mathrm{\Delta }T}\approx$−0.31), instantaneous LST (R${}_i\approx$ 0.30, daily maximum LST (R${}_{T_{max}}\approx$−0.29), daily mean LST and median LST (R${}_{\overline{T},T_{med}}\approx$−0.24), and daily minimum LST (R${}_{T_{min}}\approx$ −0.19). The strongest anticorrelation is obtained for the LST diurnal range both in morning and afternoon passes. This agrees with previous literature [56,57,58], in which a relation between SM and near-surface air temperature diurnal range was reported. However, this air temperature is usually measured at 1.5 m above the ground level, not in the top 5 cm of the soil. Results indicate that maximum LST and LST diurnal range provide complementary information of the SM–LST relationship. Since results obtained at the SMOS morning and afternoon are very similar using both the direct variables and their anomalies, hereafter only results at morning passes are presented.

Figure 6 (left) shows the correlation of SMOS SM at the morning passes and MODIS LST at different platforms and acquisition times during the study period. Correlations are less variable for satellite than for in situ observations, likely due to their different spatial scales; satellite observations represent the four 25-km pixels covering the REMEDHUS area whereas ground-based measurements are made at point scale. The correlation values obtained (in median) are: Terra day (R ≈ −0.68), Terra night (R ≈ −0.54), Aqua day (R ≈ −0.70), and Aqua night (R ≈ −0.44). The strongest anticorrrelations are obtained for MODIS LST day times, both for Terra and Aqua satellites. We hypothesize that MODIS LST from day overpasses captures better the spatial pattern induced by topography that may not be captured at night time. Note that the anticorrelation for Aqua day is slightly stronger than for Terra day. This can be explained by Aqua day passes being closer to $T_{max}$ over the REMEDHUS region than Terra day passes. These results agree with those obtained in Figure 5 (left). Further research using a geostationary satellite as MSG SEVIRI, which has a temporal resolution of 15 min, would be needed to analyze the relationship of SM with LST diurnal range.

When analyzing the anomaly time-series, Figure 6 (right) displays correlations that are weaker than when using the direct variables in absolute value, as expected, but significant in all cases: Terra day (R ≈ −0.47 in median), Terra night (R ≈ −0.13), Aqua day (R ≈ −0.46), and Aqua night (R ≈ +0.17). Note that Terra/Aqua at night overpasses produce very weak correlations, suggesting that only Terra/Aqua at day times maintain their relationship with the SM under anomalous conditions. Similar behavior is obtained at the afternoon passes (not shown). According to these results, only daily maximum LST (for the in situ analysis) and MODIS LST Aqua day (for the satellite analysis) will be shown hereafter.

Two-dimensional frequency density diagrams of in situ SM acquired at the SMOS morning passes and $T_{max}$, and SMOS SM at the morning passes and MODIS LST Aqua day during the study period are displayed in Figure 7. They inform of the frequency of occurrence of all SM–LST combinations. The normalized SM–LST space visually shows the relationship between the two variables through its distributions. A triangular shape emerges from the SM–LST relationships, with sharper edges when using in situ data (left) than when using satellite data (right). This is likely to the fact that MODIS LST Aqua day is not acquired at the maximum daily LST time, when the strongest SM–LST correlation is obtained. It is important to note that MODIS LST measurements are obtained at the skin level whereas the ground-based LST observations are acquired at 0–5 cm below the land surface. Also, daily LST fluctuations are less pronounced with depth [59]. These differences could explain the higher MODIS LST Aqua day than $T_{max}$ during very dry conditions (0–0.1 m${}^3$/m${}^3$). Similar triangles are obtained at the afternoon passes and/or using Terra day (not shown). The triangle or trapezoid was also obtained when plotting LST–NDVI space, which is used to indirectly estimate SM from TIR observations [25,26,27,28], or to disagregate SM from the synergy of microwaves and TIR [32,34].

4.3. SM–LST Relationship at the Seasonal Scale

The correlation of in situ SM at the SMOS morning passes and $T_{max}$ from the REMEDHUS network for each season is shown in Figure 8 (left). Note that the strongest anticorrelation is obtained in autumn (September–October–November, SON) with R ≈ −0.68 in median, followed by spring (March–April–May, MAM) with R ≈ −0.51, and summer (June–July–August, JJA) with R ≈ −0.48. Hence, the strongest SM–LST relationship is observed in transitional seasons (i.e., autumn and spring) because the SM is sufficiently dynamic to affect air temperature, and also LST. This result agrees with those obtained in [60], where the SM–air temperature coupling was shown to be higher in transitional zones between wet and dry climates. However, a low correlation is obtained in winter (December–January–February, DJF) with R ≈ + 0.30. This suggests that SM and LST are not coupled during winter. This could also explain the low performance of disaggregated SM products based on soil evaporation efficiency models in winter [33].

Figure 8 (right) shows the correlation of ground-based SM anomaly at the morning passes and $T_{max}$ anomaly for each season. A low correlation is also observed in winter (R ≈ +0.27 in median) and the strongest anticorrelation is obtained in summer (R $\approx -$0.49), which is the driest season, followed by spring (R $\approx -$0.43) and autumn (R $\approx -$0.31). It can be related to the fact that, under very dry conditions, an increase of SM from an individual precipitation event has an important impact on LST, due to its low thermal inertia. Also, since ET is regulated by the SM content under water-limited conditions in summer, the SM increase produces an ET increase and there is enough energy to maintain the high rate of ET until the SM content decreases. Obviously, this effect is less important in wet seasons, specially in winter, due to the lower available energy involved in the ET process and the higher thermal inertia. These results suggest that the annual cycle contains most of the SM–LST relationship information of transitional seasons, particularly in autumn. By contrast, the anomaly contains most of the SM–LST relationship information during summer. Similar behavior is obtained with emphin situ SM at the SMOS afternoon passes (not shown). Further research is needed to confirm these results in other regions with different climates.

The correlation of SMOS SM at the morning passes and MODIS LST Aqua day for each season is displayed in Figure 9 (left). In this case, all correlations are negative, including winter. The anticorrelations sorted from strong to weak are obtained in autumn (R $\approx -$0.69), spring (R $\approx -$0.66), summer (R $\approx -$0.39) and winter (R $\approx -$0.34). Thus, the transitional seasons also present the strongest anticorrelation and other seasons the weakest. We hypothesize that differences between in situ and satellite results during winter could be due to the more frequent presence of clouds than during other seasons, which produces a lower coverage and a higher LST uncertainty that could be up to 10 °C. Another reason may be a side effect of the different spatial scales represented by the remotely sensed data and ground-based data.

Figure 9 (right) shows the correlation of the detrended spaceborne time-series. The sorted anticorrelations are obtained in spring (R $\approx -$0.62), autumn (R $\approx -$0.61), summer (R $\approx -$0.39) and winter (R $\approx -$0.26). The lowest SM–LST relationship is obtained during winter, in agreement with results of Figure 8 (right). Similar behavior is obtained with the afternoon passes and/or using Terra day (not shown).

Normalized occurrence frequency density diagrams of in situ SM at the SMOS morning passes and $T_{max}$ (top row), and SMOS SM at the morning passes and MODIS LST Aqua day (bottom row) are displayed in Figure 10 for each season. Using in situ data, SM exhibits a similar dynamic range in winter than in other seasons, but there is a very low dynamic range in $T_{max}$ since there is a minimum thermal gradient of the soil in this area. This leads to the appearance of a line instead of a triangle in the SM–LST space. By contrast, satellite data have a higher dynamic range of LST Aqua day, showing a small triangle in the SM–LST space. The triangle shape is clearly visible in transition seasons (autumn and spring) both using in situ and satellite data due to the high SM and LST dynamic ranges. In summer, it is better identified using in situ than satellite data. Differences between ground-based and satellite data could be due to the different soil layer captured by in situ LST measurements (0–5 cm) and sensed by MODIS (skin), as well as to the higher thermal inertia of ground-based observations, particularly in wet soils during winter. These results are consistent with those obtained in Figure 8 and Figure 9. Similar behavior is obtained at the afternoon passes and/or using Terra day (not shown).

4.4. SM and LST Time-Series: Coupling/Decoupling and Critical SM Estimation

Figure 11 (top) shows time-series of in situ SM at the SMOS morning passes and $T_{max}$ (blue and red solid lines, respectively), and SMOS SM at the morning passes and MODIS LST Aqua day (cyan and magenta dashed lines, respectively). Daily mean precipitation from the REMEDHUS network is also shown (bottom). Winter seasons defined as December, January and February (DFJ) are shaded in light green color. Although there is a different temporal resolution for in situ and satellite data (daily versus 3-day averages, respectively), both SM estimations are consistent with precipitation and capture wetting and drying events. The SMOS SM has a dry bias with respect to in situ SM. This dry bias has been previously reported in studies [40,48,61] and remains almost constant along time. However, the MODIS LST Aqua day is very similar to $T_{max}$ in winter and is higher in other seasons, in agreement with results shown in Figure 7 and Figure 10. Note that the SM–LST relationship presents two different behaviors in both datasets. On one hand, there is SM–LST coupling (when SM decreases LST increases and vice versa), approximately corresponding to the spring, summer and autumn seasons, in which a strong anticorrelation is obtained (R $\approx -$0.5 to −0.7/R $\approx -$0.4 to −0.7 from in situ/satellite data in Figure 8 and Figure 9, respectively). This is related to a water-limited ET regime. On the other hand, there is SM–LST decoupling (when SM and LST do not show a clear relation between them) approximately corresponding to the winter seasons, in which a weak correlation is obtained (R ≈ +0.30/R ≈ −0.3). This is related to an energy-limited ET regime. The discrepancy in the sign of the correlation during winter could be due to the different spatial scales represented by the remotely sensed data and ground-based data. In addition, SM–LST decoupling periods vary from one year to another and include, in most cases, around one month after and before the defined winter. Similar behavior is obtained at the afternoon passes and/or using Terra day (not shown).

Minimum, maximum, mean and standard deviation (std) of SM and LST values obtained at the SM–LST coupling/decoupling transitions are summarized in

Table 1. Statistics obtained from SM–LST coupling/decoupling transitions along the study period.

	In Situ Data		Satellite Data
	SM	${\mathit{T}}_{\mathit{max}}$	SMOS SM	MODIS LST Aqua Day
	(m${}^{\mathbf{3}}$/m${}^{\mathbf{3}}$)	( °C)	(m${}^{\mathbf{3}}$/m${}^{\mathbf{3}}$)	( °C)
minimum	0.10	12.4	0.08	8.4
maximum	0.16	22.6	0.22	32.7
mean	0.12	16.4	0.12	15.2
std	0.06	5.9	0.05	4.4

. Note that the same mean SM is estimated in both datasets (∼0.12 m${}^3$/m${}^3$), which corresponds to the critical SM over the study area. Taking into account its variability (in situ std SM∼0.06 m${}^3$/m${}^3$; satellite std SM∼0.05 m${}^3$/m${}^3$), the estimated critical SM remains within reasonable boundaries (0.10–0.16 m${}^3$/m${}^3$ vs. 0.08–0.22 m${}^3$/m${}^3$). In addition, the estimated critical SM lies between WP and FC of the REMEDHUS network (WP = 0.08 m${}^3$/m${}^3$; FC = 0.16 m${}^3$/m${}^3$) and is between the 50% and the 80% of FC (∼75%), in agreement with [22]. The proposed methodology could be an interesting technique to measure critical SM using remote sensing. In the SM–LST transitions, the mean LST from in situ (∼16.4 °C) has a good agreement with the mean LST from satellite data (∼15.2 °C). It is observed that spaceborne LST measurements have a higher upper limit (∼32.7 °C) than ground-based measurements (∼22.6 °C). This can be explained by the fact that MODIS LST Aqua day observations have also a higher dynamic range than in situ observations, in agreement with results shown in Figure 7 and Figure 10. Considering the LST variability in both datasets (std LST ≈ 4.4 °C to 5.9 °C) the minimum LST value is similar (12.4 °C vs. 8.4 °C). Other regions with different climates should be analyzed in further research studies.

4.5. Critical SM and ET Regimes during 2014

Figure 12 shows the monthly actual evapotranspiration ($E{T}_a$) from HidroMORE, the potential evapotranspiration ($E{T}_0$) and the in situ SM from the REMEDHUS network along year 2014. It can be observed that $E{T}_0$ and $E{T}_a$ are very similar in some months of the year (November–December–January–February, NDJF) and the SM is larger than the estimated critical SM in Table 1. This result reveals an energy-limited ET regime. In this regime, the temperature regulates the ET, the $E{T}_0$ is very low since there is not enough energy and the $E{T}_a$ is almost able to reach its maximum possible value (i.e., $E{T}_0$), but there is enough SM. Thus, the ET process is independent of the SM content and has an important dependence on the sensible heat flux, which impacts on the air temperature and the LST. By contrast, $E{T}_0$ and $E{T}_a$ have a different behavior in the other months (March–April–May–June–July–August–September–October, MAMJJASO) and the SM is lower than the critical SM, indicating a water-limited ET regime. In this case, the $E{T}_0$ increases from March to July and decreases from July to September, according to the temperature variation in this area. The $E{T}_a$ is always lower than the $E{T}_0$ because the SM controls the ET process. Note that $E{T}_a$ minimum is in August, corresponding to the minimum SM, and the driest month in semi-arid climates. The energy-limited and the water-limited periods found for 2014 approximately agree with the periods between transitions in Figure 11 (top).

Table 2. Correlation (R) of $E{T}_a$ from HidroMORE and in situ SM from the REMEDHUS network, and number of observations (N) under energy-limited and water-limited conditions along 2014, considering the delayed SM with respect to the $E{T}_a$ day (delay = 0, 1, 2, 3, 4 or 5 days).

Delay in SM respect to ${\mathit{ET}}_{\mathit{a}}$	Energy-Limited		Water-Limited
	(NDJF)		(MAMJJASO)
	R	N (days)	R	N (days)
Same day	0.05 *	120	0.54	245
1 day after	−0.07 *	119	0.59	244
2 days after	−0.08 *	118	0.61	243
3 days after	−0.17 *	117	0.62	242
4 days after	−0.18	116	0.62	241
5 days after	−0.20	115	0.62	240

* Non-significant R values ($p_{value}>0.05$) are marked with an asterisk.

summarizes the correlation between daily in situ SM from the REMEDHUS network and $E{T}_a$ from HidroMORE, and the number of observations (N) for the energy-limited (NDJF) and water-limited (MAMJJASO) periods along 2014. As the persistence of SM is an important component in land–atmosphere interactions [21], the correlation has been computed considering the delayed SM with respect to the $E{T}_a$ day (with a delay between 0 to 5 days, 0 days corresponds to the same day). In the energy-limited period, most correlation values are very weak (|R| < 0.2) and others are non-significant. This result supports the argument that there is not a clear relationship between the SM and $E{T}_a$ under energy-limited conditions due to SM–LST decoupling. Instead, correlations obtained in the water-limited period (R ≈ 0.54 to 0.62) reveal an important relation of SM and $E{T}_a$ due to SM–LST coupling. Since the SM of a particular day is related to the amount of previously lost water, the correlation increases considering the SM of one or two days after the $E{T}_a$ day. This result indicates that there is a temporal lag between the ET behavior and its impact on the SM.

5. Conclusions

Nowadays, most of radiometric soil moisture (SM) retrieval and disaggregation algorithms employ the instantaneous land surface temperature (LST) in their approaches. Nevertheless, alternative LST-derived parameters, such as the daily maximum LST, may provide a better representation of the SM–LST relationship than the instantaneous LST. In order to analyze this hypothesis, we evaluated the SM and LST dynamics at daily and seasonal scales using four years of data (2011–2014) from both in situ and satellite observations acquired over the REMEDHUS network, located at the central part of the river Duero basin in Spain. The main objective was to better understand the fundamental link between the instantaneous SM (at the SMOS passes over this region) and several LST-derived parameters, and assess this relationship through the SM–LST coupling/decoupling, the evapotranspiration process and the thermal inertia. Furthermore, the SM–LST interaction was studied with the direct variables in absolute values and with their anomalies, subtracting the mean seasonal cycle of each variable (or climatology), to separate out seasonality.

Considering the entire study period, SM and LST are, as expected, anticorrelated. At the daily scale, results from in situ data show that instantaneous SM exhibits stronger anticorrelation to daily maximum LST (R $\approx -$0.7) than to the other LST-derived parameters (i.e., instantaneous LST, daily mean LST, daily median LST, daily minimum LST and LST diurnal range). When using satellite data, stronger anticorrelation is obtained between SMOS SM and MODIS LST Aqua (or Terra) day (R $\approx -$0.7) than night (R $\approx -$0.4 to −0.5). This is consistent with the in situ analysis, since the time of MODIS LST Aqua (or Terra) day is closer to the time of daily maximum LST, which is also the time of maximum potential evapotranspiration (i.e., when there is a higher atmospheric demand of water). For the anomaly time-series, results from in situ data show that LST diurnal range presents the strongest anticorrelation with SM (R $\approx -$0.3), followed by daily maximum LST (R $\approx -$0.3). This indicates that maximum LST and LST diurnal range provide complementary information of the SM-LST relationship. When using satellite data, the strongest anticorrelation is obtained with MODIS LST Aqua (or Terra) day (R $\approx -$0.5). Further research using a geostationary satellite such as MSG SEVIRI with a temporal resolution of 15 min. would be needed to further evaluate the relationship of SM with LST diurnal range. In addition, the SM–LST space of the normalized occurrence frequency density diagrams display a triangular shape in both ground-based and spaceborne observations.

The seasonal dependency of the correlation was also studied. Results from in situ data show a strong anticorrelation between instantaneous SM and daily maximum LST in spring, summer and autumn seasons (R $\approx -$0.5 to −0.7), revealing a SM–LST coupling. By contrast, a very low correlation is obtained in winter (R ≈ +0.3), which indicates a clear SM–LST decoupling. When using satellite data, results confirm a SM–LST coupling in spring, summer and autumn (R $\approx -$0.4 to −0.7), and a SM–LST decoupling in winter (R $\approx -$0.3). For the anomaly time-series, the highest correlation is obtained in summer for in situ data (R $\approx -$0.5). This can be to the related to the fact that an increase of SM from an individual precipitation event has an important impact on LST in very dry conditions due to its low thermal inertia. Differences in the correlation during winter with the direct variables and during summer with the anomalies could be due to the different spatial scales represented by the remotely sensed and ground-based data. Similar conclusions could be extracted from the normalized occurrence frequency density diagrams displayed per seasons, where the triangular shape is better defined in transitional seasons (autumn and spring) due to the higher dynamic range of SM and LST.

Collocated SM and LST time-series display SM–LST coupling periods, approximately corresponding to spring, summer and autumn, and SM–LST decoupling periods, approximately corresponding to winter. These two behaviors can be identified in both in situ and satellite measurements. A method to estimate the critical SM from the crossing points of these SM–LST coupling/decoupling transitions was proposed. A critical SM of ∼0.12 m${}^3$/m${}^3$ is estimated using both in situ and spaceborne observations. This value is consistent with field measurements of soil field capacity and wilting point.

Since the SM–LST coupling/decoupling periods are related to water/energy-limited regimes, respectively, the temporal evolution of the evapotranspiration is further analyzed for a specific year (2014). The actual evapotranspiration is simulated using the regional hydrological model HidroMORE. The correlation between the daily SM and the modeled actual evapotranspiration $E{T}_a$ confirm that SM and $E{T}_a$ only interact under water-limited conditions (with R ≈ 0.6) with a temporal lag between them of 2–3 days. A non-significant correlation between SM and $E{T}_a$ is observed under energy-limited conditions.

This study contributes to furthering present knowledge of SM and LST interactions in land surface, hydrological and climate models. Results presented could be applicable to microwave/optical synergistic SM retrieval algorithms and, particularly, for the improvement of applications based on the SM–LST link, related to evapotranspiration. Further research is needed to study the possible dependency of the different Earth’s land covers/climate types in the SM–LST relationship, and also to analyze this relationship using geostationary LST data.