A New Contextual Parameterization of Evaporative Fraction to Reduce the Reliance of the T s − VI Triangle Method on the Dry Edge

In this study, a new parameterization scheme of evaporative fraction (EF) was developed from the contextual information of remotely sensed radiative surface temperature (Ts) and vegetation index (VI). In the traditional Ts − VI triangle methods, the Priestley-Taylor parameter ∅ of each pixel was interpolated for each VI interval; in our proposed new parameterization scheme (NPS), it was performed for each isopiestic line of soil surface moisture. Specifically, ∅ of mixed pixels was determined as the weighted-average value of bare soil ∅ and full-cover vegetation ∅. The maximum Ts of bare soil (Tsmax) is the sole parameter needed as the constraint of the dry edge. This has not only bypassed the task involved in the determination of the maximum Ts of fully vegetated surface (Tcmax), but also made it possible to reduce the reliance of the Ts −VI triangle methods on the determination of the dry edge. Ground-based measurements taken during 21 days in 2004 were used to validate the EF retrievals. Results show that the accuracy achieved by the NPS is comparable to that achieved by the traditional Ts −VI triangle methods. Therefore, the simplicity of the proposed new parameterization scheme does not compromise its accuracy in monitoring EF.


Introduction
Evapotranspiration (ET), defined as the sum of water lost to the atmosphere from sources such as the soil, canopy interception and water bodies through evaporation and plant transpiration, is a primary process driving the energy and water exchange between the atmosphere, hydrosphere and biosphere.Accurate knowledge of spatially distributed ET is essential for a wide range of applications including water resources management, hydrometeorological predictions and ecological applications [1][2][3][4].As opposed to conventional observations of evapotranspiration that only represent field scale, satellite remote sensing provides an unprecedented global coverage of surface observations such as radiometric surface temperature (T s ), albedo and vegetation indices (VI), which makes it recognized as the only viable means to map regional patterns of ET at the earth's surface in a globally consistent and economically feasible manner [5][6][7].As a result, a large number of ET models have been developed from remote sensing observations over the past few decades.Overviews of these models have been provided by a number of authors [8][9][10][11].However, the determination of ET depends on so many variables that not all of them can be retrieved from remote sensing observations.As a result, most of these ET models have either adopted meteorological observations as ancillary or assumed some parameters to be constant over certain scales.This has not only limited the application of these models but also reduced the accuracy of ET estimates, especially in areas where the meteorological observations are scarce.Therefore, accurate characterization of ET distribution based on remote sensing with little or no ground observations as ancillary is still a challenge task [2,[12][13][14].
A promising method to meet this challenge is to estimate evapotranspiration based on the evaporative fraction (EF).EF defined as the ratio of latent heat flux to available total energy [15] can be estimated reasonably using only remote sensing data as inputs.Over the past few decades, a large number of techniques have been proposed to estimate EF [14,[16][17][18][19][20], and a direct estimation of EF has been shown to work well with various satellite remote sensing datasets [21][22][23].Among all these techniques, the surface temperature-vegetation index (T s − V I) triangle or trapezoid method developed by Jiang and Islam [21,24] and improved by Jiang and Islam [25] has widely been used throughout the world.It derives spatially distributed EF by using a contextual interpretation of a relationship between easily measured surface variables-radiometric surface temperature (T s ) and normalized difference vegetation index (NDVI) [26].Numerous studies have documented and interpreted the advantages of the T s − V I triangle method, and its applicability has been demonstrated in various parts of the world [13,23,[26][27][28][29][30].
The success of this method on the estimation of EF depends mainly on the correct choice of the boundaries of the T s − V I feature space, especially the dry edge that defines the maximum water stress.The development of this method in recent studies is mainly focused on the determination of the dry edge [13,[30][31][32][33][34].In general, the dry edge can be retrieved either from theoretical or empirical approach.The contrast of these two approaches is presented in Table 1.It is clear that both of these two approaches have their advantages and disadvantages.Although the observed dry edge determined by the regression method can be retrieved exclusively from remote sensing data, it is usually lower than the theoretical dry edge, and is not assigned the minimum ET [23].Moreover, the subjectivity involved would also reduce the accuracy of the EF estimates.The theoretical dry edge determined by land surface energy balance principle can overcome generally all the limitations of the empirical approach.However, a large number of parameters are needed in the determination of the theoretical boundaries.This to some extent contradicts with the original intention of the T s − V I method, because triangle methods were designed in general to avoid the need for ancillary surface data and complex parameterization of aerodynamic and surface resistances for water and heat transfer [18,28,30].Therefore, there are challenges in the determination of both the observed and theoretical dry edge, which would limit the application of the T s − V I triangle method in arid and semi-arid areas or areas with limit ancillary surface data.Development of a parameterization scheme of EF that can reduce the reliance of the T s − V I method on the dry edge is a promising method to overcome this challenge.
Table 1.Contrast of empirical and theoretical approaches for the determination of dry edge.

Advantages
It is simple and can be performed based entirely on remote sensing data.
It is performed through theoretical derivation and can remove the subjectivity involved; the theoretical dry edge determined represents the maximum water stress.

Disadvantages
Establishment of regression models involves subjectivity; the observed dry edge determined is not assigned the maximum water stress.
A large number of parameters are needed such as air temperature, near surface vapor pressure, net radiation, aerodynamic resistance.
Besides the T s − V I triangle method, there are some other significant studies performed to investigate the relationships between remotely sensed T s and VI.The method known as the temperature-vegetation index method (TVX) can be described as a representative, which has been applied widely to estimate near surface air temperature (T a ) from several satellite sensors [35][36][37][38][39][40].According to this method, remotely sensed T s tends to approach air temperature with increasing of vegetation cover, and the radiometric temperature of a full vegetated canopy is in equilibrium with the temperature of the air within the canopy [39,41].This provides a meaningful insight in the interpretation of the T s − V I triangle method mentioned above.
The objectives of this paper are twofold: (1) By the combination of two remotely sensed T s − V I relationships (the T s − V I triangle method and the TVX method), the T s − V I space was interpreted from a new perspective, and a new parameterization scheme of EF based entirely on remote sensing observations was developed, specifically with the goal of reducing the reliance of the traditional T s − V I method on the dry edge; and (2) the new parameterization scheme was compared with the traditional T s − V I triangle method, and was validated using in situ EF measurements.

Study Area and Field Measurements
This study was conducted on the Southern Great Plains (SGP) region of the US, which extends in longitude from 95 • 18 W to 99 • 30 W and in latitude from 34 • 30 N to 38 • 30 N. It was chosen as the first field measurement site by the Atmospheric Radiation Measurement (ARM) Program for several reasons including homogeneous geography and easy accessibility, wide variability of climate cloud type and surface flux properties, and large variation in temperature and rainfall.The region has a relatively flat terrain with heterogeneous land cover, which is characterized by mixed farming, interrupted forest, tall and short grass (Figure 1).Several studies have been conducted to estimate the evapotranspiration of the SGP region over the past few years [17,27,42,43].In this study, Energy Balance Bowen Ratio (EBBR) stations (Figure 1), maintained by the ARM are used for validation.The location, altitude and land cover of these stations are presented in Table 2.It should be noted that since the EBBR measurements are conducted every 30 min, only the EF measured at the time that approximates most the overpass time of satellite was used in validation.More information about the EBBR dataset can be found in the research by Wang et al. [17] and at the homepage of the ARM (www.arm.gov.com).
Remote Sens. 2017, 9, 26 3 of 25 the temperature of the air within the canopy [39,41].This provides a meaningful insight in the interpretation of the   −  triangle method mentioned above.The objectives of this paper are twofold: (1) By the combination of two remotely sensed   −  relationships (the   −  triangle method and the TVX method), the   −  space was interpreted from a new perspective, and a new parameterization scheme of EF based entirely on remote sensing observations was developed, specifically with the goal of reducing the reliance of the traditional   −  method on the dry edge; and (2) the new parameterization scheme was compared with the traditional   −  triangle method, and was validated using in situ EF measurements.

Study Area and Field Measurements
This study was conducted on the Southern Great Plains (SGP) region of the US, which extends in longitude from 95°18′W to 99°30′W and in latitude from 34°30′N to 38°30′N.It was chosen as the first field measurement site by the Atmospheric Radiation Measurement (ARM) Program for several reasons including homogeneous geography and easy accessibility, wide variability of climate cloud type and surface flux properties, and large variation in temperature and rainfall.The region has a relatively flat terrain with heterogeneous land cover, which is characterized by mixed farming, interrupted forest, tall and short grass (Figure 1).Several studies have been conducted to estimate the evapotranspiration of the SGP region over the past few years [17,27,42,43].In this study, Energy Balance Bowen Ratio (EBBR) stations (Figure 1), maintained by the ARM are used for validation.The location, altitude and land cover of these stations are presented in Table 2.It should be noted that since the EBBR measurements are conducted every 30 minutes, only the EF measured at the time that approximates most the overpass time of satellite was used in validation.More information about the EBBR dataset can be found in the research by Wang et al. [17] and at the homepage of the ARM (www.arm.gov.com).The optical and thermal infrared band data from Moderate Resolution Imaging Spectroradiometer (MODIS) Terra were used to parameterize EF.Specifically, four MODIS land surface and atmosphere products (version 5) were used in this study: daily T s product at 1-km resolution (MOD11A1), 16-day vegetation indices at 1-km resolution (MOD13A2), daily atmosphere profile product at 5-km resolution (MOD07_L2) and daily cloud product at 1-km resolution (MOD06_L2).In order to get rid of the atmosphere-contaminated data points and replace them through temporal interpolation, the 16-day NDVI extracted from MOD13A2 was time series smoothed by the LACC (locally adjusted cubic spline capping) method [44].Here it was linear interpolated to the daily value.As a result, mathematically smooth capping curves of NDVI were produced to fit the rapid seasonal change.Then following the squared relationship between fractional vegetation cover ( f c ) and NDVI proposed by Gillies et al. [45], f c was calculated.Two parameters defining the NDVI values of bare soil and full-cover vegetation are required in this transformation.Specifically, NDV I max , defining the NDVI value of full-cover vegetation, can be selected as the highest NDVI value within the scene [46,47], and NDV I min , defining the NDVI value of bare soil, can be inferred from the historical lowest NDVI value within the scene [48].In our research, the highest NDVI value (0.94) within the scene was used to define NDV I max , while the popular published NDV I min value of 0.05 was used to define NDV I min [48].MOD06_L2 and MOD07_L2 products were mainly used in the estimation of T a .Specifically, MOD06_L2 product was used to retrieve land surface temperature, while MOD07_L2 product was used to retrieve atmosphere temperature profiles and surface air pressure.Details about these two products were given by Zhu et al. [49].

Selection of Clear Sky Day Images
In order to monitor surface soil moisture status in the SGP region, daytime images for 21 days in the year 2004 with cloud cover less than 20% were selected by Sun et al. [50].The same 21 days were selected in our study to parameterize EF: DOY (Day of Year) 10, 22, 51, 67, 81, 92, 118, 152, 175, 202, 213, 239, 310, 332, 337, 339, 347, 348, 349, 355 and 360.A two-step procedure was conducted to further exclude the cloud-contamination pixels.Firstly, as Batra et al. did in their research, a cloud masking threshold (T s < 273 K and NDV I < 0) was applied to remove cloudy pixels [27].Then cloud fraction (CF) data was extracted from MOD06_L2 product, and only pixels with CF = 0 were determined as clear pixels for a given image.

Parameterization Scheme of EF Using the Traditional T s − V I Triangle Method
Based on an extension of the Priestley-Taylor equation [51] and the definition of EF, Jiang and Islam [21,24] proposed that the T s − V I triangle/trapezoid space could be used to parameterize EF as follows: where ∅ is a dimensionless parameter that accounts for aerodynamic and canopy resistances, ∆ (k•PaK −1 ) is the slope of saturated vapor pressure (SVP) at air temperature (K) and γ (k•PaK −1 ) is the psychrometric constant.∆ and ∅ are the only two unknown parameters that are needed for the estimation of EF.As shown by a number of studies [13,17,21], the sensitivity of ∆/(∆ + γ) on the variation of temperature is very small, so parameter ∆ is usually calculated directly from remotely sensed T s .The parameter ∅ in Equation ( 1) seems the same as α in the original version of the Priestley-Taylor equation, but there exists a significant difference in their physical meaning.α in the Priestley-Taylor equation is generally defined as the ratio of actual evaporation to the equilibrium evaporation, while ∅ in Equation ( 1) is interpreted as the actual surface resistance to evapotranspiration.Parameterization of ∅ is realized though the T s − V I triangle method (Figure 2).As mentioned above, the dry and wet edges of the trapezoid space represent two limiting cases of EF for each fractional vegetation cover ( f c ): the dry edge represents minimum EF in each f c class, and the wet edge represents maximum EF in each f c class.However, it should be noted that there are two kinds of dry edge in the traditional T s − V I triangle/trapezoid method [30].The first one is called the theoretical dry edge, which is determined using the surface energy balance principle.The latent heat flux along the theoretical dry edge is assumed to be zero.As a result, the Priestley-Taylor parameter ∅ along the whole theoretical dry edge is equal to zero, just as it is presented in Figure 2. Another kind of dry edge is called the observed dry edge, which is determined using statistical regression method (Table 1).The original parameterization scheme of EF proposed by Jiang and Islam [21,24] is developed under the constraint of the observed dry edge, in which the parameter ∅ along the dry edge is assumed to be varying linearly from 0 to 1.26 with the increase of f c .The actual EF for each pixel in this T s − V I space is derived using a two-step linear interpolation scheme, which is described as follows: (1) When significant advection and convection are absent, latent heat flux cannot excess available total energy, and hence ∅ has a limited range between 0 and (∆ + γ)/∆.Thus, the global minima and maxima of ∅ can be easily determined by ∅ min = 0 for Point A ( f c = 0, T s = T smax ) in the dry edge, and ∅ max = (∆ + γ)/∆ for Point D ( f c = 1, T s = T w ) in the wet edge.
(2) Given these two bounds ( ) and the range of f c , the minimum ∅ min,i for each vegetation class f c,i can be interpolated linearly as follows The maximum ∅ max,i for each vegetation class f c,i is assumed to be a constant, and is always equal to ∅ max .
(3) After the lower and upper bounds of ∅ values for each f c class have been determined, the ∅ value for each pixel within this f c class is interpolated as follows: In which T min,i = T w (5) where T smax , T cmax and T w are the y-axis of Point A, C, and D, respectively.T smax,i and T smin,i are the maximum and minimum T s for a given fractional vegetation cover f c,i , and T s,i is the land surface temperature for each pixel within this f c class.

Basic Framework
The proposed new parameterization scheme of EF was conducted by the combination of the   −  triangle method and the TVX method.The temperature-vegetation index (TVX) method, proposed by Nemani and Running [52] and Goward et al. [53], was mainly used to estimate near surface air temperature (  ) in previous studies.It is assumed that because of the combined impact of vegetation cover on the average surface thermal characteristics and on the evaporative control of energy portioning, remotely sensed land surface temperature tends to approach air temperature with increasing vegetation cover, and the radiometric temperature of a full vegetated canopy is in equilibrium with the temperature of the air within the canopy [39,41].Other assumptions involved in the TVX method are that both uniform atmospheric forcing and soil moisture conditions must take place [38].It should be noted that the uniformity of atmospheric forcing is also the assumption of the   −  triangle method, and is considered to be fulfilled on clear sky days when both the retrieval of   and VI is possible, while the assumption of uniform soil moisture is just presumably fulfilled in the isopleths of soil surface moisture within the trapezoid space.Thus, it is possible to

Basic Framework
The proposed new parameterization scheme of EF was conducted by the combination of the T s − V I triangle method and the TVX method.The temperature-vegetation index (TVX) method, proposed by Nemani and Running [52] and Goward et al. [53], was mainly used to estimate near surface air temperature (T a ) in previous studies.It is assumed that because of the combined impact of vegetation cover on the average surface thermal characteristics and on the evaporative control of energy portioning, remotely sensed land surface temperature tends to approach air temperature with increasing vegetation cover, and the radiometric temperature of a full vegetated canopy is in equilibrium with the temperature of the air within the canopy [39,41].Other assumptions involved in the TVX method are that both uniform atmospheric forcing and soil moisture conditions must take place [38].It should be noted that the uniformity of atmospheric forcing is also the assumption of the T s − V I triangle method, and is considered to be fulfilled on clear sky days when both the retrieval of T s and VI is possible, while the assumption of uniform soil moisture is just presumably fulfilled in the isopleths of soil surface moisture within the trapezoid space.Thus, it is possible to combine these two T s − V I relationships together.In fact, as the study conducted by Sandholt et al. [54] revealed, the isopleths of soil surface moisture within the T s − V I trapezoid space can be regarded as several superimposed TVX lines (Figure 2).There have been other significant studies performed to derive further the information contained in the T s − V I trapezoid space.Based on the studies conducted by Carlson [28] and Moran et al. [55], Long and Singh [30] concluded that pixels along the same isopleth of soil surface moisture has the same soil surface temperature (T soil,i ) and also the same surface temperature of the fully vegetated surface (T c,i ), and T s of a mixed surface is a weighted sum of vegetation and soil temperatures.On the basis of this conclusion and by the combination of the T s − V I triangle method and the TVX method, the proposed new parameterization scheme of EF was performed as follows.Figure 3 is a flowchart of this parameterization scheme.Five steps are demonstrated using different colors.
Remote Sens. 2017, 9, 26 7 of 25 combine these two   −  relationships together.In fact, as the study conducted by Sandholt et al. [54] revealed, the isopleths of soil surface moisture within the   −  trapezoid space can be regarded as several superimposed TVX lines (Figure 2).There have been other significant studies performed to derive further the information contained in the   −  trapezoid space.Based on the studies conducted by Carlson [28] and Moran et al. [55], Long and Singh [30] concluded that pixels along the same isopleth of soil surface moisture has the same soil surface temperature ( , ) and also the same surface temperature of the fully vegetated surface ( , ), and   of a mixed surface is a weighted sum of vegetation and soil temperatures.On the basis of this conclusion and by the combination of the   −  triangle method and the TVX method, the proposed new parameterization scheme of EF was performed as follows.Figure 3 is a flowchart of this parameterization scheme.Five steps are demonstrated using different colors.Step 1 is the calculation of the surface temperature of bare soil ( , ), which is marked with green color.According to the assumption of the TVX method, for each pixel ( , ,  , ) in Figure 2, the remotely sensed   tends to approach its air temperature ( , ) with increasing of vegetation cover under the same soil moisture conditions, and the radiometric temperature of a full vegetated canopy ( , ) is in equilibrium with the temperature of the air within the canopy.Thus,  , =  , .Variation in   for the same isopleth of soil surface moisture results essentially from the variation in   , and the remotely sensed  , is a weighted sum of vegetation and soil temperatures [30], which is defined as: For a mixed pixel,  , and  , is the surface temperature of full vegetated canopy and bare soil, respectively.Therefore, the surface temperature of the bare soil ( , ) corresponding to the isopiestic line which the given pixel ( , ,  , ) belongs to can be deduced as the intercept of this function.
Step 2 is the calculation of a simplified land surface dryness index (Temperature-Vegetation Dryness Index, TVDI) proposed by Sandholt et al. [54], which is marked with yellow color.According to Sandholt et al. [54], the TVDI for each pixel in the   −  triangle space can be retrieved as: Step 1 is the calculation of the surface temperature of bare soil (T soil,i ), which is marked with green color.According to the assumption of the TVX method, for each pixel (T s,i , f c,i ) in Figure 2, the remotely sensed T s tends to approach its air temperature (T a,i ) with increasing of vegetation cover under the same soil moisture conditions, and the radiometric temperature of a full vegetated canopy (T c,i ) is in equilibrium with the temperature of the air within the canopy.Thus, T c,i = T a,i .Variation in T s for the same isopleth of soil surface moisture results essentially from the variation in f c , and the remotely sensed T s,i is a weighted sum of vegetation and soil temperatures [30], which is defined as: For a mixed pixel, T c,i and T soil,i is the surface temperature of full vegetated canopy and bare soil, respectively.Therefore, the surface temperature of the bare soil (T soil,i ) corresponding to the isopiestic line which the given pixel (T s,i , f c,i ) belongs to can be deduced as the intercept of this function.
Step 2 is the calculation of a simplified land surface dryness index (Temperature-Vegetation Dryness Index, TVDI) proposed by Sandholt et al. [54], which is marked with yellow color.According to Sandholt et al. [54], the TVDI for each pixel in the T s − V I triangle space can be retrieved as: where T smax,i and T smin,i is the maximum and minimum surface temperature of each f c class, respectively (Figure 2).For bare soil ( f c = 0), TVDI can be defined: It is clear that there is negative correlation relationship between TVDI and soil surface moisture status.TVDI ranges from 0 to 1 with soil moisture availability decreasing from the wet edge to the dry edge.
Step 3 is the calculation of the Priestley-Taylor parameter of bare soil (∅ s,i ), which is marked with blue color.For unsaturated soil, a parameterization scheme of EF was proposed by Komatsu [56] from the soil moisture status as follows: where α = 1.26 is the Priestley-Taylor parameter, θ is the surface soil moisture, and θ c is a parameter that depends on soil type and wind speed.Similar to TVDI, θ/θ c is also a dimensionless parameter ranging from 0 to 1, which represents relative soil moisture status, but there is negative correlation relationship existing between TVDI and θ/θ c [57].Combing Equations ( 1) and ( 9) and replacing θ/θ c with (1 − TVDI soil ), the parameter ∅ in Equqtion (1) can be written as The subscript "s" indicates that the parameter ∅ is only valid for bare soil.
Step 4 is the calculation of the Priestley-Taylor parameter of full vegetated canopy (∅ c,i ), which is marked with orange color.At the same isopiestic line corresponding to ∅ s,i , the parameter ∅ of the full vegetated canopy ( f c = 1) is calculated as ∅ c,i = (∆ + γ)/∆, just as Jiang and Islam did in Section 3.1 [16,19].However, it should be noted that ∆ in Jiang and Islam's parameterization scheme is calculated by using the constant T w as input, while ∆ in our new parameterization scheme is calculated by using T a,i as input, and varies with isopleths of soil surface moisture.
Step 5 is the calculation of the Priestley-Taylor parameter of a mixed pixel (∅ i ), which is marked with red color.Along each isopiestic line, because of the invariance of soil moisture, the parameter ∅ is only determined by the variation of vegetation cover, and increases from ∅ s,i to ∅ c,i with the increasing of f c .After the lower and upper bounds of ∅ values for each isopiestic line have been determined, the ∅ value for each pixel within this isopiestic line is interpolated as follows: Similar to the parameterization scheme developed by Jiang and Islam [21,24], the proposed new parameterization scheme was also performed based on an extension of Priestley-Taylor equation.The parameter ∅ for each pixel within the trapezoid space was interpolated linearly from the established upper and lower bounds in both of these two schemes.What is different is that the interpolation procedure in Jiang and Islam's parameterization scheme was performed for each f c interval, while in the proposed new parameterization scheme it was performed for each isopiestic line of soil surface moisture.As for the determination of the limiting bounds, there is not much difference in the determination of the upper bound (∅ max,i in Jiang and Islam's scheme and ∅ c,i in our scheme).Both of them are calculated as (∆ + γ)/∆.However, ∆ in Jiang and Islam's parameterization scheme is calculated by using the constant T w as input, while ∆ in our new parameterization scheme is calculated by using T a,i as input.Therefore, ∅ max,i is a constant in the traditional approach, while ∅ c,i varies with isopleths of soil surface moisture in the NPS.Besides, because ∆ increases with the increase of T a,i , the pixels with more stressed soil moisture conditions have a lower value of ∅ c,i .Therefore, compared with the constant ∅ max in Jiang and Islam's scheme, the variable ∅ c,i in the new parameterization scheme is more reasonable, because it accounts for the impact of soil moisture conditions partly.In Jiang and Islam's parameterization scheme, the determination of the lower bound (∅ min,i ) for each vegetation class f c was interpolated linearly as ∅ min,i = ∅ max f c,i .According to the research conducted by Stisen et al. [23], the weakness of this parameterization of ∅ min,i , where evaporative fraction is not zero along the observed dry edge, is that it does not allow for the presence of water stressed full cover vegetation.In our proposed new scheme, the lower bound (∅ s,i ) for each isopiestic line was determined under bare soil conditions based on the TVDI proposed by Sandholt et al. [54] and the parameterization scheme proposed by Komatsu [56].This is an alternative way to bypass the water stress involved in the mixed surface (vegetation/bare soil), and thus can overcome the weakness involved in the parameterization of ∅ min,i .Another great advantage of the new parameterization scheme is that both the maximum surface temperature of the bare soil and fully vegetated canopy is required in the parameterization scheme proposed by Jiang and Islam [21,24], while only the maximum surface temperature of the bare soil (T smax ) is indispensable in the new parameterization scheme, and thus the reliance of the T s − V I triangle method on the dry edge has been reduced significantly.

Estimation of Near Surface Air Temperature
In the proposed new parameterization scheme, the isopiestic lines of soil surface moisture were assumed to be several superimposed TVX lines [54].Besides remotely sensed T s and VI, near surface air temperature (T a ) is another variable needed to determine the TVX lines.A simple and operational algorithm was developed by Zhu et al. [49] to retrieve the instantaneous T a for the Southern Great Plains (SGP) of the United States of America.They find that the systematic errors caused by near surface air temperature (T s a ) retrieved from MOD07_L2 product and T s retrieved from MOD06_L2 product are in completely different directions.This means that these errors can balance each other out by the combination of T s a and T s , especially when the absolute values of the errors caused by them are similar to each other, which is just the case in our study.Therefore, the instantaneous T a under clear sky conditions was estimated as the average of near surface air temperature (T s a ) retrieved from MOD07_L2 product and land surface temperature (T s ) retrieved from MOD06_L2 product.The validation results are shown in Figure 4. Specific statistical measures used in the validation are shown in Table 3.The coefficient of determination (R 2 ), MAE, RMSE, RRMSE and B are 0.94, 2.28 K, 3.02 K, 0.01 and 1.72 K, respectively.T a of the same 21 days in the SGP region was also estimated by Sun et al. [50] using an empirical method proposed by Zaksek and Schroedter-Homscheidt [58] with a R 2 of 0.91, and RMSE of 3.1 K. Therefore, the accuracy of the algorithm developed by Zhu et al. [49] applied in the SGP is acceptable.
linearly as ∅ , = ∅   , .According to the research conducted by Stisen et al. [23], the weakness of this parameterization of ∅ , , where evaporative fraction is not zero along the observed dry edge, is that it does not allow for the presence of water stressed full cover vegetation.In our proposed new scheme, the lower bound (∅ , ) for each isopiestic line was determined under bare soil conditions based on the TVDI proposed by Sandholt et al. [54] and the parameterization scheme proposed by Komatsu [56].This is an alternative way to bypass the water stress involved in the mixed surface (vegetation/bare soil), and thus can overcome the weakness involved in the parameterization of ∅ , .Another great advantage of the new parameterization scheme is that both the maximum surface temperature of the bare soil and fully vegetated canopy is required in the parameterization scheme proposed by Jiang and Islam [21,24], while only the maximum surface temperature of the bare soil (  ) is indispensable in the new parameterization scheme, and thus the reliance of the   −  triangle method on the dry edge has been reduced significantly.

Estimation of Near Surface Air Temperature
In the proposed new parameterization scheme, the isopiestic lines of soil surface moisture were assumed to be several superimposed TVX lines [54].Besides remotely sensed   and VI, near surface air temperature (  ) is another variable needed to determine the TVX lines.A simple and operational algorithm was developed by Zhu et al. [49] to retrieve the instantaneous   for the Southern Great Plains (SGP) of the United States of America.They find that the systematic errors caused by near surface air temperature (   ) retrieved from MOD07_L2 product and   retrieved from MOD06_L2 product are in completely different directions.This means that these errors can balance each other out by the combination of    and   , especially when the absolute values of the errors caused by them are similar to each other, which is just the case in our study.Therefore, the instantaneous   under clear sky conditions was estimated as the average of near surface air temperature (   ) retrieved from MOD07_L2 product and land surface temperature (T s ) retrieved from MOD06_L2 product.The validation results are shown in Figure 4. Specific statistical measures used in the validation are shown in Table 3.The coefficient of determination (R 2 ), MAE, RMSE, RRMSE and B are 0.94, 2.28 K, 3.02 K, 0.01 and 1.72 K, respectively.  of the same 21 days in the SGP region was also estimated by Sun et al. [50] using an empirical method proposed by Zaksek and Schroedter-Homscheidt [58] with a R 2 of 0.91, and RMSE of 3.1 K. Therefore, the accuracy of the algorithm developed by Zhu et al. [49] applied in the SGP is acceptable.

Statistical Measure Formula
Mean absolute error MAE = 1 Relative root mean square error

Determination of the Dry and Wet Edges
Scheme proposed by Jiang and Islam [21,24] and the proposed new parameterization scheme were performed within the triangle/trapezoid framework (Figure 2).Different algorithms for determining the dry and wet edges of the triangle/trapezoid were developed in previous studies [13,21,23,24,30,32].In this study, the observed dry edge was determined using a simple algorithm, which was described in detail by Tomas et al. [59].In simple terms, the T s pixels are firstly split into a number of bins using a NDVI step size of 0.01.Then the maximum T s value is found for each bin.After removing all the maximum T s values located to the left in the triangle graph (with corresponding lower NDVI value), a linear fit is performed through the remaining T s maximum values and the corresponding NDVI bin values, and the fit is adopted as the dry edge.
In general, there are mainly three methods in the determination of the wet edge.(1) In traditional triangle approaches, the minimum T s value in an image or the mean of the minimum T s values for each NDVI bin is adopted as the wet edge [54,59].(2) Some typical pixels of an image are taken as the wet edge.These pixels are often characterized by relatively low T s and high EF, such as inland water bodies and dense vegetation stands [18,24,60].(3) Numerous studies show that observed near surface air temperature can be used to represent the lowest T s over study area.For example, one principle of the TVX method is that the radiometric temperature of a fully vegetated canopy is in equilibrium with the air temperature within the canopy [39,41,45,61].Even some T s − V I triangle/trapezoid models that were performed based on theoretical boundaries also assumed that the near surface T a could be taken as the horizontal wet edge [30,32,50].Thus, in this study we adopted the minimum T a estimated over the whole SGP region in Section 3.2.2 as the wet edge, and set T smin = T cmin = T w = T amin .

Accuracy of EF Estimates
For simplification purpose, the parameterization scheme of EF proposed by Jiang and Islam [21,24] is called the traditional parameterization scheme (TPS) and that proposed in this paper is called the new parameterization scheme (NPS).As Batra et al. [27] and Venturini et al. [42] interpreted, there are no generally accepted methodologies to compare and validate the spatial distribution of EF over large areas because of the scale mismatch between the estimated EF map and point measurements from the ground.Nonetheless unfiltered point observations seem to be an appropriate mean to validate remote sensing applications [21,22,27,42,62].The estimated EF for the two parameterization schemes is illustrated in Figure 5a,b, respectively, and statistics are listed in Table 4. Overall, similar accuracy is achieved in these two parameterization schemes.Both of them could capture the seasonal variation of EF with R 2 higher than 0.55, although a higher standard deviation (STD) is found in the ground observations.Judging from the magnitude, the MAE, RMSE, RRMSE and B obtained in TPS are 0.11, 0.14, 0.33 and −0.03, respectively, which are just the same as those obtained in NPS.This means that EF is underestimated in both of these two schemes.In order to further investigate the difference in EF estimates between these two parametrization schemes, the scatterplots between EF retrieved from the TPS and that retrieved from the NPS are shown in Figure 6.It is clear that good agreements are observed with correlation R 2 as high as 0.98 and bias B as low as 0.00.The low values of the MAE and RMSE also indicate that there is little difference in the EF retrieved from these two parameterization schemes.The validation results of these two parameterization schemes for 11 sites are also presented in Table 4. Compared with the overall performance, significant differences are observed in the accuracy obtained in different sites.For most of the sites (E2, E4, E9, E12, E13, E18, E20, E22, E27), the accuracy obtained in both schemes is similar to or even better than that obtained for the whole study area.However, the errors obtained for E7 and E8 are much larger than those obtained for the remaining.There exists serious underestimation of EF for these two sites, and the bias B obtained in both schemes are −0.14 (E7) and −0.10 (E8), respectively.
in Table 4. Overall, similar accuracy is achieved in these two parameterization schemes.Both of them could capture the seasonal variation of EF with R 2 higher than 0.55, although a higher standard deviation (STD) is found in the ground observations.Judging from the magnitude, the MAE, RMSE, RRMSE and B obtained in TPS are 0.11, 0.14, 0.33 and −0.03, respectively, which are just the same as those obtained in NPS.This means that EF is underestimated in both of these two schemes.In order to further investigate the difference in EF estimates between these two parametrization schemes, the scatterplots between EF retrieved from the TPS and that retrieved from the NPS are shown in Figure 6.It is clear that good agreements are observed with correlation R 2 as high as 0.98 and bias B as low as 0.00.The low values of the MAE and RMSE also indicate that there is little difference in the EF retrieved from these two parameterization schemes.The validation results of these two parameterization schemes for 11 sites are also presented in Table 4. Compared with the overall performance, significant differences are observed in the accuracy obtained in different sites.For most of the sites (E2, E4, E9, E12, E13, E18, E20, E22, E27), the accuracy obtained in both schemes is similar to or even better than that obtained for the whole study area.However, the errors obtained for E7 and E8 are much larger than those obtained for the remaining.There exists serious underestimation of EF for these two sites, and the bias B obtained in both schemes are −0.14 (E7) and −0.10 (E8), respectively.

Temporal Variation of EF
The temporal variations of the observed EF versus EF retrieved from these two parameterization schemes for each site are shown in Figure 7.In general, there are no significant differences between EF retrieved from TPS and that retrieved from NPS.Both can reflect the temporal variations of EF effectively for most of the sites (E2, E4, E7, E9, E12, E18, E20, E27).However, when compared with the value range of observed EF, neither of these two schemes can represent effectively the amplitude of the variation.The value range of the observed EF is much broader than those retrieved from either the TPS or the NPS, which is especially obvious for sites E2, E9, E12, E13, E18, E20 and E22.The reasons involved can be found by investigating the minimum and maximum values of EF for each site.As Figure 7 shows, for most sites the observed minimum EF is much smaller than those retrieved from these two schemes, while the differences in the maximum EF are not very significant.Thus, the big differences in the value range are mainly caused by the overestimation of the minimum EF.It should be noted that for site E7 and E8, although accurate estimates of both the maximum and minimum EF are obtained in these two schemes, the comparison of mean values shows that EF of these two sites are underestimated most seriously.This is mainly caused by the large errors in TVDI estimates, because the research conducted by Sun et al. shows that the correlation coefficient (r) between TVDI and soil moisture for site E7 and E8 is 0.10 and 0.33, respectively, which means there is almost no correlation [50].
In general, the time series curves produced by the parameterization schemes are more similar to a smooth sinusoidal curve, while those produced by the observations have rougher shape, which takes on a more significant seasonal variation.However, a definitive conclusion from this comparison should be avoided because of the lack of long-term comparison.In addition, there are issues of mixed pixel and registration error between the remotely sensed pixel and ground observations [27].Random measurement error arises at some points in these stations on account of varied usage of sensors, measurement techniques and system calibration methods [21].The measurement errors should also kept in mind while doing such comparisons.

Temporal Variation of EF
The temporal variations of the observed EF versus EF retrieved from these two parameterization schemes for each site are shown in Figure 7.In general, there are no significant differences between EF retrieved from TPS and that retrieved from NPS.Both can reflect the temporal variations of EF effectively for most of the sites (E2, E4, E7, E9, E12, E18, E20, E27).However, when compared with the value range of observed EF, neither of these two schemes can represent effectively the amplitude of the variation.The value range of the observed EF is much broader than those retrieved from either the TPS or the NPS, which is especially obvious for sites E2, E9, E12, E13, E18, E20 and E22.The reasons involved can be found by investigating the minimum and maximum values of EF for each site.As Figure 7 shows, for most sites the observed minimum EF is much smaller than those retrieved from these two schemes, while the differences in the maximum EF are not very significant.Thus, the big differences in the value range are mainly caused by the overestimation of the minimum EF.It should be noted that for site E7 and E8, although accurate estimates of both the maximum and minimum EF are obtained in these two schemes, the comparison of mean values shows that EF of these two sites are underestimated most seriously.This is mainly caused by the large errors in TVDI estimates, because the research conducted by Sun et al. shows that the correlation coefficient (r) between TVDI and soil moisture for site E7 and E8 is 0.10 and 0.33, respectively, which means there is almost no correlation [50].
In general, the time series curves produced by the parameterization schemes are more similar to a smooth sinusoidal curve, while those produced by the observations have a rougher shape, which takes on a more significant seasonal variation.However, a definitive conclusion from this comparison should be avoided because of the lack of long-term comparison.In addition, there are issues of mixed pixel and registration error between the remotely sensed pixel and ground observations [27].Random measurement error arises at some points in these stations on account of varied usage of sensors, measurement techniques and system calibration methods [21].The measurement errors should also kept in mind while doing such comparisons.

Spatial Comparison of EF Retrieved from the TPS and NPS
Although the accuracy of EF estimated by the TPS and NPS as well as their temporal variations were compared in Section 4.1 and 4.2, respectively, further studies are still needed to investigate the spatial differences between these two parameterization schemes.Because it is difficult to evaluate the reliability of model output using point observations at a handful of sites, in this section EF retrieved from the two parameterization schemes was compared at pixel scale for the whole SGP region.Because of its low cloud cover, here only the day of DOY 239 was taken as an example.The spatial distribution of EF retrieved from these two parameterization schemes as well as their differences (EF NPS − EF TPS ) are shown in Figure 8.In general, similar spatial pattern of EF was achieved.Both of these two schemes indicate that there are obvious distinctions in EF between the west and east part of the SGP region.Specifically, EF retrieved from the TPS ranges from 0.16 to 0.96 with a mean of 0.68 and STD of 0.15, and EF retrieved from the NPS ranges from 0.19 to 0.99 with a mean of 0.67 and STD of 0.15.The scatterplots of EF retrieved from the TPS and NPS for all the pixels are presented in Figure 9.There is a good agreement of EF retrieved from these two parameterization schemes with R 2 as high as 0.96.Similar to the results obtained at site scale (Figure

Spatial Comparison of EF Retrieved from the TPS and NPS
Although the accuracy of EF estimated by the TPS and NPS as well as their temporal variations were compared in Sections 4.1 and 4.2, respectively, further studies are still needed to investigate the spatial differences between these two parameterization schemes.Because it is difficult to evaluate the reliability of model output using point observations at a handful of sites, in this section EF retrieved from the two parameterization schemes was compared at pixel scale for the whole SGP region.Because of its low cloud cover, here only the day of DOY 239 was taken as an example.The spatial distribution of EF retrieved from these two parameterization schemes as well as their differences (EF NPS − EF TPS ) are shown in Figure 8.In general, similar spatial pattern of EF was achieved.Both of these two schemes indicate that there are obvious distinctions in EF between the west and east part of the SGP region.Specifically, EF retrieved from the TPS ranges from 0.16 to 0.96 with a mean of 0.68 and STD of 0.15, and EF retrieved from the NPS ranges from 0.19 to 0.99 with a mean of 0.67 and STD of 0.15.The scatterplots of EF retrieved from the TPS and NPS for all the pixels are presented in Figure 9.There is a good agreement of EF retrieved from these two parameterization schemes with R 2 as high as 0.96.Similar to the results obtained at site scale (Figure 6), small values of the MAE, RMSE and B are also achieved at regional scale, which are 0.03, 0.04 and −0.02, respectively.
Remote Sens. 2017, 9, 26 14 of 25 6), small values of the MAE, RMSE and B are also achieved at regional scale, which are 0.03, 0.04 and −0.02, respectively.In order to further analyze the spatial distribution of the differences in EF retrieved from these two parameterization schemes, the scatterplots of (EF NPS − EF TPS ) against NDVI and T s are presented in Figure 10a,b, respectively.Positive correlations are observed in both of these two scatterplots.This means that there is complementary relationship exiting between the EF retrieved from these two schemes.For pixels with small values of NDVI and T s (in the lower-left corner of the T s − V I trapezoid space), the NPS tends to underestimate the EF retrieved from the TPS, while for pixels with large values of NDVI and T s (in the upper-right corner of the T s − V I triangle space), the NPS tends to overestimate the EF retrieved from the TPS.These findings make a great sense because it means that the errors produced by these two parameterization schemes have a distinct pattern.As a matter of fact, the application of the T s − V I triangle approaches has been recently questioned when the climatic controls for evapotranspiration are different.For example, the research conducted by Garcia et al. [63] showed that the application of the TPS should be restricted to water-limited sites excluding energy-limited regions.Considering the complementary relationship presented in Figure 10, further studies can be performed to evaluate the applicability of the TPS and NPS under different climatic controls.
In order to make a more comprehensive evaluation of the difference in EF retrieved by these two schemes, Figure 11 presents the spatial quantile statistics of EF retrieved from these two parameterization schemes as well as the spatial quantile statistics of their differences (EF NPS − EF TPS ) on each of the 21 days.In general, good temporal agreements are observed in the spatial quantile statistics of EF retrieved by these two schemes.R 2 calculated for the 25th percentile, median and 75th percentile is 0.95, 0.97 and 0.97, respectively.Judging from the seasonal variation, both of these two schemes show that the spatial median values of EF over the SGP region increases first and then deceases with increasing of time.The peak values of EF estimated by the TPS and NPS are 0.71 and 0.68, respectively (on DOY 239).As for the difference of EF retrieved by these two schemes, the median values shows that the NPS can either underestimate or overestimate the TPS.Specifically, the TPS is overestimated by the NPS during 11 days of the 21 days; during the rest 10 days, the TPS is underestimated.

Comparison with Previous Studies
Although the accuracy of the parameterization scheme developed in this paper has been evaluated in the above section, the comparison of its accuracy with previous studies is still necessary to understand its advantages.

Comparison with Previous Studies
Although the accuracy of the parameterization scheme developed in this paper has been evaluated in the above section, the comparison of its accuracy with previous studies is still necessary to understand its advantages.
Table 5 presents the summary of relevant EF estimation studies conducted based on the TPS.As mentioned in the introduction, the method adopted in previous studies for the retrieval of dry edge can be classified as theoretical and empirical approaches.In order to be consistent with the regression method used in this study, all the studies presented in Table 5 were performed within the

Comparison with Previous Studies
Although the accuracy of the parameterization scheme developed in this paper has been evaluated in the above section, the comparison of its accuracy with previous studies is still necessary to understand its advantages.
Table 5 presents the summary of relevant EF estimation studies conducted based on the TPS.As mentioned in the introduction, the method adopted in previous studies for the retrieval of dry edge can be classified as theoretical and empirical approaches.In order to be consistent with the regression method used in this study, all the studies presented in Table 5 were performed within the observed dry edge determined by empirical approach.It is clear that the accuracy of EF estimates in previous studies varies significantly with study location and satellite sensor used.Coincidently, there are two studies in which the traditional parameterization scheme was applied to the same study location as ours.EF over the SGP region was also retrieved from MODIS product in the year 2004 by Wang et al. [17], although the sixteen days selected by them for validation were not the same as ours.In their study, the MAE and B obtained were 0.14 and −0.03, which were then reduced to 0.12 and −0.02 by the introduction of the thermal inertia method.Similar results were obtained by Jiang and Islam [25], but EF over the SGP region in their study was retrieved from AVHRR product.The RMSE, B and R 2 reported were 0.12, −0.08 and 0.30, respectively.Besides the SGP region, the TPS has been applied to a wide range of other locations [22,23,26,59,64,65].According to Long et al. [33] and Tomas et al. [59], the performance of triangle method depends heavily on the resolution of satellite images.Therefore, just judging from the accuracy produced by MODIS product in Table 5, it is concluded that the accuracy achieved by the NPS in our work has reached a level comparable with those obtained by the TPS in previous studies.It should be noted that although both the TPS and NPS are conducted within the dry edge determined in Section 3.2.3, a careful inspection of the proposed new parameterization scheme shows that unlike the parameterization scheme proposed by Jiang and Islam [21,24], T cmax of the dry edge is not used in the NPS.Specifically, both T smax and T cmax are needed in the TPS for the retrieval of dry edge [30,32], while only T smax is indispensable in the NPS for the definition of the maximum water stress.Therefore, the reliance of the T s − V I triangle method on the dry edge has been reduced significantly by the NPS.Combining with the accuracy it produced (Table 5), it is reasonable to conclude that the simplicity of the proposed new parameterization scheme does not compromise its accuracy in monitoring EF.The independence of the NPS from T cmax also means that it can be performed independent from the determination of the dry edge.This makes sense in the following two aspects: (1) As mentioned in the introduction, because of the limitations of the observed dry edge, theoretical dry edge has been widely used in recent studies [30,32,50], in which both the values of T smax and T cmax are determined based on the surface energy balance principle.The introduction of the NPS has not only bypassed the task involved in the determination of T cmax , but also reduced the uncertainty of the T s − V I trapezoid method, because the mechanism involved in the surface energy balance principle over the bare soil is much simpler.(2) As presented in Figure 2, the dry edge of the T s − V I trapezoid can be regarded as a special isopleth of the soil surface moisture.Thus, it can also be determined as other isopiestic lines based on the TVX method, in which T smax can be retrieved.This means that by the combination of the T s − V I triangle method and the TVX method, the proposed new parameterization scheme can determine the boundaries that are needed by itself, and get rid of the constraint imposed by the T s − V I trapezoid framework.Details on this are presented in the following section.

Extension of the Proposed New Parameterization Scheme
The comparison with previous studies shows that although the new parameterization scheme has no superiority in accuracy over the traditional parameterization scheme, the alternative option it provides makes it possible to be performed independent from the determination of the dry edge.Regarding the dry edge of the T s − V I feature space as a special isopleth of the soil surface moisture, it can be retrieved just the same as the retrieval of other isopleths.Specifically, the T smax was determined by selecting the pixel with the maximum T s in an image.Given T s , f c and T a , T smax of the selected pixel was calculated using Equation ( 6).Although T cmax is not indispensable in the proposed new parameterization scheme, it can be assumed as T a of the selected pixel according to the TVX method.The comparison of T smax and T cmax determined above with those determined in Section 3.2.3 is shown in Figure 12a,b, respectively.Although good temporal agreements are observed in both of these two parameters with r equal to 0.98 and 0.96, respectively, there are significant differences in the errors achieved.For T smax , there is little difference in the values estimated by these two methods.The MAE, RMSE and B are 1.08 K, 1.57K and −0.06 K, respectively.This means that it is feasible to determine the constraint (T smax ) that are needed in the NPS by the combination of the T s − V I triangle method and the TVX method.However, significant differences are observed for the estimation of T cmax .The MAE, RMSE and B are 6.92 K, 7.67 K and −5.83 K, respectively.The near surface air temperature of the selected pixel is generally lower than T cmax determined by the regression approach.
The accuracy of EF retrieved from the TPS and NPS within the dry edge determined above are presented in Table 6.Compared with the results in Table 4, there is little difference in the accuracy obtained by the NPS.The R 2 , MAE, RMSE and B obtained are 0.56, 0.11, 0.14 and −0.02, respectively, which are almost the same as those presented in Table 4.This is reasonable considering the little difference in T smax presented in Figure 12a.The R 2 , MAE, RMSE and B obtained by the TPS are 0.58, 0.12, 0.15 and −0.05, respectively.All of the four statistics show that the accuracy achieved by the new dry edge is a little lower than that achieved in Table 4. Since there is little difference in T smax , the decline of the accuracy is mainly caused by the difference in T cmax .The comparison shows that the need of another parameter T cmax brings in more uncertainty to final parameterization of EF, which partly demonstrates the robustness of the proposed new parameterization scheme, since the NPS can be performed independent of the determination of T cmax .

Sensitivity of the NPS to Input Parameters
The sensitivity analysis includes the following two aspects.A comparison of these two parameterization schemes of EF shows that four temperature related parameters are required in both of these two schemes.For the TPS, they are land surface temperature ( ), maximum surface temperature of bare soil ( ), maximum surface temperature of full cover vegetation ( ) and surface temperature of the wet edge ( ), respectively.For the NPS, besides , and , near surface air temperature ( ) is required.Thus, firstly, the sensitivity of EF estimates to these four parameters was tested.A set of values with a step of 1 K were adopted here to define the variation of temperature.The strategy used for sensitivity analysis is to compare EF estimates with only one input variable changed to those without any change in inputs.The results are shown in Figure 13.Except , positive feedback effects between EF retrievals and the rest three variables are observed in both of these two parameterization schemes.Any increase (decrease) in the parameters ( , , and ) will induce the increase (decrease) in the final estimates of EF.In contrast, the correlation relationship between EF retrievals and is negative.Results show that EF estimates are most sensitive to the variation of and .Logarithmic relationships between the bias B and the variation of are observed in both of these two schemes.It means that the sensitivity of EF retrieval decreases with the increase of .Taken the new parameterization scheme as an example, a change of +6 K in causes an increase of 0.10 in EF, while a change of −6 K in

Sensitivity of the NPS to Input Parameters
The sensitivity analysis includes the following two aspects.A comparison of these two parameterization schemes of EF shows that four temperature related parameters are required in both of these two schemes.For the TPS, they are land surface temperature (T s ), maximum surface temperature of bare soil (T smax ), maximum surface temperature of full cover vegetation (T cmax ) and surface temperature of the wet edge (T w ), respectively.For the NPS, besides T s , T smax and T w , near surface air temperature (T a ) is required.Thus, firstly, the sensitivity of EF estimates to these four parameters was tested.A set of values with a step of 1 K were adopted here to define the variation of temperature.The strategy used for sensitivity analysis is to compare EF estimates with only one input variable changed to those without any change in inputs.The results are shown in Figure 13.Except T s , positive feedback effects between EF retrievals and the rest three variables are observed in both of these two parameterization schemes.Any increase (decrease) in the parameters (T smax , T cmax , T w and T a ) will induce the increase (decrease) in the final estimates of EF.In contrast, the correlation Remote Sens. 2017, 9, 26 20 of 25 relationship between EF retrievals and T s is negative.Results show that EF estimates are most sensitive to the variation of T smax and T s .Logarithmic relationships between the bias B and the variation of T smax are observed in both of these two schemes.It means that the sensitivity of EF retrieval decreases with the increase of T smax .Taken the new parameterization scheme as an example, a change of +6 K in T smax causes an increase of 0.10 in EF, while a change of −6 K in T smax results in a decrease of EF as high as 0.26.The comparison between the NPS and TPS indicates that the TPS is more sensitive to the variation of T smax than the NPS.It should be noted that the validation results in Figure 5 show that EF was underestimated by both the TPS and the NPS.The underestimation in EF may be partly caused by the underestimation of T smax , since the observed dry edge is usually lower than the theoretical dry edge.In contrast to T smax , the sensitivity of EF retrievals increases with the increase of T s .For the NPS, a change of −6 K in T s causes an increase of 0.13 in EF, while a change of +6 K in T s results in a decrease of EF as high as 0.30.The comparison shows that the sensitivity of these two parameterization schemes to T s is different when the variation of T s is in the opposite directions.The TPS is more sensitive to the variation of T s than the NPS when the variation of T s is negative.However, when the variation of T s is positive, the NPS is more sensitive to the variation of T s than the TPS.Compared with T smax and T s , the sensitivity of EF retrievals to the rest three parameters is not very significant.In general, a ±6 K change of these three parameters results in only ± 0.07 change of EF, but the sensitivity of the TPS to T w increases significantly with the increase of T w .It has even exceeded the sensitivity to T smax when the variation of temperature is beyond +5.5 K.The results presented in Figure 4 indicate that T a over the SGP region is estimated with the accuracy of a RMSE as 3.02 K. Figure 13 shows that a change of ±3.02 K in T a results in only ± 0.03 change in EF.Therefore, although T a is needed in the NPS, given the high accuracy achieved in T a estimates and the insensitivity of EF estimates to T a , it is reasonable to conclude that the uncertainty brought in by T a is very little.Besides these five temperate related parameters, the fractional vegetation cover ( f c ) is another indispensable parameter for the retrieval of EF.A set of values with a step of 0.02 were adopted in Figure 14 to define the variation of f c .Positive linear correlation relationships between EF retrievals and f c are observed in both of these two schemes.An increase (decrease) of 0.02 in f c will result in an increase (decrease) of 0.011 and 0.013 in EF, respectively.
Remote Sens. 2017, 9, 26 20 of 25 variation of is in the opposite directions.The TPS is more sensitive to the variation of than the NPS when the variation of is negative.However, when the variation of is positive, the NPS is more sensitive to the variation of than the TPS.Compared with and , the sensitivity of EF retrievals to the rest three parameters is not very significant.In general, a ±6 K change of these three parameters results in only ± 0.07 change of EF, but the sensitivity of the TPS to increases significantly with the increase of .It has even exceeded the sensitivity to when the variation of temperature is beyond +5.5 K.The results presented in Figure 4 indicate that over the SGP region is estimated with the accuracy of a RMSE as 3.02 K. Figure 13 shows that a change of ±3.02 K in results in only ± 0.03 change in EF.Therefore, although is needed in the NPS, given the high accuracy achieved in estimates and the insensitivity of EF estimates to , it is reasonable to conclude that the uncertainty brought in by is very little.Besides these five temperate related parameters, the fractional vegetation cover ( ) is another indispensable parameter for the retrieval of EF.A set of values with a step of 0.02 were adopted in Figure 14 to define the variation of .Positive linear correlation relationships between EF retrievals and are observed in both of these two schemes.An increase (decrease) of 0.02 in will result in an increase (decrease) of 0.011 and 0.013 in EF, respectively.

Conclusions
In order to reduce the reliance of T s − V I triangle method on the dry edge, a new parameterization scheme of EF was developed in this paper by the combination of the traditional T s − V I triangle method and the TVX method.For validation purpose, it was demonstrated over the SGP region and compared with the traditional parameterization scheme proposed by Jiang and Islam [21,24].
Results show that similar accuracy is achieved by these two parameterization schemes.The R 2 , MAE, RMSE, RRMSE and B obtained by the NPS are 0.58, 0.11, 0.14, 0.33 and −0.03, respectively, and those obtained by the TPS are 0.59, 0.11, 0.14, 0.33 and −0.03, respectively.A further comparison with other relevant EF estimation studies also shows that the accuracy achieved by the proposed NPS is comparable to that produced by the TPS in previous studies.Although there is little difference in the accuracy produced, there are three advantages of the proposed new parameterization scheme, which means that it can be used more widely and simply.(1) The maximum surface temperature of the full cover vegetation (T cmax ), a key parameter of the dry edge in the TPS, is not needed in the NPS.Thus, the proposed new parameterization scheme can be performed independent of the determination of the dry edge, and only the maximum surface temperature of the bare soil (T smax ) is indispensable.Consequently, this has not only bypassed the task involved in the determination of dry edge, but also reduced the uncertainty brought in by T cmax .(2) Compared with the TPS, the Priestley-Taylor parameter ∅ of a mixed pixel in the NPS is retrieved through the differentiation of soil and vegetation components.The parameter ∅ the full vegetated canopy and bare soil is estimated, respectively.It means that two-source ET models can be developed from the new parameterization scheme.(3) By regarding the dry edge of the T s − V I trapezoid space as a special isopleth of the soil surface moisture, the proposed new parameterization scheme can determine the constraint T smax by itself, which can even be extended as the dry edge of the T s − V I trapezoid space.The results of validation show that for the NPS, similar accuracy is achieved by using T smax determined above as the constraint.Thus, the simplicity of the proposed new parameterization scheme does not compromise its accuracy in monitoring EF.
The disadvantages of the new parameterization scheme mainly lie in the uncertainty it involves.The NPS is developed by the combination of the T s − V I triangle method and TVX method, so the uncertainty in both of these two methods would be introduced to the final estimates of EF.The uncertainties involved in the TVX method are mainly concentrated in the linearization of the isopiestic lines of soil moisture.Specifically, the radiometric temperature of a full vegetated canopy is assumed to be in equilibrium with the temperature of the air within the canopy.Strictly speaking, this assumption only holds true under optimal soil moisture conditions.Investigation on the uncertainty caused by this assumption needs further studies based on soil moisture observations.As for the T s − V I triangle method, besides the empirical approaches used to determine the extreme values of EF along each isopiestic line, the linear interpolation of Priestley-Taylor parameter ∅ also brings in some uncertainty [66].Moreover, considering the sensitivity of EF estimates to T smax , additional work based on the surface energy balance principle may be required to calculate the T smax of the theoretical dry edge.

Figure 1 .
Figure 1.Distribution of the Bowen ratio stations over SGP region on a United States Geological Survey (USGS) land cover map.

Figure 1 .
Figure 1.Distribution of the Bowen ratio stations over SGP region on a United States Geological Survey (USGS) land cover map.

Figure 2 .
Figure 2. A conceptual sketch of the theoretical trapezoid   −   feature space.Upper red line represents the dry edge under extreme water stressed conditions.Horizontal blue lines represents the wet edge under maximum soil moisture availability.Slanting green lines enclosed by the boundaries represent superimposed isopleths of soil moisture availability.Point A (  = 0,   =   ) represents the driest bare soil surface with the Priestley-Taylor parameter ∅ equal to zero, point B (  = 0,   =   ) represents the wettest bare soil surface with ∅ equal to 1.26, point C (  = 1,   =   ) represents the driest fully vegetated surface with the Priestley-Taylor parameter ∅ equal to zero, and point D (  = 1,   =   ) represents the wettest fully vegetated surface with the Priestley-Taylor parameter ∅ equal to 1.26.

Figure 2 .
Figure 2. A conceptual sketch of the theoretical trapezoid T s − f c feature space.Upper red line represents the dry edge under extreme water stressed conditions.Horizontal blue lines represents the wet edge under maximum soil moisture availability.Slanting green lines enclosed by the boundaries represent superimposed isopleths of soil moisture availability.Point A ( f c = 0, T s = T smax ) represents the driest bare soil surface with the Priestley-Taylor parameter ∅ equal to zero, point B ( f c = 0, T s = T smin ) represents the wettest bare soil surface with ∅ equal to 1.26, point C ( f c = 1, T s = T cmax ) represents the driest fully vegetated surface with the Priestley-Taylor parameter ∅ equal to zero, and point D ( f c = 1, T s = T cmin ) represents the wettest fully vegetated surface with the Priestley-Taylor parameter ∅ equal to 1.26.

Figure 3 .
Figure 3. Flowchart of the proposed new parameterization scheme of EF.

Figure 3 .
Figure 3. Flowchart of the proposed new parameterization scheme of EF.

Figure 4 .
Figure 4. Scatterplots of estimated and measured   .Figure 4. Scatterplots of estimated and measured T a .

Figure 4 .
Figure 4. Scatterplots of estimated and measured   .Figure 4. Scatterplots of estimated and measured T a .

Figure 5 .
Figure 5.Comparison between observed EF and estimated EF generated by the TPS (a) and the NPS (b).

Figure 5 .
Figure 5.Comparison between observed EF and estimated EF generated by the TPS (a) and the NPS (b).

Figure 6 .
Figure 6.Contrast of EF retrieved from TPS and NPS for all the sites.

Figure 7 .
Figure 7. Temporal variations of observed EF and estimated EF generated by the TPS and NPS for 11 sites.

Figure 7 .
Figure 7. Temporal variations of observed EF and estimated EF generated by the TPS and NPS for 11 sites.

Figure 8 .
Figure 8. Spatial distribution of EF retrieved from the TPS (a) and NPS (b) as well as their differences (c) on DOY 239.

Figure 8 .
Figure 8. Spatial distribution of EF retrieved from the TPS (a) and NPS (b) as well as their differences (c) on DOY 239.

Figure 9 .
Figure 9. Contrast of EF retrieved from TPS and NPS over the whole SGP region on DOY 239.

25 -Figure 10 .
Figure 10.Scatterplots of the differences in EF retrieved from the TPS and NPS versus NDVI (a), and scatterplots of the differences in EF retrieved from the TPS and NPS versus   (b).

Figure 11 .
Figure 11.Temporal variations of the spatial quantile statistics of EF retrieved from these two parameterization schemes as well as the spatial quantile statistics of their differences (EF NPS − EF TPS ) at daily scale.

Figure 10 .Figure 10 .
Figure 10.Scatterplots of the differences in EF retrieved from the TPS and NPS versus NDVI (a), and scatterplots of the differences in EF retrieved from the TPS and NPS versus T s (b).

Figure 11 .
Figure 11.Temporal variations of the spatial quantile statistics of EF retrieved from these two parameterization schemes as well as the spatial quantile statistics of their differences (EF − EF ) at daily scale.

Figure 11 .
Figure 11.Temporal variations of the spatial quantile statistics of EF retrieved from these two parameterization schemes as well as the spatial quantile statistics of their differences (EF NPS − EF TPS ) at daily scale.

Figure 12 .
Figure 12.Comparison of determined by the TVX method and regression approach (a), and comparison of determined by the TVX method and regression approach (b).

Figure 12 .
Figure 12.Comparison of T smax determined by the TVX method and regression approach (a), and comparison of T cmax determined by the TVX method and regression approach (b).

Figure 13 .
Figure 13.Sensitivity of EF retrievals to the variation of temperature related parameters.Figure 13.Sensitivity of EF retrievals to the variation of temperature related parameters.

Figure 13 .
Figure 13.Sensitivity of EF retrievals to the variation of temperature related parameters.Figure 13.Sensitivity of EF retrievals to the variation of temperature related parameters.

Figure 13 .
Figure 13.Sensitivity of EF retrievals to the variation of temperature related parameters.

Figure 14 .
Figure 14.Sensitivity of EF retrievals to the variation of   .Figure 14.Sensitivity of EF retrievals to the variation of f c .

Figure 14 .
Figure 14.Sensitivity of EF retrievals to the variation of   .Figure 14.Sensitivity of EF retrievals to the variation of f c .

Table 2 .
Location, altitude and land cover of study sites.
and   are the y-axis of Point A, C, and D, respectively. , and  , are the maximum and minimum   for a given fractional vegetation cover  , , and  , is the land surface temperature for each pixel within this   class.

Table 3 .
Statistical measures.Note: P i is predicted, O i is observed and n is the number of samples.

Table 4 .
Accuracy of EF retrieved from TPS and NPS for all the sites.

Table 4 .
Accuracy of EF retrieved from TPS and NPS for all the sites.
Figure 6.Contrast of EF retrieved from TPS and NPS for all the sites.

Table 5 .
Summary of relevant EF estimation studies based on the T s − V I triangle method.

Table 6 .
Accuracy of EF retrieved from TPS and NPS within the dry edge determined by the TVX method.