Modifying SEBAL Model Based on the Trapezoidal Relationship between Land Surface Temperature and Vegetation Index for Actual Evapotranspiration Estimation

The Surface Energy Balance Algorithm for Land (SEBAL) is widely used to estimate actual evapotranspiration (ETa). One major limitation of the SEBAL model is the subjectiveness in selecting extreme cold/hot pixels. In the present study, the SEBAL model is modified by determining the extreme cold/hot status, based on the theoretical trapezoidal relationship between land surface temperature (Ts) and Enhanced Vegetation Index (EVI), which is established for each pixel. In this way, the dependence of SEBAL model on the existence of extreme cold/hot status and the subjectiveness in selecting cold/hot pixels with SEBAL model are eliminated. The performance of the classical SEBAL model and the modified version, T-SEBAL, are compared for estimating ETa for a semi-arid catchment, and the result showed that the accuracy of ETa estimation is improved by the T-SEBAL model compared with the classical SEBAL model.


Introduction
Evapotranspiration (ET) plays an important role in hydrological cycles.It is a complex process influenced by conditions of solar radiation, atmosphere, and land surface, which makes it difficult to estimate actual evapotranspiration (referred to as ET a hereafter) accurately.The early research of ET a mainly focuses on developing meteorological approaches on the local scale [1][2][3][4][5].However, such methods can rarely be extended to large areas due to land surface heterogeneity and the dynamic nature of the heat transfer processes [6].Following the development of multi-temporal and multi-spectral remote sensing databases which can provide key surface characteristic (e.g., albedo, land surface temperature, and vegetation index), models of spatially and temporally dynamic feedback mechanism based on the land surface energy balance are proposed.The models retrieving ET a mainly relying on remotely sensed data can be roughly categorized into residual methods, which calculate latent heat flux as the residual of the surface energy balance (i.e., net radiation R n minus sensible heat flux H and soil heat flux G) [7,8], and ratio methods which estimate the ratio of latent heat flux to available energy or the ratio of ET a to potential ET [9][10][11][12].
The most widely used residual method/model so far is the Surface Energy Balance Algorithm for Land (SEBAL) proposed by Bastiaanssen et al. [7,13], which estimates regional ET a with remotely sensed land surface temperature, surface emissivity, NDVI, and minimum requirement of some routine ground meteorological observation data, i.e., sunshine duration and wind speed.The core assumption in SEBAL is a linear relationship between near-surface vertical temperature gradient and land surface temperature, and the relationship is determined by two extreme hydrological pixels, i.e., extreme wet/cold and extreme dry/hot pixels.Although SEBAL has been applied in more than 30 countries worldwide with considerable accuracy [13,14], there are still some deficiencies, including the requirement of the existence of extreme "cold" and "hot" pixels throughout the image scene, the dependence of extreme pixel selection on domain size, and the subjectiveness in selecting the extreme pixels by different users, especially in the context of complicated surfaces.
To mitigate the limitation resulted from the manual selection of cold and hot points with SEBAL model, Allen et al. [15] proposed the Mapping Evapotranspiration at high Resolution with Internalized Calibration (METRIC).METRIC departs from the SEBAL model in its use of weather-based reference ET to establish energy balance conditions at a "cold" pixel, and a value of 1.05 ET r (reference ET for 0.5 m tall alfalfa) estimated by standardized ASCE Penman-Monteith equation is considered to be ET a of the "cold" pixel.In addition, METRIC considers the effects of complicated topography.However, METRIC still need to select a "representative cold pixel" from the satellite scene, which has biophysical characteristics similar to the reference crop (alfalfa).Su [9] proposed the SEBS model to estimate the evaporative fraction based on energy balance of dry/wet limiting cases (i.e., at the wet-limit, ET takes place at its potential rate, while at the dry-limit point, the sensible heat flux reaches its maximum value R n − G) for each pixel.SEBS does not require the existence and the selection of extreme pixels, but it needs many required parameters and a relatively complex solution of the turbulent heat fluxes which may introduce considerable uncertainty when a priori knowledge about the research area is not sufficient.Moran et al. (1994) [11] proposed a measured surface minus air temperatures (T s − T a ) vs. vegetation cover (V c ) trapezoidal space which defined the extreme vertices as well-watered vegetation, water-stressed vegetation, saturated bare soil, and dry bare soil, respectively, and used measured T a , R n , and G to calculate the T s − T a for four vertices.However, using measured values of R n and G for four vertices of (T s − T a )~V c may distort the trapezoidal space significantly, because the deviation of R n between bare soil and vegetation may be up to 20% [16], and deviation of G is even larger.In addition, T s − T a measurements for the four extreme vertices of the (T s − T a )~V c are not available in real application.Therefore, Wang et al. (2011) [17] suggested a pixel-based T s ~vegetation index (VI) trapezoidal model to simplify the (T s − T a )~V c trapezoidal space, which uses an iterative procedure to consider the interaction process between T s with R n and the difference of G at different vertices.The T s ~VI trapezoidal space have been used in combination with the Priestley-Taylor equation to estimate ET a [18].Following the same line, Long et al. (2012) [19] proposed a vegetation fraction V c ~Ts trapezoidal framework with different formulation, where the warm edge is determined by analytically deriving temperatures of the bare surface with the largest water stress and the fully vegetated surface with the largest water stress, and the cold edge is set to be the regional average T a .
In the present study, we would like to develop a method to define the hot and cold pixels for the SEBAL model by calculating surface temperatures at dry bare soil points and well-watered vegetation points of the pixel-based T s ~VI trapezoid formulated by Wang et al. (2011) [17], so as to eliminate the limitations with SEBAL model.The method is described in Section 2; and the application result is presented in Section 4.

Method of
R n is given by: where, R s is the solar shortwave radiation (W•m −2 ), which is estimated by the empirical parameterization of Zillman (1972) [20] in the present study; ε s and ε a are the surface emissivity and atmospheric emissivity respectively; α is albedo (dimensionless); σ is the Stefan-Bolzman constant (5.67 × 10 −8 W•m −2 •K −4 ); T s and T a are the land surface temperature (K) and air temperature (K) at satellite observation time.
G is calculated as a friction of R n .The ratio G/R n is usually defined as a function of leaf area index (LAI) or Normalized Difference Vegetation Index (NDVI).In SEBAL, G is given by: ( ) 273.15 (0.0032 0.0062 )(1 0.978 ) where, α avg is the daily average albedo which cannot be obtained from remotely sensed images directly.
As the variation in albedo under all clear sky presents a U-shape, high at the sunrise and sunset and low at the noon, therefore a linear algorithm was suggested to calculate α avg using instantaneous α, i.e., α avg = 1.02α + 0.01 [14].
Another method for estimating G is proposed by Su [9], given by: ( where, Г c and Г s denote the ratio of G/R n for full vegetation canopy and bare soil, originally set as 0.05 and 0.315, respectively; V c is the fractional canopy coverage, defined as (EVI − EVI min )/(EVI max − EVI min ).It is found that, setting Г s = 0.315 will lead to serious overestimation of G compared with ground observations in our case study area, therefore, Г s is set to be 0.28 in the present study by a trial and error approach.As, globally, EVI mostly ranges from 0.07 for deserts to 0.7 for Amazon forests at 1-km resolution (MOD12A2) [21], and in our study area there are many pixels with EVI below 0.07, we set the value of EVI max (maximum EVI) as 0.7 for full vegetation coverage, and EVI min (minimum EVI) 0.05 for bare soil.
H is calculated by: ( where, ρ is the air density (kg•m −3 ), C p is the air specific heat at constant pressure (1004 J•K −1 •kg −1 ), T 0 is the aerodynamic temperature (K), r a is the aerodynamic resistance (m•s −1 ).As T 0 cannot be observed directly, whereas T s could be retrieved from remotely sensed images, SEBAL establishes a linear relationship relating the vertical temperature gradient dT (=T 0 − T a ) to T s , i.e.: dT = a + bT s (6) where, a and b are coefficients of the linear relationship throughout a scene, estimated using T s at two anchor points selected from the satellite image, i.e., a hot point where H = R n − G and a cold point where H = 0.At the cold point, dT = 0, whereas at the hot point, dT reaches its maximum dT hot .r a is given by: where, k is the von Karman constant (=0.4); z 2 and z 1 are the upper and lower reference heights for aerodynamic resistance taken as 2 m and 0.1 m [13], respectively; u* is the friction velocity obtained at the blending layer of 100 m to 200 m; φ h is the stability correction of atmosphere for heat transport based on Monin-Obukhov length L (dimensionless, equals ρC p u* 3 T s /(kgH)); g is gravitational acceleration (=9.8 m•s −2 ).In fact, z 1 is equivalent to the roughness length for heat transfer z oh .An alternative value of 0.01 for z 1 is also used in other studies [22,23].Paul et al. (2014) [24] found that z oh = 0.1 m and z oh = 0.01 m led to the same estimate of ET a .Thus, z oh = 0.01 m is used in the present study for Equation (7), which is assumed to be more suitable to the vegetation sparse area.
The calculation of H is the core of the SEBAL model.With dT and r a calculated, H could be estimated, and latent heat flux, LE, is immediately derived according to Equation (1).

Limitations of SEBAL Model
The existence and proper selection of cold and hot pixels in the image scene are crucial for the successful application of SEBAL, as they are essential for deriving variables a and b to estimate the temperature gradient dT.
The cold point is defined as the well-watered pixel with dense vegetation.Bastiaanssen et al. [13] chose water body as the cold point assuming that T 0 and T a are approximately equal at such sites.The hot point is defined as bare soil pixels with the highest T s [25], such as the dry bare agricultural field, but it is recommended that the hot point should not be treated as hot desert area, an asphalt floor, and biomass burning pixels [26], or other extreme hot areas which cannot represent natural surface, as the slope of the dT equation in such areas may be distorted.
However, it is cumbersome to find appropriate hot and cold pixels in some situations.For instance, water bodies does not always present in an image, especially for semi-arid areas when the image resolution is coarse.In that case, the average air temperature of the region may be used as cold point temperature [19].In a catchment experienced a long-time drought or heavy rainfall, cold or hot points may do not exist at all.Furthermore, for mountainous areas, it is not possible to define a universal linear relationship between the temperature gradient dT and surface temperature T s applicable to various topographic locations.In addition, the selection of dry and wet pixels over the area of interest (AOI) makes the SEBAL-estimated ET a subject to the sizes of the AOI and the satellite image pixels [27].
On the other hand, the calculation of H is quite sensitive to the selection of extreme pixels.It was found that a 2 K increase in T s for hot point may lead to 12% decrease in H, while for the cold point, a 2 K increase may result in 15% decrease in H in some cases [28].Due to the importance and the complexity in defining the hot and cold pixels, it is suggested that experienced image processors are needed to correctly identify hot and cold pixels [29].

Modification of SEBAL Model Based on T s ~VI Trapezoidal Relationship
In the T s ~VI trapezoidal space, the four vertices represent four extreme conditions, i.e., well-watered vegetation, water-stressed vegetation, saturated bare soil, and dry bare soil (Figure 1).The dry bare soil vertex is defined for the case of bare surface with the largest water stress, which is equivalent to the definition of hot point in SEBAL; the well-watered vegetation vertex represents the situation with full vegetation coverage and no water stress, corresponding to the cold pixel in SEBAL.Therefore, the temperatures at hot and cold points can be calculated based on the T s ~VI trapezoid pixel by pixel, and the determination of hot/cold pixels in SEBAL is modified from manual selection to physically-based calculation.The T s ~VI trapezoid based modification of SEBAL model is referred to as T-SEBAL hereafter.The pixel-based procedure for estimating H with T-SEBAL is composed of the following two parts describe in Sections 2.3.1 and 2.3.2.The core of constructing T s ~VI trapezoidal space is the calculation of T s for four extreme vertices with an iterative procedure considering the interaction among T s , R n , r a , and u*.The following formulas are used for calculating T s for each vertex in the T s ~VI trapezoidal space [11,19].
• Point 1, cold point.For full vegetation-covered and well-watered surface.
, 4 where, the subscript n (1, 2, 3 and 4) of T s refers to the nth vertex in Figure 1; r a ' is the aerodynamic resistance estimated with Equation ( 9), which is different from the formula Equation (7) used in SEBAL; r c is the canopy resistance (m•s −1 ); r cm and r cx are the minimum and maximum canopy resistance; Δ is the slope of saturated vapor pressure to air temperature; γ is the air humidity constant (hPa•C −1 ); VPD is the vapor pressure deficit of the air (hPa); C v is the volumetric heat capacity of air (1295.16J•K −1 •m −3 ); G is estimated as a fraction of R n .For full-cover vegetation edge, G/R n is negligible, say 0.05, while for dry bare soil and wet bare soil, G can approximately reach up to 0.5R n and down to 0.15R n , respectively [30].R n is estimated with Equation (2).For bare surface and full-covered vegetation, ε s are set as 0.93 and 0.993, respectively, according to the result of Snyder et al. [31].The albedo is a key parameter to estimate R n .In the work of Wang et al. [17], the surface shortwave albedo derived from MODIS product MCD43A3 is used for all the four vertices, which is not physically reasonable.Hence, values of albedo at four extreme cases are defined separately in the present study.There are some suggestions about the albedo values for different land surfaces [32][33][34].In our application, considering those suggestions and by a trial-and-error procedure, albedos are set to be 0.1 and 0.25 for saturated bare soil (point 3) and dry bare soil (point 4), and 0.18 and 0.20 for well-watered vegetation (point 1) and water-stressed vegetation (point 2).Aerodynamic resistance r a ' is given by [35]: where, u is wind speed (m•s −1 ); d is displacement height (m), given by 0.67h, with the height of vegetation h estimated from z 0m , i.e., h = 8z 0m ; z 0m is the momentum roughness length, defined as a function of NDVI, i.e., exp(a + b•NDVI), where we set a = −5.2 and b = 5.3 [36]; z 0h is the heat roughness length; ψ h and ψ m are the stability correction for heat and momentum, calculated according to the value of Monin-Obukhov length L. The heat roughness length z 0h is estimated by z 0h = z 0m /exp(kB −1 ), where kB −1 is a dimensionless parameter.Many formulas have been used for estimating kB −1 with different complexity.Verhoef et al. (1997) [37] evaluated several methods of kB −1 estimation for the bare soil site with the soil classified as silt loam and showed that the equation of Owen and Thomson (1963) [38] performed well, so it is used for the situation of bare soil, given by: where, the subscript s refers to bare soil; α is a factor lying between 0.45 and 0.7 with an average of 0.52; Re * is the roughness Reynolds number and equals to z 0m,s u*/v, with z 0m,s and v being the momentum roughness length under bare soil circumstance and kinematic viscosity of the air; Pr is the Prandtl number (0.71); m and n are constants given by 0.45 and 0.8, respectively.
For full canopy coverage situations, the formula developed by Blümel (1999) [39] is used, given by: where C k = 16.4 m −1 •s 1/2 ; D f is the leaf width, set to be 0.01 m (for open shrubland) as in NASA Global Land Data Assimilation System (GLDAS); z 0m,c is the momentum roughness length of full canopy coverage; σ α is the moentum partion factor relevant to LSAI (leaf and stem area index).If LSAI is not known, the factor (σ a LSAI 3 ) −0.25 is approximately set as 0.4 provided that LSAI ≥ 4. Canopy resistance r c is estimated as r s /LAI max , where r s is the stomatal resistance (s•m −1 ) and LAI max is the maximum LAI in full-cover vegetation.In the paper of Wang et al. [17], the maximum r s (r smax ) and minimum r s (r smin ) are set as 1500 s•m −1 and 100 s•m −1 , respectively.However, it was found that setting r smax = 1500 cannot adequately consider the constraint of evapotranspiration under full water-stress condition, therefore, we set r smax as 5000 s•m −1 as suggested by Dickinson [40].There are several schemes available for setting vegetation parameters, such as NASA GLDAS [41], and European Centre for Medium-Range Weather Forecasts (ECMWF) [42].In our application, assuming that the vegetation for the semi-arid zone is closed shrubland, according to the vegetation parameter scheme of NLDAS, r smin is set to be 175 s•m −1 and LAI max 5.
Taking point 4 as an example, the quartic equation is given by: (1 ) where, r 4 is the ratio of G to R n at point 4. By solving these equations, we get T s for each vertex of the T s ~VI trapezoid.
The procedure for constructing the T s ~VI trapezoid on a pixel basis is summarized in the dash line box of Figure 2, which may be divided into three major steps.(1) Estimate the initial value of r a ' using meteorological data without considering the two stability corrections, i.e., ψ h and ψ m ; (2) Iteratively solve the quartic equation for T s (e.g., Equation ( 12)), and calculate H, L, and r a ' sequentially for each vertex, after 10 iterations when the r a ' is stable, i.e., difference between two adjacent calculations of r a ' is smaller than five percent; (3) Verify the reasonableness of T s ~VI trapezoid by comparing the location of every pixel's (EVI, T s ) point with the hot/cold edges of the T s ~VI trapezoid, considering that there are errors in T a which is extrapolated from a few observation sites and uncertainties in the parameterization for the calculation of aerodynamic resistance, which could lead to errors in T s at vertices of T s ~VI trapezoids.If an observed (EVI, T s ) point falls outside the hot/cold edges, an adjustment is made to make the point be contained within the T s ~VI trapezoid.Figure 3 gives an example for the adjustment when a (EVI, T s ) point falls above the warm edge.

Estimate Sensible Heat Flux H
With the theoretical T s ~VI trapezoid constructed for each pixel uniquely, the calculated values of T s,1 at point 1 and T s,4 at point 4 are used as the hot point temperature and cold point temperature respectively.T s,1 and T s,4 , together with calculated R n −G at point 4 (i.e., (R n −G) 4 ), are used to estimate the coefficients a and b in Equation ( 6) for each pixel, given by: = − × s, a b T (14) where, (R n −G) 4 is the calculated sensible heat flux of point 4 (H 4 ) (Figure 1); r a,hot is the aerodynamic resistance at hot point calculated using Equation (7).Then dT and r a for each pixel are estimated using Equations ( 6) and (7), respectively.Finally, H for each pixel could be estimated with Equation (5).The procedure for estimating H is described in the right part of Figure 2, following the procedure of standard SEBAL.

Case Study Area
The Walnut Gulch Experimental Watershed (WGEW) is located in southeastern Arizona, United States, with an area of 150 km 2 contained within the upper San Pedro River Basin (Figure 4).Climate at WGEW is classified as semiarid.The mean annual air temperature and precipitation are 17.7 °C and 312 mm respectively.Sixty percent of the precipitation concentrates in monsoon season between July and September.Desert shrubs dominate the lower two-thirds of WGEW and the remaining upper one third is desert grasses.There is a very dense rain gauge network with 88 weighing rain gauges within the WGEW.The meteorological dataset of air temperature, relative humidity and wind speed, etc., are measured at three sites, i.e., Lucky Hills, Kendall and ShopMet (Figure 4), recorded at the same frequency of 20 min.Air temperature is extrapolated to the whole study area by the method described in Wang et al. [17], considering the leaf area index (LAI), elevation and topography feature.Other meteorological variables at the three observation sites are averaged and applied to the whole study area.
Products of the moderate resolution imaging spectroradiometer (MODIS) are used in the present study as remotely sensed input data for the model (Table 1).Data with 1 km and 250 m resolution are resampled to 500 m using the nearest neighbor method.Fifteen dates with clear-sky and complete meteorological observations in 2004 were selected to estimate the fluxes with the T-SEBAL and SEBAL model.To ensure a better chance that the cold/hot points can be reliably selected for each date with SEBAL model and subsequently make a fair comparison between SEBAL and T-SEBAL model, images covering an area much larger than the WGEW (as shown in Figure 4) are used in the present study.
In summary, major data used in the present study include: (1) Ground based meteorological observations: air temperature T a , wind speed u, and relative humidity U; (2) MODIS products: land surface temperature T s , satellite overpass time, surface emissivity ε s , vegetation index NDVI and EVI, albedo, and latitude (radian).(3) DEM from the Shuttle Radar Topography Mission (SRTM).

Surface Net Radiation Flux and Soil Heat Flux
The calculations for surface net radiation flux and soil heat flux are the same for both T-SEBAL and SEBAL.
Net radiation flux is estimated using Equation ( 2), and the intermediate variable-solar shortwave radiation-is calculated with the empirical parameterization by Zillman (1972) [20].Ground soil heat flux G is estimated with two equations, i.e., Equations ( 3) and (4).Comparison between calculated net radiation flux and soil heat flux with observations in 15 dates are shown in Figure 5. R n estimates are also in good agreement with the observations.By fitting the simple linear regression between estimated R n and the observations at Lucky Hills and Kendall, we found that the coefficient of determination r 2 is up to 0.79.The RMSE of estimated R n at Lucky Hills and Kendall are 39.5 W•m −2 and 29 W•m −2 respectively.As shown in Figure 5, there is generally no obvious bias between the estimated and the observed, however, R n is overestimated at Lucky Hills, while underestimated at Kendall.
Compared with ground observations, G estimates with Equation (3) poorly agree with the observations (Figure 5c), with a RMSE of 44.86 W•m −2 .Comparatively, G estimates with Equation (4) are better (Figure 5d), just slightly overestimated, with RMSE = 33 W•m −2 and average relative error 21.6% (max relative error 74%).Therefore the Equation ( 4) is used in the calculation of H for both T-SEBAL and SEBAL in later procedures.

Sensible Heat Flux and Latent Heat Flux
H is calculated differently with SEBAL and T-SEBAL because of their different ways of defining extreme cold/hot points, as described in Sections 2.1 and 2.3.
For the SEBAL model, a subjective selection procedure is conducted for selecting cold/hot points.Candidates of hot points are defined as those pixels in the satellite image scene with their surface temperatures above the 95 percentile of all pixels, NDVI less than 0.2, the slope less than 10%, and none urban and built-up land-cover.Then, the pixel with the highest albedo was chosen as the hot point from the candidates, which have similar pixels around them.As no water body exists in this semi-arid region, the lower one between the minimum T s throughout the scene and the average T a of the whole area is defined as the temperature of the cold pixel.It is found that, among the 15 dates investigated, the minimum T s is mostly larger than the average T a except for DOY345.
For the T-SEBAL model, the values of T s at the four vertices of T s ~EVI trapezoid for each pixel are estimated using the method described in Section 2.3.The calculated T s at T s ~EVI trapezoid vertices for all the pixels of the study area are plotted in the form of box-and-whisker plot in Figure 6 to illustrate the reasonability of the trapezoids for four exemplary dates in different seasons.As illustrated in Figure 6, the points in DOY183 and DOY 267 are closer to the warm edge (the line between point 4 and point 2) than in DOY75 and DOY345, indicating more intense water scarcity in DOY183 and DOY 267 than in DOY75 and DOY345.
The hot and cold point temperatures (i.e., T s,hot and T s,cold ) determined by T-SEBAL and SEBAL in 15 dates for the study area are illustrated in Figure 7 and presented in Table 2.The values of T s,hot and T s,cold of T-SEBAL shown in Figure 7 are the average of the whole image, because T-SEBAL computes T s,hot and T s,cold on a pixel basis.It is shown that, T s,hot given by T-SEBAL is higher than that by SEBAL about 1.6 K on average, and T s,cold by T-SEBAL is significantly lower than by SEBAL with about 2.9 K. Consequently, T-SEBAL retrieves larger difference between hot point and cold point temperatures.Besides the differences in the estimation for hot and cold point temperatures, from Table 2, we can also see significant differences in the values of (R n − G) at hot points calculated by T-SEBAL and SEBAL (denoted by (R n − G) 4 and (R n − G) hot , respectively).The values of (R n − G) 4 given by T-SEBAL are on average smaller than (R n − G) hot given by SEBAL.The differences in not only hot and cold point temperatures, but also the values of (R n − G), together lead to the significant difference in establishing the linear relationship between dT to T s in Equation ( 6), which is reflected in different estimates of a and b by T-SEBAL and SEBAL according to Equations ( 13) and ( 14), shown in Table 2.That means, T-SEBAL and SEBAL differ not only in their ways of determining T s,hot and T s,cold , but also in the calculation of (R n − G).
With the linear relationship determined, we can further estimate H for each pixel.Figure 8 shows the difference between SEBAL estimates and T-SEBAL estimates over the whole study area, and the scatter plots of SEBAL estimates versus T-SEBAL estimates for DOY75, DOY83, DOY267, and DOY345, respectively.The scatter plots show that the H estimates by T-SEBAL are generally well correlated with those by SEBAL, and the correlation coefficients between the two are 0.85, 0.6, 0.71, and 0.73 in DOY75, DOY183, DOY267, and DOY345, respectively.However, there are some biases between the two, that is, H estimates by T-SEBAL are on average ~117 W•m −2 , ~3 W•m −2 and ~52 W•m −2 lower in DOY345, DOY267 and DOY75, and 53 W•m −2 higher in DOY 183, than those by SEBAL.Such biases can be explained by the differences in the calculation of dT.In DOY345, b estimated by T-SEBAL is significantly lower than that by SEBAL while values of T s at the cold point given by two models have slight differences, which leads to significantly higher dT estimates by SEBAL than by T-SEBAL, especially in the areas with high T s values.As a consequence, H estimated by SEBAL is significantly larger than by T-SEBAL in DOY345.The situations are more or less similar for DOY 267 and DOY75.As to DOY 183, both T s,1 and b by T-SEBAL are smaller than SEBAL, leading to overall larger dT and consequently larger H by T-SEBAL than by SEBAL.
With H estimated for each pixel, we can easily derive LE as an energy residual with Equation (1).The scatter plots of H and LE estimates with T-SEBAL and SEBAL against observations at Lucky Hills and Kendall for 15 dates are displayed in Figure 9. From Figure 9 we see that H is underestimated by both T-SEBAL and SEBAL in general.However, in terms of RMSE (Root mean square error), MAE (Mean absolute error), MBE (Mean bias error), T-SEBAL generally performs better than SEBAL in the estimation of LE, as shown in Table 3 [8,[44][45][46][47][48].
Besides the results of LE estimation with T-SEBAL and SEBAL in the present study, Table 3 also lists some published results in the literature for Walnut Gulch watershed for comparison.It is shown that, compared with other published results, the T-SEBAL results in the present study are better in general, especially considering the fact that in T-SEBAL all energy-related variables are estimated whereas observational data are used in many published results, e.g., Norman et al. (1995) [8] and Moran et al. (1996) [44] used observed R n , Anderson et al. (1997) [45], Kustas et al. (1998) [46], and Li et al. (2008) [47]     , where P and Q are the estimation and the observation respectively.

Daily Actual ET
As the residual energy, LE is calculated according to the Equation ( 1) after estimated H, and it could be converted to instantaneous ET (ET inst ) by the latent heat of vaporization.
where ET inst is the instantaneous ET a (m•s −1 ).ρ w is density of water (~1000 kg•m −3 ).l e is the latent heat of vaporization(J•kg −1 ) calculated as: ( ) where t is the time in the day (beginning at sunrise).N is the time period between sunrise and sunset, which is calculated by the equation suggested by Jackson [49].
For the 15 dates studied in the present study, in comparison with the observations, the MAE, MBE and RMSE of daily ET a estimates with T-SEBAL are 0.42 mm, 0.1 mm, and 0.52 mm, respectively, which are all smaller than the results with SEBAL, that is, 0.58 mm, 0.49 mm, and 0.68 mm, respectively.
Distribution of daily actual ET given by both T-SEBAL and SEBAL in DOY75, DOY 183, DOY 345, and DOY 267 are plotted in Figure 10 (left and middle columns).In addition, the differences between the T-SEBAL estimates and SEBAL estimates are also plotted in Figure 10 (right column).We may find that, as expected, ET a estimates with T-SEBAL and those with SEBAL have similar spatial pattern.However, there are also some quite significant differences.Firstly, ET a estimates by T-SEBAL has an overall negative bias compared with those given by SEBAL in DOY183, and overall positive biases in DOY 75 and DOY345.Secondly, there exist clear differences in ET a estimates in the Dragoon Mountain area (on the upper right side).SEBAL gives high ET a estimates (higher than the average throughout the scene) in the mountain area in DOY75, DOY183, and DOY267, whereas T-SEBAL shows lower ET a estimates.As we know, the Dragoon Mountain area is generally rocky and poorly vegetated, thereby should have lower ET than surrounding areas, we believe that ET a estimates by T-SEBAL are more acceptable than those by SEBAL.
We further compared the average ET a estimate over the WGEW by the two models with average antecedent precipitation in the previous 10 days in the 15 dates we studied, shown in Figure 11.Clearly, ET a estimates by T-SEBAL have a better relationship with antecedent precipitation than those by SEBAL, indicating that ET a estimates by T-SEBAL are physically sounder.

Sensitivity Analysis
To assess the robustness of the T-SEBAL model, sensitivity analysis is conducted to investigate model sensitivity to variations of several important parameters and variables, including: land surface temperature T s that is derived from MODIS product; air temperature T a , relative humidity U, and wind speed u that are observed on the ground; minimum canopy resistance r cm , the ratio of G/R n , and albedos of four vertices that are set empirically or by calibration.
As sensible heat flux H is the essential output of the T-SEBAL model, the sensitivity analysis is conducted based on the change of H estimates in response of input changes.The sensitivity S i of T-SEBAL to an input i is defined as [25]: where, H 0 , H + and H -are the sensible heat flux estimated by T-SEBAL when the input equals its reference value i 0 , is increased and is decreased, respectively, with reference values used for all other inputs.
In the present study, H ± and H 0 are the mean values of the whole region as shown in Figure 4.The following sensitivity analysis results are based on the estimation of H for 15 dates.A deviation of 5 K with a step of 0.5 K is used to test the sensitivity to T s and T a , whereas a 25% deviation with a step of 2.5% change is used to test the sensitivity to other parameters.The results are presented in Figure 12.

Sensitivity to T s and T a
Figure 12a shows that T s is positively correlated with H, while T a is negatively correlated, and both correlations are approximately linear.Precisely, a 5 K increase in T s and T a result in a 16.9% increase and 27.8% decrease in the estimate of H averagely, whereas a 5 K decrease in T s and T a lead to 31.4% decrease and 17.9% increase in the estimate of H.In comparison, for SEBAL in a semi-humid region [28], a 5 K decrease or increase in T s can cause 59.8% decrease or 72.8% increase in the H estimate, which indicates more sensitivity of SEBAL than T-SEBAL.
Because the calculation of the land surface temperature of extreme points (i.e., T s,n , n = 1, 2, 3, 4) are calculated by solving quartic equations of surface temperature at extreme conditions, therefore, the changes of pixel T s do not affect the calculation of T s,hot and T s,cold , but affect the distance to hot point which could increase or decrease dT of pixels, and subsequently affect the estimation of H.
In contrast, T a is positively correlated with T s,n as shown in Equation (8a-8d).Therefore, an increase in T a would lead to the increase of T s,n and (R n −G) 4 , and consequently put the pixel in the constructed trapezoid toward the cold edge, thereby decrease the estimate of H.

Sensitivity to Wind Speed and Relative Humidity
Wind speed u is positively related to the H estimate, because the increase in u could result in the increase in friction velocity and the decrease in r a and T s,n , consequently put the T s ~EVI data points to the warm edge, resulting in the increase in H estimate.As shown in Figure 12b, a 25% increase or decrease in u could lead to 8.6% increase or 11.5% decrease in the H estimate on average.
Relative humidity (U) has smaller impacts on H estimation than wind speed.A 25% increase in U would result in only 0.7% decrease in H.Such an effect happens because the increase in U leads to the increase in ε a and the decrease in VPD, which results in the increase in T s,1 and T s,4 , consequently decrease the dT and H.

Sensitivity to Albedos at Four Vertices of T s ~VI Trapezoid
As the albedos of four vertices are set differently, thereby, the sensitivity to albedo at each point is investigated separately.As shown in Figure 12c point 4 is more sensitive to the change of albedo, with a 25% increase or decrease causing -5.9% or 4.6% deviation in the H estimate.A 25% perturbation at point 1 only causes ~0.3% variation in H.Meanwhile, H is free of any influence by albedo variations at point 2 and point 3.

Sensitivity to Minimum and Maximum Canopy Resistance
The sensitivity to minimum canopy resistance (r cm ) is illustrated in Figure 10d, which indicates the negative correlation between r cm and H estimates.That is, a 25% increase or decrease in r cm could cause 1.3% decrease or 1.5% increase in the H estimate.The cause of such an effect is that, an increase in r cm could result in an increase in the T s,1 estimate according to the Equation (8a), which would reduce the temperature difference dT, thereby, reduce the estimate of H.
The maximum canopy resistance (r cx ) barely influences H estimation on average as shown in Figure 12d.r cx is a parameter used to calculate T s,2 , which slightly affects the slope of hot edge and cause marginal influence on H estimation.

Sensitivity to G/R n at Four Vertices of T s ~VI Trapezoid
The (G/R n ) n ratio r n (n = 1, 2, 3, 4, representing 4 vertices in Figure 1) are important parameters for estimating the T s,n of the trapezoidal space.Figure 12e illustrates that, for point 4 (dry bare soil), a 25% increase or decrease in (G/R n ) 4 could result in 6.3% decrease or 4.7% increase in the estimate of H.For point 1 (wet full vegetation), H estimation is insensitive to (G/R n ) 1 , only 0.1% deviation with a 25% perturbation.Meanwhile, H is not influenced by the variation of (G/R n ) 2 and (G/R n ) 3 .−1 for Canopy and Bare Soil Surface kB −1 is the crucial variable for the estimation of z 0h .In Figure 12f, it is demonstrated that values of kB −1 at point 4 (dry bare soil) are negatively related to H estimates, that is, a 25% increase in kB −1 for dry bare soil surface causes 5.8% decrease in the estimate of H.While the kB −1 for wet full vegetation surface shows little inference on H estimation.

Comparison between T-SEBAL and SEBAL
T-SEBAL is designed to mitigate the subjectiveness and uncertainty in determining the extreme pixels in SEBAL.It intended to improve the performance of the standard SEBAL in the following two major aspects: (1) Avoid the subjective selection of hot point and cold point by calculating hot/cold point temperatures theoretically for the T s ~VI trapezoid based on energy balance for the extreme dry/wet cases and replacing hot point with the dry soil vertex and cold point with the wet full vegetation vertex of the T s ~VI trapezoid.(2) Avoid the requirement of the existence of extreme cold and hot pixels throughout the image scene and consequently avoid the dependence on domain size, by constructing the T s ~VI trapezoid on a pixel basis.
According to the design of T-SEBAL model, T s,4 , T s,1 , and (R n −G) 4 are unique for each pixel with T-SEBAL, whereas T s,hot , T s,cold , and (R n −G) hot are constant throughout the whole image scene with SEBAL.Therefore, the forms of relationships between dT and T s are different between T-SEBAL and SEBAL (see Figure 13), i.e., T-SEBAL sets up a unique linear relationship for each pixel and consequently the same T s in different pixels may correspond to different dT, whereas SEBAL builds a unified linear relationship for the whole image and therefore T s and dT have a linear one-to-one relationship.The major reason that linear relationships between dT and T s of different pixels are not consistent throughout the image scene with T-SEBAL is the spatial variation in T a .If we use the average T a of whole image to simulate the linear relationship, we can see that all pixels have similar linear relationships between dT and T s , as shown in Figure 13 (dT in T-SEBAL using average T a ).In addition, the spatial variation in T a affect not only the linear relationship between dT and T s , but also the correlation between H estimates by T-SEBAL and those by SEBAL.The correlation between H estimation by T-SEBAL when average T a is used and H estimates by SEBAL is much higher than the results presented in Figure 8, with r up to 1.0, 0.99, 0.98, and 1.0 for DOY75, DOY183, DOY267, and DOY345, respectively.Therefore, considering the spatial variation in T a when defining temperatures at hot/cold points has significant impact on the estimation of H and subsequently ET a , and the capability of doing so in T-SEBAL at a pixel basis is one of the reasons that lead to the difference in H estimation by T-SEBAL and SEBAL.

Comparison with Other Modifications of SEBAL
Considerable efforts have been made in the hydrology remote sensing community on improving the SEBAL model.One important variant of SEBAL is the METRIC model which departs from the SEBAL model in its use of weather-based reference ET to establish energy balance conditions at a "cold" pixel, but its application still depends on the manual selection of appropriate hot and cold pixels, which is quite different from the approach taken in the present study.A modification of SEBAL quite similar to T-SEBAL is M-SEBAL [19].However, besides the difference in the formulation for constructing T s ~VI trapezoid, there are significant differences in the way of defining the hot point temperature and the way of constructing the cold edge.
The hot point temperature in M-SEBAL is set on the warm edge at the point where the vegetation fraction (f c ) equals that of the pixel (f c,i ).That means, instead of setting up the linear relationship relating the vertical temperature gradient dT to T s on the basis of the entire image scene as in SEBAL, the dT~T s relationship is set on the basis of a section of the T s ~fc trapezoid (equivalent to the T s ~VI trapezoid) in M-SEBAL.The underlying assumption for setting dT~T s relationship in this way is that, LE is zero on the warm edge and H is zero on the cold edge.In T-SEBAL, however, we assume that LE is zero only on the point of bare soil with the largest water stress, but not on the point of full vegetation with the largest water stress, as it is known that the stomatal transpiration accounts for ~90% of total transpiration [50], which means that even in the case of full stomatal closure there would still be some water loss.Therefore, in the T-SEBAL, the temperature at the hot point is set the same as the temperate at the point of bare soil with the largest water stress, and the dT~T s relationship is set on the basis of the whole T s ~VI trapezoid.
The cold edge in M-SEBAL is taken as the air temperature at the image time.But as shown in Figure 7, where the air temperature is used as the cold point temperature mostly (14 out of 15 dates) in SEBAL, the air temperature is consistently higher than the calculated cold point temperatures in T-SEBAL, just contrary to the statement of Long et al., who argued that the use of T a to substitute for the theoretical cold edge would underestimate the surface temperatures over sparsely vegetated surfaces without water stress [19].Such a temperature difference at the cold edge would result in the overestimation of LE for the SEBAL model, as shown in Figure 9.

Limitations of T-SEBAL
There are several aspects of T-SEBAL needed to be worked on further, including: (1) The process of extrapolating/interpolating meteorological variables based on data at a limited number of observation sites may introduce some errors.As shown in the sensitivity analysis, the model is sensitive to air temperature and wind speed.However, only the extrapolation of air temperature is treated considering the effects of topography and LAI, whereas the wind speed is simply taken as an average of the three meteorological observation sites in our cased study area.Furthermore, although the effects of topography and LAI are considered, the accuracy of air temperature extrapolation/interpolation is still not always satisfying.(2) The parameterization under extreme conditions for constructing the T s ~VI trapezoid need to be improved, such as ratios of G/R n at extreme conditions.Compared with the G estimation with Equation (3), the method suggested by Equation (4) shows better fitness, therefore, the Equation ( 4) is used to estimate G in both SEBAL and T-SEBAL in this paper, and the G/R n at full vegetation and bare soil surface are set as 0.05 and 0.28, respectively.However, those rates do not consider the effect of soil moisture, whereas many experiments have shown the differences in G/R n for dry and wet conditions (e.g., [30]), therefore, at point 3 (wet bare soil) and point 4 (dry bare soil) we set the ratio of G/R n to be 0.15 and 0.35, respectively.In addition, several other intermediate parameters (e.g., albedo, ε s , r cx and r cm ) are also set using empirical values.The parameterization of the model for other types of ecosystems need to be further evaluated.(3) There are many formulas available for calculating aerodynamic resistance r a with different assumptions and complexity, as well as for the estimation of the dimensionless parameter kB −1 which is used to define the relationship between momentum and heat roughness length.We found that different choice of the formula for either r a or kB −1 may have significant impacts on ET a estimation results.In addition, as the land surface radiometric temperature T s , instead of the aerodynamic temperature T aero , is used in the process of constructing T s ~VI trapezoid, the kB −1 has to be adjusted, and actually has little to do with the kB −1 parameter when T aero is used [51].To avoid any confusion, Troufleau et al. (1997) even suggested using the term "radiometric kB −1 " as a contrast to the theoretical aerodynamic kB −1 [51].Norman and Becker (1995) also recommended the use of a different notation to avoid unnecessary confusion between the theoretical and empirically adjusted excess resistance [52].Verhoef et al. (1997) [37] even conclude that the concept of kB −1 is questionable, but on the other hand, using the bulk transfer equation to calculate sensible heat from remotely sensed surface temperature requires the use of kB −1 [35].As a consequence of using the adjusted "radiometric kB −1 ", the calculation of aerodynamic resistance r a ' is different from the theoretical aerodynamic resistance r a .That is why we used two formulas for calculating aerodynamic resistance However, the theoretical basis of doing so need further investigation.

Conclusions
The surface energy balance algorithm for Land (SEBAL) is widely used to estimate actual evapotranspiration (ET a ).Two major disadvantages of SEBAL model are its requirement of the existence of extreme cold and hot pixels in the study area, and the subjectiveness in selecting the two extreme pixels.A modified version of SEBAL model, T-SEBAL, is proposed in the present study to determine the extreme cold/hot points using a parameterized procedure based on the theoretical trapezoidal relationship between land surface temperature (T s ) and Enhanced vegetation index (EVI), so as to reduce the uncertainty in selecting hot and cold pixels and make the SEBAL model more robust.
The T-SEBAL model avoids the subjective selection of hot point and cold point in SEBAL model, does not need the existence of extreme cold and hot pixels throughout the image scene, and therefore has no dependence on domain size.In addition, as T-SEBAL sets up a unique linear relationship between dT (difference between aerodynamic temperature and air temperature) and T s on a pixel basis, whereas SEBAL builds a unified linear relationship for the whole image, therefore T-SEBAL model has better capability to characterize the spatial variability of complex landscapes.The performances of the classical SEBAL model and T-SEBAL are compared for estimating ET a for a semi-arid area, the Walnut Gulch Experimental Watershed, in USA, for 15  with SEBAL, respectively.Additionally, for the 15 dates studied in the present study, the MAE, MBE, and RMSE of daily ET a estimates with T-SEBAL are 0.42 mm, 0.1 mm, and 0.52 mm, respectively, which are all smaller than the results of SEBAL, that is, 0.58 mm, 0.49 mm, and 0.68 mm, respectively.In addition, ET a estimates with T-SEBAL have better physical soundness than with SEBAL in the terms of the relationship between ET a values with geomorphologic locations and the antecedent precipitation.The results indicate that, with the advantage of calculating hot/cold point temperatures theoretically on a pixel basis, T-SEBAL has better applicability for complex landscapes than standard SEBAL model.
There are several aspects of T-SEBAL needed to be worked on further, including the way of extrapolating meteorological variables based on data at a limited number of observation sites, the parameterization under extreme conditions for constructing the T s ~VI trapezoid, and the formulation of kB −1 and the aerodynamic resistance r a .
Modifying SEBAL Using the Trapezoidal Relationship between Land Surface Temperature and Vegetation Index 2.1.SEBAL Model SEBAL model estimates the surface net radiation flux R n (W•m −2 ), the soil heat flux G (W•m −2 ), and the sensible heat flux H (W•m −2 ) firstly, then estimates the latent heat flux LE (W•m −2 ) as an energy residual by:

Figure 1 .
Figure 1.Structure of the T s ~VI trapezoidal space.

Figure 2 .
Figure 2. Procedure for estimating sensible heat flux H with T-SEBAL.

Figure 3 .
Figure 3. Adjustment of T s when the (EVI, T s ) point is above the warm-edge of the calculated T s ~VI trapezoid.

Figure 4 .
Figure 4. Location of the study area, flux stations and meteorological stations.

Figure 6 .
Figure 6.Box-and-whisker plots of calculated T s at four vertices of T s ~EVI trapezoids of four dates in 2004: (a) DOY75; (b) DOY183; (c) DOY267; (d) DOY345.(For each boxplot, the bottom and top of the box are the first and third quartiles, the band inside the box is the median, the lower whisker is the lowest datum still within 1.5 IQR (interquartile range) of the lower quartile, and the upper whisker is the highest datum still within 1.5 IQR of the upper quartile).

Figure 7 .
Figure 7.Comparison of the land surface temperature at hot points and cold points for SEBAL and T-SEBAL.
used observed R s to estimate LE.

Figure 8 .
Figure 8. Differences between H estimated by T-SEBAL and H estimated by SEBAL (T-SEBAL estimates minus SEBAL estimates) (left panels) and the scatter plots T-SEBAL estimates versus SEBAL estimates (right panels) in four dates.

Figure 9 .
Figure 9.Comparison of T-SEBAL and SEBAL estimated (a) H and (b) LE versus observations.
Jackson et al. (1983) [49] related diurnal trend of the ratio of ET inst to daily value (ET d ) to a sine function:

Figure 10 .
Figure 10.Spatial distribution of daily ET a (mm) estimated by T-SEBAL (left column) and SEBAL (middle column) and their differences (T-SEBAL estimates minus SEBAL estimates) (right column) in four dates.

Figure 11 .
Figure 11.Comparison of the average ET a estimated by T-SEBAL and SEBAL vs. 10-day antecedent precipitation over the WGEW.
, in comparison with other points, the H estimation at Relative change of G/R n (%) Relative change of H (%)
dates in 2004.The results of latent heat flux LE estimation demonstrate that, T-SEBAL model outperformed the classical SEBAL model in terms of root mean square error (RMSE), mean absolute error (MAE) and mean bias error (MBE), that is, 49.6 W•m −2 , 40.4 W•m −2 and 14.4 W•m −2 with T-SEBAL vs. 57.3W•m−2 , 48.9 W•m −2 and 45.8 W•m −2

Table 1 .
MODIS product data used.

Table 2 .
Variables and parameters at extreme points given by Surface Energy Balance Algorithm for Land (SEBAL) and T-SEBAL for 15 dates.

Table 3 .
Comparison of different LE estimation accuracies for Walnut Gulch watershed (W•m −2 ).
n , G, VPD, and T a , etc. Note: LH and KD represent the two flux stations Lucky Hills and Kendall in Walnut Gulch watershed, respectively.a. MAE is given by ∑ |P Q | ; b.MBE is given by