Soil Water Content Assessment: Critical Issues Concerning the Operational Application of the Triangle Method

Knowledge of soil water content plays a key role in water management efforts to improve irrigation efficiency. Among the indirect estimation methods of soil water content via Earth Observation data is the triangle method, used to analyze optical and thermal features because these are primarily controlled by water content within the near-surface evaporation layer and root zone in bare and vegetated soils. Although the soil-vegetation-atmosphere transfer theory describes the ongoing processes, theoretical models reveal limits for operational use. When applying simplified empirical formulations, meteorological forcing could be replaced with alternative variables when the above-canopy temperature is unknown, to mitigate the effects of calibration inaccuracies or to account for the temporal admittance of the soil. However, if applied over a limited area, a characterization of both dry and wet edges could not be properly achieved; thus, a multi-temporal analysis can be exploited to include outer extremes in soil water content. A diachronic empirical approach introduces the need to assume a constancy of other meteorological forcing variables that control thermal features. Airborne images were acquired on a Sicilian vineyard during most of an entire irrigation period (fruit-set to ripening stages, vintage 2008), during which in situ soil water content was measured to set up the triangle method. Within this framework, we tested the triangle method by employing alternative thermal forcing. The results were inaccurate when air temperature at airborne acquisition was employed. Sonic and aerodynamic air temperatures confirmed and partially explained the limits of simultaneous meteorological forcing, and the use of proxy variables improved model accuracy. The analysis indicates that high spatial resolution does not necessarily imply higher accuracies.

Two such methods are the "thermal inertia" and "triangle" methods. The thermal inertia method (e.g., [10][11][12]) determines θ using a combination of thermal and visible images. Because the ability of a material to accumulate heat and release it within a defined time interval strongly depends on the water content of the material itself, θ can be indirectly inferred by measuring the variation of the body temperature over time. In particular, the short wave albedo, compared to the day-night temperature difference, can be used to estimate the soil surface water content because the amplitude of the temperature variation is a direct consequence of the amount of energy reaching the soil in the shortwave and of the soil properties, including bulk density, thermal conductivity, and heat capacity. For bare ground, soil thermal properties strongly depend on water content and soil density. Thermal inertia is therefore higher for wetter soils because water in the soil pores absorbs heat, and thus soil temperature changes more slowly than in dry soils.
Because the variation in surface temperature depends on the surface water content, the triangle method exploits this relationship [13] through evaporation over bare soil and transpiration over dense vegetation. Vegetation cover determines applicability of these methods: (i) the triangle method depends on the temperature difference between canopy and air, a relationship based on canopy resistance that varies with soil water availability; thus, the method is more suitable for vegetated than bare soils; (ii) in contrast, operational thermal inertia is based on the hypothesis that some fluxes are linearly related to surface temperature; thus, the method is more suitable for bare soil than vegetated areas.
Our research focuses on a method based on the feature space defined by optical and thermal data. In particular we discuss advantages and limitations of the triangle method and a simplified version of the thermal inertia method. According to Petropoulos et al. [14], the wide base of the triangular envelope depicts the relatively higher sensitivity of the bare soil temperature to its water content changes compared to the much lower sensitivity of progressively denser vegetation, which is highlighted by the triangle's much narrower vertex.
We implemented a model based on the empirical approach of the triangle method [15,16], named for the triangle enveloping the typical trapezoidal shape of the scatterplot of the land surface temperature (LST) vs. a vegetation index (VI), such as the Normalized Difference Vegetation Index (NDVI). The minor trapezoid base corresponds to densely vegetated areas and highlights the low temperature variability caused by the thermoregulation mechanism of the plants, which are able to take water (if available) from soil layers not directly linked to soil surface evaporation. In contrast, the major base of the trapezoid corresponds to low or absent vegetation coverage and confirms the strong relationship between surface θ and surface temperature of bare soils.
Each interval of vegetation coverage can be classified by the direct correspondence between the "warm edge" characterized by dry conditions and thus also referred to as "dry edge", and the "cold edge" characterized by wet conditions and alternatively named "wet edge". The soil water content index can then be defined as a function of the relative position to the triangle edges of a generic point of the scatterplot. One assumption is that the time series extent be sufficient to include all possible soil water content conditions (from residual to saturation) and soil vegetation coverage to assess θ based on the hypotheses of linear variation of isopleths from the wet to the dry edge.

Theoretical Background
The soil water content is evaluated using both VIs and surface temperatures via the triangle method. A fundamental assumption is that, given a large number of pixels describing a full range of soil surface water content and fractional vegetation cover, sharp boundaries in the data describe real physical limits: bare soil, pixels fully covered by vegetation, and lower and upper limits of the surface θ (completely dry and at field capacity) [15].
Under these hypotheses, Jackson et al. [17] showed that air temperature minus canopy temperature can be used as an index of crop water stress. Their results are based on the surface energy balance equation for a crop canopy, taking into account the net radiation and soil heat flux (Rn and G0, respectively) and the sensible and latent heat fluxes (H and λE, respectively). Under atmospheric neutral conditions, sensible and latent heat fluxes can be expressed as where ρ is the air density, cp is the air heat capacity, Tc and Ta are the canopy and air temperatures, ec * and ea are the saturated and vapor pressures of the air, γ is the psychrometric constant, and ra and rc are the aerodynamic and canopy resistances to water transport. The temperature difference between canopy and above air can be written as [17]: where Δ is the slope of the saturated vapor pressure vs. temperature relationship (ec * − ea * )/(Tc − Ta), and ea * − ea is the vapor pressure deficit. Theoretically, Equation (1) can be used to define two boundary limits for well-watered and dry crops. For well-watered crops and full cover, the canopy resistance reaches potential value, rc = rcp, whereas in bare soil, rc is approximated as zero. The potential value of canopy resistance is obtained per unit of Leaf Area Index (LAI) from minimum stomata resistance: rcp = rSmin/LAI [18]. These resistances determine the theoretical lower limit (wet edge).
For no water availability and full vegetation cover, the canopy resistance rises to the stomata closure value rc = rcx = rSmax/LAI as a function of the leaf stomata closure value [18], whereas for dry bare soil rc is set to infinity, rc = ∞. These resistances determine the theoretical upper limit (dry edge).
Following Ortega-Farias et al. [19], rcp over vineyards could be set to 25 s·m −1 (cultivar Savignon), although values ranged between 10 and 100 s·m −1 [20]. By following Giordani et al. [21], rcx over vineyards associated with nearly complete stomata closure could be set to 2000 s·m −1 . Thus, defining the value that describes the actual behavior of the vegetation is difficult. Regarding LAI influence, if none of the controlling variables is limiting, rs could be roughly assumed to be 40 s·m −1 [22]. Net radiation is obtained by summing shortwave and longwave net radiations: The equation is implicit for surface temperature and thus must be solved iteratively.
Canopy resistance (rc) can be expressed as a function of the bulk stomata resistance (rs) per unit of active (or green) LAI that actively contributes to the surface heat and vapor transfer: The bulk stomata resistance is crop-specific and differs among crop varieties (and crop management). Because it usually increases as the crop matures and ripens, it is difficult to accurately evaluate in operational applications.
Stomata resistance depends on several environmental variables, including photosyntetically active radiation (PAR); atmospheric water vapor deficit (AMD); atmospheric carbon dioxide concentration (CO2), air temperature (Ta); and θ. Conventionally, stomata resistance (rs, defined as the inverse of conductance, G = rs −1 ), is parameterized as a function of reduction coefficients accounting for these actual environmental variables, based on their minimum values observed under optimal conditions: Within the diachronic empirical approach, the mutually compensating effect of meteorological forcing (PAR, AMD, Ta and CO2) on rs must be verified. This evaluation must occur at least at acquisition time, and preferably carried out with the same sun elevation angle (airborne platform) or at the same time of day (sun-synchronous satellites).
The empirical method bypasses efforts in assessing rc by adjusting dry and wet edges to the experimental scatterplot. This means that improvements can be achieved by considering a VI time series and the difference (ΔT) between LST and a reference temperature T * . We tested different T * and propose a procedure to select the best model parameters.

Methods
The triangular shape [23] characterizing the scatterplot, relies on the assumption that the temperature of an elementary heterogeneous surface is given by the linear combination of the soil and vegetation radiative temperatures, and these are in linear relationship with the evaporation and transpiration processes, respectively.
Given the VI, the soil water content of a pixel is often assumed to be proportional to the ratio between its thermal difference with the dry edge compared to the total thermal excursion between dry and wet edges. Two straight lines identify dry and wet edges. Within this research, a θ index for the generic k element (NDVIk, ΔTk) is given by the ratio of two angles [24], βk and α. Angle βk is given by the wet edge and the straight line joining the generic pixel of the scatterplot and the triangle vertex. Angle α is between the dry and wet edges and is thus proportional to the whole variability ( Figure 1). This index is then converted into absolute values using the residual (θres) and saturation (θsat) outer values. This index strongly depends on the positioning of the dry and wet edges, which are determined by means of linear regression of minima and maxima percentiles for a given NDVI class. The residual and saturated θ values, θres and θsat, were determined by laboratory analyses on samples collected within the study area.

Edge Determination through a Diachronic High Spatial Resolution Dataset Acquired on a Small Area
This analysis aimed to determine how an operational triangle method approach can be developed to assess θ over a small, vegetated area using a diachronic Earth Observation (EO) dataset of thermal and VIS/NIR images characterized by high spatial resolution. Once the method is applied on an empirical basis, wet and dry edges can be determined directly from the optical-thermal feature space by setting up few parameters ( Figure 2, tuning parameters). Over small areas, a diachronic approach is valuable because time-series describe, for given vegetation cover, a wider range of variability of θ; thus, including all the NDVI-ΔT pairs within the scatterplot could be opportune (and in some cases needed). The use of high spatial resolution does not necessarily imply best accuracies in θ retrieval, and thus the appropriate resolution of θ maps requires a further discussion. Wet and dry edges must be assessed at the finest spatial resolution (usually the resolution of thermal images). Spatial aggregation reduces both the range of variability and total number of pairs, threatening the characterization of θ minima and maxima values.

Parameter Tuning
Parameters were defined to optimize the match between remote sensing θ (θRs) and θ measured in situ ( Figure 2). Parameters included minima and maxima percentiles (Pmin and Pmax) and minimum NDVI of the wet and dry edges (NDVImin,D and NDVImin,W); whereas the cluster shape allows fixing maximum NDVI (NDVImax,D and NDVImax,W) with confidence.
The wet and dry edges were determined by applying a linear regression to the minimum and maximum percentiles of NDVI classes; this choice gives a certain degree of freedom to the operator because it directly influences the final result. The experiment shows that percentiles cannot be set a priori, and thus, the percentiles that best fit the in situ data must be determined.
Increasing Pmin and decreasing Pmax bring the edges closer to the cluster. The choice of different percentiles determines the value of α (Figure 1), representing the range of variability of θ (Δθ). The more Pmin increases, the more βk decreases ( Figure 1); consequently, θ rises. As a first attempt (rough tuning) the following percentiles were tested: Pmin = [2,5,10,20,30] and Pmax = [98, 95, 90, 80, 70]. Resulting maps represent θ driving evaporation and transpiration of the soil vegetation system; however, pixels characterizing other elements of the study area could still remain within the scatterplot. Those pixels could distort the expected trapezoidal shape of the cluster. Hence, it is possible, sometime necessary, to remove them. Although input images can be masked, both NDVImin,D and NDVImin,W have to be set to this aim. During the rough tuning, minimum NDVI arrays of NDVImin,W = NDVImin,D = [0.1, 0.2, 0.3, 0.4] were tested, and NDVImax,W and NDVImax,D were fixed to 0.90 and 0.99, respectively. Parameters led to 450 model runs for each T * .

Reference Temperature Tests
Although in principle Ta, measured simultaneously with LST should be employed, a weak agreement with in situ θ suggests testing alternative T * (Table 1). In particular, two tests were conducted to analyze (a) whether Ta does or does not work ( * test); (b) if operational T * directly retrieved from thermal images could be used ( * test); and (c) if including thermal admittance could lead to more accurate results ( * test). The * test includes: (i) Ta measured by a thermoigrometer installed 2.75 m above ground level (a.g.l.), ≈1 m above vegetation top; (ii) the sonic air temperature (Tsonic), retrieved by sound velocity measured by a 3D sonic anemometer installed ≈3 m a.g.l.; its value closely approximates Ta. This test verifies if the weak results obtaining using Ta are eventually due to thermoigrometer malfunctioning; (iii) the aerodynamic temperature (Taero) retrieved by Two-Source Energy Balance model based on intercalibration (TSEB-IC [25]). Although Taero is not suitable for an operational use, it determines turbulent fluxes that underlie the theoretical basis of the method, thus identifying the limits of Ta in describing the fluxes.
The * test assesses the operational use of temperatures of particular targets, eventually included within the scene: (i) the water surface temperature of irrigation ponds (Tlake); (ii) the temperature of dense-well watered vegetation (Tveg).
The * test verifies if accuracy is improved when thermal admittance is accounted for; to this aim, two temperatures were tested: (i) the minimum soil surface radiometric temperature (Tsoil_min); (ii) the minimum air temperature (Ta_min), which is easier to obtain because it is usually measured by standard meteorological stations.

Parameter Configuration and Reference Temperature Selection
The choice of the reference temperature and parameter configuration characterizes θRs coherently to the measured θ. The choice is based on the determination coefficient (r 2 ), Student test value (T), Fischer test value (F), and the difference between maximum and minimum percentiles used to define the dry and wet edges (ΔP). Parameters must be configured to maximize r 2 , T, F, and ΔP; thus, a synthetic parameter configuration maximizing the product, I = r 2 ·T·F·ΔP, was chosen. The mean absolute error (MAE) is qualitatively verified to assume a low value. The slope and intercept (m and q, respectively) between estimated and measured data are qualitatively verified to approach the unity (m), and the null value (q). After the first approximation of the parameters, fine-tuning yields the optimal configuration.

Optimal Spatial Aggregation
Parameters are configured at the spatial resolution of the thermal acquisition using the whole diachronic dataset; the whole dataset eventually accounts for a sample size large enough to characterize dry and wet edges for the whole range of the vegetation index (outer θ conditions). Thus, the output RS of θ is originally the one of thermal acquisition. Nevertheless, not necessarily the better RS (more detailed) implies the better accuracy of θ; this is could be due to several factors including the spatial variability of the soil vegetation system, the bidirectional reflectance distribution function (BRDF) of both thermal and optical (VIS/NIR) bands, and the row orientation. Thus, the spatial aggregation of θ maps could increase their accuracy. Different aggregation scales were tested (2, 3, 4, 5 and 6 times RS), and the aggregation RS representing the best agreement with the in situ data was selected by maximizing r 2 , minimizing q and MAE, and approximating m to unity (Figure 3).

Study Area
The study area includes the experimental fields of the "Tenute Rapitalà" farm within the Camporeale district (Sicily), a specific zone that grants the Alcamo Denomination of Controlled Origin (DOC) status. The fields are characterized by an average slope of 10%, mainly S-SW aspect (<100°), and elevations ranging from 290 to 320 m a.s.l. The soil texture (USDA) was classified as loam (20% clay, 29% silt) to silty clay loam (37% clay, 45% silt). The organic content ranged between 1.4% and 1.7%. The residual and saturation volumetric water contents (θR and θS) were 0.040 and 0.453 m 3 ·m −3 , respectively, whereas volumetric water content at field capacity (θFC) was 0.239 m 3 ·m −3 . According to Centeno et al. [26], permanent wilting point (θWP) was assumed to be 0.096 m 3 ·m −3 (cultivar "Tempranillo", sandy loam soil).
The experimental field included four vineyards separated by cross-shaped section-breakers. The vineyards were cultivated adopting a vertical trellis system (Figure 4, lower-left box), with the cultivar "Nero d'Avola" at NW and NE; the SE field was partly cultivated with "Cabernet Franc" cultivar and partly with "Nero d'Avola", and finally the SW field was partly cultivated with "Syrah" cultivar and partly with "Nero d'Avola" (Figure 4). Other cultivars included "Insolia", "Cabernet Sauvignon", "Sauvignon Blanc", "Merlot", "Grillo", "Catarratto Lucido" and "Viogner". Plant spacing was 240 cm between rows and 95 cm between the plants within the rows. Plant density was 4386 plants per ha; rows were oriented N-NE to S-SW.

Data
Remote sensing images in the short and long waves were acquired by an airborne platform, SKY ARROW 650 TC/TCNS, at a height of ≈1000 m a.g.l. The VIS and NIR images (three bands) were acquired by means of a multispectral MS4100 camera (Duncantech, Auburn, CA, USA, 767-832, 650-690, and 530-570 nm), and the thermal infrared band (TIR) images were acquired by means of a SC500/A40M camera (Flir, Wilsonville, OR, USA, 7.5-13 μm). The nominal pixel resolution was approximately 0.6 m for the VIS/NIR bands and 1.7 m for the TIR channel. The VIS/NIR images were aggregated at a coarser spatial resolution. Details about the field campaign and remote sensing data can be found in Maltese et al. [24]. The installed micro-meteorological station is described in [27].
Volumetric water content was measured with a CS616 Water Content Reflectometer (Campbell Scientific, Inc., Logan, UT, USA). The probe consisted of two 30 cm long stainless steel rods introduced vertically in the soil between 10 and 40 cm below ground level. An ARG100 Tipping Bucket Raingauge (Campbell Scientific, Inc.) installed 1 m above the vegetation top measured rainfall, tipping once for each 0.2 mm of rain. LAI was measured using an optoelectronics instrument (LAI-2000 Plant Canopy Analyzer, by LI-COR ® , Lincoln, NE, USA). Radiometric temperature of the bare soil was measured by an IRTS-P Precision InfraRed Temperature Sensor pyranometer (Campbell Scientific, Inc.).
Vineyards were characterized by low active LAI in the first ten days of June (≈1.2 m 2 ·m −2 , the average value at pixel scale during the first two acquisitions). LAI increased up to ≈1.7 m 2 ·m −2 in the last ten days of July (third image acquisition) and decreased again to a minimum of ≈0.95 m 2 ·m −2 at the beginning of September (last acquisition).
Simultaneous to the image acquisitions, in situ measurements were carried out, including θ, measured using both the traditional thermo-gravimetric method and a handheld TDR (FieldScout TDR300, Spectrum Technologies Inc., Aurora, IL, USA). During each airborne acquisition date (11 June, 3 and 22 July, and 3 September 2008 corresponding to 163, 185, 204, 247 days of year, DOYs, respectively) 36 TDR measurements of bulk dielectric dielectric permittivity (ε) of the soil matrix were performed for each of the 6 selected plots (Figure 4): three (A, D and F) within the Nero D'Avola cultivar; two within a Cabernet Franc (C and D), and the last (E) in the Syrah plot. Half of these measurements were acquired in the upper soil layer (0-20 cm) and the remaining half in the upper part of the root zone (~25-35 cm), totaling about 1100 observations. Soil water content measured in the upper part of the root zone was compared with the remote sensing assessment. Additionally, soil samples were collected to apply the thermo-gravimetric method and to assess the bulk density, needed to convert reflectometric measurements into θ. The Soil Hydrology Laboratory of the Department of Engineering and Agro-Forest Technology (ITAF) of the University of Palermo performed the soil characterization (including grain size curve, residual and saturation θ, organic content, etc.).
The short wave images were calibrated to spectral reflectance and corrected for atmospheric influence, applying the empirical line method [28] requiring for ground reflectance measurements over targets taken simultaneously with the images. The thermal images were converted into surface radiometric temperatures by means of a linear regression with in situ temperature measurements using the emissivity maps retrieved as a function of the NDVI, as proposed by Sobrino et al. [29]. To do this, a number of spectroradiometric and radiometric measurements were performed on artificial surfaces (black and white panels) and on natural surfaces with homogeneous radiometric characteristics at the spatial acquisition scale (bare soils, roads, small irrigation reservoirs).
The vineyards were irrigated using a controlled water deficit technique, a drop irrigation system able to supply 4 L·h −1 to each plant. The farm irrigation strategy included nine drop irrigations ranging between 4 and 7 h in duration each from 28 May to 18 August. Watering was provided twice between the end of May and June and 4 times in the first twenty days of July during flowering, fruit setting, and veraison. Water management during this period has a direct effect on the subsequent berry average weight, and irrigation between the end of July and August aims to optimize growth as well as berry sugar and acidity contents.
Rainfall events during the study period were few and characterized by low accumulations. Values ranged between 6 and 8 mm in the first 4 events during 10-28 May, and the subsequent 3 months were dry until 14 September, when a significant rainfall occurred (196 mm).
During this research, the θ extremes chosen to characterize the dry and wet edges were 0.040 m 3 ·m −3 (θR) and 0.453 m 3 ·m −3 (θS), determined by in situ soil sampling and laboratory measurements.
A micro-meteorological station measured several variables, including PAR, Tair, CO2, and θair, during the whole period. Recorded data were used to evaluate the product of the stomata conductance reduction factors (FPAR FTa FCO 2 Fθair); their product was almost constant with time (variation coefficient ≈ 14%).

Parameter Tuning
As reported in the Methods section, the positioning of wet and dry edges was determined by applying a linear regression to the minimum and maximum ΔT percentiles of NDVI classes; thus, their positions depend on the choice of percentiles Pmin and Pmax. Statistics of r 2 , T-test, F-test, ΔP, and MAE (the case of T * = Ta_min, was chosen as an example) for varying Pmin m, and q ( Figure 5) show clear uncertainties in the choice of a suitable Pmin. High ranges of variability of these indices confirm how the parameter settings affect the result. Thus, the use of a synthetic index (see index I, Section 3) could be valuable for the parameter selection. No noticeable behavior of these statistics has been found with Pmax. For each T * , statistics characterizing in situ vs. remote sensing θ of the 450 model runs were analyzed. The process for selecting the best parameter set was explained in Methods. Both the best set of parameters and statistical indices characterizing in situ θ vs. θRs were synthesized (Table 2). Results of the * test revealed that the best parameter setup using Ta as reference temperature produces low accuracy maps (I ≈ 0.03), and Tsonic revealed that these weak performances are not related to measurement inaccuracies or failures (providing a small increment of I, ΔI ≈ +0.06). Taero should drive actual H and, subsequently, λET; thus, better results are expected if Taero is used instead of Ta because it determines turbulent fluxes. Results (Table 2) confirm these hypotheses (ΔI ≈ +0.22); thus, we conclude that Ta is not suitable as reference temperature.
Results of the * test showed that some reference temperatures directly retrievable on thermal images provided an accurate assessment of θ. In particular, the temperature of a small irrigation pond, Tlake, achieved a similar level of agreement with in situ θ (ΔI ≈ +0.21) if compared to that obtained using Taereo (lower accuracy is obtained using Tveg, ΔI ≈ +0.18). The * test revealed that accurate results can be achieved without the use of ancillary meteorological data.
The * test confirmed that accuracy could be improved by accounting for thermal admittance. Indeed, agreement with in situ θ notably increased (ΔI ≈ +0.38). However, measures of surface soil temperature are not commonly available; thus, employing Tsoil_min has limited operational utility. From an operational standpoint we tested Ta_min (easily available) as a surrogate of Tsoil_min. Ta_min lead to results (ΔI ≈ +0.37) comparable to those achieved by Tsoil_min. The level of agreement with in situ θ was even higher if compared to that obtained using Taereo.
Because the best setup configuration using Ta_min provided the higher ΔI during the rough tuning procedure, T * = Ta_min was chosen for a subsequent fine-tuning of parameters to optimize the parameter setup (Table 3). Results revealed that all statistical indices are improved, in particular ΔI remarkably increased up to ≈+0.70. A scatterplot of NDVI vs. ΔT pairs for the whole diachronic dataset is plotted in Figure 6 (red, orange, cyan and blue color ramp indicates increasing θ). Individual clusters highlight that an EO acquired on a small area during a single DOY was not able to characterize both dry and wet edges, and a diachronic dataset is required. Minima θ were reached on DOY 163 (pairs drawn in red) and determined the dry edge, whereas maxima θ were reached on DOY 204 (pairs represented in blue) and determined the wet edge.

Spatial Resolution Analysis
As reported in Methods, higher RS do not necessarily lead to higher accuracy of θ maps; thus, different aggregation scales were tested (Figure 7). Six times RS, 6.8 m (approximately four times the distance between rows, 7.2 m), improved the retrieval accuracy; r 2 increased from ≈0.89 to ≈0.91; m approached the unity, from ≈0.96 to ≈0.97; and q and MAE were minimized, from ≈0.43 and ≈0.83 to ≈0.35 and ≈0.79, respectively. Results are significant mainly because comparable or even better accuracy can be achieved by acquiring images at lower RS, thus reducing the imaging cost.

Soil Water Content Spatial Distribution
The θ map characterizing DOY 204 (Figure 8, left panel; lakes and buildings are masked out) was chosen because it exhibits the higher spatial variability. The temporal variability (Figure 8, right panel) shows that θ approaches the dry edge on the whole area at the beginning of irrigation period (DOY 163); on DOY 185, average θ increases, and its distribution shows a statistical tendency toward wet pixels; the skewness is confirmed on DOY 204 when θ values are maximized. Finally, at the end of the irrigation period (DOY 247), soil dries and θ average value and distribution are similar to the beginning of the irrigation period.
In irrigation practices, only a percentage of available water capacity (AWC) is depleted because plants start to experience water stress even before soil water reaches θwp. Therefore, a maximum allowed depletion (MAD, %) of the AWC must be specified. For example, MAD values of 50%-55% were used by Ortega-Farias and Acevedo [30] to schedule the irrigation of a Cabernet Sauvignon vineyard located in Chile. The actual percentage of available water (CAW) to the plant roots (CAW = θ − θWP/θFC − θWP) at the acquisition time at RS = 6.8 m is plotted in Figure 8 (upper right panel). The isopleths allow identifying areas above a chosen MAD.

Conclusions and Outlook
This research highlights the limits of the triangle method if applied to small areas with high spatial resolution images on a strict theoretical basis. Additional efforts are required to parameterize the theoretical method. The empirical method bypasses these efforts by setting up some parameters and adjusting the outer edges bounding the feature space defined by optical and thermal images. An operational procedure was proposed that exploited a diachronic approach to reduce operator influence on setting up model parameters. Results confirmed that the use of a diachronic approach includes a wider range of variability of soil water content (θ) given the vegetation cover, thus making the method applicable on small areas using high spatial resolution data.
Different alternative reference temperatures were tested, and some achieved high θ accuracy. Results revealed that if measures of simultaneous air temperature (Ta) are used as reference temperature, weak agreement with in situ θ is achieved (r 2 ~ 0.1, MAE ~ 3.3). The unsuitability of Ta was confirmed by the weak results achieved also using the sonic temperature (Tsonic) (r 2 ~ 0.3, MAE ~ 2.7) and by the accurate results obtained using the aerodynamic temperature (Taero) (r 2 ~ 0.9, MAE ~ 1.4). The retrieval of reference temperatures directly from thermal images (if suitable targets are available) promises advantages. The water surface temperature of an irrigation pond (Tlake) and the temperature of a well-watered dense vegetation (Tveg) were tested. According to the in situ θ, Tlake was better suited as a reference temperature (r 2 ~ 0.9, MAE ~ 1.3) than Tveg (r 2 ~ 0.6, MAE ~ 1.8).
The * test closely agreed with in situ θ, confirming that minima temperatures of soil (Tsoil_min) allows accounting for the soil admittance (r 2 ~ 0.85, MAE ~ 1.1); if the Tsoil_min is unknown, the use of minimum air temperature (Ta_min) as surrogate also yields accurate assessments (r 2 ~ 0.9, MAE ~ 1.1). In addition, these data are available by standard meteorological stations.
A final test was employed to verify if the higher spatial resolution implies more accurate results compared with in situ observations. This test, performed by aggregating the remotely assessed θ on the best T * and model parameter set, confirmed that the optimum spatial resolution (r 2 ~ 0.9, MAE ~ 0.8 at RS ~ 7m) is a multiple of the spatial land fragmentation of the observed soil-vegetation system (~three times the distance between rows, in this case).
There are concerns about the assumption of linear behavior of the isopleths within the temperature-vegetation index space. Stisen et al. [31] suggested a nonlinear interpretation of the surface temperature-vegetation index domain. Krapez et al. [32] simulated the isopleths of the root zone between residual and saturation values within a T-NDVI space using the SEtHyS model [33]. Simulated isopleths were nonlinear, showing curvatures for low and high vegetation covers. Curvature for low vegetation cover may originate from aerodynamic parameterization for bare or scarcely vegetated soils; nonlinear behavior for dense vegetation is reduced for high root zone water content.
Future research must take into account the nonlinear behavior of the isopleths within the temperature-vegetation index space and the different sensitivities of bare soil temperature to θ changes if compared with the vegetated surfaces.

rcx and rcp
Stomata and potential closure resistances, resp. rSmin and rSmax Minimum and maximum stomata resistances, resp.