Combining the trapezoidal relationship between land surface temperature and vegetation index with the Priestley-Taylor equation to estimate evapotranspiration

Evapotranspiration (ET) plays an important role in the hydrological cycle. A method of combining the Priestley-Taylor (P-T) equation with a trapezoidal space between land surface temperature (Ts) and enhanced vegetation index (EVI) is proposed based on the principle of energy balance. Generally, this method is divided into three major parts: (1) construct the Ts versus EVI (Ts-VI) trapezoidal space for calculating the Ts at four extreme conditions (i.e. well-watered vegetation, water-stressed vegetation, saturated bare soil and dry bare soil); (2) calculate the P-T coefficient for each pixel according to the position of the observed (EVI, Ts) point in the trapezoid space; (3) calculate actual ET of the pixel using the P-T equation. The method is validated using Landsat-8 images and ground-observed data for a semi-humid area in China. The result shows that the ET estimates match the observations well, which indicates the effectiveness the proposed method here.


INTRODUCTION
Evapotranspiration (ET) plays an important role in hydrological cycle.Numerous methods have been proposed to estimate actual ET.Among those models, the land surface temperature−vegetation index feature space related methods were applied widely (Jiang and Islam, 1999;Venturini et al., 2004).Price (1990) found that the scatterplot of remotely sensed surface temperature (Ts) was related to vegetation index (VI) approximately with a triangular distribution and such a spatial context could be used to infer regional scale ET.Basically, the triangular feature space was plotted by utilizing substantial pixels representing complete dry to wet status, and full vegetation coverage containing an upper declining envelope (dry edge) and a lower nearly horizontal envelope (wet edge).The method is insensitive to initial atmospheric and surface conditions, net radiation and atmospheric correction and has been used widely.However, the triangular space method has a major limitation in constructing the feature space, as the warm and cold edges are mainly determined by fitting a linear equation or set empirically, which requires a large number of pixels over an area with a wide range of soil wetness and fractional vegetation coverage.Moran et al. (1994) estimated Water Deficit Index (WDI) for each pixel given the vegetation cover and dry/wet condition based on the principle of energy balance and consideration of a trapezoidal relationship (referred to VIT trapezoid hereafter) between the difference of Ts and air temperature (Ta), Ts−Ta and fractional vegetation cover.Generally, VIT trapezoidal space resolved the two previous problems of triangular space.However, Ts−Ta cannot be observed by remote sensing techniques.Meanwhile, in the construction of the VIT trapezoid, the dynamic interaction between Ts−Ta and Rn (i.e. the observed Rn was used directly for Ts estimation of four vertexes) was not accounted, which may bring large bias because the Rn deviation could reach up to 20% between bare land and vegetation cover surface (Daughtry, 1990).Therefore, Wang et al. (2011) developed an algorithm to simplify the VIT trapezoid to a Ts-VI trapezoidal space with an iterative procedure to consider the dynamic interactions between Ts−Ta and Rn.
The objectives of this paper are to propose a method of combining the land surface temperature versus enhanced vegetation index (Ts-VI) trapezoid space for each pixel and the Priestley-Taylor (P-T) equation to estimate regional actual ET.

TS-VI TRAPEZOIDAL SPACE
The core equation of VIT space is given by: where, Ts is land surface temperature (K); Ta is air temperature (K); Rn is net radiation (W m -2 ), G is soil heat flux (W m -2 ) estimated as a fraction of Rn; VPD is the vapour pressure deficit of the air (hPa); Δ is the slope of saturated vapour pressure to air temperature; γ is the air humidity constant (hPa•C −1 ); rc is the canopy resistance (m•s −1 ); ra is the aerodynamic resistance (m s -1 ).Wang et al. (2011) simplified the VIT space as Ts-VI trapezoidal space (Fig. 1) by moving Ta in equation ( 1) to the right side and constructing an iterative process to describe dynamic interaction between Ts and Rn, ra, H, etc.That is, this process repeatedly adjusts the parameters (Rn, ra and H etc.) to match the change of Ts.Then, ra, Rn, H and Ts under extreme conditions were calculated.Ts at four Ts-VI trapezoid vertexes are calculated from the equations proposed by Wang et al. (2011).ra in equation ( 1) is calculated by: ( ) ( ) where, k is the von Karman constant (=0.4); u (m s -1 ) is wind speed observed at reference height, z (2 m); zom and zoh are the momentum roughness length and heat roughness length; ψh and ψm are the stability correction for heat and momentum; d is the displacement length (m).
The procedure for constructing Ts-VI trapezoidal space on a pixel basis is summarized in Fig. 2 (see the process in dash box), which is divided into three major steps: (a) Estimate the initial value of ra using meteorological data without considering ψh and ψm; (b) Iteratively solve the quartic equation of Ts), calculate H, Monin-Obukhov length L and ra sequentially for each vertex until the ra is stable, i.e. difference between two adjacent calculations of ra is smaller than 5%; (c) verify the reasonableness of Ts-VI trapezoid.If an observed (EVI, Ts) point falls outside the hot/cold edges, an adjustment is made to make the point be contained within the Ts-VI trapezoid.
After the process to construct the Ts-VI trapezoidal space for each pixel shown in Fig. 2, the Priestley-Taylor (P-T) coefficients (α) at four vertexes are estimated using the calculated Rn, G and H at each vertex (see the process in solid box in Fig. 2), given by: ( ) ( )( ) where the subscript i represents the parameters related to four vertexes.
Two linear interpolations were conducted to estimate a pixel's α at a given observed Ts and EVI position (EVI, Ts).First, on the basis of the assumption that α and Ts are linearly related to EVI at warm and cold edges, the extreme α (αmax, αmin) and Ts ( Ts,min, Ts,max) under a given EVI value are estimated, i.e., α at point A and point B in Fig. 1, respectively.Secondly, on the linear relation assumption between α and Ts at a given EVI, α for an observed (EVI, Ts) point is where the subscript real represents the value at point C in Fig. 1.With αreal estimated, LE of the pixel with a given (EVI, Ts) can be estimated with P-T Equation ( )( ) Compared with the VIT by Moran et al. (1994), the Ts-VI method builds the dynamic interactive relationship between Ts and Rn, ra, H, ψh and ψm under four extremes.Therefore, the calculated Rn, H, ra, Ts are closely related to extreme conditions, which could reflect the interaction of parameters within atmosphere−vegetation−soil surface system.
Meanwhile, it should be noted that ET at the warm edge is not set as 0, which is different from the triangular method and VIT method.Many studies have shown that, even in the case of full stomatal closure, there is still some water loss because the stomatal transpiration accounts for ~90% of total transpiration (Levitt, 1972) while the remainder is derived from the cuticle.

STUDY AREA AND DATA
The Wudaogou Hydrological Experimental Station (117°21′E, 33°9′N) is located in the middle of the Huaihe River basin in China.Wudaogou has a semi-humid climate, with an annual average temperature around 14.6° and annual average precipitation 840mm.
One tower equipped with an eddy covariance system and five-layer meteorological gradient observations was installed in Wudaogou station in June 2012.Shortwave and longwave radiation, soil heat flux, latent heat flux and sensible heat flux are observed at half-hour intervals; air temperature, relative humidity and wind speed are observed at 2, 5, 10, 15, and 20 m height.
Landsat-8 images were used.Land surface radiance temperature (Ts) is estimated using the single channel method proposed by Jiménez-Muñoz et al. (2014) due to stray light coming into the TIRS detector from outside the planned field of view arousing inaccuracy of Band 11.

APPLICATION OF TS-VI TRAPEZOID METHOD AND VALIDATION
The calculated Rn and G were compared to observations on three dates (Table 1).The root mean squared error (RMSE) and mean absolute error (MAE) of Rn are 43.2 and 42 W m -2 , respectively, The iterative process integrate P-T equation with Ts-VI trapezoidal space calculate α at four vertexes according to Eq. ( 3) .and the overestimate is 10 W m -2 .As for G, estimated value indicates good agreement with observations showing no obvious bias, with the RMSE and MAE of 20.7 W m -2 and 16.7 W m -2 .
According the procedure described in Fig. 2, we built Ts-VI trapezoidal space for each pixel.The calculated Ts at the Ts-EVI trapezoid vertices for all the pixels of the study area are plotted in the form of box-and-whisker plots in Fig. 3 to illustrate the reasonability of the trapezoids.As shown in Fig. 3, the pixels for all the three dates are closer to the cold edge (the line between point 3 and point 1) than the warm edge, indicating more energy is used to evaporate water than to heat the air.
Intermediate variables, Rn, H and G at four extremes are calculated in the process of building Ts-VI feature space, which are used to estimated P-T coefficients (α) by Eq. ( 3) and energy balance equation.At point 1, average of α of 3 dates is 1.61, higher than 1.26 due to the value of 1.85 estimated in DOY326 in year 2013.Brutsaert (1982) pointed out that α values at potential ET conditions may increase to approximately 1.6-1.8when advective conditions are experienced.
α ranges from 0.09 to 0.24 with an average of 0.16 at point 2 consistent with the cognition that little transpiration happens even under the full water-stressed condition.α are averaged at 1.31 and 0 at point 3 and point 4, respectively.
After estimating α at four vertexes, the α of the (Ts, EVI) pair for each given pixel is calculated with the linear interpolation equation (4) described in Section 2. Consequently, LE is calculated by equation (3).Comparing with the observations, the RMSE and MAE are 38.7 and 33.2 W m -2 , respectively.The value of the Mean Bias Error (MBE) indicates that LE is overestimated by around 13.6 W m -2 .
As Fig. 4 shows, the differences between water body and land are more obvious in DOY246 and DOY121 than in DOY326.For DOY246, mean LE of water body is 463W m -2 ranging from 400 to 550 W m -2 with a standard deviation of 32.0 W m -2 , while the average and standard deviation of LE at the land surface are 265.6 and 57.0 W m -2 , respectively, ranging from 58 to 400 W m -2 , corresponding to the two peaks in the histogram (middle column), i.e. one around 270 W m -2 (mainly land crop) and the other around 490 W m -2 (water body).The estimated LE at Wudaogou observation tower is close to that of water body, although the land cover is dominated by wheat in winter and spring, and maize in summer and autumn.Comparing the spatial distribution of LE and Ts shown in Fig. 4, we can see that they have strong consistency.

CONCLUSION
A method of estimating regional actual ET is proposed by combining the land surface temperature versus enhanced vegetation index (Ts-VI) trapezoid space and the Priestley-Taylor equation on a pixel-basis.The method was applied to a semi-humid area with satisfying results.The prominent feature of the method is the use of a physically-based trapezoidal framework to account for the effect of variations in Ts and EVI on ET estimation, which reduces the dependence on the existence of the full range of dry to wet status and vegetation coverage.Meanwhile, minimum meteorological parameters are required, which reduces the uncertainties on spatial interpolation.However, some aspects need to be improved and be worked on further, including the accuracy of extrapolating/interpolating meteorological variables based on data at a limited number of observation sites, and the rationality of Ts-VI trapezoidal space.

Fig. 1
Fig. 1 Structure of the Ts-VI trapezoidal space with four vertexes representing four extreme conditions (i.e.well-watered vegetation, water-stressed vegetation, saturated bare soil and dry bare soil).

Fig. 2
Fig. 2 Procedure for estimating sensible heat flux H with T-SEBAL.calculated, i.e. point C in Fig. 1.The integrated equation for calculating α of a given (EVI, Ts) is given by: , , , , for Ts calculate Rn, H initial ra, assuming ψh = ψm = 0 calculate Monin-Obukhov length L calculate ψh, ψm, and ra Adjust Ts Ts, Rn, G, H, at four vertexes N r stable?Space reasonable ?

Fig. 3
Fig. 3 Box-and-whisker plots of calculated Ts at four vertices of Ts-EVI trapezoids of three dates: (a) DOY246 in 2013; (b) DOY326 in 2013; (c) DOY121 in 2014.(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).

Fig. 4
Fig. 4 Spatial distribution of LE (left column), histogram of LE (middle column) and spatial distribution of ts (right column) in three dates: DOY246 in 2013, 326 in 2013, 121 in 2014; the pentagram represents the location of the observation tower.

Table 1
Statistics on estimations versus ground observations.