Next Article in Journal
The 2013 FLEX—US Airborne Campaign at the Parker Tract Loblolly Pine Plantation in North Carolina, USA
Previous Article in Journal
Remotely Monitoring Ecosystem Water Use Efficiency of Grassland and Cropland in China’s Arid and Semi-Arid Regions with MODIS Data
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Civil Engineering, Universidad Politécnica de Cartagena, Paseo Alfonso XIII, 52, 30203 Cartagena, Spain
2
Department of Food Engineering and Agricultural Equipment, Universidad Politécnica de Cartagena, P. Alfonso XIII, 48, 30203 Cartagena, Spain
*
Author to whom correspondence should be addressed.
Submission received: 20 February 2017 / Revised: 30 May 2017 / Accepted: 10 June 2017 / Published: 16 June 2017

Abstract

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

Graphical Abstract

1. 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].
λE = RnGH,
where the different fluxes, expressed in W·m−2, are net radiation (Rn), 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” (Ew), defined as the evaporation of a well-moist surface under conditions of minimal regional advection [19]:
λ E w = α [ Δ Δ + γ ] ( R n G ) ,
where Δ (Pa·K−1) is the slope of the saturated vapor pressure curve at the prevailing Tair, γ 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:
EF = λE/A
with A = RnG. Equation (3) indicates that, besides EF, the estimation of λE requires the assessment of two other fluxes, Rn 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 (rs) 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.

2. Materials and Methods

2.1. 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].
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].

2.2. 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.
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 (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.
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 = RnG) 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.

2.3. Remote Sensing Data and Other Data

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

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

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

2.4. Data Processing

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

2.4.2. ET-Retrieval Algorithm

Estimation of Biophysical Variables

The surface albedo was estimated as follows [43]:
a = 0.2212 ρ T M 1 + 0.2569 ρ T M 2 + 0.1787 ρ T M 3 + 0.2295 ρ T M 4 + 0.0815 ρ T M 5 + 0.0322 ρ T M 7
where ρ T M 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 fc is necessary for the assessment of other variables involved in the images correction, such as emissivity and thermal-infrared atmospheric corrections. We assumed that fc was related to NDVI following the relationship [47]:
f c = ( N D V I N D V I s N D V I v N D V I s ) 2 i f N D V I < N D V I s f c = 0 i f N D V I > N D V I v f c = 1
where NDVIs and NDVIv refer to soil and vegetation NDVI, respectively. Following [48], we fixed values of NDVIv = 0.5 and NDVIs = 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]:
ε = 0.979 0.035 ρ T M 3 D N V I < 0.2 ε = 0.0968 + 0.004 f c 0.2 N D V I 0.5 ε = 0.99 N D V I > 0.5
where ρTM3 is the band 3 reflectance of Landsat (visible-red), and fc 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
L S T = γ [ ε 1 ( ψ 1 L s e n + ψ 2 ) + ψ 3 ] + δ
where Lsen 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 − Tair 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:
λ E T = φ [ Δ Δ + γ ] ( R n G )
with φ given by:
φ = φ max DT max D T o b s D T max D T min
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” (DTmaxDTobs) and “n” (DTmaxDTmin). Equation (9) implies that the evaporative fraction is equal to:
E F = φ [ Δ Δ + γ ]
Once φ, and thus EF (Equation (3)) was obtained, the next step was the estimation of the available energy A = RnG, at each overpass time. The net radiation at the surface was derived from the following equation [54]:
R n = ( 1 a ) S + L d L u p
where S is downward shortwave radiation, a is surface albedo, and Ld and Lup are downward and upward longwave radiation, respectively.
The net longwave radiation L n = L d L u p was calculated as [55]:
L n = L d L u p = σ ε s ε a T a i r 4 σ ε s L S T
where σ is Stefan–Boltzmann constant (5.67 × 10−8 Wm−2·K−4), Tair (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].
ε a = 9.2 × 10 6 T a i r 2

Ground Heat Flux

The usual approach to estimate G is to consider that as a constant fraction of Rn (β = G/Rn) [57]. In the present work, G was obtained from a direct relationship linking β to the evaporative fraction (EF) applying the approach proposed by [24].
β = a + b E F
with a = 0.23 and b = −0.22.

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

2.5. Performance Analysis

According the recommendations of various authors [35,58,59,60], the root mean square error (RMSE) (Equation (15)), mean absolute error (MAE) (Equation (16)) and mean absolute percent error (MAPE) (Equation (17)) were used to assess the model performance during the evaluation process.
RMSE = i = 1 n ( E i M i ) 2 n
M A E = 1 n i = 1 n | E i M i |
M A P E = 100 M ( 1 n i = 1 n | E i M i | )
where Ei and Mi 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.

3. Results

3.1. Predicted vs. Ground Truth Values

Table 2 presents the values of NDVI and fc 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 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.

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

3.3. 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 Tair 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 Tair. 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: (EFTR − EFRC)/EFTR. For the whole irrigated area, the relative differences were not higher than 20% (Figure 7).

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

4. Discussion of Results

4.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 DTmin o DTmax 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.

4.2. Uncertainty Due to Rn

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 Rn, three variables have important influence: albedo, emissivity, and LST. For this reason, it is important to highlight the bias in Rn 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 Rn, and available energy between observed and estimated datasets (Figure 10).

4.3. 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 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, rs. He showed that the α value at overpass time of the satellite (10:00–11:00 a.m.) is close to 1 when rs is in the range of 150 to 250 s/m, which is the order of magnitude of rs 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 rs, a mean value of α could be derived [22] and used as a more realistic value of the PT coefficient for ET retrieval.

4.4. Overall Accuracy

The results of the validation were satisfactory, with a mean RMSE of ~30 W/m2 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%.

5. 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 (rs). The selected α value is the same for the two field and the rest of citrus and vegetable crops with fairly high rs 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 Rn 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 (rs and aerodynamic resistance, ra) 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 rs and ra 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.

Acknowledgments

The first author acknowledges the several scholarships received in the last years from the projects: “REDSIM Remote-Sensing based DSS for Sustainable Drought-Adapted Irrigation Management” REDSIM and “SIRRIMED Sustainable Use of Irrigation Water in the Mediterranean Region”. In addition, the support from HYDROCLIM (ref. CGL2012-39895-C02-01), funded by National Plan of R&D of National Secretary of Research, Science and Technology of MINECO (Spain) and related with the R&D Group of Water Resources Management-UPCT, is acknowledged.

Author Contributions

Martínez-Pérez, J. and García-Galiano, S. conceived and designed the experiments. Martínez-Pérez, J. performed the experiments. Martín-Gorriz B. collected the data from the flux towers. Martínez-Pérez, J., García Galiano, S. and Baille, A. wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Corwin, D.L.; Plant, R.E. Applications of apparent soil electrical conductivity in precision agriculture. Comput. Electron. Agric. 2005, 46, 1–10. [Google Scholar] [CrossRef]
  2. Hetzroni, A.; Peters, A.; Ben-Gal, A. Towards precision management of orchards: Using automated monitoring to build a GIS-based spatial decision support system. In Proceedings of the International Conference on Agricultural Engineering, Valencia, Spain, 8–12 July 2012. [Google Scholar]
  3. Eberbach, P.; Pala, M. Crop row spacing and its influence on the partitioning of evapotranspiration by winter-grown wheat in Northern Syria. Plant Soil 2005, 268, 195–268. [Google Scholar] [CrossRef]
  4. Falkenmark, M.; Rockstrom, J. Balancing Water for Humans and Nature: The New Approach in Ecohydrology; Earthscan Press: London, UK, 2004; p. 272. [Google Scholar]
  5. Tang, R.; Li, Z.L.; Tang, B. An application of the Ts-VI triangle method with enhanced edges determination for evapotranspiration estimation from MODIS data in arid and semi-arid regions: Implementation and validation. Remote Sens. Environ. 2010, 114, 540–551. [Google Scholar] [CrossRef]
  6. Urrea Mallebrera, M.; Mérida Abril, A.; García Galiano, S.G. Segura River Basin: Spanish Pilot River Basin Regarding Water Scarcity and Droughts. In Agricultural Drought Indices; Sivakumar, M.V.K., Motha, R.P., Wilhite, D.A., Wood, D.A., Eds.; World Meteorological Organization: Geneva, Switzerland, 2010; pp. 2–12. [Google Scholar]
  7. Choudhury, B.J.; Ahmed, N.U.; Idso, S.B.; Reginato, R.J.; Daughtry, C.S. Relations between evaporation coefficients and vegetation indices studied by model simulations. Remote Sens. Environ. 1994, 50, 1–17. [Google Scholar] [CrossRef]
  8. Carlson, T.N.; Capehart, W.J.; Gillies, R.R. A new look at the simplified method for remote sensing of daily evapotranspiration. Remote Sens. Environ. 1995, 54, 161–167. [Google Scholar] [CrossRef]
  9. Carlson, T.N.; Gillies, R.R.; Schmugge, T.J. An interpretation of methodologies for indirect measurement of soil water content. Agric. For. Meteorol. 1995, 77, 191–205. [Google Scholar] [CrossRef]
  10. Bastiaanssen, W.G.M.; Menenti, M.; Feddes, R.A.; Holtslag, A.A.M. A remote sensing surface energy balance algorithm for land. Part I: Formulation. J. Hydrol. 1998, 212–213, 198–212. [Google Scholar] [CrossRef]
  11. Bastiaanssen, W.G.M.; Pelgrum, H.; Wang, J.; Ma, Y.; Moreno, J.F.; Roerink, G.J. A remote sensing surface energy balance algorithm for land (SEBAL) Part 2: Validation. J. Hydrol. 1998, 212–213, 213–229. [Google Scholar] [CrossRef]
  12. Norman, J.M.; Kustas, W.P.; Humes, K.S. Source approach for estimating soil and vegetation energy fluxes from observations of directional radiometric surface temperature. Agric. For. Meteorol. 1995, 77, 263–293. [Google Scholar] [CrossRef]
  13. Kustas, W.P.; Norman, J.; Anderson, M.C.; French, A. Estimating subpixel surface temperatures and energy fluxes from the vegetation index–radiometric temperature relationship. Remote Sens. Environ. 2003, 85, 429–440. [Google Scholar] [CrossRef]
  14. Melesse, A.; Nangia, V. Spatially distributed surface energy flux estimation using remotely-sensed data from agricultural fields. Hydrol. Process. 2005, 19, 2653–2670. [Google Scholar] [CrossRef]
  15. Beven, K.J.; Fisher, J. Remote sensing and scaling in hydrology. In Scaling up in Hydrology Using Remote Sensing; Stewart, J.B., Engman, E.T., Feddes, R.A., Kerr, Y., Eds.; Wiley: Chichester, UK, 1996; pp. 1–18. [Google Scholar]
  16. Price, J.C. Using spatial context in satellite data to infer regional scale evapotranspiration. IEEE Trans. Geosci. Remote Sens. 1990, 28, 940–948. [Google Scholar] [CrossRef]
  17. Jiang, L.; Islam, S. Estimation of surface evaporation map over Southern Great Plains using remote sensing data. Water Resour. Res. 2001, 37, 329–340. [Google Scholar] [CrossRef]
  18. Sandholt, I.; Rasmussen, K.; Andersen, J. A simple interpretation of the surface temperature/vegetation index space for assessment of surface moisture status. Remote Sens. Environ. 2002, 79, 213–224. [Google Scholar] [CrossRef]
  19. Priestley, C.H.B.; Taylor, R.J. On the assessment of surface heat flux and evaporation using large-scale parameters. Mon. Weather Rev. 1972, 100, 81–92. [Google Scholar] [CrossRef]
  20. Long, D.; Singh, V.P. Assessing the impact of end-member selection on the accuracy of satellite-based spatial variability models for actual evapotranspiration estimation. Water Resour. Res. 2013, 49, 2601–2618. [Google Scholar] [CrossRef]
  21. García, M.; Fernández, N.; Villagarcía, F.; Domingo, F.; Puigdefábregas, J.; Sandholt, I. Accuracy of the Temperature–Vegetation Dryness Index using MODIS under water-limited vs. energy-limited evapotranspiration conditions. Remote Sens. Environ. 2014, 149, 100–117. [Google Scholar] [CrossRef]
  22. De Bruin, H. A Model for the Priestley-Taylor Parameter α. J. Clim. Appl. Meteorol. 1983, 22, 572–578. [Google Scholar] [CrossRef]
  23. Kim, J.; Hogue, T.S. Evaluation of a MODIS triangle—Based evapotranspiration algorithm for semi-arid regions. J. Appl. Remote Sens. 2013, 7. [Google Scholar] [CrossRef]
  24. Tanguy, M.; Baille, A.; González-Real, M.; Lloyd, C.; Cappelaere, B.; Kergoat, L.; Cohard, J. A new parameterisation scheme of ground heat flux for land surface flux retrieval from MODIS products. J. Hydrol. 2012, 454–455, 113–122. [Google Scholar] [CrossRef]
  25. Pereira, A.R. The Priestley-Taylor parameter and the decoupling factor for estimating reference evapotranspiration. Agric. For. Meteorol. 2004, 125, 305–313. [Google Scholar] [CrossRef]
  26. Eichinger, W.E.; Parlange, M.B.; Stricker, H. On the concept of equilibrium evaporation and the value of the Priestley Taylor coeficient. Water Resour. Res. 1996, 32, 161–164. [Google Scholar] [CrossRef]
  27. Davies, J.A.; Allen, C.D. Equilibrium, potential and actual evaporation from cropped surface in southern Ontario. J. Appl. Meteorol. 1973, 12, 649–657. [Google Scholar] [CrossRef]
  28. Komatsu, T.S. Toward a robust phenomenological expression of evaporation efficiency for unsaturated soil surface. J. Appl. Meteorol. 2003, 42, 1330–1334. [Google Scholar] [CrossRef]
  29. Romero Diaz, M.A.; Lopez Bermudez, F.; Cabezas, F. Erosion and fluvial sedimentation in the River Segura basin (Spain). Catena 1992, 19, 379–392. [Google Scholar] [CrossRef]
  30. Jiménez-Martínez, J.; García-Aróstegui, J.L.; Hunink, J.E.; Contreras, S.; Baudron, P.; Candela, L. The role of groundwater in highly human-modified hydrosystems: A review of impacts and mitigation options in the Campo de Cartagena-Mar Menor coastal plain (SE Spain). Environ. Rev. 2016, 24, 377–392. [Google Scholar] [CrossRef]
  31. Allen, R.G.; Pereira, L.S.; Raes, D.; Smith, M. Crop Evapotranspiration: Guidelines for Computing Crop Requirements; FAO: Rome, Italy, 1998. [Google Scholar]
  32. Twine, T.E.; Kustas, W.P.; Norman, J.M.; Cook, D.R.; Houser, P.R.; Meyers, T.P.; Prueger, J.H.; Starks, P.J.; Wesely, M.L. Correcting eddy-covariance flux underestimates over a grassland. Agric. For. Meteorol. 2000, 103, 279–300. [Google Scholar] [CrossRef]
  33. Schuepp, P.H.; Leclerc, M.Y.; MacPherson, J.I.; Desjardins, R.L. Footprint prediction of scalar fluxes from analytical solutions of the diffusion equation. Bound. Layer Meteorol. 1990, 50, 355–373. [Google Scholar] [CrossRef]
  34. Chuvieco, E.; Hantson, S. Plan Nacional de Teledetección. Procesamiento Estándar de Imágenes Landsat. Documento Técnico de Algoritmos a Aplicar; Instituto Geográfico Nacional: Madrid, Spain, 2010. [Google Scholar]
  35. Timmermans, W.J.; Kustas, W.P.; Anderson, M.C.; French, A.N. An intercomparison of the Surface Energy Balance Algorithm for Land (SEBAL) and the Two-Source Energy Balance (TSEB) modeling schemes. Remote Sens. Environ. 2007, 108, 369–384. [Google Scholar] [CrossRef]
  36. Cristóbal, J.; Ninyerola, M.; Pons, X. Modeling air temperature through a combination of remote sensing and GIS data. J. Geophys. Res. 2008, 113, 2156–2202. [Google Scholar] [CrossRef]
  37. McVicar, T.R.; Jupp, D.L.B. Using covariates to spatially interpolate moisture availability in the Murray -Darling Basin: A novel use of remote sensed data. Remote Sens. Environ. 2002, 79, 199–212. [Google Scholar] [CrossRef]
  38. Parton, W.J.; Logan, J.A. A model for diurnal variation in soil and air temperature. Agric. Meteorol. 1981, 23, 205–216. [Google Scholar] [CrossRef]
  39. Remund, J.; Wald, L.; Lefèvre, M.; Ranchin, T.; Page, J. Worldwide Linke Turbidity Information. Available online: http://hal.archives-ouvertes.fr/docs/00/46/57/91/PDF/ises2003_linke.pdf (accessed on 23 December 2016).
  40. Wald, L.; Albuisson, M.; Best, C.; Delamare, C.; Dumortier, D.; Kift, R.; Ratto, C.; Reise, C.; Kunz, S.; Pade, J. SoDa: A project for the integration and exploitation of networked solar radiation databases. In Environmental Communication in the Information Society; Pillmann, W., Tochtermann, K., Eds.; International Society for Environmental Protection: Vienna, Austria, 2002. [Google Scholar]
  41. Muneer, T. Solar Radiation and Daylight Models for Energy Efficient Design of Buildings; Butterworth-Heinemann: Oxford, UK, 1997. [Google Scholar]
  42. Hofierka, J.; Suri, M. The solar radiation model for Open source GIS: Implementation and applications. In Proceedings of the Open Source GIS-GRASS Users Conference, Trento, Italy, 11–13 September 2002. [Google Scholar]
  43. Menenti, M. Physical Aspects and Determination of Evaporation in Deserts Applying Remote Sensing Techniques; Institute for Land and Water Management Research (IWC): Wageningen, Holanda, 1984. [Google Scholar]
  44. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W.; Harlan, J.C. Monitoring the Vernal Advancement of Retrogradation of Natural Vegetation; NASA/GSFC: Greenbelt, MD, USA, 1974.
  45. Asrar, G.; Fuchs, M.; Kanemasu, E.T.; Hatfield, J.L. Estimating of absorbed photosynthetic radiation and leaf area index from spectral reflectance in wheat. Agron. J. 1984, 76, 300–306. [Google Scholar] [CrossRef]
  46. Myneni, R.B.; Hall, F.G.; Sellers, P.J.; Marshak, A.L. The interpretation of spectral vegetation indexes. IEEE Trans. Geosci. Remote Sens. 1995, 33, 481–486. [Google Scholar] [CrossRef]
  47. Sobrino, J.A.; Jimenez-Muñoz, J.C.; Paolini, L. Land surface temperature retrieval from LANDSAT TM 5. Remote Sens. Environ. 2004, 90, 434–440. [Google Scholar] [CrossRef]
  48. Carlson, T.N.; Ripley, D.A. On the relation between NDVI, fractional vegetation cover, and leaf area index. Remote Sens. Environ. 1997, 62, 241–252. [Google Scholar] [CrossRef]
  49. Sobrino, J.A.; Jiménez Muñoz, J.C.; Soria, G.; Romaguera, M.; Guanter, L.; Moreno, J.; Plaza, A.; Martínez, P. Land surface emissivity retrieval from different VNIR and TIR sensors. IEEE Trans. Geosci. Remote Sens. 2008, 46, 316–327. [Google Scholar] [CrossRef]
  50. Jiménez-Muñoz, J.C.; Cristóbal, J.; Soria, G.; Ninyerola, M.; Pons, X. Revision of the Single-Channel Algorithm for Land Surface Temperature Retrieval from Landsat Thermal-Infrared Data. IEEE Trans. Geosci. Remote Sens. 2009, 47, 339–349. [Google Scholar] [CrossRef]
  51. Kalma, J.D.; McVicar, T.R.; McCabe, M.F. Estimating land surface evaporation: A review of methods using remotely sensed surface temperature data. Surv. Geophys. 2008, 29, 421–469. [Google Scholar] [CrossRef]
  52. Long, D.; Singh, V.P. A modified surface energy balance algorithm for land (M-SEBAL) based on a trapezoidal framework. Water Resour. Res. 2012, 48. [Google Scholar] [CrossRef]
  53. Batra, N.; Islam, S.; Venturini, V.; Bisht, G.; Jiang, L. Estimation and comparison of evapotranspiration from MODIS and AVHRR sensors for clear sky days over the Southern Great Plains. Remote Sens. Environ. 2006, 103, 1–15. [Google Scholar] [CrossRef]
  54. Bisht, G.; Venturini, V.; Islam, S.; Jiang, L. Estimation of the net radiation using MODIS (Moderate Resolution Imaging Spectroradiometer) data for clear sky days. Remote Sens. Environ. 2005, 97, 52–67. [Google Scholar] [CrossRef]
  55. Xiong, Y.J.; Qiu, G.Y. Simplifying the revised three-temperature model for remotely estimating regional evapotranspiration and its application to a semi-arid steppe. Int. J. Remote Sens. 2014, 35, 2003–2027. [Google Scholar]
  56. Swinbank, W. Long-wave radiation from clear skies. Q. J. R. Meteorol. Soc. 1963, 89, 339–348. [Google Scholar] [CrossRef]
  57. Brutsaert, W. Evaporation into the Atmosphere. In An Introduction to Environmental Biophysics, 2nd ed.; Campbell, G.S., Norman, J.M., Eds.; Springer: New York, NY, USA, 1998. [Google Scholar]
  58. Willmott, C.J. Some Comments on the Evaluation of Model Performance. Am. Meteorol. Soc. 1982, 63, 1309–1313. [Google Scholar] [CrossRef]
  59. Anderson, M.C.; Norman, J.M.; Mechikalski, J.R.; Torn, R.D.; Kustas, W.P.; Basara, J.B. A multiscale remote sensing model for disaggregating regional fluxes to micrometeorological scales. J. Hydrometeorol. 2004, 5, 343–363. [Google Scholar] [CrossRef]
  60. Li, F.; Kustas, W.P.; Prueger, J.H.; Neale, C.M.U.; Jackson, J.J. Utility of remote sensing-based two-source energy balance model under low- and high-vegetation cover conditions. J. Hydrometeorol. 2005, 6, 878–891. [Google Scholar] [CrossRef]
  61. Malbéteau, Y.; Merlin, O.; Gascoin, S.; Gastellu, J.P.; Mattar, C.; Olivera-Guerra, L.; Khabba, S.; Jarlan, L. Normalizing land surface temperature data for elevation and illumination effects in mountainous areas: A case study using ASTER data over a steep-sided valley in Morocco. Remote Sens. Environ. 2017, 189, 25–39. [Google Scholar] [CrossRef]
  62. Wang, K.; Wang, P.; Li, Z.; Cribb, M.; Sparrow, M. A simple method to estimate actual evapotranspiration from a combination of net radiation, vegetation index, and temperature. J. Geophys. Res. 2007, 112. [Google Scholar] [CrossRef]
  63. Hengl, T.; Heuvelink, G.B.; Tadić, M.P.; Pebesma, E.J. Spatio-temporal prediction of daily temperatures using time-series of MODIS LST images. Theor. Appl. Climatol. 2012, 107, 265–277. [Google Scholar] [CrossRef]
  64. Cherif, I.; Alexandridis, T.K.; Jauch, E.; Chambel-Leitao, P.; Almeida, C. Improving remotely sensed actual evapotranspiration estimation with raster meteorological data. Int. J. Remote Sens. 2015, 36, 4606–4620. [Google Scholar] [CrossRef]
  65. Jiang, L.; Islam, S.; Guo, W.; Jutla, A.S.; Senarath, S.U.S.; Ramsay, B.H.; Eltahir, E.A.B. A satellite-based Daily Actual Evapotranspiration estimation algorithm over South Florida. Glob. Planet. Chang. 2009, 67, 62–77. [Google Scholar] [CrossRef]
  66. Nemani, R.R.; Keeling, C.D.; Hashimoto, H.; Jolly, W.M.; Piper, S.C.; Tucker, R.B.; Myneni, R.B.; Running, S.W. Climate-Driven Increases in Global Terrestrial Net Primary Production from 1982 to 1999. Science 2003, 300, 1560–1563. [Google Scholar] [CrossRef] [PubMed]
  67. Mattar, C.; Franch, B.; Sobrino, J.A.; Corbari, C.; Jiménez-Muñoz, J.C.; Olivera-Guerra, L.; Skokovic, D.; Sória, G.; Oltra-Carri, R.; Julien, Y.; et al. Impacts of the broadband albedo on actual evapotranspiration estimated by S-SEBI model over an agricultural area. Remote Sens. Environ. 2014, 147, 32–42. [Google Scholar] [CrossRef]
  68. Villalobos, F.J.; Testi, L.; Moreno-Perez, M.F. Evaporation and canopy conductance of citrus orchards. Agric. Water Manag. 2009, 96, 565–573. [Google Scholar] [CrossRef]
  69. Gillies, R.R.; Carlson, T.N. Thermal remote sensing of surface soil water content with partial vegetation cover for incorporation into climate models. J. Appl. Meteorol. 1995, 34, 745–756. [Google Scholar] [CrossRef]
  70. Moran, M.S.; Clarke, T.R.; Inoue, Y.; Vidal, A. Estimating crop water deficit using the relation between surface-air temperature and espectral vegetation index. Remote Sens. Environ. 1994, 49, 246–263. [Google Scholar] [CrossRef]
Figure 1. Study area: (a) Segura River Basin false-color mosaic of Landsat, natural pseudo-color (June 2011), extension of study (boundary of scene 199-34) and area of analyses (Campo de Cartagena); (b) SPOT image of farms (mosaic of 2012 year); and (c) eddy-covariance flux tower. Coordinate system: ETRS89 UTM Zone 30N.
Figure 1. Study area: (a) Segura River Basin false-color mosaic of Landsat, natural pseudo-color (June 2011), extension of study (boundary of scene 199-34) and area of analyses (Campo de Cartagena); (b) SPOT image of farms (mosaic of 2012 year); and (c) eddy-covariance flux tower. Coordinate system: ETRS89 UTM Zone 30N.
Remotesensing 09 00611 g001
Figure 2. Land cover map of CRCC. Coordinate system: ETRS89 UTM Zone 30N.
Figure 2. Land cover map of CRCC. Coordinate system: ETRS89 UTM Zone 30N.
Remotesensing 09 00611 g002
Scheme 1. Scheme of implemented processes.
Scheme 1. Scheme of implemented processes.
Remotesensing 09 00611 sch001
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. 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.
Remotesensing 09 00611 sch002
Figure 3. Comparison of remote sensing estimated variables and ground truth at overpass time satellite: (a) LST; and (b) Rn.
Figure 3. Comparison of remote sensing estimated variables and ground truth at overpass time satellite: (a) LST; and (b) Rn.
Remotesensing 09 00611 g003
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. 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.
Remotesensing 09 00611 g004
Figure 5. Boxplot of EF considering α = 1.3 and α = 1, two shapes, on dates 11/07/2010, 29/09/2010, 16/11/2010 and 28/06/2011: (a,c,e,g) trapezoidal shape (TR); and (b,d,f,h) rectangular shape (RC).
Figure 5. Boxplot of EF considering α = 1.3 and α = 1, two shapes, on dates 11/07/2010, 29/09/2010, 16/11/2010 and 28/06/2011: (a,c,e,g) trapezoidal shape (TR); and (b,d,f,h) rectangular shape (RC).
Remotesensing 09 00611 g005
Figure 6. Frequency histograms of EF for α = 1, two shapes, on dates 11/07/2010, 29/09/2010, 16/11/2010 and 28/06/2011: (a,c,e,g) trapezoidal shape (TR); and (b,d,f,h) rectangular shape (RC).
Figure 6. Frequency histograms of EF for α = 1, two shapes, on dates 11/07/2010, 29/09/2010, 16/11/2010 and 28/06/2011: (a,c,e,g) trapezoidal shape (TR); and (b,d,f,h) rectangular shape (RC).
Remotesensing 09 00611 g006
Figure 7. Spatial distribution of relative EF differences for rectangular to trapezoidal shape, on dates: (a) 11/07/2010; (b) 29/09/2010; (c) 16/11/2010; and (d) 28/06/2011. Coordinate system: ETRS89 UTM Zone 30N.
Figure 7. Spatial distribution of relative EF differences for rectangular to trapezoidal shape, on dates: (a) 11/07/2010; (b) 29/09/2010; (c) 16/11/2010; and (d) 28/06/2011. Coordinate system: ETRS89 UTM Zone 30N.
Remotesensing 09 00611 g007
Figure 8. Relationship between observed and estimated values of ET, for α = 1.3 ((a) trapezoidal shape; and (b) rectangular shape); and α = 1 ((c) trapezoidal shape; and (d) rectangular shape).
Figure 8. Relationship between observed and estimated values of ET, for α = 1.3 ((a) trapezoidal shape; and (b) rectangular shape); and α = 1 ((c) trapezoidal shape; and (d) rectangular shape).
Remotesensing 09 00611 g008
Figure 9. Contrast of observed EF versus estimated EF, for the two farms.
Figure 9. Contrast of observed EF versus estimated EF, for the two farms.
Remotesensing 09 00611 g009
Figure 10. Relationship between observed and estimated values of available energy (A = Rn − G) for the trapezoidal space.
Figure 10. Relationship between observed and estimated values of available energy (A = Rn − G) for the trapezoidal space.
Remotesensing 09 00611 g010
Table 1. Details of the used images and irrigation conditions.
Table 1. Details of the used images and irrigation conditions.
DateSeasonIrrigation Conditions
Farm AFarm B
wfv (m3·m−3)wfv (m3·m−3)
Mean ± StdMean ± Std
10/09/2009Summer0.361 ± 0.0130.453 ± 0.002
08/05/2010Spring0.440 ± 0.0360.385 ± 0.002
24/05/2010Spring0.432 ± 0.0370.390 ± 0.003
11/07/2010Summer0.400 ± 0.0210.425 ± 0.002
27/07/2010Summer0.427 ± 0.0230.465 ± 0.002
29/09/2010Autumn0.450 ± 0.0360.477 ± 0.001
16/11/2010Autumn0.474 ± 0.0130.454 ± 0.001
02/12/2010Autumn0.447 ± 0.0370.471 ± 0.002
04/02/2011Winter0.436 ± 0.0540.463 ± 0.002
28/06/2011Summer0.321 ± 0.0220.307 ± 0.019
wfv: water fraction by volume (m3·m−3).
Table 2. Values of NDVI and fc assessed from the satellite images for the pixel of ET tower flux measurement in Farms A and B (F), respectively.
Table 2. Values of NDVI and fc assessed from the satellite images for the pixel of ET tower flux measurement in Farms A and B (F), respectively.
DateFNDVIfc
10/09/2009A0.340.24
B0.541.00
08/05/2010A0.440.64
B0.601.00
24/05/2010A0.460.74
B0.541.00
11/07/2010A0.410.48
B0.561.00
27/07/2010A0.430.53
B0.631.00
29/09/2010A0.480.74
B0.591.00
16/11/2010A0.420.55
B0.551.00
02/12/2010A0.541.00
B0.611.00
04/02/2011A0.470.82
B0.521.00
28/06/2011A0.480.89
B0.621.00
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 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.
DateShape Dry Edge Wet Edge
a LSTb LSTa DTb DTa’ LSTa’ DT
10/09/2009TR47.7−11.223.46−15.915.24−14.6
RC44.45 18.6 15.24−14.6
08/05/2010TR42.6−16.023.3−11.910.4−15.1
RC38.2 20.7 10.4−15.1
24/05/2010TR48.2−20.528.5−22.918−5.9
RC45.7 23.9 18−5.9
11/07/2010TR51.3−17.524.2−17.927.3−5.3
RC47 19.8 27.3−5.3
27/07/2010TR51.3−17.521−17.515.8−12.6
RC44.4 16.9 15.8−12.6
29/09/2010TR38.1−11.223.3−22.711.5−9.7
RC35.2 17.6 11.5−9.7
16/11/2010TR40.3−29.029.3−33.29.7−6.0
RC23.8 21.6 9.7−6.0
02/12/2010TR41.6−38.625.8−34.14.1−8.6
RC24.6 16.5 4.1−8.6
04/02/2011TR41.5−38.636.4−47.77.8−6.6
RC30.8 22.2 7.8−6.6
28/06/2011TR60.1−17.235.0−22.329.3−7.5
RC54.94 28.12 29.3−7.5
Table 4. Influence of Tair in the Triangle method for α = 1.27.
Table 4. Influence of Tair in the Triangle method for α = 1.27.
LSTDT
FR2RMSE (W/m2)MAE (W/m2)MAPE %R2RMSE (W/m2)MAE (W/m2)MAPE %
TRA0.8184.8770.1844.270.8045.2636.6519.82
B0.9191.4983.0551.270.8558.9853.3729.90
RCA0.8286.9977.1546.470.8054.1944.1822.98
B0.91103.1394.0054.180.8578.5266.4533.79
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 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.
Date α = 1.3α = 1
ETETAEAPEETAEAPE
W/m2W/m2W/m2%W/m2W/m2%
FTowerTRRCTRRCTRRCTRRCTRRCTRRC
10/09/09A175.7171.2164.24.511.62.66.6127.8122.64853.127.330.2
B197.0224.8243.127.746.114.123.4167.8141.229.33014.917.5
08/05/10A171.2173188.21.9171.19.9130236.641.112.2245.4
B118.0171.8202.553.884.445.671.5129.2238.411.227.89.510.4
24/05/10A224.4281.7318.857.394.525.542.1209.7214.214.73.76.61.8
B238.0290.9358.153120.122.350.5216.619921.428.6916.8
11/07/10A266.2305.1321.638.955.414.620.8226.5103.239.754.414.934.5
B232.0322354.290122.238.852.7238.987.66.928.6324.6
27/07/10A210.4282.9288.372.577.934.437210.294.20.212.40.111.6
B222.02843216298.927.944.6211217.9112.351.1
29/09/10A170.4256.6269.886.299.450.658.4189.5181.119.215.911.28.1
B215.8276.230660.490.22841.8203.9151.811.833.85.528.6
16/11/10A157.6130.5139.327.118.217.211.696.8265.160.827.138.611.4
B155.4120135.135.420.322.813.189.2262.166.230.142.613
02/12/10A116.299.1117.5171.414.71.274.1237.84215.836.27.1
B123.781107.342.716.434.513.360.8225.262.99.450.84.4
04/02/11A106.7115.7126.69208.518.786.3100.120.455.319.135.6
B110.1102.9118.17.286.57.3778033.143.730.135.3
28/06/11A220.3279.1293.358.873.126.733.2207.68812.622.15.720.1
B251.0324.4362.273.3111.229.244.3240.726810.3174.16.8
MAE MAPEA 37.346.819.624 29.825.318.315.4
B 50.671.82736.2 26.42717.417
RMSEA 4758.7 34.930.7
(W/m2)B 55.383.4 33.630.2
Table 6. Influence of α selection in EF values.
Table 6. Influence of α selection in EF values.
α = 1.3α = 1
FμΣMax.Min.µσMax.Min.
TRA0.570.100.730.390.440.070.560.30
B0.550.110.740.370.430.080.570.29
RCA0.600.100.760.420.460.070.590.32
B0.620.110.800.430.480.080.620.33

Share and Cite

MDPI and ACS Style

Martínez Pérez, J.Á.; García-Galiano, S.G.; Martin-Gorriz, B.; Baille, A. Satellite-Based Method for Estimating the Spatial Distribution of Crop Evapotranspiration: Sensitivity to the Priestley-Taylor Coefficient. Remote Sens. 2017, 9, 611. https://0-doi-org.brum.beds.ac.uk/10.3390/rs9060611

AMA Style

Martínez Pérez JÁ, García-Galiano SG, Martin-Gorriz B, Baille A. Satellite-Based Method for Estimating the Spatial Distribution of Crop Evapotranspiration: Sensitivity to the Priestley-Taylor Coefficient. Remote Sensing. 2017; 9(6):611. https://0-doi-org.brum.beds.ac.uk/10.3390/rs9060611

Chicago/Turabian Style

Martínez Pérez, José Ángel, Sandra G. García-Galiano, Bernardo Martin-Gorriz, and Alain Baille. 2017. "Satellite-Based Method for Estimating the Spatial Distribution of Crop Evapotranspiration: Sensitivity to the Priestley-Taylor Coefficient" Remote Sensing 9, no. 6: 611. https://0-doi-org.brum.beds.ac.uk/10.3390/rs9060611

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop