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

: 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 inﬂuence 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 e ﬀ ects 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 speciﬁc environmental and climatic conditions (strong seasonality and rainfall out of phase with vegetation growth), Normalized Di ﬀ erence 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 2 ~0.8) over the range of frequencies of the main harmonics.


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. Commons, the free media repository) and the catchment river network (on the right, blue lines) and discharge and meteorological gauge stations (red and green dots, respectively).

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 bidirectional surface reflectance and is delivered every 16 days at a 250-meter spatial resolution [21].
Night-time and day-time LST composite images (LSTD and LSTN) were collected from the 1 st of January 2001 to the 31 st of December 2012.
A reference dataset for the soil water content dynamic was derived from the LISFLOOD rainfallrunoff 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, ULIS, 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 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).

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 2 ≈ 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 2 decreases to ≈0.45, showing that LISFLOOD performances increase as the temporal scale increases.
Geosciences 2019, 9, x FOR PEER REVIEW 4 of 18 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 2 ≈ 0.62) (Figure 2, top-right corner), also the soil water content modeled by LISFLOOD, ULIS, 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 2 decreases to ≈0.45, showing that LISFLOOD performances increase as the temporal scale increases.

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] LSTD and LSTN images allowed to apply the thermal admittance version of the triangle method (i.e., using ΔLST = LSTD−LSTN). 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.

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, Geosciences 2020, 10, 23 5 of 18 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]: 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). 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]: where, and 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).

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.

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 2 ) on canopy resistance is neglected under the hypothesis that the atmospheric CO 2 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.

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.

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.

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).

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.

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. 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).

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) Δ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.    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 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 ∆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.
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. ΔLST temporal behavior is less stable; maxima values occur between DOYs ≈210 and ≈230 (redyellow 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. ΔLST temporal behavior is less stable; maxima values occur between DOYs ≈210 and ≈230 (redyellow 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 2 , 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).
Geosciences 2019, 9, x FOR PEER REVIEW 11 of 18 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 (ULIS 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 2 , of the linear trend line between ΔLSTt and NDVIt+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 2 of the linear trend line between TVXt and ΔLSTt rises up to 0.69 (negative correlation, Figure 9c); the r 2 of the linear trend line between TVXt and NDVIt+88 rises from 0.01 to 0.65 (positive correlation, Figure 7d). 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, ULIS, NDVI, and ΔLST time-series (a, b, c, and d panels, respectively). The resulting TVX behavior (Figure 9, black line) is in phase opposition to ∆LST, and out of phase with NDVI (88 days): The r 2 of the linear trend line between TVX t and ∆LST t rises up to 0.69 (negative correlation, Figure 9c); the r 2 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 (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. Noticeably, ULIS 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 (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. 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.

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.

Stress Factor Analysis
The TVX was compared to ULIS, 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.
Time (dd/m) Figure 11. Soil water content modeled via LISFLOOD (ULIS, black lines) vs. soil water content index estimated via the TVX (red lines).
Using a moving average window (dashed lines) for TVX (Figure 11, red line) and ULIS (black line), it is possible to observe a general sinusoidal behavior of the two variables. Generally, the timeseries are characterized by ~2 months (64-72 days) phase shift; with TVX and ULIS 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. ULIS 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 ULIS 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). 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).

Stress Factor Analysis
The TVX was compared to ULIS, in terms of magnitude, temporal behavior, and phase shift for the whole study period. Figure 11 shows results for a limited period (2005)(2006)  Using a moving average window (dashed lines) for TVX ( Figure 11, red line) and ULIS (black line), it is possible to observe a general sinusoidal behavior of the two variables. Generally, the timeseries are characterized by ~2 months (64-72 days) phase shift; with TVX and ULIS 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. ULIS 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 ULIS 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 ).

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 −1 (Figure 14, upper and lower panels, respectively).
Remarkably, the best matching is observed between TVX and P normalized amplitudes of power spectra (r 2 ≈ 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 2 ≈ 0.73, panel b, top-right corner). Although r 2 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 −1 ) corresponding to the top right point of the top-right plot r 2 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 −1 instead of ±0.01 d −1 (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. We further compared TVX and ULIS to the measured total daily rainfall, P, as normalized amplitudes of the power spectrum in the range of frequencies ±0.01 d −1 (Figure 14, upper and lower panels, respectively). Remarkably, the best matching is observed between TVX and P normalized amplitudes of power spectra (r 2 ≈ 0.80, panel a, top-right corner), while the comparison between ULIS and P normalized amplitudes of power spectra result in a lower determination coefficient (r 2 ≈ 0.73, panel b, top-right corner). Although r 2 indicate strong correlations, a systematic bias is shown for both of P vs. TVX and ULIS (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 ULIS, 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 −1 ) corresponding to the top right point of the top-right plot r 2 decreases to ≈ 0.23 and 0.03 for TVX vs. P and ULIS vs. P, respectively. Moreover, determination coefficients decrease to ≈0.49 and ≈0.43, for TVX and ULIS, respectively, if we consider the whole range of frequencies: ±0.0625 d −1 instead of ±0.01 d −1 (not shown here), meaning that most

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 2~0 .8) than for U LIS (r 2~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.