Soil Water Content Diachronic Mapping: An FFT Frequency Analysis of a Temperature–Vegetation Index

Abstract

Among the indirect estimation approaches of soil water content in the upper layer of the soil, the “triangle method” is one of the most common that relies on the simple relationship between the optical and thermal features sensed via Earth Observation. These features are controlled by water content at the surface and within the root zone but also by meteorological forcing including air temperature and humidity, as well as solar radiation. Night- and day-time MODIS composites of land-surface temperature (LST) allowed applying a version of the triangle method that takes into account the temporal admittance of the soil. In this study, it has been applied to a long time-series of pair images to analyze the seasonal influence of the meteorological forcing on a triangle method index (or temperature–vegetation index, TVX), as well as to discuss extra challenges of the diachronic approach including seasonality effects and the variability of environmental forcing. The Imera Meridionale basin (Sicily, Italy) has been chosen to analyze the method over a time-series of 12 years. The analysis reveals that, under these specific environmental and climatic conditions (strong seasonality and rainfall out of phase with vegetation growth), Normalized Difference Vegetation Index (NDVI) and LST pairs move circularly in time within the optical vs. thermal feature space. Concordantly, the boundaries of the triangle move during the seasons. Results showed a strong correlation between TVX and rainfall normalized amplitudes of the power spectra (r² ~0.8) over the range of frequencies of the main harmonics.

1. Introduction

Soil water content plays an important role in the vegetated natural and agricultural ecosystems. The accurate assessment of the spatiotemporal traits of this fundamental variable is crucial for agricultural practices (e.g., irrigation) and, in general, for the detection of stress conditions over vegetated lands.

Among the indirect estimations of soil water content, θ, in the upper soil layer, the “triangle method” is a methodology based solely on the joint analysis of the optical and the thermal features (i.e., vegetation indices, VIs, and Land Surface Temperature, LST, respectively) sensed via Earth Observation, EO [1].

The traditional single acquisition formulation of the triangle method is based on the hypothesis that meteorological conditions are fairly constant across the area. Differently to other methods (i.e., the thermal inertia method [2,3]), it requires the use of a single pair of VI/LST images over a scene where a wide range of θ, covering both dry and wet extremes under a full range of vegetation coverage. Maltese et al. [4] highlighted the limits of this method if applied over small areas with high spatial resolution images, thus, the authors proposed a diachronic based empirical method to determine the outer edges. Besides, Maltese et al. [4] tested alternative reference temperatures, achieving a high correlation with the observed θ (r² up to ~0.9). Similarly, other researchers proposed to use alternative variables as the reference temperature. Wang et al. [5] proposed to adopt the day-night LST difference, using daily MODIS (MODerate-resolution Imaging Spectroradiometer) images to estimate the evaporative fraction. Minacapilli et al. [6] explored a diachronic approach of the triangle method (named the Time-Domain Triangle Method) to estimate the actual evapotranspiration using MODIS images. Stisen et al. [7] suggested employing the surface temperature increases between 8.00 local time and midday, using (Spinning Enhanced Visible and InfraRed Imager) SEVIRI thermal images to estimate evapotranspiration with a triangle method–thermal inertia combined approach over a four-month-long period.

However, the multi-temporal approach carries some extra challenges in the correct parameterization of the two bounds, mainly related to the effect of seasonality on the involved variables. Indeed, since the periodic fluctuations of the main atmosphere meteorological forcing (i.e., air temperature, humidity, and solar radiation) directly influence the thermal regime of the land surface, it is necessary to isolate the contribution of θ fluctuations on LST from that of other sources.

Similarly, the phenological development of the vegetation (of which the remotely sensed VIs, as the Normalized Difference Vegetation Index—NDVI, are a proxy) is only partially controlled by the vegetation and soil water availability, but it may be also controlled by other meteorological factors.

Optical vegetation indices, such as NDVI, may be a conservative indicator of soil water content. As vegetation remains green after initial water stress [8,9], NDVI is a slow-responding indicator because it shows some inertia to change its values under water stress forcing, whereas, LST increases forthwith water stress [10,11,12]; thus, defining a temporal relationship between NDVI, LST, and θ could be not straightforward. This issue is undoubtedly important in semi-arid regions where θ is affected by the alternation of rainy winter and warm and dry summer seasons, which results in a cycle of θ asynchronous with solar radiation and dependent variables.

The primary objective of this study is to analyze a temperature-vegetation index-based method to retrieve soil water content temporal dynamics and to assess whether the diachronic approach of the triangle method is suitable to mimic the behavior of soil water content in the first centimeters of soil and for long periods. Secondarily, the connections among remotely sensed LST and NDVI signals and meteorological forcing (e.g., Higuchi et al. [13]) are also analyzed, in relation to the modeled triangle method-based index and soil water content dynamics.

With this aim in mind, EO data from the MODIS satellite were collected over a hydrological basin located in southern Sicily (Italy), whereas θ values were obtained from the outputs of the LISFLOOD rainfall-runoff model [14]. As in situ measurements of θ were not available, the hydrological modeled θ were considered as a proxy of the wetting status of the soil for comparison.

2. Materials and Methods

2.1. Study Area

The method has been tested on the Imera Meridionale catchment located in southern Sicily, Italy (Figure 1). The catchment extends from the Madonie mountains (North) to the Mediterranean Sea (South) for about 2000 km². The land use includes water bodies, wooded areas, homogeneous cropland (mainly cereals) and sparse cropland (vineyards, orchards, and olive groves).

The area is characterized by a typical Mediterranean climate, with moderate rainfall during autumn and winter and very high air temperature and scarce precipitation during summer. In a previous study [15], the rainfall regime over the catchment was characterized, which exhibits a mean annual precipitation that ranges from about 450 mm in coastal areas up to 800 mm in mountain areas with an average annual value of 520 mm. Almost no precipitation occurs in July–August, while maximum total precipitation occurs in December–February (≈140–170 mm). Occasionally, severe rainfall events occur in spring (from March to May), with alternation of dry periods and rainy days. Cammalleri et al. [16] and Scordo et al. [17], instead, characterized the climate of the basin by means of a Péguy climograph [18]. In particular, the coastal zone is classified as arid climate from April–May to September–October (≈5–7 months), while the central part experiences arid climate for a shorter period (May to September, 5 months), higher altitudes have arid climate from May–June to August–September (≈3–5 months).

The catchment is gauged with a hydrometric station operated by Servizio Osservatorio delle Acque (Sicilian regional government, Regione Siciliana) [19], namely the Drasi hydrometric station (red dot in right panel of Figure 1) and 5 remote controlled rain gauge stations operated by SIAS (Servizio Informativo Agrometeorologico Siciliano, Regione Siciliana) [20]: Licata, Riesi, Delia, Caltanissetta Enna, and Petralia Sottana (green dots in the right panel of Figure 1).

2.2. Data

Total daily rainfall, P, records were acquired for the whole period by SIAS, together with daily average air temperature and humidity and total daily incoming solar radiation. Daily total discharges were acquired by the hydrometric gauge station at Drasi (37.12.30N 13.59.34E; 130 m a.s.l.) [19].

Eight- and sixteen-day composite LST and NDVI MODIS products (MOD11A2 and MOD13Q1) were processed. The MOD11A2 product is composed of the daily (both day and night) 1-kilometer LST product (MOD11A1) and stored on a 1-kilometer grid as the average values of clear-sky LSTs during an 8-day period. The MOD13Q1 product is computed from atmospherically corrected bi-directional surface reflectance and is delivered every 16 days at a 250-meter spatial resolution [21].

Night-time and day-time LST composite images (LST_D and LST_N) were collected from the 1st of January 2001 to the 31st of December 2012.

A reference dataset for the soil water content dynamic was derived from the LISFLOOD rainfall-runoff model [14]; this model, developed by the Joint Research Centre (JRC) of the European Commission, is designed to reproduce the hydrology of large and trans-national European river catchments at 5 × 5 km of spatial resolution and it currently runs operationally within the European Flood Awareness System (EFAS) [22]. The model has been extensively calibrated and validated over Europe in terms of river discharge [23,24]. LISFLOOD also provides valuable information on soil water content status in two layers, i.e., topsoil and subsoil [25]. Within this study, the LISFLOOD soil water content index, U_LIS, of the topsoil (root zone) was considered; please note that the root zone refers to a variable depth depending on the land cover.

To assess indirectly the LISFLOOD performance in the investigated area, the daily-simulated river discharge was compared with in situ discharge data collected at the Drasi gauge station. In addition, to assess the model performance in catching the temporal dynamics of discharge, simulated and observed discharge time-series were compared in the power spectra domain by applying a fast Fourier-transformation (FFT).

As the power spectra of simulated and observed discharges show a strong determination coefficient (r² ≈ 0.62) (Figure 2, top-right corner), also the soil water content modeled by LISFLOOD, U_LIS, is assumed to match the actual soil water content at least in terms of frequencies of peaks occurrence. However, by restricting the comparison to sub-yearly harmonics only (by removing the top-left point) r² decreases to ≈0.45, showing that LISFLOOD performances increase as the temporal scale increases.

2.3. Methods

2.3.1. Thermal/Optical Features Space Characterization

The soil water content indirect estimation was performed by applying the so-called “triangle method.” The soil water content of a generic pixel is determined by looking at two variables: One representing a vegetation index (VI, e.g., NDVI or the vegetation fractional cover) and the land surface temperature. The model is an empirical implementation of the triangular shape of the VI–LST scatterplot [11]. The method requires the detection of dry (warm) and wet (cold) edges of the triangle; then, for each pixel, a wetness index is derived on the base of its relative position, with respect to the above-mentioned edges, within the VI–LST scatterplot. Several versions of the triangle method were proposed in the literature; in this research, accordingly to Maltese et al., [4] LST_D and LST_N images allowed to apply the thermal admittance version of the triangle method (i.e., using ΔLST = LST_D − LST_N). Indeed, within this implementation, the scatterplot is made up of the difference between the day- and night-time land surface temperatures here assumed as reference temperature, vs. NDVI. The implementation is based on Maltese et al. [2] whose θ estimates closely agreed with in situ measures by accounting for the soil admittance. The thermal admittance version of the triangle method was tested on high spatial resolution (~2 m) airborne images acquired on a small Sicilian vineyard (tens of hectares) in the period from mid-June to the beginning of September. The irrigation farm strategy was a controlled water deficit technique, by means of a drip irrigation system. Few rainfall events, characterized by low values, occurred within the study period. The climate was similar to that of the Imera Meridionale catchment at comparable altitudes.

According to Carlson [26], even at NDVI maximum values, LST varies, resulting from variation in soil admittance with changing θ in the soil depth, thus conferring a trapezoidal shape (instead of a perfect triangular shape) to the scatterplot distribution.

On a single temporal acquisition, it is expected that a wider base (hereinafter base) of the trapezoid highlights a high thermal variability of bare soils due to surface θ evaporation processes and soil admittance properties [2]; a narrow base, instead, corresponds to densely vegetated areas and highlights the low-temperature variability due to the plants thermoregulation through transpiration and to thermal properties of the underneath soil [3].

Under the assumption that the time-series for each interval of vegetation coverage includes all possible soil water content conditions (from residual to saturation), both dry (“dry or warm” edge) and wet (“wet or cold” edge) conditions can be detected. The TVX is then defined as a function of the relative position to the edges of a generic pixel within the scatterplot, based on the hypothesis of linear variation of isopleths from the wet to the dry edge. The relationship between LST and θ over a vegetated surface can be explained through the canopy resistance and its role in the surface energy balance. The canopy resistance depends on the Leaf Area Index, LAI, and stomata resistance; this latter is controlled by several other meteorological factors, as exemplified by the empirical approach proposed by Jarvis [27]. In the empirical triangle methods, it is convenient to use NDVI instead of LAI as it produces a nearly linear relation between temperature and vegetation index [28]. The relationship between LAI and NDVI is non-linear (moreover, it varies with vegetation species) and expresses the Beers-type light attenuation through the canopy.

The triangle method edges and, subsequently, the TVX is adapted from [2]:

$$\{\begin{array}{c}TVX=1-\frac{\alpha }{\beta }\\ \alpha =\arctan \left|\frac{m_D-m_W}{1+m_D{m}_W}\right|\\ \beta =\arctan \left|\frac{m_D-m}{1+m_{D}m}\right|\end{array},$$

(1)

where, $m_W$ and $m_D$ refer to the slope of the dry and wet edges in the ΔLST–NDVI space and m is the slope of the equation of the line passing through a generic pixel of coordinates (ΔLST, NDVI) and the vertex of the triangle that can be obtained by intersecting the dry and the wet edges.

The methodology is made of different analyses (grey boxes in Figure 3); specific outputs of each analysis are reported in bulleted lists (white boxes).

2.3.2. Time-Domain Analysis

Seasonality effects: By arranging NDVI and ΔLST images in a matrix having the day of year (DOY) on lines and the year sequence on columns, the seasonality effects can be highlighted by identifying cold and warm periods to understand also the phase shift between ΔLST and NDVI; this visualization also allows us to identify changes occurring year by year. Dry and wet edges are determined on the basis of NDVI–ΔLST values collected in summer and winter seasons, respectively.

The clusters positioning in the ΔLST–NDVI feature space is compared to the Péguy climograph to analyze similarities among climate and land surface EO characteristic behavior during the year.

2.3.3. Stress Factors Analysis

The distribution of stomatal conductance in the NDVI–ΔLST space is evaluated to highlight observable stress due to other environmental factors only (i.e., the θ stress factor is not taken into account). To this aim, according to Noilhan and Planton [28], the optimal canopy resistance, r_c, for a vegetated surface is computed by the ratio of the stomata resistance and the effective LAI (big leaf hypothesis). LAI is determined as a function of the NDVI, whereas the actual canopy resistance is derived on the basis of some set of environmental stress functions [29]; these latter are assumed to be non-interacting functions (e.g., incident solar radiation, air humidity and temperature, and carbon dioxide concentration dependent), which allow characterizing the ambient conditions that influence the canopy resistance, r_c, as well as θ. In this case study, the stress effect of carbon dioxide concentration (CO₂) on canopy resistance is neglected under the hypothesis that the atmospheric CO₂ is fairly constant [30].

Stomata Conductance Spatial Distribution: The influence of the environmental stress factors is evaluated by computing the canopy conductance related to the solar radiation, C_Rs, the air temperature, C_Tair, and the air humidity, C_Uair. In this study, canopy resistances have been evaluated in the spatial distributed form within the NDVI vs. ΔLST space according to the empirical approach of Jarvis [27]. The limiting value of the solar radiation is calculated according to Noilhan and Planton [28]. Prino et al. [30] calibrated the function C_Tair taking into account the reduced plant activity for temperatures different from the optimum air temperature. Sepúlveda and Kliewer [31] measured the stomata response of some grapevine cultivars to high temperature and they found the range of variability of the optimum air temperature. The reduction factor for atmospheric water vapor deficit is given by the atmospheric specific humidity and the saturated specific humidity at a given air temperature as proposed by Stewart and Gay [32].

In this framework, as a diachronic approach is considered additionally to the soil water content dynamics, variations of the above-mentioned environment stress factors could play a role. Need to specify that, as these equations were calibrated for specific cultivars, whereas within this research they were applied over the whole catchment, only spatial patterns are considered for discussion. An interpretation of the temporal-spatial distributions of C_Rs, C_Tair, and C_Uair, is provided in the “results and discussion” section. The distribution of the conductance reduction factor, C, in the NDVI–ΔLST space is obtained as in the following: (i) for each MODIS acquisition time a C value is computed (using meteorological data); (ii) the NDVI–ΔLST space is divided into sectors (50 intervals for both ΔLST and NDVI); (iii) for each sector and DOY only sectors accounting for at least the 20% of ΔLST–NDVI pairs were selected; (iv) the average of C is associated to the sectors by assigning the pair corresponding to the common acquisition time. The time average spatial distribution is evaluated for the conductance reduction factor due to solar radiation, air temperature and humidity (C_RS, C_Tair, and C_Uair, respectively). Other factors, such as soil nutrients availability and moderated transpiration (e.g., [33]) are not taken into account in this analysis.

2.3.4. Frequency Domain Analysis

Once clarified the influence of environmental factors on TVX, the index is compared to the LISFLOOD θ, U_LIS, in terms of magnitude, temporal behavior, and phase shift for the whole study period. To analyze TVX response to rapid changes in θ due to rainfall events, TVX and U_LIS were rescaled and compared within selected years where frequent rainfall events occurred (2005–2006). All correlations are classified according to Evans [34].

Soil water content indexes (i.e., TVX and U_LIS) can differ in magnitudes as outcomes of different methods and dependently on the depth of the soil they refer to. However, we expect a similar frequency of occurrence to the main meteorological forcing, the rainfall; thus, variables were analyzed in the frequency domain via a FFT. The FFT allows one to compute the discrete Fourier-transform (DFT) by decomposing a sequence of values into components of different frequencies [35,36]. A Shift zero-frequency component was added to the center of the spectrum to describe each time-series by a series of harmonics. Increasing harmonics number, from 0 to 2, describes the inter-annual, annual, and semi-annual components. The sum of the main harmonics describes well the whole seasonal variability [37].

Additionally, to the forward FFT transformation, the phases of the more relevant frequencies were investigated. In particular, once the FFT is applied to each time-series, its more relevant frequencies are selected and bounded by means of band pass filters; signal harmonics are then obtained by applying the inverse FFT, iFFT. This allowed characterizing each more relevant harmonic in terms of peak frequency, periodicity, and peak time occurrence. From an operational point of view, FFT and iFFT were applied through Matlab functions [38]: Fast Fourier-transform (FFT), shift zero-frequency component to the center of the spectrum (FFTshift), and inverse fast Fourier-transform (iFFT), inverse zero-frequency shift (iFFTshift), respectively, for the forward and backward transformations.

2.3.5. Comparisons

U_LIS and TXV time-series were qualitatively compared over the whole period. Over a restricted temporal period U_LIS and TXV were compared with the rainfall time-series. Finally, TVX and U_LIS were compared with the measured total daily rainfall, P.

3. Results and Discussion

3.1. Thermal/Optical Features Space Characterization

The distribution of ΔLST–NDVI pairs during the whole time-series (Figure 4a) highlights two main clusters (green pixels in the scatterplot) characterized by low NDVI and high ΔLST and vice versa by low ΔLST and high NDVI.

A lower number of pairs link these clusters conferring the trapezoidal shape to the whole pattern. The trapezoid exhibits a “hole” in the middle, at least at low spatial and temporal resolutions, likely due to the periodic behavior of meteorological forcing on surface temperature and vegetation indexes. Two x-profiles passing through the maxima ((~8.2, ~0.55) and (~21.5, ~0.25)) of the main clusters highlight normalized density variability (Figure 4b) close to the wet edge lower than that close to the dry edge; this results in a higher sensitivity of the position of the points determining the wet edge especially for low NDVI, confirming the retrieval of Maltese et al. [4].

Dry and wet edges (red and blue lines, respectively) were calculated by extracting minima and maxima percentiles (0.02 and 0.98, respectively). The periodic behavior of meteorological forcing should be also responsible for a less noticeable difference between the wider and the narrow bases (corresponding to low and high NDVI in Figure 4a).

3.2. Time-Domain Analysis

By plotting ΔLST–NDVI pairs in time (Figure 5), starting at the beginning of the year in the winter season, it is noticeable that pairs start to move clockwise with the highest NDVI and the lowest ΔLST.

Figure 5 shows a whole sequence of ΔLST–NDVI pairs (time-domain mosaics) spanned over a year (from April 2002 to March 2003). In April (upper left panel in Figure 5) ΔLST is relatively high ~20 °C with NDVI~0.4. In July–August (upper right panel in Figure 5), the cloud moves to the highest ΔLST and lowest NDVI values. In November (lower left panel), the NDVI slightly increases, whereas, ΔLST considerably decreases (~12 °C). The lowest ΔLST is reached in December–January (~5 °C, lower right panel in Figure 5) with high NDVI values (0.6). NDVI remains fairly constant until March (NDVI~0.65) whereas ΔLST increases to ~15 °C.

The distribution of ΔLST–NDVI pairs over the time is confirmed also by analyzing the other years of the study period, with some peculiarities (Figure 6); in particular, (i) the cloud spread changes their extension due to both ΔLST and NDVI; (ii) the two main clusters in some years (e.g., 2003–2004, 2005–2007, and 2009–2012) are more distinguishable than in others; (iii) the hole in the middle is some years (e.g., 2001, 2003, 2008, and 2012) is more pronounced than in the remaining years.

The typical trapezoidal shape, over small-medium areas, is only appreciable if NDVI and ΔLST are collected during at least one year, meaning that dry and wet edges positions are strongly dependent on the choice of the period.

Figure 7 shows the NDVI and ΔLST time-series (a and b panels, respectively). DOYs are arranged in-column while years are arranged in-line.

The higher NDVI values are found between DOYs (1–129) and (305–365) (green pixels). Several exceptions are found; i.e., years 2002 and 2005–2009 show lower NDVI values between DOYs 1 and 50. Minima NDVI are found between DOYs ≈190 and ≈260 (brown); two intermediates occur around DOY 140–160 and 260–305 (white).

ΔLST temporal behavior is less stable; maxima values occur between DOYs ≈210 and ≈230 (red-yellow pixels), even if several years exhibit maxima ΔLST shifted forward or backward. Winter is characterized by lower ΔLST (yellow-blue) because of increasing in both fractions of vegetation and average θ.

By selecting dry and wet pixels in the NDVI-ΔLST scatterplot (Figure 4a) in a region of interest (ROY) ± 2 °C buffering the wet and dry edges, it is emphasized that dry pixels (Figure 8, red areas) and wet ones (blue areas) were mainly remotely sensed during the spring-summer seasons and autumn-winter seasons, respectively.

In particular, dry pixels spread widely approximately between March and August, whereas wet pixels spread between November and January. In addition, pixels defining dry and wet edges change in time over different years.

Hereinafter, soil water content variables (U_LIS and TVX), input remote sensing variables (ΔLST and NDVI) and stomata conductance reduction factors were compared over a pixel at the station of Drasi.

The ΔLST (Figure 9a, red line) and NDVI (green line) time-series exhibit an average phase shift of ≈ 1.5 months (48 days): The determination coefficient, r², of the linear trend line between ΔLST_t and NDVI_t+48 rises from 0.36 (correlation was negative) to 0.71 (negative correlation, Figure 9b).

The resulting TVX behavior (Figure 9, black line) is in phase opposition to ΔLST, and out of phase with NDVI (88 days): The r² of the linear trend line between TVX_t and ΔLST_t rises up to 0.69 (negative correlation, Figure 9c); the r² of the linear trend line between TVX_t and NDVI_t+88 rises from 0.01 to 0.65 (positive correlation, Figure 9d). Due to the parallelism of dry and wet edges (at least in this case), TVX inherits the quick response behavior of ΔLST in comparison to the smooth temporal variability characterizing the NDVI.

The FFT and iFFT allowed comparing the main harmonics (Figure 10), and the sum of the main harmonics (black line) with the TVX, U_LIS, NDVI, and ΔLST time-series (a, b, c, and d panels, respectively).

Noticeably, U_LIS is characterized by the inter-annual (0 harmonic, of a periodicity ranging between 1528 to 1800 days) and annual harmonics (1 harmonic); whereas, NDVI and ΔLST and thus TVX, are characterized also by a 6-months harmonic (2 harmonic). Peaks of the different harmonics are shifted to each other as summarized in

Table 1. Variables’ main harmonics: frequencies and times. For the harmonic “0” no peak frequency was reported since a zero-shift FFT algorithm was applied.

Variable(s)	Harmonic	Peak Frequency (d−1 × 10−3)	Band-Pass Range (d−1 × 10−3)	Periodicity (d)	Peak Time (DOY)
ULIS	0		0–0.6	1700	1337/3041
	1	2.737	2.50–3.00	365	57
TVX	0		0–0.6	1800	1305/3105
	1	2.737	2.50–3.00	365	348
	2	5.457	5.25–5.70	182	168/350
ΔLST	0		0–0.6	1630	1433/3065
	1	2.737	2.50–3.00	365	200
	2	5.457	5.25–5.70	182	60/240
NDVI	0		0–0.6	1528	1537/3065
	1	2.737	2.50–3.00	365	70
	2	5.457	4.80–5.70	182	70/250

(Peak time column). The peak of harmonic zero is not reported as the shift zero-frequency component to the center of the spectrum was applied, while the peak time is computed and reported through the inverse zero-frequency shift.

3.3. Stress Factor Analysis

The TVX was compared to U_LIS, in terms of magnitude, temporal behavior, and phase shift for the whole study period. Figure 11 shows results for a limited period (2005–2006) characterized by the highest rainfall events.

Using a moving average window (dashed lines) for TVX (Figure 11, red line) and U_LIS (black line), it is possible to observe a general sinusoidal behavior of the two variables. Generally, the time-series are characterized by ~2 months (64–72 days) phase shift; with TVX and U_LIS peaks occurring in December and February, respectively; whereas minima occur in May and August, respectively; it is noteworthy that sinusoidal behaviors are typical of variables affected by a seasonal pattern. U_LIS is usually characterized by a steep upward sloping ramp (see the period between August 2015 and February 2016) and gentle downward sloping ramp (see the period between January and August 2005) if compared to TVX ramps. The behavior is typical of soil water content dynamic, where the wetting is faster than the drying and drying rates are expected to be most negative at the beginning of a dry down and trend toward zero [12].

The phase shift between TVX and U_LIS series is to ascribe to the role of the canopy conductance reduction factors. Figure 12 shows the spatial distributions of the average conductance reduction factors accounting for not optimal air temperature (left panel), air humidity (central panel), and solar radiation (right panel).

Environmental stress factors exhibit an inhomogeneous distribution: C_Tair causes the maximum reduction when pixels are bare or close to the dry edge; C_Rs exhibits a similar distribution but characterized by lower reductions; C_Uair tends to compensate this behavior assuming maxima when pixels are densely vegetated or close to the wet edge; although quantitatively, conductance reduction factors do not totally compensate each other. Results are coherent to Goward et al. [29] reporting that the changes in ground temperature due to incident solar radiation have similar magnitude as those caused by changes in modeled θ. They also found that wind speed plays a significant but modest role in determining both vegetation and soil temperatures. The authors also tested the influence of surface water vapor deficit that was found to have a smaller influence than wind speed; confirming that other stress factors can affect the estimations of the soil water content through the diachronic application of the triangle method, even though the role played by these factors is not quantified.

Even exhibiting a phase shift, as mentioned before, actual values of TVX and U_LIS range around their own sinusoidal carrier waves as a consequence of a quick response to changes in θ. Indeed, showing U_LIS and TVX time-series through two different y-axes (U_LIS in the primary axis, TVX in the secondary axis), it is highlighted that these variables behave synchronously to rainfall events (Figure 13 upper panel principal and secondary y-axis, respectively over selected years 2005–2006); even if the evaluation is not straightforward due to different time resolutions (8 days for TVX and 1 day for U_LIS).

3.4. Frequency-Domain Analysis

Based on this evidence, although the triangle method allows retrieving a stress index related to environmental factors, short-term changes in TVX appear to be related to θ through changes in soil admittance and evaporative control on the available energy.

However, the correlation between the residuals (U_LIS and TVX subtracted respectively by their main harmonics) if looking at the whole time-series (2001–2012) is weak (r ~ 0.21, rmse ~ 0.13), even though r is significant; indeed it is greater than the critical value, ≈0.146 (for level α= 0.001, degree of freedom equal to 500, two-tails test).

We further compared TVX and U_LIS to the measured total daily rainfall, P, as normalized amplitudes of the power spectrum in the range of frequencies ±0.01 d⁻¹ (Figure 14, upper and lower panels, respectively).

Remarkably, the best matching is observed between TVX and P normalized amplitudes of power spectra (r² ≈ 0.80, panel a, top-right corner), while the comparison between U_LIS and P normalized amplitudes of power spectra result in a lower determination coefficient (r² ≈ 0.73, panel b, top-right corner). Although r² indicate strong correlations, a systematic bias is shown for both of P vs. TVX and U_LIS (characterized by slopes of 0.82 and 0.5, respectively, and almost negligible intercepts). In addition, high rmse values (0.03 and 0.02 for P vs. TVX and U_LIS, respectively) are not negligible if compared to the range of variation the variables for harmonics “2” and higher (i.e., ~0–0.15). It should be highlighted that removing the main harmonic (≈ -0.0027 d⁻¹) corresponding to the top right point of the top-right plot r² decreases to ≈ 0.23 and 0.03 for TVX vs. P and U_LIS vs. P, respectively. Moreover, determination coefficients decrease to ≈0.49 and ≈0.43, for TVX and U_LIS, respectively, if we consider the whole range of frequencies: ±0.0625 d⁻¹ instead of ±0.01 d⁻¹ (not shown here), meaning that most of the correlation between P vs. TXV or U_LIS is due the harmonic “1” which represent the annual periodicity.

4. Conclusions

The application of the TVX methods using a multi-temporal approach ensured to incorporate both dry and wet extreme conditions for a given vegetation cover class. However, it has been demonstrated that, in this diachronic formulation, the assumption that the temperature vs. vegetation index features space is mainly controlled by soil water availability is not physically realistic. Indeed, the position of ΔLST–NDVI pairs and wet and dry edges as well, depend not only on root zone water stress but also on environmental stress factors (e.g., distance from the optimal air temperature, excessive incident solar radiation, and air humidity). This result is in accordance with other authors (e.g., Zhang et al. [39]) introducing an iterative method for locating the actual dry edge that eliminates the effects of four controlling factors (albedo, air temperature, surface resistance to evapotranspiration, and actual vapor pressure near the surface).

The temporal analysis reveals that the NDVI–ΔLST pairs move during the seasons following a circular behavior. Moreover, two main clusters, corresponding to the two main periods detected according to Pèguy climographs for the study area, are retrieved in the scatterplot. The dry period, generally from April–May to August–October, corresponds to low NDVI and high ΔLST; whereas, the remaining months are characterized by a temperate climate as well as higher NDVI values and the lower difference between day and night temperatures. Empirical evidence supports the hypothesis that TVX should be interpreted as a vegetation stress index rather than a soil water content index, at least in regions characterized by such a climate. However, the effect of the soil water content can be seen in TVX short term dynamics, as quick increases in TVX occur after rainfall events. This result confirms the potential value of the diachronic triangle method approach in surface hydrology.

The FFT analysis allowed to characterize TVX (and thus NDVI and ΔLST) and U_LIS time-series in terms of carrier waves. More specifically, it was observed that only the carrier wave of U_LIS is described by two harmonics (0- and 1-harmonics). The comparisons between TVX and U_LIS power spectra normalized amplitudes of the main harmonics with the rainfall one showed a better match for TVX (r² ~ 0.8) than for U_LIS (r² ~ 0.7). In this regard, the outcomes of this study suggest that the diachronic application of the triangle method has potential in retrieving time-series of the surface soil water content.