Satellite-Based Method for Estimating the Spatial Distribution of Crop Evapotranspiration : Sensitivity to the Priestley-Taylor Coefficient

This work discusses an operational method for actual evapotranspiration (ET) retrieval from remote sensing, considering a minimum quantity of ancillary data. The method consists in a graphical approach based on the Priestley-Taylor (PT) equation, where the dry soil and non-limiting water conditions are defined by land surface temperature (LST) and vegetation index (VI) space, both retrieved from remote sensing. Using ET tower flux measurements and Landsat 5 TM images of an irrigation scheme in southeast Spain, a sensitivity analysis of ET spatial distribution was performed for the period 2009–2011 with respect to: (i) the shape (trapezoidal or rectangular) of the LST-VI space; and (ii) the value of the PT coefficient, α. The results from ground truth validation were satisfactory, both shapes providing similar performances in estimating ET, with root mean square error ~30 W/m2 and relative difference ~10% with respect to tower-based measurements. Importantly, the best fit with ground data was found for α close to 1, a somewhat different value from the commonly used value of 1.27, indicating that substantial error might arise when using the latter value. Overall, our study underlines the importance of a more precise knowledge of the actual value of α coefficient when using ET retrieval methods based on the LST-VI space.


Introduction
Remote sensing (RS) techniques, information technology and geospatial methods are potentially useful to analyze the spatial distribution of land evapotranspiration (ET) in agricultural regions [1].This information is highly valuable for implementing water management practices, such as precision irrigation and fertilization that will optimize inputs, improve yields and promote sustainable agricultural practices [2].Such RS-based tools are especially welcome in semiarid areas where water shortages are a major obstacle to agricultural production, economic welfare and sustainable development [3][4][5].This is the case of the Segura River Basin (SRB), a highly productive agricultural region of southeast Spain, with the highest agricultural water demand and the lowest percentage of renewable water resources of all Spanish basins [6].Cultivated lands represent 43% of its area, with one-third of the surface under irrigation (more than 269,000 ha).The water from non-conventional resources (water transfer, desalination and reuse) accounts for more than 65% of the total available [6] with the consequent higher costs than conventional resources.In this context, accurate estimations of crop water consumption (i.e., ET) are crucial to optimize water allocation, save water and reduce irrigation-related costs.
The ET-retrieval methods can be classified into two main groups.The first one is based on the so-called "residual" method, which determines the latent heat flux as the residual term of the surface energy balance (SEB) [7].
where the different fluxes, expressed in W•m −2 , are net radiation (R n ), latent heat flux (λE, with λ = latent heat of vaporization in J•g −1 and E = actual evaporation in g•m −2 •s −1 ), sensible heat flux (H) and soil heat flux (G).Several variants of the residual method are: (i) the simplified method SM [8,9]; (ii) the SEBAL model (Surface Energy Balance Algorithm for Land [10,11]); and (iii) the TSEB model (Two Source Energy Balance [12][13][14]).Such algorithms are currently used because of their relative simplicity, but require ancillary ground data (such as air temperature and wind speed), and their extrapolation from the meteorological stations to other pixels [15].In semiarid regions, the high values of the sensible heat flux (H) make the method very sensitive to small errors in the estimation of H.
The second group is formed by the so-called "graphical methods", based on the interpretation of the land surface temperature (LST) vs. vegetation index (VI) scatterplot.This method was first suggested by [16], and later successfully applied to retrieve ET from remote sensing data [17].The LST-VI space usually presents a geometrical shape (typically, triangular, trapezoidal or rectangular), whose boundaries can be interpreted in terms of evaporative limits [18].The upper limit ("dry" or "warm" edge) corresponds to null ET, while the lower limit ("wet" or "cold" edge) represents the maximum ET in the region.A commonly used assumption is that the pixels of the "wet" edge of the LST-VI space correspond to areas with ET equal to the "wet areal evaporation" (E w ), defined as the evaporation of a well-moist surface under conditions of minimal regional advection [19]: where ∆ (Pa•K −1 ) is the slope of the saturated vapor pressure curve at the prevailing T air , γ is the psychrometric constant (Pa•K −1 ) and α is the Priestley-Taylor coefficient (PT, dimensionless).
The position of the pixel in the LST-VI space with respect to the dry and wet edge allows the determination of the evaporative fraction, EF, defined as the ratio of latent heat flux λE to available energy, A, as follows: with A = R n − G. Equation (3) indicates that, besides EF, the estimation of λE requires the assessment of two other fluxes, R n and G. Compared to the residual method, the graphical method presents the important advantage of including the valuable information contained in the vegetation index, an indicator of vegetation cover and health status.Nevertheless, the use of vegetation indices is considered for most SEB approaches.Different drawbacks of the graphical methods were analyzed [20,21].Therefore, several shortcomings affect the accuracy and reliability of the method.One of them is the uncertainty on the true value of the Priestley-Taylor coefficient (α).A generally accepted value for α is 1.27 [19], although it is recognized that α depends on the characteristics and vegetation cover of the surface.In particular, the authors of [22] found that α could range from 0.9 to 1.5, depending among others on the resistance to water vapor transfer (r s ) of the vegetated surface.In semiarid regions, values of α were assumed substantially higher than the standard value of 1.27.The authors of [21] considered values of α between 1.6 and 1.8, in two shrublands of south Spain due to local conditions of strong advection.Globally, most authors considered α equal 1.27, in their studies in four sites in North America (e.g., [23]); in a semi-arid area of East Asia [5]; and in West African countries [24].According to [25], the main factors related with α are the atmospheric vapor pressure deficit, soil moisture content, wind run, atmospheric stability condition and the turbulent sensible heat flux.The α value could present systematic variations with the time of the day and season of the year [22,26], and with the soil moisture content [27,28].To the best of our knowledge, a sensitivity analysis of the method results to α value selection has not been yet performed.
A second shortcoming is the inherent uncertainty in the determination of the dry and wet edges, deriving from the more or less arbitrary selection of the geometrical shape and limits.With the aim to reduce the dose of subjectivity and ambiguity in the application of the method, [20] analyzed the impact of end-member selection on the performance and mechanism for error propagation in the spatial variability of fluxes, highlighting that the choice of the shape lead to significant distortions in the spatial distributions of the surface fluxes.
One of the main objectives of the present work is to validate a spatio-temporal distributed method of crop evapotranspiration retrieval from remote sensing, considering a minimum quantity of ancillary data.The selected method for ET assessment corresponds to a graphical approach based on the Priestley-Taylor equation, where the dry soil and non-limiting water conditions are defined from the space conformed by land surface temperature (LST) and vegetation index (VI) both retrieved from remote sensing.
In the graphical approach, we have to fix the Priestley-Taylor coefficient (α).The α parameter corresponds to well-watered conditions, and it could differ among crops/vegetation types because of a differentiated response to climatic stress.Citrus crops are well-known for their strong stomatal closure in response to vapor pressure deficit (VPD), therefore diminishing substantially the advective term of the Penman-Monteith (PM) equation (which includes the aerodynamic resistance and VPD), even when the crop is well-watered.As Equation ( 2) is a simplified form of the PM equation, this suggests a value of α lower than 1.27 for crops such as citrus which exert a strong stomatal control on water loss.Therefore, with this focus, the assessment of sensitivity to the Priestley-Taylor coefficient is a crucial objective to derive accurate spatial distributions of ET.The evaluation of influence of the LST-IV space edges on the ET spatial distribution is the other main objective of the work.To reach these aims, several Landsat 5 TM images were considered for the period 2009-2011, as well as ET tower flux measurements as ground truth in an agricultural irrigated region of southeast of Spain.

Study Area
The study area is located in the Segura River Basin, southeast of Spain (Figure 1a), with an extension approximately to half of one satellite Landsat image (199-34 scene).Nevertheless, our analysis focuses on the irrigation scheme of Campo de Cartagena (CRCC).Figure 1a corresponds to a false-color mosaic of Landsat images for the study area (red color for vegetated surfaces, and intense red color for irrigated lands).Figure 1b presents the location of the two citrus farms (SPOT image, mosaic 2012), where surface fluxes were monitored by means of the eddy-covariance technique (EG flux tower, Figure 1c).
The Segura River Basin is subject to a Mediterranean climate with strong contrasts, frequent droughts, torrential rainfall and high air temperatures in summer [29].The extreme heat, insolation, and irregularity of precipitation resemble a sub-desert climate due to the prolonged anticyclonic situation, which gives a stable hot and dry sunny weather.It can be considered that the climate of the study area is semiarid, but somewhat tempered by the presence of irrigated crops.
The Campo de Cartagena basin is a flat region with a gentle slope of around 1%.It is surrounded by small mountain ranges, and in the east it is open to the Mediterranean Sea through a hypersaline coastal lagoon, the Mar Menor.The main climatic characteristics correspond to a mean annual temperature of 18 • C, an average annual precipitation of 300 mm (distributed into a few intensive events which take place mainly in spring and autumn seasons), and an average potential evapotranspiration of 1275 mm/year [30].The Campo de Cartagena basin is a flat region with a gentle slope of around 1%.It is surrounded by small mountain ranges, and in the east it is open to the Mediterranean Sea through a hypersaline coastal lagoon, the Mar Menor.The main climatic characteristics correspond to a mean annual temperature of 18 °C, an average annual precipitation of 300 mm (distributed into a few intensive events which take place mainly in spring and autumn seasons), and an average potential evapotranspiration of 1275 mm/year [30].
Land use is dominated by agriculture (see Figure 2).Water requirements are satisfied by groundwater pumped from aquifers, surface water resources from the inter-basin Tajo-Segura Water Transfer (since 1979 provides more than one-third of the total water demand), with water from seawater desalination plants and brackish groundwater [30].

Ground Data
To validate the estimates of the ET remote sensing retrieval algorithm, in situ measurements of λET, H, Rn and G were performed simultaneously in the two commercial citrus orchards located in the CRCC.The farms (Figure 1b) correspond to:

•
Farm A (37 ha) with 6-year old orange trees (Citrus sinensis var.Navelate) planted at a tree spacing of 4.5 m × 3 m, and Leaf Area Index (LAI) ≈ 3.This value corresponds to field LAI.

•
Farm B (16 ha) with 30-year old orange trees (Citrus sinensis var.Navelate) planted at a tree spacing of 6 × 4 m, and high LAI ≈ 4 to 5. The field LAI was higher at Farm B than at Farm A, because of the large crown area of the adult trees, compared to the much smaller crown area of the young trees.Therefore, fraction cover was about 0.70-0.75 at Farm B and 0.35-0.45 at Farm A. Land use is dominated by agriculture (see Figure 2).Water requirements are satisfied by groundwater pumped from aquifers, surface water resources from the inter-basin Tajo-Segura Water Transfer (since 1979 provides more than one-third of the total water demand), with water from seawater desalination plants and brackish groundwater [30].

Ground Data
To validate the estimates of the ET remote sensing retrieval algorithm, in situ measurements of λET, H, R n and G were performed simultaneously in the two commercial citrus orchards located in the CRCC.The farms (Figure 1b) correspond to:

•
Farm A (37 ha) with 6-year old orange trees (Citrus sinensis var.Navelate) planted at a tree spacing of 4.5 m × 3 m, and Leaf Area Index (LAI) ≈ 3.This value corresponds to field LAI.

•
Farm B (16 ha) with 30-year old orange trees (Citrus sinensis var.Navelate) planted at a tree spacing of 6 × 4 m, and high LAI ≈ 4 to 5. The field LAI was higher at Farm B than at Farm A, because of the large crown area of the adult trees, compared to the much smaller crown area of the young trees.Therefore, fraction cover was about 0.70-0.75 at Farm B and 0.35-0.45 at Farm A. The two orchards were irrigated to satisfy 100% of the standard crop evapotranspiration, ETc, throughout the whole year.ETc was estimated as the product of reference evapotranspiration, ETo, and the crop coefficient for orange trees from the FAO Penman-Monteith method [31].
The monitored variables (Figure 1c) were: The lack of closure in the surface energy balance equation was checked by comparing for each 30 min measurement the difference between available energy (A = Rn − G) and the sum of the EC fluxes (λE + H).On a daily average, the sum was found 15% to 20% lower than A, the difference being small in the morning (~5% to 10% at overpass time) and somewhat higher (up to 20-25%) in the mid afternoon.To account for the lack of closure, we applied the correction proposed by [32], assuming that the Bowen ratio was kept unchanged.
A simplified footprint analysis [33] was used to estimate the relative contribution of a given area within the plot to the total measured flux.The analysis indicated that more than 90% of the measured λE came from an upwind area less than 160 m away from the measurement tower.As the minimum fetch was close to 200 m (south sector), we did not apply footprint corrections according to the wind direction.The two orchards were irrigated to satisfy 100% of the standard crop evapotranspiration, ETc, throughout the whole year.ETc was estimated as the product of reference evapotranspiration, ETo, and the crop coefficient for orange trees from the FAO Penman-Monteith method [31].
The monitored variables (Figure 1c) were: • at 1.5 m above the tree crowns air temperature (T air ) and relative humidity, solar radiation (RG, Kipp & Zonen pyranometer CMP3, Delft, the Netherlands) and crown temperature (Tc, Apogee infrared thermometer (IRT), Logan, UT, USA); The lack of closure in the surface energy balance equation was checked by comparing for each 30 min measurement the difference between available energy (A = R n − G) and the sum of the EC fluxes (λE + H).On a daily average, the sum was found 15% to 20% lower than A, the difference being small in the morning (~5% to 10% at overpass time) and somewhat higher (up to 20-25%) in the mid afternoon.To account for the lack of closure, we applied the correction proposed by [32], assuming that the Bowen ratio was kept unchanged.
A simplified footprint analysis [33] was used to estimate the relative contribution of a given area within the plot to the total measured flux.The analysis indicated that more than 90% of the measured λE came from an upwind area less than 160 m away from the measurement tower.As the minimum fetch was close to 200 m (south sector), we did not apply footprint corrections according to the wind direction.

RS Data
Landsat 5 TM images present an adequate spatial resolution to address water management issues at farm-scale.The images were downloaded from the server of the Spanish National Remote Sensing Programme (PNT), coordinated by the National Geographic Institute (IGN) with re-sample of 25 m by cubic convolution method (or nearest neighbor algorithm in some cases) [34].A total of 10 images from September 2009 to June 2011 were processed and used to validate with ground truth.The images were geometrically corrected to the ETRS-89 geodetic reference system.The Landsat 5 TM images present six multispectral bands with a spatial resolution of 30 m, with the exception of thermal infrared band (120 m), and spatial resolution of 16 days.The corresponding images of atmospheric water vapor content, provided by the product MOD05_L2 from the MODIS Terra platform, were used to correct the Landsat Thermal-Infrared data from atmospheric effects.Table 1 presents the dates and seasons of the used images, and the irrigation conditions.

Land Cover
A detailed spatial distribution of land cover was used to eliminate outliers.The data were provided by the Spanish Land Cover Information System (SIOSE system) at a scale of 1:25,000, which categorizes the Earth surface according to their biophysical properties (e.g., urban use, crops, forest, etc.) and characterizes the territory (e.g., industrial use, commercial use, etc.).
As it is presented in Figure 2, agriculture is the main land use, in particular irrigated intensive agriculture.The main irrigated crops are horticultural crops (lettuce, broccoli, melon, and others) and citrus trees (oranges and lemons), where the drip is the primary irrigation method (90%) due to water scarcity and the requirement of water conservation.On the other hand, the most representative rainfed crops (covers only ≈ 6%) are almond, winter cereals, and olive.
A method to filter the SIOSE data was applied to identify the areas without vegetation (e.g., beaches, bare soils, fired areas, etc.), artificial covers (e.g., parking areas, urbanized areas, etc.), and water (lakes, sea, etc.).We decide to filter the water category for these reasons: (i) water bodies present NDVI <0, in our algorithm we look for the wet edge for NDVI >0.5; and (ii) in the studied semi-arid area, the water bodies are very small reservoir only used to store water for irrigation.
The authors of [35] indicated that tuning end-members for land cover could minimize error between observed and estimated values.

Other Data
Additional processed information were:

•
Air temperature data (T air ) from meteorological stations (maximum and minimum daily) to assess T air at satellite overpass time (10:30 a.m.).A multivariate regression technique based on meteorological data and spatial covariates was applied to the spatial distribution of maximal and minimal temperatures according to the methodology proposed by [36,37].The explanatory variables were elevation, latitude, longitude, and the distance to the sea.Then, the T air at the satellite time overpass was obtained by a model proposed by [38] in which its diurnal rhythm is given by a sinusoidal progression during daytime, and a decreasing exponential curve during the night.

•
DEM and spatial distribution of monthly average of Linke turbidity [39] generated in the SoDa project [40] to estimate downward longwave radiation (S) at satellite overpass time according to a solar radiation model proposed by [41], integrated into the GIS-Open source [42] which simulates its three components: direct solar radiation, diffuse sky radiation and reflected radiations.Spatial variation of solar radiation due to terrain and terrain-shadowing effects are also included in the model.The required inputs of the model are map elevations and topographic attributes (slope and aspect maps), Linke atmospheric turbidity, surface albedo, day of year, solar declination, and satellite overpass time.
With these variables mentioned above, we attempt to model the different meteorological conditions over such large area.

Radiometric and Topographical Corrections
The radiometric, optical, topographical, emissivity, and thermal-infrared atmospheric corrections of Landsat images were performed following the specifications in [34].Cloud masks were also applied.
From these corrected data, the values of the inputs required for the ET-retrieval algorithm were derived: LST, NDVI, vegetation fraction cover (f c ), albedo and emissivity (ε).For clarification, Scheme 1 presents a simplified diagram of implemented processes.

Other Data
Additional processed information were: • Air temperature data (Tair) from meteorological stations (maximum and minimum daily) to assess Tair at satellite overpass time (10:30 a.m.).A multivariate regression technique based on meteorological data and spatial covariates was applied to the spatial distribution of maximal and minimal temperatures according to the methodology proposed by [36,37].The explanatory variables were elevation, latitude, longitude, and the distance to the sea.Then, the Tair at the satellite time overpass was obtained by a model proposed by [38] in which its diurnal rhythm is given by a sinusoidal progression during daytime, and a decreasing exponential curve during the night.

•
DEM and spatial distribution of monthly average of Linke turbidity [39] generated in the SoDa project [40] to estimate downward longwave radiation (S) at satellite overpass time according to a solar radiation model proposed by [41], integrated into the GIS-Open source [42] which simulates its three components: direct solar radiation, diffuse sky radiation and reflected radiations.Spatial variation of solar radiation due to terrain and terrain-shadowing effects are also included in the model.The required inputs of the model are map elevations and topographic attributes (slope and aspect maps), Linke atmospheric turbidity, surface albedo, day of year, solar declination, and satellite overpass time.
With these variables mentioned above, we attempt to model the different meteorological conditions over such large area.

Radiometric and Topographical Corrections
The radiometric, optical, topographical, emissivity, and thermal-infrared atmospheric corrections of Landsat images were performed following the specifications in [34].Cloud masks were also applied.
From these corrected data, the values of the inputs required for the ET-retrieval algorithm were derived: LST, NDVI, vegetation fraction cover (fc), albedo and emissivity (ε).For clarification, Scheme 1 presents a simplified diagram of implemented processes.
Scheme 1. Scheme of implemented processes.Scheme 1. Scheme of implemented processes.

ET-Retrieval Algorithm Estimation of Biophysical Variables
The surface albedo was estimated as follows [43]: where ρ TM i (i = 1 to 7) are the reflectance values in the bands 1 to 7 of Landsat TM.NDVI [44] is a well-known indicator, strongly related to the fraction of photosynthetically active radiation (fPAR), and hence is closely associated with vegetation activity or greenness [45,46].
The estimation of the f c is necessary for the assessment of other variables involved in the images correction, such as emissivity and thermal-infrared atmospheric corrections.We assumed that f c was related to NDVI following the relationship [47]: where NDVI s and NDVI v refer to soil and vegetation NDVI, respectively.Following [48], we fixed values of NDVI v = 0.5 and NDVI s = 0.2, according to the technical document for processing high resolution images of PNT [34] Surface emissivity was calculated by considering NDVI in three different cases [49]: where ρ TM 3 is the band 3 reflectance of Landsat (visible-red), and f c is the vegetation fraction cover, following the specifications of [34].
The algorithm used is a revision single-channel (SC) algorithm developed by [50], which was particularized to the thermal-infrared (TIR) channel (band 6) to Landsat 5 TM sensor where L sen is the spectral radiance at the sensor's aperture (W•m −2 •sr −1 •µm −1 ), ε is the surface emissivity; γ and δ are two parameters dependent on the Planck's function; and ψ 1 , ψ 2 , ψ 3 are atmospheric momentum transport coefficients.Coefficients to estimate the atmospheric functions are identified from the database of TIGR61 [50].The product MOD05_L2 provided by the TERRA-MODIS sensor is used to estimate the atmospheric water vapor content.

Interpretation of LST vs. NDVI Space
As suggested by [51], we used the difference DT = LST − T air instead of the absolute value of LST.The geometrical envelope of the pixels in the DT-NDVI space is generally considered to fall within the three following cases (Scheme 2):

•
Case 1: Triangular shape, as it was used by [17,18]; • Case 2: Trapezoidal shape (e.g., [21,52]); and • Case 3: Rectangular shape (e.g., [53]).For all images of the CRCC area, the DT-NDVI space was far from matching a triangular shape.Therefore, the automatically detection of end members (dry and wet edges) was performed considering only the trapezoidal and rectangular shapes.

ET Estimation
The graphical method based on an interpolation of the Priestley-Taylor formula (Equation ( 2)), was used to assess the spatio-temporal distribution of ET.The pixels of the DT-NDVI space were assumed to have ET equal to: with φ given by: where max φ = α (fully wet surface).According to [17], the maximum and minimum values of DT (DTmax and DTmin) corresponded to φ = 0 and φ = max φ , respectively; and DT values were assumed to scale linearly with the evaporative fraction, EF.In Scheme 2 for a pixel "i" with an observed DTobs, φ was obtained by the ratio between "l" (DTmax − DTobs) and "n" (DTmax − DTmin).Equation ( 9) implies that the evaporative fraction is equal to: For all images of the CRCC area, the DT-NDVI space was far from matching a triangular shape.Therefore, the automatically detection of end members (dry and wet edges) was performed considering only the trapezoidal and rectangular shapes.

ET Estimation
The graphical method based on an interpolation of the Priestley-Taylor formula (Equation ( 2)), was used to assess the spatio-temporal distribution of ET.The pixels of the DT-NDVI space were assumed to have ET equal to: with ϕ given by: where ϕ max = α (fully wet surface).According to [17], the maximum and minimum values of DT (DT max and DT min ) corresponded to ϕ = 0 and ϕ = ϕ max , respectively; and DT values were assumed to scale linearly with the evaporative fraction, EF.In Scheme 2 for a pixel "i" with an observed DT obs , ϕ was obtained by the ratio between "l" (DT max − DT obs ) and "n" (DT max − DT min ).Equation ( 9) implies that the evaporative fraction is equal to: Once ϕ, and thus EF (Equation ( 3)) was obtained, the next step was the estimation of the available energy A = R n − G, at each overpass time.The net radiation at the surface was derived from the following equation [54]: where S is downward shortwave radiation, a is surface albedo, and L d and L up are downward and upward longwave radiation, respectively.The net longwave radiation L n = L d − L up was calculated as [55]: where σ is Stefan-Boltzmann constant (5.67 × 10 −8 Wm −2 •K −4 ), T air (Kelvin) at screen level, ε s is surface emissivity, and ε a is the apparent emissivity of the sky, calculated from the empirical equation proposed by [56].
Ground Heat Flux The usual approach to estimate G is to consider that as a constant fraction of R n (β = G/R n ) [57].In the present work, G was obtained from a direct relationship linking β to the evaporative fraction (EF) applying the approach proposed by [24].

Automatic Determination of End-Members
In the present work, a modified method of the original one proposed by [5] was applied to automatically identify the wet and dry edges.The edges were identified by filtering out the spurious pixels, and applying an iterative method which minimizes the root mean square error (RMSE).An automatic algorithm was implemented to identify the maximum and minimum subset pixels to identify the end-members.In summary, the algorithm presents the following steps: (i) The range of NDVI is divided in M intervals (M ≤ 20) and each interval into N subintervals (N ≥ 5).(ii) The DT average (DTaver) and standard deviation (σ) of each subinterval is calculated as an initial state.Then, for each subinterval the maximum and minimum DT is identified.(iii) The subinterval is omitted if DT is less than (DTaver-σ), recalculating the new average value and σ. (iv) DTaver is taken as the maximum or minimum DT for this given interval.(v) A linear regression of DT and NDVI values is performed, assessing the Root Mean Square Error (RMSE).The interval will be omitted if DT is less than DTc (calculated with coefficients a' + b' obtained from the linear regression) minus 2 times RMSE.(vi) A subset of average DT (end-members maximum o minimum) for each interval is obtained.
Once the subset of end-members were identified, the wet and dry edges of the trapezoidal and rectangular domains were defined applying the following methodology:

•
Pixels with a dominant soil fraction, NDVI < 0.3 were omitted to define the dry edge.In the case of wet edge, only pixels with complete vegetation cover values (NDVI > 0.5) were considered [21], as presented in Scheme 2.

•
In the case of trapezoidal shapes, the dry edge is defined by linear adjustment to maximum end-members, and the wet edge is considered a constant value equal to average of minimum end-members (wet edge and dry edge 2 in Scheme 2).

•
For rectangular shapes, the dry edge is estimated as a constant equal to higher value of maximum of maximum subset, with NDVI > 0.3.The dry edge is equally considered as a constant equal to the lower value of minimum end-members with NDVI >0.5 (wet edge and dry edge 3 in Scheme 2).
where E i and M i are paired model estimates/predictions and measured/observed variables, respectively, while n is the number of data points, and M is the mean value of the observed variable.

Predicted vs. Ground Truth Values
Table 2 presents the values of NDVI and f c assessed from satellite images for the ten overpass dates, for the pixel of ET tower flux measurement in Farms A and B, respectively.Higher values of NDVI were observed in Farm B in comparison with Farm A, a logical result considering the greater LAI of Farm B (LAI ≈ 5) with respect to Farm A (LAI ≈ 3).If the total surface of the farms is considered.The mean value of NDVI ranged from 0.48 to 0.57 in Farm B, and from 0.30 to 0.50 in Farm A. The comparison between estimated and observed values of LST and R n at overpass time indicated that the retrieval method for these variables provided accurate estimates (Figure 3a,b).In particular, R n , which is the basic variable required for the quantification of ET (Equation ( 9)) was estimated with a MAPE value of 15% (Figure 3b), a rather satisfactory result.According to [61], the slight underestimation of LST values could be justified by the procedure for the estimation of emissivity from NDVI.
The comparison between estimated and observed values of LST and Rn at overpass time indicated that the retrieval method for these variables provided accurate estimates (Figure 3a,b).In particular, Rn, which is the basic variable required for the quantification of ET (Equation ( 9)) was estimated with a MAPE value of 15% (Figure 3b), a rather satisfactory result.According to [61], the slight underestimation of LST values could be justified by the procedure for the estimation of emissivity from NDVI.

Sensitivity of EF to α
In order to identify the proper α value for the study area, we searched for the value of α which minimizes the RMSE between the observed and estimated values of EF.It was done considering the sensitivity to DTmin and DTmax.The actual fully wet conditions are observed due to irrigation, while DTmax represents the areas at "barely" fully stress.The edges were identified from the error of EF to α range (0.7, 1.3) selection, with an increment of 0.1 for both shapes-trapezoidal (TR) and rectangular (RC).The results of the analysis are presented in Figure 4.The results were similar for both farms, indicating that RMSE was minimized for α close to 1.1 in the case of the TR shape, and to 1.0 for the RC shape.
In Figure 4a, the wet and dry edges are estimated for the whole range of NDVI (by all DTmax and DTmin values).In Figure 4b, the DTmin values are considered for NDVI > 0.5 (from Scheme 2, it corresponds to wet edge 2 TR shape or wet edge 3 RC shape), and the dry edge is estimated across the whole range of NDVI (all DTmax values are considered) Then, in Figure 4c, DTmax values are omitted with dominant soil fraction NDVI < 0.3 (from Scheme 2, it corresponds to dry edge 2 TR shape and dry edge 3 RC shape), and the wet edge is estimated across the whole range of NDVI (all DTmin values are considered).Finally, Figure 4d represents the situation when the restrictions defined in Scheme 2 are considered.
Considering all range of NDVI values, the smallest errors are close to 1.1-1.2. Figure 4c show a light dismiss compared to Figure 4a, with smaller errors close to 1-1.1.Nevertheless, the more relevant restriction is related to the wet edge (Figure 4b), showing values close to 1. Finally, the best results indicated that RMSE was reduced for α close to 1.1 in the case of the TR shape, and to 1.0 for RC shape (Figure 4d).
We consider the latest approach over others because we obtain the best error, in concordance with [21].On the other hand, factors such as the spatial resolution and the small size of farms lead us to seek more reliability in our results omitting pixels contaminated by other land covers.

Sensitivity of EF to α
In order to identify the proper α value for the study area, we searched for the value of α which minimizes the RMSE between the observed and estimated values of EF.It was done considering the sensitivity to DT min and DT max .The actual fully wet conditions are observed due to irrigation, while DT max represents the areas at "barely" fully stress.The edges were identified from the error of EF to α range (0.7, 1.3) selection, with an increment of 0.1 for both shapes-trapezoidal (TR) and rectangular (RC).The results of the analysis are presented in Figure 4.The results were similar for both farms, indicating that RMSE was minimized for α close to 1.1 in the case of the TR shape, and to 1.0 for the RC shape.
In Figure 4a, the wet and dry edges are estimated for the whole range of NDVI (by all DT max and DTmin values).In Figure 4b, the DT min values are considered for NDVI > 0.5 (from Scheme 2, it corresponds to wet edge 2 TR shape or wet edge 3 RC shape), and the dry edge is estimated across the whole range of NDVI (all DT max values are considered) Then, in Figure 4c, DT max values are omitted with dominant soil fraction NDVI < 0.3 (from Scheme 2, it corresponds to dry edge 2 TR shape and dry edge 3 RC shape), and the wet edge is estimated across the whole range of NDVI (all DT min values are considered).Finally, Figure 4d represents the situation when the restrictions defined in Scheme 2 are considered.
Considering all range of NDVI values, the smallest errors are close to 1.1-1.2. Figure 4c show a light dismiss compared to Figure 4a, with smaller errors close to 1-1.1.Nevertheless, the more relevant restriction is related to the wet edge (Figure 4b), showing values close to 1. Finally, the best results indicated that RMSE was reduced for α close to 1.1 in the case of the TR shape, and to 1.0 for RC shape (Figure 4d).
We consider the latest approach over others because we obtain the best error, in concordance with [21].On the other hand, factors such as the spatial resolution and the small size of farms lead us to seek more reliability in our results omitting pixels contaminated by other land covers.The boxplots of EF obtained with α = 1 and α = 1.3 in the case of both trapezoidal and rectangular shapes indicate that the mean value of EF for both shapes were rather close (Figure 5a,b).The reason for the presence of clear outliers could be that, even though clouds and cirrus were filtered, the imperfect cleaning of shadows could produce anomalous values of reflectances.When comparing the boxplots of the two shapes, a lower number of outliers for the rectangular shape (Figure 5b) than for the trapezoidal one (Figure 5a) can be observed.The boxplots of EF obtained with α = 1 and α = 1.3 in the case of both trapezoidal and rectangular shapes indicate that the mean value of EF for both shapes were rather close (Figure 5a,b).The reason for the presence of clear outliers could be that, even though clouds and cirrus were filtered, the imperfect cleaning of shadows could produce anomalous values of reflectances.When comparing the boxplots of the two shapes, a lower number of outliers for the rectangular shape (Figure 5b) than for the trapezoidal one (Figure 5a) can be observed.

Sensitivity of EF Spatial Pattern to Space Shape
Table 3 presents the edges identified by the developed algorithm for trapezoidal and rectangular shapes, considering DT and LST to define the edges.A decrease in the shape of approximately 20 • C (between LST and DT) is identified, especially in summer months, due to the effect of maximum temperature in these months.However, for winter months, the decrease is less relevant.To evaluate the influence of T air in the results, the errors were quantified for each method (Table 4) considering LST and DT, respectively.In Table 4, considering α = 1.27, the errors were reduced around 20% with the inclusion of T air .These results accord with those obtained by [51].Analyzing the edges defined by DT, only one image has a similar shape for both configurations, with a slope of 10%, and the trapezoidal shape is more close to a rectangle shape.The other images have slopes of 15-35% where the shape is more close to a trapezoid (degenerated triangle).Finally, only one image can be considered as a triangular shape with a slope of approximately 50%.In most cases, degenerated triangles, such as was defined by [20], are observed.
As noted above, α close to 1 was better performance; with this in mind, we analyze EF histograms, represented in Figure 6.These were more concentrated for the trapezoidal shape in comparison with the rectangular shape.In general, the rectangular shape (Figure 6b,d,f,h) presents higher CV in comparison with trapezoidal shape (Figure 6a,c,e,g).The relative differences between the spatial distributions of EF estimated for trapezoidal and rectangular shape were calculated as: (EF TR − EF RC )/EF TR .For the whole irrigated area, the relative differences were not higher than 20% (Figure 7).
The relative differences between the spatial distributions of EF estimated for trapezoidal and rectangular shape were calculated as: (EFTR − EFRC)/EFTR.For the whole irrigated area, the relative differences were not higher than 20% (Figure 7).

Results of ET
Finally, we compared the values of ET estimated for α = 1.0 and α = 1.3 to the observed ET (Table 5).The value of 1.3 led to a systematic overestimation of ET (Figure 8a,b), while α = 1 provided a much better fit (Figure 8c,d).We also found that the trapezoidal shape globally provided smaller errors (MAE, MAPE and RMSE) in comparison with the rectangular shape, and that ET estimated assuming a trapezoidal shape was less sensitive to the selection of the end members.

Results of ET
Finally, we compared the values of ET estimated for α = 1.0 and α = 1.3 to the observed ET (Table 5).The value of 1.3 led to a systematic overestimation of ET (Figure 8a,b), while α = 1 provided a much better fit (Figure 8c,d).We also found that the trapezoidal shape globally provided smaller errors (MAE, MAPE and RMSE) in comparison with the rectangular shape, and that ET estimated assuming a trapezoidal shape was less sensitive to the selection of the end members.In the case of comparing EF for both farms, the variations are small because these are well irrigated crops.The low correlation with respect to the observed data, are justified because the method is not precise enough to discriminate small variations of EF.For both farms, Table 6 presents the EF values for each α selected value, and TR and RC shapes.The identified pattern in EF from Figure 5 is reproduced by the results of Table 6 for α = 1 in comparison with α = 1.3: a decrease of the mean, and extreme values as well as their variation.
Removing two winter days with values of observed EF (EFobs) too high with respect to the remaining values-due to high experimental errors when evaporation and available energy were low-the estimated values of EF (EFest) were in relatively good agreement with the observed ones (Figure 9): for Farm A, the mean values for EFest and EFobs were 0.46 and 0.45, respectively, and for Farm B, they were 0.49 and 0.44, respectively, with a relative RMSE of 15% in both cases, for α = 1.In the case of comparing EF for both farms, the variations are small because these are well irrigated crops.The low correlation with respect to the observed data, are justified because the method is not precise enough to discriminate small variations of EF.For both farms, Table 6 presents the EF values for each α selected value, and TR and RC shapes.The identified pattern in EF from Figure 5 is reproduced by the results of Table 6 for α = 1 in comparison with α = 1.3: a decrease of the mean, and extreme values as well as their variation.
Removing two winter days with values of observed EF (EF obs ) too high with respect to the remaining values-due to high experimental errors when evaporation and available energy were low-the estimated values of EF (EF est ) were in relatively good agreement with the observed ones (Figure 9): for Farm A, the mean values for EF est and EF obs were 0.46 and 0.45, respectively, and for Farm B, they were 0.49 and 0.44, respectively, with a relative RMSE of 15% in both cases, for α = 1.

Differences Induced by Space Shape
In the methods based on the analysis of the LST-VI space, the shape selection plays a critical role in defining the position of the end-members, and therefore on the resultant value of ϕ, and consequently EF and ET [20].For this reason, our algorithm uses a statistic method based on iterative process to avoid a misinterpretation of maximum and minimum values (it eliminates subjectivities); in addition, we try to correct the overestimation of DT min o DT max taking into account the average values in each interval of NDVI such as it was explained in Scheme 2. However, it is necessary to take into account certain uncertainties in DT due to several sources of error: (i) LST bias [62], i.e., due to a lack of correction of LST data by elevation [61], or by a noise component [63]; (ii) errors in the interpolation of maximum and minimum temperature [35,36,64] and their residuals distribution [63]; (iii) due to the assessment method of instantaneous temperature; and (iv) DT bias due in relation with the suitability of image extension and resolution [20,21].
Nevertheless, as shown in Section 3.3, the obtained differences in our case of study are mainly induced by the dry edge (dry edge 2 or dry edge 3 in Scheme 2), which is parallel to the NDVI axis in the case of the rectangular shape, and decreased with NDVI in the trapezoidal shape.However, all spatial patterns are degenerate triangles such as it was defined by [20] with the following differences: • a lower number of outliers with the rectangular shape, because this shape allows to include within its bounds a larger number of "valid" pixels than the trapezoidal shape; and • for the same reason, the coefficient of variation, considering only the valid pixels, was higher in the case of the rectangular shape.
Such differences due to the shape did not appear to affect to a large extent the resulting values of EF and ET, as the values of RMSE and MAE were rather similar for both shapes, although the trapezoidal shape was better adapted in the case of Farm A (Table 2) and therefore could be preferred.The authors of [20] also found that the two shapes provided similar accuracy, advising for the rectangular form whenever there was no clear indication that the space presents a trapezoidal trend.According to [65], the rectangular form could be appropriate for large areas and coarse spatial resolutions.

Uncertainty Due to R n
According to [66], the temperature and the influence of radiation are more important than the scarcity water, being the energy available that controls the surface evaporation [61].In the balance of R n , three variables have important influence: albedo, emissivity, and LST.For this reason, it is important to highlight the bias in R n due: (i) LST and hence emissivity [54,61]; and (ii) the use of surface reflectance instead to calculate surface albedo according to [67].Nevertheless, the spatial resolution of LST (120 m) controls the accuracy of this method and hence the errors are balanced, so we obtain a good correlation of R n , and available energy between observed and estimated datasets (Figure 10).

Uncertainty Due to the Value of α
The value chosen for the parameter α has an important influence on the value of EF, hence of ET.According to our results, a value of α = 1 resulted in the lowest errors for EF, and therefore provided the best ET estimates.This may be explained by the fact that most of the wet edge of the LST-NDVI space correspond to irrigated crops, which in the area are citrus and vegetable crops with a fairly high surface resistance to water vapor transfer (rs).The authors of [22] clearly demonstrated

Uncertainty Due to the Value of α
The value chosen for the parameter α has an important influence on the value of EF, hence of ET.According to our results, a value of α = 1 resulted in the lowest errors for EF, and therefore provided the best ET estimates.This may be explained by the fact that most of the wet edge of the LST-NDVI space correspond to irrigated crops, which in the area are citrus and vegetable crops with a fairly high surface resistance to water vapor transfer (r s ).The authors of [22] clearly demonstrated that the control of evaporation in vegetated lands depends strongly on the physiological characteristics of the vegetation, in particular of the value of the surface resistance, r s .He showed that the α value at overpass time of the satellite (10:00-11:00 a.m.) is close to 1 when r s is in the range of 150 to 250 s/m, which is the order of magnitude of r s for citrus orchards similar to those growing in the studied area [68].
Obviously, such a difference of approximately 25% in the value of α resulted in a similar order of error in the quantitative estimation of EF and ET (Figure 5a-d).Such a substantial error might cast doubt on the overall robustness of the method, unless a reasonable estimate of α could be made for the area under study.A solution could be to collect information-for instance from land use maps-of the type of crop or vegetation mostly growing in the area, and of their mean surface resistance.From the knowledge of the mean r s , a mean value of α could be derived [22] and used as a more realistic value of the PT coefficient for ET retrieval.

Overall Accuracy
The results of the validation were satisfactory, with a mean RMSE of ~30 W/m 2 for the two selected shapes of the LST-NDVI space.The overall accuracy of our method was also similar to that obtained in more recent studies, such as [21], or [20,52].In particular, the latter authors stressed that the error committed in the end-members selection could be high because of outliers.In our method, we intended to minimize the number of outliers through two types of pretreatments:

•
In a first step, filtering by land use types, eliminating by this way possible conflicts between land use class and EF values [35] and avoiding erroneous positioning of the wet and dry edges • In a second step, considering that the dry edge was determined by all pixels with NDVI > 0.3, and the wet edge by all pixels with NDVI > 0.5 [21].
In addition, two other elements might have resulted beneficial to the overall precision

•
The determination of G, by means of a robust formulation in function of EF, which implicitly includes the effect of soil moisture and soil properties on soil heat flux [24].

•
As it was demonstrated by other authors, the consideration of air temperature variation at overpass satellite time [69,70] and the correction of by air temperature of [51], reduce the errors in 15-30%.

Conclusions
The present study underlines the importance of a more precise knowledge of the actual value of the PT coefficient when using ET retrieval methods based on the LST-VI space.Both ET tower flux measurements and Landsat 5 TM images of an irrigation scheme in the southeast Spain, were considered.
The results of the sensitivity analysis of ET retrieval from remote sensing with respect to the shape (trapezoidal or rectangular) of the LST-VI space, does not appear to affect substantially the overall accuracy of the method.However, the relevant role of the variation of NDVI lead us to prefer TR shape which presents less variability, hence more reliability than RC shape.
More critical is the choice of the value of the PT coefficient, which is likely to vary in a large range around the standard value of 1.27.The α value has an important influence on the value of EF, hence of ET.According to our results, a value of α = 1 resulted in the lowest errors for EF, and therefore provided the best ET estimates.This may be explained by the fact that most of the wet edge of the LST-NDVI space correspond to irrigated crops, which in the area are citrus and vegetable crops with a fairly high surface resistance to water vapor transfer (r s ).The selected α value is the same for the two field and the rest of citrus and vegetable crops with fairly high r s at overpass time of the satellite (10:00-11:00 a.m.)The irrigated areas in such an arid region (with a precipitation of 250 to 300 mm per year in the Campo de Cartagena) are those which present the maximum evaporation rate and therefore the lowest land surface temperature.
The good results in ET retrieval from remote sensing, is mainly justified by the high correlation between observed and estimated R n datasets, and hence, available energy.However, the method is not enough precise to discriminate small variations of EF, as it was demonstrated with evidences.
In conclusion, the proposed methodology using the PT equation presents several advantages in comparison with the application of PM equation: • two parameters should be fix in PM method (r s and aerodynamic resistance, r a ) instead of one parameter in PT equation; • the order of magnitude of α (=ϕ max ) is more or less identified and known, while much more uncertain are the corresponding values for r s and r a of PM equation; and • the air VPD should be known to calculate ET with the PM equation, therefore complicating the retrieval process as it is very difficult to get correct estimates of the air humidity near the ground from satellite data.
Finally, our study presents a simplified methodology to obtain a more accurate ET retrieval from remote sensing, demonstrating its usefulness in comparison with other conventional methodologies.

• at 1 . 5 m
above the tree crowns air temperature (Tair) and relative humidity, solar radiation (RG, Kipp & Zonen pyranometer CMP3, Delft, the Netherlands) and crown temperature (Tc, Apogee infrared thermometer (IRT), Logan, UT, USA); • Soil heat flux was measured by means of heat flux plates (REBS, model HFT-3.1,Seattle, WA, USA) buried 5 mm below the surface, near to drippers (wet bulbs) and in the middle of the rows (dry soil)]; and • The turbulent energy fluxes, H and λE, were measured at 1.5 m above the trees by an eddy-covariance system comprising a Campbell Scientific Inc. (Logan, UT, USA) CSAT-3 sonic anemometer measuring high-frequency (10 Hz) three-dimensional wind speed and a LICOR (Lincoln, NE, USA) LI-7500 open path infrared gas analyzer measuring CO2 and H2O mixing ratios in absolute mode at 10 Hz.

Scheme 2 .
Scheme 2. Conceptual representation of DT-NDVI relationship.Points represent the area of each shape: (Case 1) grey pixels represent triangular area; (Case 2) green and yellow pixels represent the trapezoidal area; and (Case 3) red and blue pixel represent the rectangular area.

Scheme 2 .
Scheme 2. Conceptual representation of DT-NDVI relationship.Points represent the area of each shape: (Case 1) grey pixels represent triangular area; (Case 2) green and yellow pixels represent the trapezoidal area; and (Case 3) red and blue pixel represent the rectangular area.

Figure 3 .
Figure 3.Comparison of remote sensing estimated variables and ground truth at overpass time satellite: (a) LST; and (b) Rn.

Figure 3 .
Figure 3.Comparison of remote sensing estimated variables and ground truth at overpass time satellite: (a) LST; and (b) Rn.

Figure 4 .
Figure 4. Sensitivity of EF to α range (0.7, 1.3) with an increment of 0.1 for trapezoidal shape (TR) and rectangular shape (RC).Black dots and triangles represent MAPE and RMSE, respectively, in the following cases: (a) across the complete range of NDVI values; (b) considering pixels with NDVI values >0.5 to wet edge (wet edge 2 or 3 in Scheme 2); (c) the complete range of NDVI to dry edge omitting pixel with dominant soli fraction NDVI < 0.3 to dry edge (dry edge 2 or 3 in Scheme 2); and (d) all range of NDVI in wet edge dry and wet edge 2 or 3, depending on pattern to study.

Figure 4 .
Figure 4. Sensitivity of EF to α range (0.7, 1.3) with an increment of 0.1 for trapezoidal shape (TR) and rectangular shape (RC).Black dots and triangles represent MAPE and RMSE, respectively, in the following cases: (a) across the complete range of NDVI values; (b) considering pixels with NDVI values >0.5 to wet edge (wet edge 2 or 3 in Scheme 2); (c) the complete range of NDVI to dry edge omitting pixel with dominant soli fraction NDVI < 0.3 to dry edge (dry edge 2 or 3 in Scheme 2); and (d) all range of NDVI in wet edge dry and wet edge 2 or 3, depending on pattern to study.

Figure 9 .
Figure 9. Contrast of observed EF versus estimated EF, for the two farms.

Figure 9 .
Figure 9. Contrast of observed EF versus estimated EF, for the two farms.

Figure 10 .
Figure 10.Relationship between observed and estimated values of available energy (A = Rn − G) for the trapezoidal space.

Figure 10 .
Figure 10.Relationship between observed and estimated values of available energy (A = R n − G) for the trapezoidal space.

Table 1 .
Details of the used images and irrigation conditions.

Table 2 .
Values of NDVI and f c assessed from the satellite images for the pixel of ET tower flux measurement in Farms A and B (F), respectively.

Table 3 .
Dry and wet edge for trapezoidal (TR) and rectangular (RC) shapes, where a and b are the regression coefficients of the dry edge for LST and DT, respectively, of dry edge, while a' is the regression coefficient for DT of wet edge.

Table 4 .
Influence of T air in the Triangle method for α = 1.27.

Table 5 .
Errors in the assessment of ET for the available satellite images: absolute error (AE), absolute percentage error (APE), mean absolute (MAE), mean absolute percentage error (MAPE) and root mean square error (RMSE).TR = trapezoidal space, RC = rectangular space.

Table 6 .
Influence of α selection in EF values.

Table 5 .
Errors in the assessment of ET for the available satellite images: absolute error (AE), absolute percentage error (APE), mean absolute (MAE), mean absolute percentage error (MAPE) and root mean square error (RMSE).TR = trapezoidal space, RC = rectangular space.

Table 6 .
Influence of α selection in EF values.