Evapotranspiration Partitioning at Field Scales Using TSEB and Multi-Satellite Data Fusion in The Middle Reaches of Heihe River Basin, Northwest China

Daily evapotranspiration (ET) and its components of evaporation (E) and transpiration (T) at field scale are often required for improving agricultural water management and maintaining ecosystem health, especially in semiarid and arid regions. In this study, multi-year daily ET, E, and T at a spatial resolution of 100 m in the middle reaches of Heihe River Basin were computed based on an ET partitioning method developed by combing remote sensing-based ET model and multi-satellite data fusion methodology. Evaluations using flux tower measurements over irrigated cropland and natural desert sites indicate that this method can provide reliable estimates of surface flux partitioning and daily ET. Modeled daily ET yielded root mean square error (RMSE) values of 0.85 mm for cropland site and 0.84 mm for desert site, respectively. The E and T partitioning capabilities of this proposed method was further assessed by using ratios E/ET and T/ET derived from isotopic technology at the irrigated cropland site. Results show that apart from early in the growing season when the actual E was reduced by plastic film mulching, the modeled E/ET and T/ET agree well with observations in terms of both magnitude and temporal dynamics. The multi-year seasonal patterns of modeled ET, E, and T at field scale from this ET partitioning method shows reasonable seasonal variation and spatial variability, which can be used for monitoring plant water consumption in both agricultural and natural ecosystems.


Introduction
Approximately one-third of Earth's land surface is occupied by semiarid and arid regions, where limited water resources makes it very challenging to provide industrial, agricultural, and ecological water use requirements. In many of these regions, most of the pressures on water resources is driven by evapotranspiration (ET) or water use by irrigated agriculture, which consumes nearly 80% of water resources from open water bodies and aquifers [1,2]. Increasing water consumption from agricultural irrigation results in serious degradation of natural vegetation in many of inland river basins in regions of water scarcity [3]. There is a pressing need to improve agricultural irrigation management to enhance water use efficiency in semiarid and arid regions.
For ET partitioning, transpiration (T) from crops is beneficial water use associated with crop production while evaporation (E) from bare soil is regarded as wasted water from an agronomic point of view [4,5]. Thus, the water efficient agricultural irrigation management strategy is to maximize T and minimize E. To utilize and manage water resources rationally and maintain ecosystem health in semiarid and arid regions, accurate spatiotemporal patterns of ET and its partitioning between E and T at field scale are required to assess and monitor the spatial and temporal distribution of regional water use and stress [6,7].
Over the last few decades, much research has been undertaken to quantify regional ET using remote sensing data over a wide range of temporal and spatial scales and climate conditions [8][9][10]. However, partitioning ET to E and T has mainly been achieved by measurement methods, such as micro-lysimeter, sap-flow, and isotropic analyzer [4,11,12]. Micrometeorological methods using eddy covariance measurements of water vapor and CO 2 , however, have also emerged (e.g., Scanlon and Kustas, [13]) and continue to be refined and improved [14]. Nonetheless, these measurement methods not only require complex instrumentation and data processing, they are also difficult to extrapolate from local patch scale to larger spatial scales [5,15].
A limited number of remote sensing-based models have been developed to estimate ET and its components E and T at large spatial scales by using land surface parameters (e.g., land surface temperature (LST), leaf area index (LAI), and fractional vegetation cover ( f c )) derived from satellite observations. Such remote sensing-based models include the two source energy balance (TSEB) model [16,17], the Priestley-Taylor Jet Propulsion Lab (PT-JPL) model [18], the Global Land Evaporation Amsterdam (GLEAM) model [19], and the Penman-Monteith moderate resolution imaging spectroradiometer (PM-MODIS) [20]. Most of these models have been applied to accurately estimate ET at regional and global scales, but E and T partitioning capabilities of these models are rarely evaluated [15]. The TSEB model originally proposed by Norman et al. [16] and subsequent versions have been applied widely to partitioning ET to E and T under wide range of spatial scales and environmental conditions [5,17,[21][22][23][24].
For agricultural water management and ecosystem health assessment, time-continuous ET and components E and T at field or finer scales (≤ 100 m) are required to provide detailed information about water use [22,25]. Unfortunately, time-continuous values at field scale cannot be directly computed using a single satellite system, because no single satellite system is capable of providing optical and thermal-infrared observations at both high spatial and high temporal resolutions.
Recently, efforts have been made to estimate time-continuous daily ET at field scales based on multi-satellite data fusion combined with remote sensing-based ET models [26,27]. The existing studies for estimating time-continuous daily ET at field scales can be divided into two main categories: (1) ET fusion scheme and (2) input parameter fusion scheme. The first scheme achieves daily ET at field scales via fusing modeled multi-scale ET products using remote sensing data from satellites with high spatial resolution (e.g., Landsat) and satellites with high temporal resolution (e.g., MODIS) [7,22,25,26,[28][29][30][31][32][33]. In the second scheme, daily field-scale ET is directly estimated based on remote sensing-based ET models (e.g., TSEB and surface energy balance system (SEBS)) using input parameters with high spatiotemporal resolution, which is fused from multiple satellite sensors [27,[34][35][36]. Both schemes have been successfully used to estimate daily field-scale ET over rainfed and irrigated agricultural areas [22,25,28,29,33,35,36], forested landscapes [31,32,35,36], and vineyards [7,30]. Ma et al. [36] reported that both schemes perform similarly for homogeneous surfaces while the input parameter fusion scheme may have slightly better performance for heterogeneous surfaces. In these previous studies, the spatial and temporal adaptive reflectance fusion model (STARFM) [37] has been widely applied for producing ET and input parameters of ET models with high spatiotemporal resolution [7,22,25,26,[29][30][31][32][33]35].
Although daily field-scale ET has been achieved over different landscapes and environmental conditions in these previous studies, there are only a few studies estimating daily field-scale E and T based on multi-satellite data fusion. Following the basic scheme of estimates of daily field-scale ET, an ET partitioning method based on TSEB and STARFM is proposed for estimating daily field-scale E and T in this study. The ET partitioning method is applied to estimate daily field-scale ET, E, and T using MODIS and Landsat data collected over an arid region in China, which is mainly covered by irrigated cropland and sparsely vegetated desert. Moreover, the performance of this method is assessed by comparing modeled and measured instantaneous turbulence fluxes and daily ET for irrigated cropland and desert steppe, and by comparing modeled and measured ratio of E and T to ET over irrigated cropland. Furthermore, spatiotemporal water consumption of study area is analyzed using ET, E, and T maps at field scale produced using this ET partitioning method.

Study Area
The evapotranspiration partitioning method was applied in a desert-oasis zone, which is located in the middle reaches of Heihe River Basin in northwest China, as shown in Figure 1. This study area is 30 km × 30 km characterized by sparsely vegetated desert (36%), and an artificial oasis (>60%) dominated by irrigated maize, vegetables, residential areas, orchards, and poplar shelterbelts. This study area has an arid continental monsoon climate characterized by rainy summers and dry winters with over 60% of precipitation occurring in the summer [38]. The annual mean air temperature is approximately 6-8 • C [39] and the annual precipitation is about 100-250 mm, while the annual potential evaporation ranges from 1200 mm to 1800 mm [40]. The water resources for this area are the Heihe River and groundwater [39]. There is considerable amount of water consumed by irrigated cropland surrounded by desert creating an oasis type setting. Given that a fraction of the irrigated water is lost via soil evaporation and therefore not used in crop biomass accumulation, accurately partitioning of evapotranspiration into evaporation and transpiration is crucial to evaluate strategies for reducing E that will lead to water conservation.
In 2012, the Heihe Watershed Allied Telemetry Experimental Research-Multi-Scale Observation Experiment on Evapotranspiration (HiWATER-MUSOEXE) was launched for studying the spatial-temporal variations as well as the impact factors of ET within this desert-oasis ecosystem [41][42][43][44]. In HiWATER-MUSOEXE experiment, the intensive observations of ET and atmospheric and biophysical data were collected over this study area between 3 May and 21 September in 2012 [42], while long term observations of ET in cropland and desert steppe have continued since 2012 [41].

Field Measurements
Two study sites, Daman superstation (DM) and Huazhaizi desert station (HZZ), for obtaining long-term observations in HiWATER-MUSOEXE were selected for applying and evaluating estimates of ET and ET partitioning in this study, as shown in Figure 1.
The Daman superstation (Figure 1c) is located in cropland covered by seeded maize, which is usually planted in late April and harvested in late September. In the early stage of growth, the maize grows quickly and reaches its maximum canopy height (~2 m) and LAI (~4.5) at the end of June [23]. The soil moisture in this cropland is relatively high (ranges from 0.2 to 0.5 m 3 /m 3 ) due to flood irrigation and precipitation during the whole growing season [23]. To reduce the loss of evaporated water from the soil surface for improving water use efficiency, plastic films are used in this cropping system. The plastic films cover approximately 60% of the soil surface. The Huazhaizi desert site (Figure 1d) is located in a natural desert steppe dominated by perennial short shrub-Reaumuria soongaria with constant canopy height of 0.3 m, maximum LAI of~0.7, and soil moisture ranging between 0.07 and 0.2 m 3 /m 3 [45]. At each site, a tower-based observation system is installed and is still in operation, which includes a four-component radiometer, automated weather stations (AWSs), eddy covariance (EC) systems, and soil heat flux plates, and soil moisture and temperature measurement systems (SMTMS) buried around the EC tower. The meteorological observations for model input and surface fluxes for model evaluation span from 1 May to 31 October for the years 2012-2016 were extracted from the HiWATER-MUSOEXE data archive at the Cold and Arid Regions Science Data Center at Lanzhou (http://card.westgis.ac.cn/).
A detailed description of all weather and surface energy balance measurements is provided in Xu et al. [44] and Liu et al. [41]. A brief description of the measurements used in this study is provided here. Meteorological data including air temperature, air pressure, wind speed and direction, and humidity at height of 5 m above ground level (agl) for DM and 3 m agl for HZZ were obtained from AWSs at 10 min averaging intervals. The four components of radiation (incoming and outgoing shortwave and longwave radiation) at height of 12 m agl for DM and 2.65 m agl for HZZ were measured by a 4-way radiometer at 10 min averaging intervals. The surface soil heat flux was calculated using heat flux plate measurements at a depth of 0.06 m combined with soil moisture and temperature at a depth of 0.02 m and 0.04 m measured by SMTMS at 10 min averaging intervals based on the one-dimensional equation of soil thermal diffusion [46]. The surface soil heat flux, radiation observations, and meteorological data were processed to 30 min averages.
The sensible (H) and latent (LE) heat fluxes used to evaluate the accuracy of modeled turbulent fluxes and daily ET were measured from EC systems installed at 4.5 m agl for DM and 2.85 m agl for HZZ. The EC system consists of an open-path infrared gas analyzer and a 3-D sonic anemometer with a sampling frequency of 10 Hz. The high frequency data collected from EC system were processed to produce 30 min average turbulent fluxes of sensible and latent heat and net carbon flux exchange using Edire software (http://www.geos.ed.ac.uk/abs/-research/micromet/EdiRE/). For more details of the EC instrumentation and processing of the 10 Hz data, please refer to Liu et al. [43] and Xu et al. [44].
Generally, the sum of H and LE from EC system is less than the difference between net radiation (R n ) measured from radiometer and the ground soil heat flux (G) calculated from soil heat plates, as known as surface energy imbalance or lack of energy balance closure. The average closure ratios (i.e., (H + LE)/(R n − G)) during study periods for DM and HZZ were 0.80 and 0.79, respectively. To compare with the modeled fluxes, which requires energy conservation, energy balance closure was imposed on the observed fluxes. Given that the surface energy imbalance is mainly associated with underestimation of LE from EC systems [47], the forced closure fluxes were achieved by using the residual-LE closure method (assigning the residual (R n − G − H − LE) to observed LE) [48].
Measurements of δ 18 O of ecosystem water pools and their flux ratios from an isotropic analyzer were obtained to evaluate the accuracy of partitioned E and T at Daman superstation from 27 May to 21 September in 2012 during HiWATER-MUSOEXE experiment [12]. In this study, the ratio of E from soil and T from maize canopy to the total ET from cropland (E/ET and T/ET) were directly computed by Wen et al. [12] using the isotope method, which is based on the isotopic compositions of water fluxes from soil, foliage, and composite land surface. According to Wen et al. [12], the isotopic compositions of ET (δ ET ) was determined by flux-gradient technique, the isotopic composition of E (δ E ) was calculated using Craig-Gordon model, and the composition of T (δ T ) was estimated based on xylem water under isotopic steady state assumption. Finally, E/ET and T/ET were estimated by conservation principles based on δ ET , δ E , and δ T .

Model Inputs
The ET partitioning method based on TSEB and data fusion can be driven by meteorological, radiation, and remote sensing data at coarse and fine spatial resolution. To estimate regional ET, E, and T, the regional meteorological forcing including air temperature (T a ), wind speed (u), relative humidity (RH) and air pressure, and radiation forcing including downward solar radiation (DSR) and longwave radiation (DLR) were determined by interpolating the observations at two study sites based on the inverse distance weighted method.
The remote sensing data used in this study included the MODIS surface reflectance product (MOD09GA) and LST product (MOD11A1) at 500m-1km spatial resolution and level-1 data of Landsat 7 Enhanced Thematic Mapper Plus (ETM+) and Landsat 8 Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS) at 30m-100m spatial resolution. To obtain more continuous MODIS LST, a robust LST reconstruction method proposed by Sun et al. [49] was used to recover clear sky LST for cloud-contaminated pixels in MODIS scenes with <25% cloud cover. Both MODIS reflectance and LST products were projected to UTM projection and resampled using nearest neighbor algorithm to spatial resolution of 100 m, which is consistent with Landsat 8 TIRS images. All the Landsat 7 data used in this study have striped gaps due to the scan-line corrector (SLC) failure of ETM+ since May 2003. Therefore, a simple gap-filled method proposed by Scaramuzza et al. [50] was used to create filled scenes of Landsat 7. To derive Landsat reflectance images from visible near infrared and shortwave bands of Landsat 7/8, the Fast line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) module of ENVI were applied for atmospheric correction. The Landsat LST were derived by using the radiative transfer based method [51] to thermal infrared bands of Landsat 7/8. Landsat reflectance and LST images were resampled to spatial resolution of 100 m. Excluding the cloud-contaminated data, a total of 416 days of MODIS products and 34 days of Landsat products were used in this study from 1 May to 31 October for years 2012-2016, as described in Table 1. The key land surface inputs required for TSEB including LST, f c , LAI, and canopy height (h c ) were derived based on STARFM using MODIS and Landsat data. The fused LST were directly used in TSEB model algorithms. The vegetation inputs f c , LAI, and h c were calculated based on the normalized difference vegetation index (NDVI), which was estimated using fused red and near-infrared land reflectance. The value of f c was estimated based on a NDVI-based formulation (i.e., f c = 1 − NDVI max −NDVI NDVI max −NDVI min p , where NDVI max = 0.8 and NDVI min = 0.1 are the maximum and minimum NDVI values in study area, coefficient p = 0.6 and p = 0.7 were used for vegetated area and desert, respectively) from Gutman and Ignatov [52] and LAI was determined from f c using a logarithmic relationship (i.e., LAI = ) proposed by Choudhury [53]. For vegetated desert and forest, h c was considered as constant (h c = 0.3 m for vegetated desert [45] and h c = 8 m for forest [35]). For the oasis region growing mainly irrigated maize, h c was calculated using a relationship between h c and NDVI (i.e., h c = 2.9184 × NDVI 1.65 ) [35].
To specify surface roughness parameters based on land use types in TSEB, the monthly land cover map at 30 m resolution in study area were collected from dataset of land cover map of Heihe River Basin derived by using multiple classifiers and multisource remote sensing data [54]. The land cover map was aggregated to spatial resolution of 100 m by majority class assignment to match up to the fused reflectance and LST inputs.

ET Partitioning Method
The ET partitioning method based on TSEB and multi-satellite data fusion is illustrated in Figure 2. For clear-sky days, the fused land surface reflectance at visible near infrared and shortwave bands and LST with MODIS temporal resolution and Landsat spatial resolution were derived based on STARFM. Then, the instantaneous land surface fluxes and the ratio of latent flux to insolation for soil ( f SUNs ), canopy ( f SUNc ), and composite land surface ( f SUN ) were estimated based on TSEB by using land surface parameters calculated from fused data combined with regional meteorological inputs interpolated from observations. For cloudy days, the instantaneous f SUNs , f SUNc , and f SUN were obtained by linear interpolation based on an assumption that f SUNs , f SUNc , and f SUN changes linearly between the two adjacent clear days. Finally, the time continuous regional daily E, T, and ET were calculated based on the assumption of f SUNs , f SUNc , and f SUN is constant during the daylight hours [3].

TSEB Model
To estimate surface energy fluxes from soil and canopy components, Norman et al. [16] originally proposed the TSEB model, which has undergone several revisions [17,[55][56][57][58][59] to accommodate more diverse land cover conditions. In TSEB, the surface energy balance for soil, canopy and the composite land surface are expressed as where is net radiation, is latent heat flux, is sensible heat flux, is ground soil heat flux (all in W/m 2 ). Subscript and represent soil and canopy components, respectively. Apart from , all the land surface fluxes in Equation (3) can be calculated as the sum of soil and canopy components in Equation (1) and (2).
is usually calculated as the fraction of where is an empirical coefficients, which varies with time, soil type, and soil moisture. For satellite overpass time around solar noon, is assumed to be constant value of 0.35 for many soil type and moisture conditions [60].
The net radiation for soil ( ) and canopy ( ) are given as [17], where ↓ and ↓ are incoming longwave and shortwave radiation (W/m 2 ), respectively. and are the transmittances of longwave and shortwave radiation through the canopy,

TSEB Model
To estimate surface energy fluxes from soil and canopy components, Norman et al. [16] originally proposed the TSEB model, which has undergone several revisions [17,[55][56][57][58][59] to accommodate more diverse land cover conditions. In TSEB, the surface energy balance for soil, canopy and the composite land surface are expressed as where R n is net radiation, LE is latent heat flux, H is sensible heat flux, G is ground soil heat flux (all in W/m 2 ). Subscript s and c represent soil and canopy components, respectively. Apart from G, all the land surface fluxes in Equation (3) can be calculated as the sum of soil and canopy components in Equation (1) and (2). G is usually calculated as the fraction of R ns [16], where c G is an empirical coefficients, which varies with time, soil type, and soil moisture. For satellite overpass time around solar noon, c G is assumed to be constant value of 0.35 for many soil type and moisture conditions [60]. The net radiation for soil (R ns ) and canopy (R nc ) are given as [17], where L ↓ and S ↓ are incoming longwave and shortwave radiation (W/m 2 ), respectively. τ longwave and τ solar are the transmittances of longwave and shortwave radiation through the canopy, respectively. ε and α are the emissivity and albedo, respectively. Subscript s and c represent soil and canopy components, respectively. T s and T c are soil and canopy components temperature (K), respectively. By permitting the interaction between soil and canopy fluxes, Norman et al. [16] proposed a "series" version of TSEB, which is more robust to compute the sensible heat fluxes [61,62]. Following Norman et al. [16], H, H s , and H c are expressed as where ρ is the density of air (kg/m 3 ), C P is the specific heat of air (J/(kg·K)), T ac is the air temperature in the canopy air layer (K), T a is the air temperature in the soil-canopy combined surface (K). r a is the bulk aerodynamic resistance (s/m), which is computed by using Monin-Obukhov surface layer similarity theory. The soil resistance r s (s/m) is derived according to Kustas and Norman [17]. r x is the boundary layer resistance of the canopy elements (s/m).
In TSEB, T s and T c in Equations (5), (6), (8), and (9) are related to satellite derived directional surface temperature (T R (θ)) and the fraction vegetation cover ( f c (θ)), as expressed in [63], where θ is the view angle of thermal sensor. T s and T c can be calculated by iteration for Equations (2)- (10) with an initial value of LE c , which is estimated based on Priestley-Taylor equation, where f g is the fraction of LAI that is green, ∆ is the slope of the saturation vapor pressure versus temperature curve (kPa/ • C), γ is the psychrometer constant (~0.066 kPa/ • C). The α PT term is the Priestley-Taylor parameter, which is assumed to be constant value of 1.26 for unstressed transpiring canopy. For stressed conditions which can cause LE s < 0, TSEB iteratively reduces the value of α PT from 1.26 until LE s > 0 [63]. In addition, for the sparsely vegetated desert, the modification for initial value of LE c proposed by Zhuang and Wu [59] was adopted since it reduced bias between measured and modeled H. This modification was not as reliable in partitioning ET between E and T as the original TSEB for the cropland and therefore was not applied to this land cover type. The instantaneous surface energy fluxes for soil, canopy, and composite land surface at MODIS overpass time are solved for iteratively with Equations (1)- (11) and are then upscaled to daily ET, E, and T based on the ratio of instantaneous to daily insolation approach [3,22,[29][30][31], where R si are instantaneous insolation. Subscript i represents soil, canopy, or composite land surface. LE id and R sid are daily latent heat flux and insolation for soil, canopy, and composite land surface, respectively. Then LE, LE s , and LE c are converted to daily ET, E, and T by dividing the latent heat of evaporation.

STARFM
The STARFM, originally proposed by Gao et al. [37] for fusing MODIS and Landsat surface reflectance, has been widely used to obtained high spatiotemporal information of LST [35,36], land use types [64], and ET [22,26,[29][30][31][32][33]. In this study, STARFM was used to produce LST and land surface reflectance at spatial resolution of 100 m with temporal resolution of MODIS (hereinafter referred to as Landsat-like). When Landsat data are not available, the predicted Landsat-like data is determined by using existing Landsat/MODIS data pair/pairs based on a weighting function, which is developed using a moving search window method. In STARFM, the weighting function for central pixel of moving window can be expressed as where L and M denote the land surface reflectance/LST of Landsat and MODIS, respectively. w is the size of searching window. (x w/2 , y w/2 ) and (x i , y i ) are the location of central pixel and the neighboring pixel in searching window, respectively. The time variable t 0 and t k are the predication date and acquisition date of Landsat/MODIS data pair/pairs, respectively. k is the number of Landsat/MODIS data pair/pairs. W ijk is the weight that represents the contribution of each neighboring pixel to the predicted central pixel. The weight W ijk is determined by spectral difference between MODIS and Landsat pair at t k , temporal difference between MODIS data at t 0 and t k , and the distance between central pixel and neighboring pixel. For more details of STARFM, please refer to Gao et al. [37].

Landsat-Only Interpolation
To evaluate the value added by MODIS in the fusion process, the daily ET estimated based on fused data and the daily ET estimated based on Landsat data using a simple Landsat-only interpolation scheme were compared with the observed daily ET from the EC flux, radiation and meteorological observations. The Landsat-only interpolation scheme has been widely used to obtain time continuous daily ET in Mapping Evapotranspiration with Internalized Calibration (METRIC) model based on the ratio of ET to reference ET ( f RET ) [65], where ET re f is the reference ET estimated based on the Penman-Monteith formulation of the Food and Agriculture Organization standard using meteorological data [66]. In this scheme, the time continuous f RET is determined by linear interpolation based on estimated f RET on Landsat overpass dates. Then, the time continuous daily ET is derived as the product of f RET and ET re f .

Evaluation of Land Surface Fluxes at The Study Sites
To assess the performance of the ET estimating method, the modeled and measured land surface fluxes including R n , G, H, and LE for the satellite overpass times were compared at the two study sites. For matching with the source areas of observations, the modeled land surface fluxes based on TSEB and fusion data were averaged over upwind source area (two pixels) at the two study sites [35,67]. The scatter plots in Figure 3 show the comparisons of modeled and measured land surface fluxes at DM and HZZ sites during 1 May to 31 October from 2012 to 2016. In addition, the statistical difference between measured and modeled land surface fluxes are given in Table 2, including the bias error computed as the average difference between measured and modeled fluxes, the root mean square error (RMSE), and the mean absolute percent difference (MAPD).  From Figure 3 and Table 2, it appears that modeled , , , and agree well with the forced closure observations, with the scatter of points falling generally along the 1:1 line. RMSE values of ~25 W/m 2 for , ~40 W/m 2 for , ~55 W/m 2 for , and ~60 W/m 2 for , indicate reasonable surface flux partitioning at both DM and HZZ sites. Compared with , , and , modeled has the best agreement with measured yielding a MAPD of 3% and 4% for DM site and HZZ site, respectively. The difference between modeled and measured is mainly caused by the considerable mismatch in source areas between soil heat flux plates (0.1-1 m) and the resolution of modeled fluxes (200 m) [68]. Compared with HZZ site covered by sparse vegetation, the modeled at DM site covered by irrigated maize with wider range of soil moisture and vegetation cover during the whole growing season, yields greater scatter with RMSE of 41 W/m 2 and MAPD of 47%, but shows a smaller bias of 12 W/m 2 .
For , the difference between modeled and observations from EC system illustrated in Figure 3c and 3g indicates a slightly underestimate with bias of −20 W/m 2 for DM site and −24 W/m 2 for HZZ site, which is consistent with the results from earlier studies [45,68]. For the residual of the surface energy balance , the values estimated using TSEB and fusion data yield good agreement with Bias of ~10 W/m 2 for both DM and HZZ sites when compared with flux tower observations using the residual-LE closure method. The MAPD value for at DM site (12%) is less than the uncertainty of the EC measured (18%) reported by Wang et al. [69]. However, the MAPD value is as larger as 32% at HZZ site due in part to the fact that the overall range in values are relatively small under sparse arid vegetation cover and water-limited conditions.  From Figure 3 and Table 2, it appears that modeled R n , G, H, and LE agree well with the forced closure observations, with the scatter of points falling generally along the 1:1 line. RMSE values of 25 W/m 2 for R n ,~40 W/m 2 for G,~55 W/m 2 for H, and~60 W/m 2 for LE, indicate reasonable surface flux partitioning at both DM and HZZ sites. Compared with G, H, and LE, modeled R n has the best agreement with measured yielding a MAPD of 3% and 4% for DM site and HZZ site, respectively. The difference between modeled and measured G is mainly caused by the considerable mismatch in source areas between soil heat flux plates (0.1-1 m) and the resolution of modeled fluxes (200 m) [68]. Compared with HZZ site covered by sparse vegetation, the modeled G at DM site covered by irrigated maize with wider range of soil moisture and vegetation cover during the whole growing season, yields greater scatter with RMSE of 41 W/m 2 and MAPD of 47%, but shows a smaller bias of 12 W/m 2 .
For H, the difference between modeled H and observations from EC system illustrated in Figure 3c,g indicates a slightly underestimate with bias of −20 W/m 2 for DM site and −24 W/m 2 for HZZ site, which is consistent with the results from earlier studies [45,68]. For the residual of the surface energy balance LE, the values estimated using TSEB and fusion data yield good agreement with Bias of~10 W/m 2 for both DM and HZZ sites when compared with flux tower observations using the residual-LE closure method. The MAPD value for LE at DM site (12%) is less than the uncertainty of the EC measured LE (18%) reported by Wang et al. [69]. However, the MAPD value is as larger as 32% at HZZ site due in part to the fact that the overall range in LE values are relatively small under sparse arid vegetation cover and water-limited conditions.

Evaluation of Daily ET at The Study Sites
To further evaluate the accuracy of the proposed ET estimating method as well as to assess the value added by MODIS in the fusion process, the modeled daily ET based on the proposed method and Landsat-only interpolation scheme were compared with the residual LE observations at the two flux tower sites (Figure 4). The error statistics for modeled daily ET based on both two methods are listed in Table 3. Additionally, the time series of observed daily ET and modeled ET based on both two methods for the period 1 May to 31 October from 2012 to 2016 are shown in Figure 5.
Sens. 2020, 12, x FOR PEER REVIEW 1 valuation of Daily ET at The Study Sites o further evaluate the accuracy of the proposed ET estimating method as well as to asses added by MODIS in the fusion process, the modeled daily ET based on the proposed me andsat-only interpolation scheme were compared with the residual LE observations at the ower sites ( Figure 4). The error statistics for modeled daily ET based on both two method in Table 3. Additionally, the time series of observed daily ET and modeled ET based on ethods for the period 1 May to 31 October from 2012 to 2016 are shown in Figure 5.      For the DM site, the modeled daily ET using both methods agree well with observations, with similar scatter falling along the 1:1 line (Figure 4). The modeled daily ET using Landsat-only interpolation method show acceptable accuracy with bias, RMSE, and MAPD values of 0.04 mm, 0.91 mm, and 17%, respectively. By contrast, the modeled daily ET using TSEB and fusion data have slightly lower bias, RMSE and MAPD values of 0.01 mm, 0.85 mm, and 16%, indicating slightly better performance using this daily ET interpolation method. Given that there is sufficient irrigation and precipitation for the irrigated maize field to meet atmospheric demand at the DM site, the variation of ET is mainly affected by radiation. From Figure 5, we can see that the similar ET trends with low values correspond to decreased insolation on cloudy days were produced by the two daily ET interpolation methods during the whole study period. However, neither model captures the observed ET values falling below 3 mm. This is likely due to neither interpolation method providing accurate estimates under these low radiation conditions. Given that the lack of available Landsat data during DOY 125-170 in 2012 and 2014, and DOY 140-240 in 2015, the daily ET modeled from Landsat-only interpolation scheme differ greatly from both daily ET observed from EC and modeled from TSEB using the fusion data. This is an issue using Landsat-only interpolation, especially in regions with persistent cloud cover.
For HZZ site, better agreement between modeled and observed daily ET is achieved by using TSEB and data fusion. The modeled values of daily ET from the Landsat-only interpolation scheme are clearly underestimated when the observed values of daily ET are greater than 3 mm (Figure 4d). Daily ET modeled based on TSEB using fusion data yields bias, RMSE, and MAPD values of −0.25 mm, 0.84 mm, and 29%, respectively. By comparison, the underestimated daily ET based on Landsat-only interpolation scheme yields bias of −0.45 mm, RMSE of 1.00 mm, and MAPD of 33%. The HZZ site is located in a natural desert steppe with sparsely clumped vegetation and low soil moisture. As a result, the significant fluctuation of time series of daily ET is mainly caused by rainfall events. From Figure 5, it appears that the trend of daily ET modeled based on Landsat-only interpolation scheme is hard to capture the oscillations in ET caused by rainfall events, which significantly impacts soil evaporation and the subsequent rapid drying that again strongly impacts E. Since there is MODIS imagery often available soon after rainfall events, the more accurate temporal trend of daily ET is captured with relatively large ET increase after a rainfall event followed by drying of the soil surface, hence decreasing E and ET. Therefore, using TSEB and fusing MODIS and Landsat significantly improves daily ET over rainfed ecosystems, particularly with sparse vegetation cover where rapid changes in E significantly impact ET.

Evaluation of ET Partitioning Over The Cropland
The assessment of the ET partitioning method based on TSEB and fusion data was conducted by comparing with the ratios of E/ET and T/ET obtained from isotopic technology at DM site from 27 May (DOY 148) to 21 September (DOY 265) in 2012 [12], as illustrated in Figure 6 and Table 4. For the period when the plastic film has deteriorated (DOY >170), the proposed method performs well for partitioning ET to E and T at the cropland site, with both values scattering around 1:1 line. However, outliers in E/ET and T/ET are under the early growing season when there is a plastic film significantly affecting the partitioning. For the whole growing season, this partitioning method based on TSEB and fusion data yields a bias of 0.1, a RMSE of 0.2, and a MAPD of 107% for E/ET while it produces a bias of −0.1, a RMSE of 0.2, and a MAPD of 17% for T/ET. The large MAPD value for E/ET is due to two factors: (1) the impact of the plastic film on E and T partitioning early in the growing season and (2) with most E/ET values below 0.3, any slight deviation in E/ET estimate significantly increases the value of the MAPD statistic. The time series of observed and modeled E/ET and T/ET are illustrated in Figure 7. From Figure  7, it is confirmed that the large discrepancy between measured and modeled E/ET and T/ET for the group of points illustrated in Figure 6 occurs in the beginning of the growing season. As a previous study [23] concluded, this is due to the fact that the actual E/ET is reduced by the plastic film barrier  The time series of observed and modeled E/ET and T/ET are illustrated in Figure 7. From Figure 7, it is confirmed that the large discrepancy between measured and modeled E/ET and T/ET for the group of points illustrated in Figure 6 occurs in the beginning of the growing season. As a previous study [23] concluded, this is due to the fact that the actual E/ET is reduced by the plastic film barrier deployed at the start of the growing season. This thin and disposal plastic film, covering approximately 60% of the soil surface, degrades from radiation and farming operations during the early growth stage of the maize. As a result, the plastic film is not an effective vapor barrier affecting observed E/ET and T/ET by the middle and later part of the growing season. The modeled E/ET and T/ET using TSEB and fusion data agree well with the observations in terms of both magnitude and temporal dynamics by the middle of the growing season. This suggests that the ET partitioning method predicted by TSEB can partition ET into E and T accurately over irrigated cropland as long as there are no artificial barriers reducing soil E. For the middle and late part of measurement period (DOY >170), the modeled E/ET from TSEB and fusion data yields bias, RMSE, and MAPD values of 0, 0.11, and 53%, respectively, and the modeled T/ET results in bias, RMSE, and MAPD values of 0, 0.11, and 10%, respectively. The time series of observed and modeled E/ET and T/ET are illustrated in Figure 7. From Figure  7, it is confirmed that the large discrepancy between measured and modeled E/ET and T/ET for the group of points illustrated in Figure 6 occurs in the beginning of the growing season. As a previous study [23] concluded, this is due to the fact that the actual E/ET is reduced by the plastic film barrier deployed at the start of the growing season. This thin and disposal plastic film, covering approximately 60% of the soil surface, degrades from radiation and farming operations during the early growth stage of the maize. As a result, the plastic film is not an effective vapor barrier affecting observed E/ET and T/ET by the middle and later part of the growing season. The modeled E/ET and T/ET using TSEB and fusion data agree well with the observations in terms of both magnitude and temporal dynamics by the middle of the growing season. This suggests that the ET partitioning method predicted by TSEB can partition ET into E and T accurately over irrigated cropland as long as there are no artificial barriers reducing soil E. For the middle and late part of measurement period (DOY >170), the modeled E/ET from TSEB and fusion data yields bias, RMSE, and MAPD values of 0, 0.11, and 53%, respectively, and the modeled T/ET results in bias, RMSE, and MAPD values of 0, 0.11, and 10%, respectively.

Spatiotemporal Patterns of E, T, and ET
Overall model validation in Section 3 indicates that the proposed ET partitioning method produced satisfactory estimates of ET and its partitioning, with statistical errors within the expected and reported errors in literature. The daily ET modeled based on TSEB using fusion data at study sites is comparable or had slightly better accuracy than previous studies [3,35,36,67,68]. The magnitude and temporal dynamics of modeled E/ET and T/ET in this study are consistent with results

Spatiotemporal Patterns of E, T, and ET
Overall model validation in Section 3 indicates that the proposed ET partitioning method produced satisfactory estimates of ET and its partitioning, with statistical errors within the expected and reported errors in literature. The daily ET modeled based on TSEB using fusion data at study sites is comparable or had slightly better accuracy than previous studies [3,35,36,67,68]. The magnitude and temporal dynamics of modeled E/ET and T/ET in this study are consistent with results from model output based on the dual temperature difference (DTD) approach using MODIS data [3] and output from the two-source variational data assimilation (TVDA) system using in situ observations [23]. This indicates the current modeling approach incorporating data fusion can provide reliable E, T, and ET maps using Landsat and MODIS data.
The monthly E, T, and ET maps derived from modeled daily values from May to October from 2012 to 2016 in the study domain were produced for evaluating the E versus T as a function of land cover. From Figure 8, it appears that the seasonal variation and spatial diversity in monthly E fluctuates slightly from month to month and is affected by different land cover conditions. Given that the important factor influencing the magnitude of E-soil moisture, which is almost exclusively affected by the frequency of precipitation in rainfed environments [70], the monthly precipitation derived from averaged observation of two study sites is listed in Table 5. For sparsely vegetated desert site, the seasonal fluctuation of monthly E coincides with that of precipitation illustrated in Table 5, due to the fact that the rate of E is correlated positively with soil moisture under water-limited condition [71,72]. The value of E for the desert site is relatively high in the middle of the summer season (i.e., from June to August) since the highest amount of precipitation normally occurs during this period, with a monthly mean precipitation of~30 mm. On the other hand, precipitation is relatively low prior to June and after August with monthly mean precipitation of <15 mm.
from model output based on the dual temperature difference (DTD) approach using MODIS data [3] and output from the two-source variational data assimilation (TVDA) system using in situ observations [23]. This indicates the current modeling approach incorporating data fusion can provide reliable E, T, and ET maps using Landsat and MODIS data.
The monthly E, T, and ET maps derived from modeled daily values from May to October from 2012 to 2016 in the study domain were produced for evaluating the E versus T as a function of land cover. From Figure 8, it appears that the seasonal variation and spatial diversity in monthly E fluctuates slightly from month to month and is affected by different land cover conditions. Given that the important factor influencing the magnitude of E-soil moisture, which is almost exclusively affected by the frequency of precipitation in rainfed environments [70], the monthly precipitation derived from averaged observation of two study sites is listed in Table 5. For sparsely vegetated desert site, the seasonal fluctuation of monthly E coincides with that of precipitation illustrated in Table 5, due to the fact that the rate of E is correlated positively with soil moisture under waterlimited condition [71,72]. The value of E for the desert site is relatively high in the middle of the summer season (i.e., from June to August) since the highest amount of precipitation normally occurs during this period, with a monthly mean precipitation of ~30 mm. On the other hand, precipitation is relatively low prior to June and after August with monthly mean precipitation of <15 mm.   The monthly E for the irrigated agricultural area (oasis) shows complex fluctuations during growing season that is influenced by solar irradiance, soil moisture, as well as vegetation cover. In the early part of growing season, the soil surface is partially covered by vegetation and soil moisture is relatively high due to the flood irrigation implemented once prior to the growing season and twice in the early growing season [23,73]. As a result, the rate of E in May and June is relatively high and fluctuates from year to year with changes in precipitation and timing of the irrigation. In the middle part of the growing season, the vegetation cover is relatively high, which rapidly reduces the amount of the solar irradiation and wind that reaches the soil surface both causing a significant reduction in E even if there are high soil moisture conditions. This results in the monthly E in July and August being smaller than that in May and June, even though the soil moisture is still relatively high due to frequent precipitation events and additional irrigation events. Compared to August, the rate of E in September is slightly increased mainly caused by the leaf senescence and harvesting crop that leads to the decrease of interception of solar irradiation and attenuation of wind. Lowest solar irradiation as well as relatively low soil moisture conditions due to low precipitation leads to the monthly E in October being the smallest during the study period.
The monthly T maps are illustrated in Figure 9. Considering that the seasonal variation trend of rate of T generally follows the temporal variability in LAI [3,5,72,73], the spatial pattern of monthly mean LAI derived from fusion data during the study period is shown in Figure 10. From Figures 9 and 10, we can see that the seasonal variation, spatial variability, and differences in monthly T in different years is almost identical with the discrepancies in monthly mean LAI between different years. The spatio-temporal pattern for T distinguishes the oasis containing irrigated crops from the surrounding desert. During the whole growing season, the monthly T can be regarded as constant at about 20 mm for the desert vegetation and residential areas, while the T maps show seasonal variation for the oasis area where the value of T is driven by vegetation growth and biomass accumulation. The monthly T for oasis area increases from May to June with vegetation growth and reaches maximum in July with highest biomass accumulation of vegetation and then decreases from August to October associated with the leaf senescence.
As shown in Figure 11, the monthly ET, the sum of E and T, has significant seasonal variation over the six month period and spatial patterns for the different land cover types in the study domain. The ET values for cropland is clearly greater than surrounding desert during the whole growing season. Although the value of ET is dominated by rate of E for sparsely vegetated desert while dominated by rate of T for the irrigated croplands/oasis, the seasonal variation trend of ET in the oasis is similar to the surrounding desert. Comparison of ET maps in different years (cf. Figure 11), the trends in the seasonal variation in ET remain relatively consistent from 2012 to 2016. Naturally, the monthly ET in the middle growing season is larger than that in early and late growing season due to maximum plant biomass and precipitation/irrigation inputs. However, differences in amount and timing of the precipitation and irrigation events, plant growth rates and peak biomass differ from year to year. As a result, there are some apparent discrepancies in monthly cumulative ET for different years.  As shown in Figure 11, the monthly ET, the sum of E and T, has significant seasonal variation over the six month period and spatial patterns for the different land cover types in the study domain. The ET values for cropland is clearly greater than surrounding desert during the whole growing season. Although the value of ET is dominated by rate of E for sparsely vegetated desert while dominated by rate of T for the irrigated croplands/oasis, the seasonal variation trend of ET in the oasis is similar to the surrounding desert. Comparison of ET maps in different years (cf. Figure 11), the trends in the seasonal variation in ET remain relatively consistent from 2012 to 2016. Naturally, the monthly ET in the middle growing season is larger than that in early and late growing season due to maximum plant biomass and precipitation/irrigation inputs. However, differences in amount and timing of the precipitation and irrigation events, plant growth rates and peak biomass differ from year to year. As a result, there are some apparent discrepancies in monthly cumulative ET for different years.

Water Consumption by Land Type
To characterize water consumption by different land covers in the study domain, the monthly value of E and T per unit area for cropland, woodland, wetland, water, and desert in 2013 are demonstrated in Figure 12. The densely vegetated area including cropland, woodland, and wetland significantly consume more water than the sparsely vegetated desert. The monthly E and T shows that there are differences in water loss from soil surface and water consumption for vegetation growth of densely vegetated land types caused by different soil moisture and phenology, growth, and development characteristics of the vegetation. The woodland consumes greatest amount of water for vegetation growth in early and late of growing season, while the cropland yielded the most plant water use in the middle of growing season.

Water Consumption by Land Type
To characterize water consumption by different land covers in the study domain, the monthly value of E and T per unit area for cropland, woodland, wetland, water, and desert in 2013 are demonstrated in Figure 12. The densely vegetated area including cropland, woodland, and wetland significantly consume more water than the sparsely vegetated desert. The monthly E and T shows that there are differences in water loss from soil surface and water consumption for vegetation growth of densely vegetated land types caused by different soil moisture and phenology, growth, and development characteristics of the vegetation. The woodland consumes greatest amount of water for vegetation growth in early and late of growing season, while the cropland yielded the most plant water use in the middle of growing season.
The cropland, covering 53% of study area, consumes the greatest amount of surface water and groundwater for the study domain mainly from irrigation. To distinguish water consumption of the cropland for the years 2012 to 2016, the monthly E, T, and ET per unit area are illustrated in Figure 13. In general, the total water consumption of the cropland area remains relatively constant from 2012 to 2016 with cumulative ET~675 mm during May to October but the water loss through evaporation from soil surface increases gradually, especially in July and August. The ratio of cumulative E to cumulative The cropland, covering 53% of study area, consumes the greatest amount of surface water and groundwater for the study domain mainly from irrigation. To distinguish water consumption of the cropland for the years 2012 to 2016, the monthly E, T, and ET per unit area are illustrated in Figure  13. In general, the total water consumption of the cropland area remains relatively constant from 2012 to 2016 with cumulative ET ~675 mm during May to October but the water loss through evaporation from soil surface increases gradually, especially in July and August. The ratio of cumulative E to cumulative ET from May through October yields values of 0. 31

Conclusions
To provide detailed information of water use for agricultural water management at field and farm to landscape scales in semiarid and arid regions, a daily field-scale ET partitioning method developed by combining STARFM and TSEB was utilized and verified in the middle reaches of Heihe River Basin in this study. This method was developed by following the input parameter fusion scheme for estimating daily field-scale ET, which estimates ET, E, and T based on TSEB model formulations using land surface input parameters derived from data with spatiotemporal characteristics of MODIS and Landsat fused using STARFM.  The cropland, covering 53% of study area, consumes the greatest amount of surface water and groundwater for the study domain mainly from irrigation. To distinguish water consumption of the cropland for the years 2012 to 2016, the monthly E, T, and ET per unit area are illustrated in Figure  13. In general, the total water consumption of the cropland area remains relatively constant from 2012 to 2016 with cumulative ET ~675 mm during May to October but the water loss through evaporation from soil surface increases gradually, especially in July and August. The ratio of cumulative E to cumulative ET from May through October yields values of 0.31 for 2012, 0.41 for 2013, 0.47 for 2014, 0.46 for 2015, and 0.53 for 2016, respectively. With this increasing trend in E/ET, better irrigation strategies and water management decisions are urgently needed to reduce the water loss from the soil in the irrigated cropland area.

Conclusions
To provide detailed information of water use for agricultural water management at field and farm to landscape scales in semiarid and arid regions, a daily field-scale ET partitioning method developed by combining STARFM and TSEB was utilized and verified in the middle reaches of Heihe River Basin in this study. This method was developed by following the input parameter fusion scheme for estimating daily field-scale ET, which estimates ET, E, and T based on TSEB model formulations using land surface input parameters derived from data with spatiotemporal characteristics of MODIS and Landsat fused using STARFM.

Conclusions
To provide detailed information of water use for agricultural water management at field and farm to landscape scales in semiarid and arid regions, a daily field-scale ET partitioning method developed by combining STARFM and TSEB was utilized and verified in the middle reaches of Heihe River Basin in this study. This method was developed by following the input parameter fusion scheme for estimating daily field-scale ET, which estimates ET, E, and T based on TSEB model formulations using land surface input parameters derived from data with spatiotemporal characteristics of MODIS and Landsat fused using STARFM.
In this study, daily ET, E, and T at spatial resolution of 100 m from May through October for the years 2012 to 2016 for the study area were estimated by using the ET partitioning method. A comparison of modeled instantaneous land surface fluxes and residual closure observations derived from the EC systems and AWSs at cropland and desert sites indicate that reliable surface fluxes were estimated using this method. In addition, the value added by incorporating MODIS in the fusion process for estimating daily ET, compared to Landsat-only interpolation scheme was assessed by comparing with ET observations for the irrigated cropland and desert sites. The results demonstrate that more accurate daily ET is produced by using the data fusion-based TSEB, with RMSE values of less than 1 mm/day for cropland and desert. We further tested the E and T partitioning capabilities of this proposed method for the cropland using E/ET and T/ET observed from isotopic technology. Apart from early growing season when the actual E was reduced by plastic film covering, the modeled E/ET and T/ET agreed well with observations in terms of both magnitude and temporal dynamics.
Maps of monthly E, T, and ET from May through October for the years 2012 to 2016 in the study area show reasonable temporal and spatial variability. The monthly E varies from month to month due to compounding effects from different land cover types being affected by precipitation, irrigation events and rate and amount of accumulated vegetation cover/biomass over the growing season. Meanwhile monthly T shows a consistent seasonal trend associated with vegetation growth and biomass production. Moreover, water consumption by different land cover types was conducted based on spatiotemporal distributions of E and T, reflecting the influence of soil moisture, fractional vegetation cover and phenology on E and T partitioning. The general relative increase in water loss through evaporation from bare soil in irrigated cropland during the growing season from 2012 to 2016 indicates that irrigation and water management strategies are urgently needed for conserving water resources in the Heihe River Basin.
As far as modeling the effect of plastic film layers for reducing E, previous studies indicate that the accuracy in estimating ET and its partitioning using Priestly-Taylor [74] and Shuttleworth-Wallace [75] modeling approaches with plastic film-mulching can be improved by incorporating the effect of mulching fraction on E. Based on these previous studies, we plan to develop and evaluate a refinement in the TSEB algorithm for computing E for fields containing plastic films. Due to the limitation of observations of ET partitioning into E and T, model validation could only be conducted for the oasis-desert ecosystem. Future plans include verifying the utility of this ET partitioning method over a wider range of land cover types using eddy covariance-based techniques [13,14].