Intercomparison of Three Two-Source Energy Balance Models for Partitioning Evaporation and Transpiration in Semiarid Climates

Evaporation (E) and transpiration (T) information is crucial for precise water resources planning and management in arid and semiarid areas. Two-source energy balance (TSEB) methods based on remotely-sensed land surface temperature provide an important modeling approach for estimating evapotranspiration (ET) and its components of E and T. Approaches for accurate decomposition of the component temperature and E/T partitioning from ET based on TSEB requires careful investigation. In this study, three TSEB models are used: (i) the TSEB model with the Priestley-Taylor equation, i.e., TSEB-PT; (ii) the TSEB model using the Penman-Monteith equation, i.e., TSEB-PM, and (iii) the TSEB using component temperatures derived from vegetation fractional cover and land surface temperature (VFC/LST) space, i.e., TSEB-TC-TS. These models are employed to investigate the impact of component temperature decomposition on E/T partitioning accuracy. Validation was conducted in the large-scale campaign of Heihe Watershed Allied Telemetry Experimental Research-Multi-Scale Observation Experiment on Evapotranspiration (HiWATER-MUSOEXE) in the northwest of China, and results showed that root mean square errors (RMSEs) of latent and sensible heat fluxes were respectively lower than 76 W/m2 and 50 W/m2 for all three approaches. Based on the measurements from the stable oxygen and hydrogen isotopes system at the Daman superstation, it was found that all three models slightly overestimated the ratio of E/ET. In addition, discrepancies in E/T partitioning among the three models were observed in the kernel experimental area of MUSOEXE. Further intercomparison indicated that different temperature decomposition methods were responsible for the observed discrepancies in E/T partitioning. The iterative procedure adopted by TSEB-PT and TSEB-PM produced higher LEC and lower TC when compared to TSEB-TC-TS. Overall, this work provides valuable insights into understanding the performances of TSEB models with different temperature decomposition mechanisms over semiarid regions.


Introduction
Evapotranspiration (ET) observations and modeling are crucial in water cycle studies [1][2][3][4]. Water scarcity is one of the main factors constraining agricultural development in arid or semiarid areas. The knowledge of ET, as well as the mechanism of ET partitioning into evaporation (E) and transpiration (T), is very important for precise quantification of the water balance in water resources planning and management, optimizing crop production, identifying crop stress and drought impacts, and evaluating the effects of climate change on water yields [5][6][7][8][9]. As satellite remote sensing is widely used to obtain information on the regional water and heat balance, it has been used to derive several global multi-decadal ET datasets that arouse extensive concern [10][11][12][13][14][15].
Over the last few decades, several remote sensing-based ET models have been proposed to estimate regional surface heat fluxes via satellite-derived land surface temperature (LST) [5,6,[16][17][18][19][20]. Specifically, a one-source model, such as the Surface Energy Balance Algorithm for Land (SEBAL) [16,17], Simplified Surface Energy Balances Index (S-SEBI) [21], Surface Energy Balance System (SEBS) [20,22], and Mapping EvapoTranspiration at high Resolution with Internalized Calibration (METRIC) [5,6], often use LST and empirical resistance corrections to estimate surface heat fluxes. Two-source models, such as the two-source energy balance model (TSEB) [19], the two-source time-integrated model (TSTIM) [23], Pixel Component Arranging and Comparing Algorithm (PCACA) [24], and the enhanced two-source evapotranspiration model for land (ETEML) [25] were developed to make use of LST to estimate a sensible heat flux (H) and latent heat flux (LE) for the soil and canopy separately. Other alternative practical approaches have been proposed based on vegetation fractional cover and LST (VFC/LST) space [26,27]. Related research includes the work by Moran et al. [28], Jiang and Islam [29,30], Stisen et al. [31], and Shu et al. [32]. Extensive reviews of remote sensing-based methodologies on surface heat flux estimation can be found in the works of Courault et al. [33], Kalma et al. [34] and Li et al. [35].
As an appealing modeling scheme, two-source energy balance models can estimate evapotranspiration (ET), evaporation (E), and transpiration (T) of vegetated surfaces, which has important applications in water resources management for irrigated crops in arid and semiarid areas [36]. To estimate E and T separately, the soil temperature (T S ) and canopy temperature (T C ) need to be accurately measured or derived, making temperature decomposition the core of the two-source modeling approach. Generally, four categories of methods have been developed to decompose remotely-sensed temperature into soil and canopy temperatures. (1) The first category calculates soil and vegetation component temperatures using dual-or multi-angular thermal infrared measurements [37]. For instance, a dual-source model using bi-angular thermal infrared measurements was developed by Jia [38] and Jia et al. [39]. However, currently limited sensors have bi-or multi-angular thermal infrared channels and thus constrained the availability of remote sensing data source in this category of method. (2) The second category assumed that LST is the sum of vegetation and soil temperatures weighted by the vegetation fractional cover (f c ), i.e., LST ≈ f c T C + T S (1 − f c ). This kind of method assumed that the LST of a highly vegetated pixel in the neighborhood would be a reasonable approximation of T C , and a selected neighborhood area of the target pixel based on certain thresholds in the corresponding NDVI image [40,41]. As T C is approximately determined, T S can be derived accordingly using T S = (LST − f c T C )/(1 − f c ). (3) The third category retrieves component temperatures using an iterative procedure based on an energy balance and resistance network. Typically, the TSEB model proposed by Norman et al. [19] uses a system of temperature gradient-resistance equations that are solved by an iterative procedure, in which an initial estimate of plant transpiration is determined by the Priestley-Taylor equation (TSEB-PT). The procedure terminates when soil evaporation (LEs) exceeds zero, and Tc and Ts are recalculated until energy balance closure is reached. Recently, TSEB was revised by Colaizzi et al. [42] and Colaizzi et al. [43] using the Penman-Monteith equation (TSEB-PM) instead of the Priestley-Taylor equation in order to characterize vegetation transpiration, and this revised version was found to be more applicable for advective semiarid climates. space [24,28,44,45]. Moran et al. [28] introduced a water deficit index (WDI) based on the VFC/LST trapezoid space and extended the application of a crop water stress index (CWSI) over fully-to partially-vegetated surface areas. Zhang et al. [24] proposed a method to retrieve land surface component temperatures based on the interpretation of soil wetness isolines within a VFC/LST trapezoidal space. The main assumption of this method is that the isolines in the VFC/LST trapezoid space is mainly controlled by soil water availability and the isolines within a VI/LST trapezoidal space are used to decompose the composite temperature into component temperatures. This method was recently revised and adopted by Sun et al. [46], Long and Singh [47], Yang and Shang [48], Yang et al. [25], and Sun [49] to develop two-source models.
Combined with satellite remote sensing, the TSEB models have been extensively used to estimate regional ET. However, how the accuracy of ET estimations from TSEB models is affected by its component temperature estimation and E/T partitioning are rarely reported, and requires careful study. The objective of this study is to evaluate the capabilities of three TSEB models in predicting surface fluxes under various ranges of soil moisture contents, and then analyze how the different performances of three models are attributed to the different adopted schemes in component temperature decomposition and E/T partitioning. Section 2 presents the description of the three two-source models: TSEB-PT, TSEB-PM, and TSEB-T C -T S . Section 3 introduces the pertinent experiment campaign in the study area, i.e., Heihe Watershed Allied Telemetry Experimental Research-The Multi-Scale Observation Experiment on Evapotranspiration (HiWATER-MUSOEXE) Campaign, remotely-sensed data used for driving the TSEB models, and the ground measurements for the models' assessment. Section 4 first reports the fluxes estimation accuracies from three TSEB models with retrievals from ASTER imagery. In Section 4.2, E/T partitioning is intercompared within three models and evaluated against the benchmark obtained using the stable oxygen and hydrogen isotopes technique. In the remainder of Section 4, intercomparison using component temperatures generated from three models was conducted to further investigate the uncertainty in E/T partitioning. Section 5 discusses the advantages and limitations in this work, and the final section provides a conclusion.

TSEB-PT Model
The TSEB model was originally developed by Norman et al. [19] to make use of remotely-sensed radiometric surface temperatures to estimate soil evaporation and canopy transpiration. This model has been modified by Kustas and Norman [50,51] through improving the soil surface resistance formulation and net radiation partitioning between the soil and canopy components. The net radiation is partitioned between the vegetated canopy and soil, and can be expressed as: where R n is net radiation (W/m 2 ), R ns and R nc are the net radiation for soil and vegetation canopy (W/m 2 ), respectively; H and LE are the sensible and latent heat fluxes (W/m 2 ), respectively, and G is the soil heat flux (W/m 2 ). The energy balance for the soil and vegetated canopy can be expressed as: where H s and H c are the sensible heat fluxes for the soil and canopy respectively (W/m 2 ), LE s and LE c are the latent heat fluxes for the soil and canopy, respectively (W/m 2 ). G is parameterized with the phase-difference approach proposed by Santanello and Friedl [52]: where t is the solar time angle (s), a is the amplitude parameter (dimensionless), b is the period length (s), and c is the shift (s). In this study, parameters a, b, and c take the values of 0.3, 86,400, and 10,800 following Colaizzi et al. [43] and Song et al. [53]. In this study, the series resistance network form was applied, in which H c , H s , and the sum of both terms are calculated as: where ρ is the air density (kg/cm 3 ), C P is the specific heat of air (J/kg/K), T S is the soil temperature (K), T C is the canopy temperature (K), T AC and T A are the air temperature within the canopy boundary layer and air temperature (K), respectively, r A is the aerodynamic resistance (s/m), r x is the resistance in the boundary layer near the canopy (s/m), and r s is the resistance to heat flux in the boundary layer above the soil surface (s/m). r A , r x , and r s are calculated according to Norman et al. [19] and Kustas and Norman [50].
The TSEB-PT model uses a modified Priestley-Taylor formulation to parameterize the canopy transpiration: where α PT is the Priestley-Taylor parameter (dimensionless), f G is the fraction of green vegetation (dimensionless), ∆ is the slope of the saturation vapor pressure versus temperature curve (kPa/ • C), and γ is the psychrometric constant (kPa/ • C). An initial estimate of T C can be derived as follows: Accordingly, T S is calculated with an initial estimate of T S , and then r s can be estimated with the temperature gradient between the soil and canopy described in Kustas and Norman [51]. From Equations (5) to (8), the component sensible heat flux H S can be calculated and latent heat fluxes from canopy LE C and soil surface LE S are solved as residual terms. In order to obtain a realistic estimation of surface heat fluxes under water stressed conditions, the α PT is iteratively decreased until LEs exceeds zero. The detailed description of the TSEB model and the parameterization of the resistance network can be found in Norman et al. [19] and Kustas and Norman [51].

TSEB-PM Model
The TSEB model was revised by Colaizzi et al. [42] and Colaizzi et al. [43] using the Penman-Monteith equation instead of the Priestley-Taylor formulation to account for the impact of advection over semiarid climates. This revised version of the TSEB model is termed TSEB-PM. The effects of varying the vapor pressure deficit can thus be taken into account in the TSEB-PM model. The canopy transpiration is characterized using the Penman-Monteith equation: and T c is initialized as follows: where γ * = γ(1 + r c /r A ), r c is the bulk canopy resistance (s/m), r A is the aerodynamic resistance between the canopy and the air above the canopy (s/m), and e a and e s are the actual and saturation vapor pressure of the air (kPa), respectively. Similar to TSEB-PT, the TSEB-PM model was iteratively implemented as described in Section 2.1. During the iterative procedure, r c increases from 10 s/m with an increment of 20 s/m and terminates at 5000 s/m, or until LEs exceeds zero. The comprehensive introduction of the TSEB-PM can be found in Colaizzi et al. [42] and Colaizzi et al. [43].

TSEB-T C -T S Model
TSEB was modified by Kustas and Norman [54] to calculate turbulent heat fluxes using canopy and soil component temperatures that were measured or derived from other methods. This modified TSEB model is called TSEB-T C -T S in this study. The Priestley-Taylor iteration procedure is not applied in TSEB-T C -T S , but the remainder of the physical framework of TSEB-T C -T S is identical to that of TSEB-PT. The within-canopy temperature (T AC ) is estimated from derived component temperatures as follows: Consequently, the component sensible heat fluxes H S and H C are directly calculated from Equations (5) and (6), and the component latent heat fluxes LE C and LE S are calculated as residual terms from Equations (2) and (3).
The VFC/LST space is adopted to retrieve component temperatures by Zhang et al. [24], Zhang et al. [55], Sun et al. [40], Merlin et al. [44], and Merlin et al. [45]. Based on the theoretical determination of the dry and wet edges of the VFC/LST space, this method was further revised and adopted by Long and Singh [47] and Yang and Shang [48] to develop two-source models. In the traditional method, the VFC/LST trapezoid is mainly determined manually based on the 2-D VFC/LST scatter plot, which would cause great uncertainty. Recently, Yang et al. [25] improved this method and proposed to make use of the VFC/LST trapezoid space for each pixel. Different from traditional method, the four theoretical points for each pixel are determined based on an energy balance model and the Penman-Monteith equations. In this study, this pixel-wise surface temperature decomposition method is adopted. The employed VFC/LST space is determined with four theoretical points using the following equations from Moran et al. [28].
For the vertex of dry bare soil, the difference between LST and Ta can be derived as follows: For the vertex of saturated bare soil, For the vertex of well-watered fully-cover vegetation, For the vertex of water-stressed fully-cover vegetation, where VPD (kPa) is the vapor pressure deficit of the air at temperature T a , γ is the psychrometric constant (kPa/ • C), and r a is aerodynamic resistance (s/m) estimated with the equations proposed by Thom [56]. r cp is the canopy resistance at potential evapotranspiration (s/m), and r cx is the maximum canopy resistance (s/m). As such, a VFC/LST trapezoid space is determined for each pixel and the measurements of (LST-Ta) and the fractional vegetation cover would be located within the trapezoid space. Based on the VFC/LST trapezoid space determined above, the composite radiometric surface temperature of each pixel can be decomposed to soil and canopy component temperatures [24,25,55]: where K s is the slope of the isoline that passes through the point located in VFC/LST space, and can be derived by interpolating the slope of the warm edge and that of the cold edge. The detailed description of this method can be found in Yang et al. [25]. Once T S and T C are determined, the component fluxes in TSEB-T C -T S can be estimated.
In addition to T S and T C , TSEB models can estimate evaporation (E) and transpiration (T) of vegetated surface, such information can be used to calculate component water stress of vegetation and soil. Following the equations from Yang et al. [25], the crop water stress index for the canopy component (CWSIc) and soil water deficit index for soil component (SWDI S ) are calculated as follows: EP c and EP s are the potential transpiration rate and potential soil evaporation rate, respectively.

HiWATER-MUSOEXE Campaign and Ground-Based Measurements
Heihe Watershed Allied Telemetry Experimental Research (HiWATER) was performed at the Heihe River Basin of northwestern China with airborne, satellite-borne, and ground-based remote sensing experiments at various scales during 2012-2015 [57]. As one of most important thematic experiments of the HiWATER, the Multi-Scale Observation Experiment on Evapotranspiration (MUSOEXE) over the Zhangye oasis provided multiscale data sets of meteorological elements and land surface parameters that facilitate the development and validation of ET models over heterogeneous surfaces [58,59]. Figure 1 shows the distribution of flux towers and the land use classifications in the MUSOEXE.
MUSOEXE involved a multi-scale observation campaign over heterogeneous surfaces by using an observation matrix composed of 21 stations. Each station was equipped with a set of eddy covariance (EC) system and automatic weather system (AWS). Turbulent heat fluxes are measured with EC system and the raw data were processed using EdiRe software and averaged over 30 min. Wind speed, wind direction, air temperature, vapor pressure, net radiation, and atmospheric pressure were measured in AWS with 10-min intervals. The soil heat fluxes were measured using three heat flux plates located 6 cm below the ground's surface at each site. In order to better represent the surface soil heat flux, the Thermal Diffusion Equation and Correction (TDEC) method proposed by Yang and Wang [60] was applied to correct the observed soil heat flux with the observed soil moisture and temperature profile. The residual method [61] was adopted to adjust sensible and latent heat fluxes by forcing the energy balance closure.
The isotopic composition of atmospheric water vapor provided rich information on the hydrological cycle and gaseous exchange processes between the terrestrial vegetation and atmosphere. Due to technical and instrumental limitations, the measurements of the isotopic composition of water vapor were limited to discrete campaigns and discrete samples. During the HiWATER program, the isotopic composition of water vapor from the surface air was measured in a corn cropland (100.3722 • E, 38.8555 • N) at Daman superstation using a flux gradient method, using a cavity ring-down spectroscopy (CRDS) water vapor isotope analyzer. Considering water vapor as a mixture of ET from an ecosystem that carried unique isotopic signals from plant transpiration and soil evaporation separately, the measured isotopic composition of water was used to partition ET into evaporation and transpiration. Details of the isotope experiment and its calibration procedure are given by Huang and Wen [62] and Wen et al. [63]. In this study, the ratio of T over ET was collected to evaluate the reliability of three TSEB models on evaporation and transpiration partitioning. and are the potential transpiration rate and potential soil evaporation rate, respectively.

HiWATER-MUSOEXE Campaign and Ground-Based Measurements
Heihe Watershed Allied Telemetry Experimental Research (HiWATER) was performed at the Heihe River Basin of northwestern China with airborne, satellite-borne, and ground-based remote sensing experiments at various scales during 2012-2015 [57]. As one of most important thematic experiments of the HiWATER, the Multi-Scale Observation Experiment on Evapotranspiration (MUSOEXE) over the Zhangye oasis provided multiscale data sets of meteorological elements and land surface parameters that facilitate the development and validation of ET models over heterogeneous surfaces [58,59]. Figure 1 shows the distribution of flux towers and the land use classifications in the MUSOEXE.

Remote Sensing Data and Derivation of Related Variables
In this study, nine scenes from the Advance Spaceborne Thermal Emission and Reflection Radiometer (ASTER) on board Terra during the experiment period were collected. These nine scenes were acquired on 15 June, 24 June, 10 July, 2 August, 11 August, 18 August, 27 August, 3 September, and 12 September of 2012 (DOYs: 167, 176, 192, 215, 224, 231, 240, 247, 256, respectively). The LST and the land surface albedo were provided by "Heihe Plan Science Data Center, National Natural Science Foundation of China" (http://www.heihedata.org). The LST data were retrieved by Li et al. [64] using a temperature/emissivity separation (TES) algorithm proposed by Gillespie et al. [65], combined with the water vapor scaling atmospheric correction method [66]. Land surface albedo was retrieved from a Charge Coupled Device (CCD) camera on HJ-1 satellite by Sun et al. [67].
Three data sets including the visible, near-infrared (NIR) bands of ASTER and albedo from HJ-1, were all re-sampled to 90 m to be consistent with the thermal infrared band in spatial resolution. In addition, the leaf area index (LAI), crop height, and fractional vegetation cover for HiWATER-MUSOEXE were derived based on the empirical relationships proposed by Yang et al. [68].

Validation of Three TSEB Models over MUSOEXE
TSEB-PT, TSEB-PM, and TSEB-T C -T S were applied to the Zhangye oasis using ground-based and satellite-derived observations introduced in Section 3, and the model performances were evaluated using flux measurements from the MUSOEXE observation matrix. As a first step of model validation, the flux estimates were averaged over the upwind source area for each flux tower [69], and flux Remote Sens. 2018, 10, 1149 8 of 20 measurements were linearly interpolated to temporally match the time of satellite overpass. Figure 2 shows the validation scatter plot for each energy balance component (R n , G, LE, and H) from TSEB-PT, TSEB-PM, and TSEB-T C -T S estimations. Validation statistics comparing the three models' performances are summarized in Table 1. Results show that estimated R n from all three models were in good agreement with tower observations, and the absolute mean biases in estimated R n for three models are all below 10 W/m 2 . Besides, all models slightly overestimated G, with the RMSEs all slightly exceeding 37 W/m 2 . Despite the different parameterization schemes used in three models (for instance, both TSEB-PT and TSEB-PM use the iterative procedure to calculate component surface temperatures, whereas TSEB-T C -T S employs the theoretical VFC/LST trapezoid for component temperature decomposition), they all exhibited comparable skills in estimation of H and LE, which can be indicated by the similar RMSEs (≈44.9-47.9 W/m 2 for H, and ≈61. 8 The spatial distributions of H and LE over the Zhangye oasis based on TSEB-PT, TSEB-PM, and TSEB-TC-TS for the satellite overpass on July 10 of 2012 are shown in Figure 3. Generally, the spatial patterns of surface fluxes were similar between the three models. In addition, the contrasting features over the oasis and the surrounding sandy and Gobi desert were clearly observed from all three models. Specifically, the Zhangye oasis, which mainly comprises of irrigated farmland, exhibited an average LE over 400 W/m 2 . On the contrary, across the sandy and Gobi desert, LE was typically below 300 W/m 2 and H was over 100 W/m 2 .    The spatial distributions of H and LE over the Zhangye oasis based on TSEB-PT, TSEB-PM, and TSEB-T C -T S for the satellite overpass on July 10 of 2012 are shown in Figure 3. Generally, the spatial patterns of surface fluxes were similar between the three models. In addition, the contrasting features Remote Sens. 2018, 10, 1149 9 of 20 over the oasis and the surrounding sandy and Gobi desert were clearly observed from all three models. Specifically, the Zhangye oasis, which mainly comprises of irrigated farmland, exhibited an average LE over 400 W/m 2 . On the contrary, across the sandy and Gobi desert, LE was typically below 300 W/m 2 and H was over 100 W/m 2 . In summary, all three employed models performed similarly in estimating H and LE despite using different schemes for modeling canopy transpiration and deriving component surface temperatures. With substantial ground observations from the tower-based network, the performances of the three models over MUSOEXE were proven reliable. The difference in derived component temperatures and partitioned evaporation and transpiration among the three models are analyzed in the following two sections.

Intercomparison of E/T Partitioning from Three TSEB Models
The spatial distributions of canopy transpiration (LEc) and soil evaporation (LEs) based on TSEB-PT, TSEB-PM, and TSEB-TC-TS on 10 July of 2012 are shown in Figure 4. The spatial patterns of LEc and LEs derived from the three models are similar, and the range of LEc for irrigated farmland was about ≈350-500 W/m 2 . With respect to LEs, the difference derived from the three models can be visually discerned from Figure 4. In comparison with TSEB-PT and TSEB-PM, TSEB-TC-TS tended to produce higher LEs, especially over the sparsely vegetated area around the residential area and the sandy Gobi desert.
In this study, both CWSIc and SWDIs were further derived to compare the performances of the three TSEB models on detecting vegetation and soil water stresses. The spatial distributions of CWSIc and SWDIs on 10 July 2012 based on three models are shown in Figure 5. The sandy and Gobi desert pixels are masked in CWSIc images in order to more clearly reveal the difference of the three TSEB models in detecting vegetation water stress over the oasis. TSEB-PM and TSEB-PT show similar performances regarding the detection of vegetation stress with both CWSIc values close to zero. On the contrary, TSEB-TC-TS seemed to detect a higher level of vegetation water stress, with a CWSIc peak at 0.45. In addition, the spatial heterogeneity was more prominent in TSEB-TC-TS, since the CWSIc gradient in the oasis and urban areas are clearly seen from the upper right subplot of Figure 5. With respect to SWDIs, the spatial distributions were quite similar between the three TSEB models. The SWDIs values for the oasis was smaller than the surrounding sandy and Gobi desert, and the contrasting features over the oasis and the surrounding sandy and Gobi desert were clearly observed from SWDIs subplots from all three models. In summary, all three employed models performed similarly in estimating H and LE despite using different schemes for modeling canopy transpiration and deriving component surface temperatures.
With substantial ground observations from the tower-based network, the performances of the three models over MUSOEXE were proven reliable. The difference in derived component temperatures and partitioned evaporation and transpiration among the three models are analyzed in the following two sections.

Intercomparison of E/T Partitioning from Three TSEB Models
The spatial distributions of canopy transpiration (LEc) and soil evaporation (LEs) based on TSEB-PT, TSEB-PM, and TSEB-T C -T S on 10 July of 2012 are shown in Figure 4. The spatial patterns of LEc and LEs derived from the three models are similar, and the range of LEc for irrigated farmland was about ≈350-500 W/m 2 . With respect to LEs, the difference derived from the three models can be visually discerned from Figure 4. In comparison with TSEB-PT and TSEB-PM, TSEB-T C -T S tended to produce higher LEs, especially over the sparsely vegetated area around the residential area and the sandy Gobi desert.
In this study, both CWSIc and SWDIs were further derived to compare the performances of the three TSEB models on detecting vegetation and soil water stresses. The spatial distributions of CWSIc and SWDIs on 10 July 2012 based on three models are shown in Figure 5. The sandy and Gobi desert pixels are masked in CWSIc images in order to more clearly reveal the difference of the three TSEB models in detecting vegetation water stress over the oasis. TSEB-PM and TSEB-PT show similar performances regarding the detection of vegetation stress with both CWSIc values close to zero. On the contrary, TSEB-T C -T S seemed to detect a higher level of vegetation water stress, with a CWSIc peak at 0.45. In addition, the spatial heterogeneity was more prominent in TSEB-T C -T S , since the CWSIc gradient in the oasis and urban areas are clearly seen from the upper right subplot of Figure 5.
With respect to SWDIs, the spatial distributions were quite similar between the three TSEB models. The SWDIs values for the oasis was smaller than the surrounding sandy and Gobi desert, and the contrasting features over the oasis and the surrounding sandy and Gobi desert were clearly observed from SWDIs subplots from all three models.  The ratios of LEC/LE (i.e., T/ET) measured using the stable oxygen and hydrogen isotopes technique during HiWATER-MUSOEXE are utilized to evaluate the performances of three TSEB models on partitioning E and T. Figure 6     The ratios of LEC/LE (i.e., T/ET) measured using the stable oxygen and hydrogen isotopes technique during HiWATER-MUSOEXE are utilized to evaluate the performances of three TSEB models on partitioning E and T. Figure 6   The ratios of LE C /LE (i.e., T/ET) measured using the stable oxygen and hydrogen isotopes technique during HiWATER-MUSOEXE are utilized to evaluate the performances of three TSEB models on partitioning E and T. Figure 6 shows the intercomparison of LE C /LE between the three models and ground measurements. The ratios of LE C /LE were underestimated on DOYs of 176, 192, 215, 224, 231, and 240, while slight overestimation of LE C /LE occurred on DOYs of 247 and 256 in all three models. The LE C /LE ratios derived from three models were very close in most cases at this site. The mean observed LE C /LE ratio was 84.7%, while the mean estimated LE C /LE ratios were 76.7%, 76.9% and 77.0% for TSEB-PT, TSEB-PM, TSEB-T C -T S . All three models seemed to slightly underestimate the LE C /LE ratio. However, the observed LE C /LE ratios exhibited a decline during September (DOYs of 247 and 256), mainly due to leaf senescence, which was not characterized by all three models. The mean observed LEC/LE ratio was 84.7%, while the mean estimated LEC/LE ratios were 76.7%, 76.9% and 77.0% for TSEB-PT, TSEB-PM, TSEB-TC-TS. All three models seemed to slightly underestimate the LEC/LE ratio. However, the observed LEC/LE ratios exhibited a decline during September (DOYs of 247 and 256), mainly due to leaf senescence, which was not characterized by all three models. In order to further investigate the difference regarding E and T partitioning between the three models, a pixel-based comparison of LEC and LES in the kernel experimental area of MUSOEXE was conducted and shown in Figure 7, with color shading indicating pixel density. The statistics for a pixel-based comparison of LEC and LES in the kernel experimental area are summarized in Table 2. Most points are under the 1:1 line in Figure 7a,b, which suggests that TSEB-PM tended to produce higher LEC than TSEB-PT and TSEB-TC-TS. The mean differences (MDs) for the pairs of TSEB-PT/TSEB-PM and TSEB-TC-TS/TSEB-PM were 2.9 W/m 2 and 18.1 W/m 2 , respectively (mean differences were calculated by subtracting the former model of the pair from the latter model). LEC estimated from TSEB-PT was comparatively higher than that from TSEB-TC-TS, with MD being 15.  In order to further investigate the difference regarding E and T partitioning between the three models, a pixel-based comparison of LE C and LE S in the kernel experimental area of MUSOEXE was conducted and shown in Figure 7, with color shading indicating pixel density. The statistics for a pixel-based comparison of LE C and LE S in the kernel experimental area are summarized in Table 2. Most points are under the 1:1 line in Figure 7a,b, which suggests that TSEB-PM tended to produce higher LE C than TSEB-PT and TSEB-T C -T S . The mean differences (MDs) for the pairs of TSEB-PT/TSEB-PM and TSEB-T C -T S /TSEB-PM were 2.9 W/m 2 and 18.1 W/m 2 , respectively (mean differences were calculated by subtracting the former model of the pair from the latter model). LE C estimated from TSEB-PT was comparatively higher than that from TSEB-T C -T S , with MD being 15.2 W/m 2 for the TSEB-T C -T S /TSEB-PT pair. Correspondingly, comparing to TSEB-PT and TSEB-T C -T S , TSEB-PM tended to underestimate LEs as shown in Figure 7e,f, with MDs of −4.2 W/m 2 and −2.8 W/m 2 for the pairs of TSEB-PT/TSEB-PM and TSEB-T C -T S /TSEB-PM, respectively. The LE S estimations from TSEB-PT and TSEB-T C -T S for the kernel experimental area were similar, with a correlation coefficient R being approximately 1.0, MD being 1.4 W/m 2 , and the mean absolute difference (MAD) being 6.2 W/m 2 . Since the component temperature decomposition has a major impact on E/T partitioning, to explore the mechanism underlying the observed difference on LE C (LE S ) from three models, a further intercomparison on the derived component temperature is conducted in the next section.

Intercomparison of Tc and Ts Derived from Three TSEB Models
Intercomparison of the component temperatures derived from the three models is shown in Figure 8. The TC derived from TSEB-PM showed a relatively homogenous pattern over the Zhangye oasis and the average of TC approximated 301 K. On the contrary, TC derived from TSEB-PT and TSEB-TC-TS showed a much larger spatial variability, and the contrast between farmland and residential area is more discernable in both models. However, the spatial contrast between farmland and residential area was less significant for the TSEB-TC-TS derived Ts (Figure 8).  Table 2. Statistics summarizing the intercomparison of LE C and LE S derived from TSEB-PT, TSEB-PM, and TSEB-T C -T S in the kernel experimental area.

Intercomparison of Tc and Ts Derived from Three TSEB Models
Intercomparison of the component temperatures derived from the three models is shown in Figure 8. The T C derived from TSEB-PM showed a relatively homogenous pattern over the Zhangye oasis and the average of T C approximated 301 K. On the contrary, T C derived from TSEB-PT and TSEB-T C -T S showed a much larger spatial variability, and the contrast between farmland and residential area is more discernable in both models. However, the spatial contrast between farmland and residential area was less significant for the TSEB-T C -T S derived Ts (Figure 8). A pixel-based comparison of decomposed TC and TS in the kernel experimental area was conducted and is shown in Figure 9. It is noticed that the component temperatures were much more scattered compared to the component fluxes in the scatter plot of Figure 7. The statistics for the pixel-based comparison of TC and TS in the kernel experimental area are listed in Table 3. Generally, TSEB-TC-TS estimates higher TC compared to TSEB-PT and TSEB-PM, and TSEB-PT tended to overestimate TC in relation to TSEB-PM, with the MD of TSEB-PT/TSEB-PM pair being −0.2 K. A pixel-based comparison of decomposed T C and T S in the kernel experimental area was conducted and is shown in Figure 9. It is noticed that the component temperatures were much more scattered compared to the component fluxes in the scatter plot of Figure 7. The statistics for the pixel-based comparison of T C and T S in the kernel experimental area are listed in Table 3. Generally, TSEB-T C -T S estimates higher T C compared to TSEB-PT and TSEB-PM, and TSEB-PT tended to overestimate T C in relation to TSEB-PM, with the MD of TSEB-PT/TSEB-PM pair being −0.2 K. Table 3. Statistics summarizing the intercomparison of T C and T S derived from TSEB-PT, TSEB-PM, and TSEB-T C -T S in the kernel experimental area. Different strategies for deriving T C and T S were responsible for the observed differences in surface fluxes between the three modeling approaches. As previously stated, both TSEB-PT and TSEB-PM applied an iterative approach to derive component temperatures. Although the identical temperature decomposition method was applied to TSEB-PT and TSEB-PM, the component temperatures derived from both models showed a significant difference, with correlation coefficient R being only 0.13 for T C . Different from TSEB-PT and TSEB-PM, TSEB-T C -T S adopted the VFC/LST trapezoid space to estimate the component temperatures, which would explain the distinct characteristics shown in Figure 8. Specifically, for the irrigated farmland, the component temperatures derived from TSEB-T C -T S were comparable to those derived from the other two models. Due to the lack of constraint on adjusting the canopy transpiration, the vegetation component temperature from TSEB-PT and TSEB-PM would be very close to pixels in well-watered and fully-covered vegetation areas. Different strategies for deriving TC and TS were responsible for the observed differences in surface fluxes between the three modeling approaches. As previously stated, both TSEB-PT and TSEB-PM applied an iterative approach to derive component temperatures. Although the identical temperature decomposition method was applied to TSEB-PT and TSEB-PM, the component temperatures derived from both models showed a significant difference, with correlation coefficient R being only 0.13 for TC. Different from TSEB-PT and TSEB-PM, TSEB-TC-TS adopted the VFC/LST trapezoid space to estimate the component temperatures, which would explain the distinct characteristics shown in Figure 8. Specifically, for the irrigated farmland, the component temperatures derived from TSEB-TC-TS were comparable to those derived from the other two models. Due to the lack of constraint on adjusting the canopy transpiration, the vegetation component temperature from TSEB-PT and TSEB-PM would be very close to pixels in well-watered and fully-covered vegetation areas.

TSEB-PT vs. TSEB-PM TSEB-Tc-Ts vs. TSEB-PM TSEB-Tc-Ts vs. TSEB-PT
The scientific rationale can be explained as follows. Based on the VFC/LST space, Figure 10 illustrates the temperature decomposition methods adopted in TSEB-PT, TSEB-PM, and TSEB-TC-TS. Both TSEB-PT and TSEB-PM applied an iterative approach to derive the component temperatures and this approach assumed vegetation transpiration at the potential rate as an initial value. Due to the lack of constraint on adjusting the canopy transpiration, the vegetation component temperature from TSEB-PT and TSEB-PM would be very close to Point C ( Figure 10) in well-watered and fully vegetated areas. Among TSEB-PM and TSEB-PT, the former employs the Penman-Monteith equation to consider the varying VPD under advective semiarid climates, and this would lead to higher LEC compared to TSEB-PT. Consequently, the soil component temperature from TSEB-PM was higher than that from TSEB-PT. Different from TSEB-PM and TSEB-PT, TSEB-TC-TS adopted the VFC/LST trapezoid space to estimate component temperatures. TSEB-TC-TS assumed that the vegetation and soil share the same water pool, and the slope of each isoline in the VFC/LST space could be derived by interpolating the slopes of both dry and cold edges. The temperature decomposition methods adopted in the three TSEB models is illustrated in Figure 10, the soil surface temperatures derived from TSEB-PM, TSEB-PT, and TSEB-TC-TS are denoted by Ts1, Ts2, and Ts3, and the vegetation canopy temperatures derived from the same three models are denoted by Tv1, Tv2, and The scientific rationale can be explained as follows. Based on the VFC/LST space, Figure 10 illustrates the temperature decomposition methods adopted in TSEB-PT, TSEB-PM, and TSEB-T C -T S . Both TSEB-PT and TSEB-PM applied an iterative approach to derive the component temperatures and this approach assumed vegetation transpiration at the potential rate as an initial value. Due to the lack of constraint on adjusting the canopy transpiration, the vegetation component temperature from TSEB-PT and TSEB-PM would be very close to Point C ( Figure 10) in well-watered and fully vegetated areas. Among TSEB-PM and TSEB-PT, the former employs the Penman-Monteith equation to consider the varying VPD under advective semiarid climates, and this would lead to higher LE C compared to TSEB-PT. Consequently, the soil component temperature from TSEB-PM was higher than that from TSEB-PT. Different from TSEB-PM and TSEB-PT, TSEB-T C -T S adopted the VFC/LST trapezoid space to estimate component temperatures. TSEB-T C -T S assumed that the vegetation and soil share the same water pool, and the slope of each isoline in the VFC/LST space could be derived by interpolating the slopes of both dry and cold edges. The temperature decomposition methods adopted in the three TSEB models is illustrated in Figure 10, the soil surface temperatures derived from TSEB-PM, TSEB-PT, and TSEB-T C -T S are denoted by T s1 , T s2 , and T s3 , and the vegetation canopy temperatures derived from the same three models are denoted by T v1 , T v2 , and T v3 . It is clear that T v1 and T v2 were underestimated compare to T v3 , and this led to a higher T/ET partition for TSEB-PM and TSEB-PT compared to TSEB-T C -T S , which consequently overestimated LE C when compared to TSEB-T C -T S . A component temperatures decomposition method based on isolines in the VFC/LST space provided a further constraint for vegetation transpiration and would be a good substitution for an iterative approach.
Tv3. It is clear that Tv1 and Tv2 were underestimated compare to Tv3, and this led to a higher T/ET partition for TSEB-PM and TSEB-PT compared to TSEB-TC-TS, which consequently overestimated LEC when compared to TSEB-TC-TS. A component temperatures decomposition method based on isolines in the VFC/LST space provided a further constraint for vegetation transpiration and would be a good substitution for an iterative approach. Figure 10. Illustration of the temperature decomposition methods adopted in the three TSEB models. Ts1, Ts2, and Ts3 denote the soil surface temperatures derived from TSEB-PM, TSEB-PT, and TSEB-TC-TS, respectively, and Tv1, Tv2, and Tv3 denote the vegetation canopy temperatures derived from the same three models.

Reliability of the Employed TSEB Models in Estimating Surface Fluxes
The HiWATER-MUSOEXE data have been used extensively to validate the land surface flux models, including one-source and two-source models [70][71][72]. Using the same set of ground-based observations, Ma et al. [70] applied a revised SEBS model to estimate regional heat fluxes in the Heihe River Basin, and their assessment indicates that the RMSEs of the modeled H and LE are 56.9 W/m 2 and 74.8 W/m 2 , respectively. Huang et al. [71] integrated a Normalized Difference Water Index (NDWI) as a water stress index into SEBS through the modification of the parameter kB -1 , and showed RMSEs of 79.8 W/m 2 in H and 84.1 W/m 2 in LE for the revised SEBS. In this study, three two-source models were applied to the middle reach of the Heihe River Basin. The RMSEs of H and LE from TSEB-PT are 47.5 and 75.3 W/m 2 , respectively, the RMSEs of H and LE from TSEB-PM are 44.9 and 70.6 W/m 2 , respectively, and the values from TSEB-TC-TS are 47.9 and 61.8 W/m 2 , respectively. Overall, the performances of two-source modeling approaches were reliable in relation to previously published studies.

Discrepancies in E/T Partitioning between the Three TSEB Models
Despite the comparable skills of ET estimations within the three TSEB models, different assumptions and formulations were adopted by the three different models, and discrepancies in E/T partitioning among the three models were observed in the kernel experimental area of MUSOEXE. In this study, both CWSIC and SWDIS were derived to compare the performances of three models regarding the detection of vegetation and soil water stresses. The results indicated that different from TSEB-PT and TSEB-PM, TSEB-TC-TS had the potential to detect vegetation water stress. In addition, the E/T partitioning efficacies of the three TSEB modeling approaches were evaluated Figure 10. Illustration of the temperature decomposition methods adopted in the three TSEB models. T s1 , T s2 , and T s3 denote the soil surface temperatures derived from TSEB-PM, TSEB-PT, and TSEB-T C -T S , respectively, and T v1 , T v2 , and T v3 denote the vegetation canopy temperatures derived from the same three models.

Reliability of the Employed TSEB Models in Estimating Surface Fluxes
The HiWATER-MUSOEXE data have been used extensively to validate the land surface flux models, including one-source and two-source models [70][71][72]. Using the same set of ground-based observations, Ma et al. [70] applied a revised SEBS model to estimate regional heat fluxes in the Heihe River Basin, and their assessment indicates that the RMSEs of the modeled H and LE are 56.9 W/m 2 and 74.8 W/m 2 , respectively. Huang et al. [71] integrated a Normalized Difference Water Index (NDWI) as a water stress index into SEBS through the modification of the parameter kB −1 , and showed RMSEs of 79.8 W/m 2 in H and 84.1 W/m 2 in LE for the revised SEBS. In this study, three two-source models were applied to the middle reach of the Heihe River Basin. The RMSEs of H and LE from TSEB-PT are 47.5 and 75.3 W/m 2 , respectively, the RMSEs of H and LE from TSEB-PM are 44.9 and 70.6 W/m 2 , respectively, and the values from TSEB-T C -T S are 47.9 and 61.8 W/m 2 , respectively. Overall, the performances of two-source modeling approaches were reliable in relation to previously published studies.

Discrepancies in E/T Partitioning between the Three TSEB Models
Despite the comparable skills of ET estimations within the three TSEB models, different assumptions and formulations were adopted by the three different models, and discrepancies in E/T partitioning among the three models were observed in the kernel experimental area of MUSOEXE. In this study, both CWSI C and SWDI S were derived to compare the performances of three models regarding the detection of vegetation and soil water stresses. The results indicated that different from TSEB-PT and TSEB-PM, TSEB-T C -T S had the potential to detect vegetation water stress. In addition, the E/T partitioning efficacies of the three TSEB modeling approaches were evaluated using the measurements from the stable oxygen and hydrogen isotopes system. It was found that all three models tended to slightly underestimate the ratio of T/ET at the Daman site. Aside from the intercomparison between the three two-source models, it was found that TSEB-PM tended to generate a higher LEc estimation than the other two models, especially for the partially vegetated areas. LE S derived from TSEB-PT and TSEB-T C -T S were quite similar. The differences in component temperatures derived from the three models, as well as the aerodynamic resistance, were responsible for this divergence. Besides, the declined LE C /LE ratios observed due to leaf senescence by the end of the growing season was not well characterized by all three models. This indicated the necessity of adding a green vegetation fraction in the future work, as suggested by Kustas et al. [73].

Impact of Temperature Decomposition Accuracies on ET Estimations
To separately estimate E and T, component temperatures are indispensable in the TSEB modeling approach and the temperature decomposition method is the core in E/T partitioning process. In this study, two categories of temperature decomposition methods were intercompared: (1) an iterative procedure based on an energy balance resistance network; and (2) a VFC/LST space method with the assumption that the isolines could be used to decompose composite temperatures. The iterative procedure is commonly adopted in the TSEB model and this approach is verified by different researchers [69,74,75]. The VFC/LST-based approach and related contextual-based method are applied by different researchers, such as Zhang et al. [24], Zhang et al. [55], Merlin et al. [44], Merlin et al. [45], Long and Singh [47], Yang et al. [48], Song et al. [53], and Sun [49].
Further intercomparison indicated that the differences in component temperatures derived from two categories of models were responsible for the discrepancy in ET estimates. As there was no constraint for canopy transpiration to terminate the procedure when LE S exceeded zero in the iterative procedure, this procedure adopted by TSEB-PT and TSEB-PM may have derived a higher LEc and lower Tc compared to TSEB-T C -T S . This is consistent with the findings of Anderson et al. [76], who pointed out the same case under conditions of moderate stomatal closure for TSEB model. In this study, we found the temperature decomposition based on VFC/LST added a further constraint on vegetation transpiration, and this category of method could be a substitute for the iterative method. Besides, it was found that TSEB-T C -T S performed better at detecting the vegetation stress than the other two models. Although a similar temperature decomposition approach was applied to TSEB-PT and TSEB-PM, a significant difference was observed between them. For instance, the R between the two models was only 0.13 for the derived T C . Two reasons may be responsible for this: (1) TSEB-PM applied the Penman-Monteith equation to characterize the canopy transpiration and thus took into account the effect of varying VPD over different underlying surfaces, which was a different case than in TSEB-PT. (2) Resistances were involved in the iterative procedure and the resistance would impact the temperature decomposition. In this study, the theoretically defined VFC/LST trapezoid space was adopted to estimate the component temperatures in TSEB-T C -T S . A small amount of points may have been located outside the trapezoid, and this was mainly caused by ill-parameterized aerodynamic resistance. However, this portion was limited to 1% of the total pixels. Future studies are required to further investigate the influence of resistance on the component temperature decomposition.

Conclusions
In this study, three two-source modeling approaches were evaluated using the ground-based observations from the HiWATER-MUSOEXE campaign. Validated with tower observations, the RMSEs of H and LE were lower than 50 W/m 2 and 76 W/m 2 , respectively, and the results demonstrated that the three TSEB models were capable of predicting reliable surface heat fluxes. The measurements from the stable oxygen and hydrogen isotopes system were used to evaluate the capabilities of three models in E and T partitioning, and it was found that all three models appeared to slightly underestimate the ratio of T/ET.
A further intercomparison of the component temperature decomposition among the three models was conducted to explore the underlying mechanism for the observed differences. Results indicated that the interactive methods applied in TSEB and TSEB-PM may have produced higher LE C and lower T C compared to TSEB-T C -T S due to lack of constraint on vegetation transpiration. Based on the soil moisture isoline in the VFC/LST space, the VFC/LST-based temperature decomposition method added a further constraint on vegetation transpiration, and could be used as a substitution for the interactive procedure adopted in the original TSEB model.