A coupled remote sensing and the Surface Energy Balance with Topography Algorithm (SEBTA) to estimate actual evapotranspiration over heterogeneous terrain

Evapotranspiration (ET) may be used as an ecological indicator to address the ecosystem complexity. The accurate measurement of ET is of great significance for studying environmental sustainability, global climate changes, and biodiversity. Remote sensing technologies are capable of monitoring both energy and water fluxes on the surface of the Earth. With this advancement, existing models, such as SEBAL, SEBI and SEBS, enable us to estimate the regional ET with limited temporal and spatial coverage in the study areas. This paper extends the existing modeling efforts with the inclusion of new components for ET estimation at different temporal and spatial scales under heterogeneous terrain with varying elevations, slopes and aspects. Following a coupled remote sensing and surface energy balance approach, this study emphasizes the structure and function of the Surface Energy Balance with Topography Algorithm (SEBTA). With the aid of the elevation and landscape information, such as slope and aspect parameters derived from the digital elevation model (DEM), and the vegetation cover derived from satellite images, the SEBTA can account for the dynamic impacts of heterogeneous terrain and changing land cover with some varying kinetic parameters (i.e., roughness and zero-plane displacement). Besides, the dry and wet pixels can be recognized automatically and dynamically in image processing thereby making the SEBTA more sensitive to derive the sensible heat flux for ET estimation. To prove the application potential, the SEBTA was carried out Correspondence to: C. S. Liu (csliu@re.ecnu.edu.cn) to present the robust estimates of 24 h solar radiation over time, which leads to the smooth simulation of the ET over seasons in northern China where the regional climate and vegetation cover in different seasons compound the ET calculations. The SEBTA was validated by the measured data at the ground level. During validation, it shows that the consistency index reached 0.92 and the correlation coefficient was 0.87.


Introduction
Stream flow, plant transpiration and photosynthesis, evaporation, soil moisture changes, and groundwater recharge are intimately related with each other to form water balance dynamics on the surface of the Earth. Plants transpire water through stomata to maintain the favorable temperatures necessary for metabolic processes while simultaneously absorbing carbon dioxide (CO 2 ) from atmosphere. Actual evapotranspiration (ET) is comprised of not only evaporation during photosynthesis, but also ET from the root zone storage, ET from the percolation to water table, and ET supplied by capillary rise from the water table. As ET is affected by both water and energy balances in the soil-plant-atmosphere continuum (SPAC), it involves many complex processes in the nexus of water and thermal cycles at the surface of the Earth. The surface soil moisture driven by precipitation and irrigation actually sustains the ET phenomenon that mainly occurs from the tree crowns with a fully closed canopy. In addition to the water content in the soil, ET is also driven by the water vapor gradient in the air in relation to the net radiation at the Published by Copernicus Publications on behalf of the European Geosciences Union. ground level. The partitioning of net radiation between sensible and latent heat fluxes is in turn dependent on the amount of available soil moisture on the surface of the Earth.
Large-scale ET estimation is of great concern in numerous practices from regional water resources management to local irrigation scheduling (Bastiaanssen et al., 1998a, b;Kite and Droogers, 2000;Schuurmans et al., 2003). Basically, point measurement of ET can be made possibly using a few ground-based methods, such as: (i) crop coefficient and climatic parameters, (ii) soil moisture monitoring and vapor flux measurement, and (iii) energy balance and water vapor transfer using the eddy-covariance method. These insitu measurements can only provide point measurements of ET and they do not account for spatial variability of ET in a large scale. Satellite remote sensing technologies therefore offer comparative advantages to monitoring water/heat fluxes at the regional scale (Anderson et al., 2005;Boegha et al., 2002;Chen et al., 2005;Ottle et al., 1989). Recent advances in satellite remote sensing and ET estimation algorithms have collectively enabled us to develop regional maps of ET to present the spatiotemporal variations within varying weather conditions. Two of the popular approaches being used for ET estimation are the reflectance-based crop coefficient method (e.g., potential evapotranspiration -PET) and the energy balance method (e.g., daily actual evapotranspiration -DAET) (Kogan, 1995;Jiang and Islam, 1999;Mecikalski et al., 1999;Norman et al., 1995Norman et al., , 2000Norman et al., , 2003Olioso et al., 2005;Jia et al., 2003).
In early stage, Seguin et al. (1983) and  calculated difference between surface radiative temperature and air temperature to estimate the large-scale regional ET with the aid of satellite thermal infrared information. Menenti et al. (1993Menenti et al. ( , 2001) also presented calculation of the surface ET based on the surface energy balance index (SEBI). Roerink et al. (2000) proposed the method of simplified surface energy balance index (S-SEBI), by which the evaporative fraction was calculated before estimating the ET. Su et al. (2000Su et al. ( , 2001Su et al. ( , 2002Su et al. ( , 2003 developed the surface energy balance system (SEBS) for estimating surface ET. By determining a series of physical parameters, the established algorithm linking surface roughness of heat transfer with dynamic characteristics of the surface leads the SEBS model to produce evaporative fraction, and then drive ET in the SEBS. Bastiaanssen et al. (1998aBastiaanssen et al. ( , b, 2000a had also systematically published a series of essays in regard to estimating regional surface ET with the Surface Energy Balance Algorithm for Land (SEBAL) model based on surface energy balance equation using Landsat TM data. The SEBAL model was widely used in several countries and areas under different weather conditions (Bastiaanssen et al., 1998a(Bastiaanssen et al., , b, 2000a(Bastiaanssen et al., , b, 2002(Bastiaanssen et al., , 2005Allen, et al., 2005). Besides, Kustas et al. (1990Kustas et al. ( , 1996aKustas et al. ( , b, 1999Kustas et al. ( , 2006) estimated ET within a two-source energy balance model of heat transfer using remote sensing data.
After the availability of the Moderate Resolution Imaging Spectroradiometer (MODIS), Hafeez et al. (2002) estimated the regional distribution of ET in upstream region of Pumapanga River in the Philippines with SEBAL model using MODIS data. These estimates based on the Penman-Monteith method were compared with the measured data collected from the meteorological stations and they were in a good agreement. Nishida et al. (2003a, b) calculated the evaporative fraction of the vegetation and bare soil (e.g., EF v and EF s ) at the pixel level, and then calculated the total evaporative fraction in pixels with weight of vegetation coverage. With the difference between of net radiation and soil heat flux, at last the surface ET in pixel can be calculated. Nagler et al. (2005a, b) estimated ET with empirical relationship between MODIS Vegetation Index products and ET via the use of a regression model.
According to the classification of Courault et al. (2005), the ET algorithms can be divided into four categories: (i) empirical direct methods (Choudhury et al., 1984(Choudhury et al., , 1994(Choudhury et al., , 1998, (ii) residual methods of the energy budget (SEBAL, S-SEBI, SEBS) (Bastiaanssen et al., 1998a, b;Roerink et al., 2000;Su et al., 2002), (iii) deterministic methods that is the Soil-Vegetation-Atmosphere Transfer (SVAT) model, and (iv) vegetation index methods (Kogan et al.,1995;Allen et al., 2005;Neale et al., 2005;Garatuza-Payan and Watts, 2005). The models with residual methods of the energy budget are better models to be operated in regional ET estimate with remote sensing data. The Surface Energy Balance with Topography Algorithm (SEBTA) proposed in this paper belongs to the second category. Comparing SEBTA model with other models (SEBAL, S-SEBI, SEBS) in the same category, Table 1 summarizes some of the salient characteristics for the purpose of differentiation. Advances in SEBTA were made possible through improving the ET model at the practical level on the basis of all previous work. With the same energy balance principle and boundary layer similarity theory, SEBTA establishes its regional water and heat fluxes via using MODIS satellite remote sensing data and introduces the impact factors of topography and land cover into the context of ET calculations. To achieve such a goal, SEBTA requires retrieving the relevant surface parameters (land surface temperature, surface albedo, etc.) and vegetation structure parameters (vegetation cover, normalized difference vegetation index (NDVI), etc. for use in conjunction with conventional ground truth data and meteorological data. Overall, SEBTA demonstrates three advantages: (1) the improvement with the topography information -the improvement of surface kinetic parameters that mainly include the surface roughness (z om ) and the zero-plane displacement (d) both of which can be calculated dynamically in combination with digital elevation model (DEM) and land cover; (2) roughness estimation, and (3) automatic calculation for separating wet and dry pixels -the discrimination of dry and wet pixels automatically with spatial characteristics expressed by land surface temperature (LST) and modified Table 1. The comparison of characteristics of different ET models associated with residual methods of the energy budget (Zhan et al., 1996;Timmermans et al., 2007)  soil adjusted vegetation index (MSAVI). All these three improvements were integrated with continuous calculations of water/heat fluxes in support of the ET calculations under heterogeneous terrains in a large area and long time period over years. With these advantages, large-scale ET maps with higher resolution than the interpolated data of ground-based measurements on a long-term basis can then be produced to provide practical support for applications such as agricultural irrigation and drought monitoring especially in those areas with heterogeneous terrains. Case study in northern China was drawn for the purpose of demonstration of the SEBTA modeling practice.

The fundamental theory of the SEBTA
To estimate the DAET, all of the heat fluxes available at the surface of the Earth are deemed as a response to the streams of solar and thermal radiation. They can be disaggregated for heating and cooling of the air (sensible heat), evaporating water from soil and vegetation (latent heat), and heating or cooling of the soil (soil conduction). Monitoring and modeling the energy balance at the Earth's surface has a critical role in improving the estimation of ET within the coupled energy and water cycles. The most common way to estimate ET is to solve for the latent heat flux, LE, as a residual in the energy balance equation for the land surface. However, the net radiation and soil heat flux are the driving factors of ET which deeply affect the closure in the energy balance equation. The residual method can be expressed as: where R n is the net radiation (Wm −2 ), G is the soil heat flux (Wm −2 ), and H is the sensible heat flux (Wm −2 ). The quantity R n -G is commonly called the "available energy"; remote sensing methods for estimating these components are described in Kustas and Norman (1999). With remotely sensed estimates of solar radiation, differences between these estimates and ground-observed R n -G are typically within 10% of one another . As one of the residual methods of the energy budget, the SEBTA was developed based on the energy balance principle and aerodynamics turbulence theory. Figure 1 holistically delineates the working flowchart of the SEBTA that is particularly organized to show how the terrain complexity can be taken into account as an integral part of the energy balance equations. Whereas the left hand side describes the use of MODIS images to support albedo, surface emissivity, land surface temperature, and vegetation index, the middle portion and right hand side signifies the refined adjustment of the ET algorithm reflecting the impact of heterogeneous terrain at varying temporal and spatial scales. The following equations elucidate the rationale embedded in those building blocks of Fig. 1.
In Eq. (1), net radiation at the bottom of Fig. 1 can be calculated on the basis of the land surface radiation as follows: where R s↓ is the incident solar short wave radiation (Wm −2 ), also known as total solar radiation; α is the surface albedo; ε s is the surface emissivity; σ is the Stefan-Boltzmann constant (5.6696 × 10 −8 W m −2 K −4 ); T s is the surface or canopy temperature (K), retrieved from remote-sensing data, such as MODIS data; T a is the air temperature (K) of reference height (Z2); ε a is the atmospheric emissivity, and it can be calculated by the empirical formula (Bastiaanssen et al., 1998a, b).
The R s↓ can be expressed as:  (Bastiaanssen et al., 1998a(Bastiaanssen et al., , 1998b  where, I 0 is the solar constant (1.367 W m −2 ); R 0 is the vertical incidence of solar radiation on top of atmosphere (W m −2 ); θ is the sun zenith angle on slope; τ b is the direct atmospheric transmittance; (1/ρ) 2 is the revised coefficient of Sun-Earth distance or Earth orbit. Following the same logic as shown in the literature, the SEBTA produces the instantaneous soil heat flux that is defined as a function in terms of surface albedo, vegetation index, and surface temperature (Bastiaanssen et al., 2000a, b).
where NDVI is the Normalized Difference Vegetation Index; and T s is the surface temperature (K). In particular, G water = 0.5R n is employed for water body in the study area. The sensible heat flux (H ) is in a form of the heat exchange between the surface and a level in the atmosphere, which can be expressed as (Zibognon et al., 2002): where ρ a is the air density (kg m −3 ); c p is the air heat capacity at constant pressure (1.004.07 J/(kg K)); T s is the surface or canopy temperature (K); T a is the air temperature (K) of reference height (Z2); dT is the temperature difference (K) over the two heights of Z2 and Z1; r ah is the aerodynamic resistance (m s −1 ) between Z2 and Z1. The calculation of H of each pixel has to be carried out with an iterative procedure to minimize the discrepancy, which is deemed as a methodological advancement in this study. Figure 2 shows the flowchart of iterative procedure for the H calculation. According to Eq. (6), there two unknown parameters (dT and r ah ) to obtain the value of H . It is known that two types of pixels with extreme digital number (DN) (e.g., very dry, very wet) in satellite images are oftentimes present, and can be referred to as "dry" or "wet" points in image processing. To calculate the aerodynamic resistance, we assume that the atmospheric temperature and land surface temperature at the reference height has a linear relationship.
Thus, the H under unchanging or stable aerodynamic resistance may be estimated by an iterative procedure in calculation with respect to Monin-Obukhov similarity theory. In our analysis, the wind speed at 200 m height whose value is relatively stable was selected as the wind speed in the mixing layer, which would not be affected by surface roughness. Such a wind speed, u 200, can therefore be estimated based on the measurements at the ground meteorological stations based on the logarithmic wind profile. Yet, this may not be always true in view of terrain complexity. 9 points in image processing. To calculate the aerodynamic resistance, we assume that the 4 atmospheric temperature and land surface temperature at the reference height has a linear 5 relationship.
where u 200 and u 10 are the wind speed at 200 m and 10 m height, respectively (m s −1 ); z om station is the roughness length for momentum transfer around the meteorological stations and the value is assigned with Table 2. The aerodynamic resistance under the neutral atmospheric condition can be expressed as: where k is the von Karman constant (=0.41, dimensionless); Z 2 is the reference height (= 200 m); Z 1 = z oh , is the heat transfer roughness (m); u * is the wind-friction velocity (m s −1 ), which can be expressed by the logarithmic wind profile equations under neutral atmospheric condition as expressed below: where u(z) is the average wind speed (m s −1 ); z is the height for wind speed measurement (m); d is zero-plane displacement height (m), and z om is the surface roughness (m). Under neutral atmospheric condition, the wind-friction velocity, u * , can be obtained by Eq. (9) and then aerodynamic resistance, r ah , can be estimated by using Eq. (8). The calculation of H can then be carried out via the following steps: where a and b are the coefficients in the linear equation above; T s is the land surface temperature (K). With this consideration, SEBTA may be empowered to address come delicate phenomena such as the local advection effects, especially for the regions of image where irrigated vegetation is in vicinity of bare soil. The pixel covered with water or with high coverage of vegetation where the temperature could be very low may be selected as a "wet" point. Assuming that all the available heat fluxes at this pixel location is used for ET (LE Wet = R n,Wet G Wet ), which leads to H Wet = 0 and dT Wet = 0. On the other hand, the pixel covered with low degree of vegetation where the temperature could be very high may be selected as a "dry" point. Assuming that there is no evaporation at this point (LE dry = 0), which leads to calculate H dry = R n,dry G dry and dT dry = (R n,dry G dry )r ah,dry /ρ a c p. Then we can get:   (Brutaert, 2005).

Surface description z om station
Large water surfaces ("average"), Snow, mud flats, Smooth runways 0.0001-0.0005 Short grass 0.008-0.02 Long grass, prairie 0.02-0.06 Short agricultural crops 0.05-0.10 Tall agricultural crops 0.10-0.20 Prairie or short crops with scattered bushes and 0.2-0.40 tree clumps, Continuous bushland Bushland in rugged and hilly (50-100 m) terrain 1.0-2.0 Mature pine forest 0.80-1.5 Tropical forest 1.5-2.5 Fore-Alpine terrain (200-300 m) with scattered 3.0-4.0 tree stands Thus, the dT can be re-organized as follows: in which (Conrad et al., 2007): where ρ is the air density; λET r is the potential ET; dT wet , T wet , H wet , R nwet , G wet (or dT dry , T dry , H dry , R ndry , G dry ) are the temperature difference, surface temperature, sensible heat flux, net radiation flux, and soil heat flux at wet point (or dry point), respectively. In this practice, the surface temperature of the study area was adjusted to the same temperature of reference height (average height of the region) to eliminate impact of surface temperature difference on the calculation of the latent heat flux caused by the elevation change. These coefficients in these equations are set by empirical values based on our local observations. The surface temperature of the study area was adjusted to the same temperature of reference height (average height of the region) with the T s dem equation as follows: where h is the altitude above the sea level (m); h mean is the average altitude of the study area (m); and T s is the surface temperature (K). 0.0065 is derived based on a constant lapse rate (−0.65 K/100 m) (Stull, 2001). The dT (temperature difference between Z2 and Z1) can also be determined by the following equation: When the polar-orbital satellite overpasses the study area, the air temperature can be determined as below: By identifying dT and r ah , the sensible heat flux H can finally be obtained according to Eq. (6). Due to the instability of the planetary boundary layer, Monin-Obukhov similarity theory was adopted, from which the Monin-Obukhov length, L, can be linked with the parameters of r ah , u * , H , dT and air density ρ a via the following equation where g is the acceleration constant of gravity (9.807 m s −2 ). This can be used to support the iterative procedure for the calculation of H when the value of r ah becomes stable (Fig. 2). Further, in this paper, wind-friction velocity (u * ) and aerodynamic resistance (r ah ) were revised by using integrated stability correction functions for momentum and heat transfer based on the similarity theory (Su, , 2002, respectively. where m (x) and h (x) are integrated stability correction functions for momentum and heat transfer; d is zero-plane displacement height (m); k is the von Karman constant (=0.41, dimensionless); L is Monin-Obukhov length (m), which reflects the atmospheric turbulence conditions of the near-surface layer used to determine the stability of atmospheric stratification condition. The adjustment function of dynamic roughness ( m (x)) and heat transfer stability ( h (x)) under unstable atmospheric conditions (L < 0), respectively, are defined as follows: The adjustment function of dynamic roughness ( m (x)) and heat transfer stability ( h (x)) under stable atmospheric conditions (L > 0), respectively, are defined as follows: The adjustment function of dynamic roughness ( m (x)) and heat transfer stability ( h (x)) under the neutral atmospheric conditions (| L | → ∞, x → 0), respectively, are defined as follows: where x = z/L. The above equation system can be finally summarized by the dashed line at the bottom of the Fig. 1 showing the streamlines of required iterative calculation in Fig. 2 (Gieske, 2003). Finally, based on the energy balance equation, the instantaneous latent heat flux can be computed as shown in Eq. (1). Then, instantaneous ET (mm) can be calculated as: in which λ= [2.501-0.00236 (T a -273.15)] × 10 6 (Xie, 1991). The ET estimated by such a remote sensing method represents instantaneous conditions. To achieve the estimation of the daily ET, the instantaneous ET must be extended temporally. Several commonly used methods, including statistical empirical method, sine relations method and the ratio method, can be applied. In this paper, we followed the ratio method to fulfill such estimation from instantaneous ET to daily ET (Jackson et al., 1981. Polar-orbital satellites generally overpass our study area at the local time of 10:00 to 11:00 a.m. Based on the daily variation features of the evaporative fraction, the instantaneous ET estimated by the remote sensing method can be extended to daily ET throughout the day and 24 h ET may be accumulated. Yet sometimes, 24 h evaporative fraction is not appropriate and is not equivalent to daily evaporative fraction as at night time this evaporative fraction is different and normally there is no transpiration by crops in case agricultural land is the predominant one. In this study, we assume: where 24 is for 24 h average evaporative fraction; d is for daily evaporative fraction; is for the instantaneous evaporative fraction; LE d is the cumulative ET during the day; F d is an accumulated value during the day for a energy component in the energy balance equation, which could be a net or effective energy of radiation; LE and F are instantaneous values during the day; ER is evaporative ratio; If F take the LE + H (effective energy R n −G), then ER is the evaporative fraction (EF) . Thus, after the evaporative fraction can be obtained, we can produce the daily ET based on the formula as below: where R n24 is the daily net radiation flux (W m −2 ); G 24 is the daily soil heat flux (W m −2 ); and the other parameters are same as defined above. The integral of sunshine time can produce the total solar radiation energy on a daily basis over the unit area. Ignore the changes of the distance and the solar declination in day time, the total solar radiation of 24 h that accounts for the effects of terrain on R s can be defined as: where t is time of the sunshine hours starting from midnight as 0; t 1 and t 2 are sunrise and sunset time, respectively; I 0 is the solar constant (1.367 W m −2 ); τ b is direct atmospheric transmittance. The −ω sr and ω ss are solar time angles on a hourly basis corresponding to sunrise and sunset time, respectively, and T is the total day time (h), by which t = T (ω + π )/2π, dt = (T /2π )dω. Hence, the computation of effects of terrain on R s is achieved by the following formula.
Combined with the computation of effects of terrain on H (Eq. 6) by the effects on dT (Eq. 13) and T s (Eq. 14) achieved, the SEBTA can be functional for the following ET estimation.
Since the soil heat flux remains the roughly same from daytime upward to night time downward, the daily cumulative value of soil heat flux is assumed as approximately zero and the soil heat flux is therefore negligible in general. Combination of Eq. (21), the formula (22), the 24 h ET can be simplified as follows: According to the above parameter settings and modeling mechanisms, computer program can be designed by using Z. Q. Gao et al.: A coupled remote sensing and the Surface Energy Balance Arc/Info ® 9.0 Macro Language (AML) and Compaq Visual FORTRAN ® 6.5 mixed-language programming to generate the ultimate SEBTA computational code. The SEBTA computer package can be operated in a Microsoft Windows system using the ESRI GRID module as the major data format.

The realization of heterogeneous terrain in the SEBTA
To fully elucidate the advantage of the SEBTA, there is a need to realize how the SEBTA assimilates the heterogeneous terrain data to generate the solar radiation at the pixel level. The impacts of the topographic factors such as elevation, slope, and aspect on net radiation (R n ) and sensible heat (H ) can be addressed by the SEBTA reflecting the changing total solar energy per unit area received from the Sun over time, which is manipulated through the integral of total sunshine time in the day. Because the heterogeneous terrain may cause mutual shade of sunshine and the actual terrain undulating is irregular, therefore, the sunshine time in the day per unit of slope cannot be expressed with available mathematical formula. Albeit the total amount of radiation in a day per unit of slope can only be obtained by running the integral over a sequence of temporal subsections (Hasager et al., 1999). The following detailed illustration shows how the SEBTA achieves this sophisticated manipulation. Firstly, the information of slope, aspect, latitude, longitude, and elevation of each grid can be obtained using the DEM data at the time step (e.g., 0.5 h in our case). The time interval for integration over solar time angle can be selected subjectively. In a range of [−ω sr , ω ss ], the time interval may be set up as ω, and the sunshine time period during a day is divided into n time periods where n may be obtained by using n = int((ω ss + ω sr )/ ω) + 1 (i.e., int is the integer function). Over the n +1 time intervals, the solar time angle expressed from the start time to the end time within the n time periods are as below: where ω is solar time angle; −ω sr is angle of sunrise time; ω ss is angle of sunset time; ω is time step in a range of [−ω sr, ω ss ]. At the same time, the topographic shading coefficient of corresponding time interval can be calculated with the aid of the Hillshade function embedded in the GRID module of the Arc/Info system and elevation of grid can be provided by using DEM: where h is the elevation with a slope unit; SHADOW is keyword in this Hillshade function; A s is the sun azimuth; Z s is the sun zenith angle; and b i is the topographic shading coefficient (i = 0,1,2,...,n). Whether the Sun can shine on the top of each grid has to be determined by the following equation in each time period: where g 0 is start time of topographic shading; and g i is defined as a topographic shading parameter.
In the period of (ω i−1 ,ω i ), the condition of sunshine per unit of slope depends entirely on the sunshine condition of two consecutive time points, if the Sun can not shade over the whole time period. Otherwise, half of time period can be taken into account as the Sun can shin (or shade) partially. As the g i is defined as a topographic shading parameter, it would support the calculation of solar radiation in Eq. (2), leading to pass the information of heterogeneous terrain into the regular calculation of the surface net radiation. Our observational evidence confirms that the same reference height may be applied to convert the remotely sensed land surface temperature to the air temperature so as to minimize any possible inconsistency of elevation basis in the calculation of the sensible and soil heat fluxes due to the elevation changes. This logic is carried out by properly handling the T s dem in Eq. (14). With the aid of the DEM data, the impacts of surface characteristics on the SEBTA can be fully reflected in the calculation of R n and H .

Automatic determination of the dry and wet pixels in the SEBTA
One of the key factors to improving the ET estimation is to discriminate the wet and dry pixels automatically so as to determine the sensible heat flux (H ) more accurately by getting a better understanding of the vertical temperature difference (dT = (T s − T a )) and the aerodynamic resistance (r ah ). This can be made possible by examining the extreme conditions present on the satellite images. Two conditions may be distinguished, namely dry and wet areas. To determine which pixel belongs to a dry area or a wet area, a linear theory proposed by Menenti and Bastiaanssen (1989) has to be applied. In this study, however, the problem of mixed pixels was not considered to ease the processing of automatic computation Different vegetation index such as NDVI, Leaf Area Index (LAI), Soil adjusted Vegetation Index (SAVI), etc., and land surface temperature (LST) provide some evidence to help discriminate the wet versus dry pixels on the satellite images. In previous studies, Moran et al. (1991Moran et al. ( , 1996 suggested the use of 0.1 of SAVI to represent complete bare soil, and 0.8 of SAVI to represent complete vegetation cover. Verstraeten et al. (2005) suggested using combinations of NDVI and LST. Allen et al. (2005) chose the dry and wet pixels with the help from combined LAI and LST. In this study, we adopted the combination of modified soil adjusted vegetation Index (MSAVI) and LST to discriminate dry and wet pixel.  (27), (28) and (29). The roughness z oh can ally obtained by equation (29) e H eff is the effective height of objects in the pixel (m), which can be estimated by As a consequence, Fig. 3 applies the Vegetation Index/Temperature Trapezoid (VITT) concept to differentiate the dry versus wet points with respect to combined vegetation index and surface temperature. It uniquely defines four extreme conditions at the four corners within the Trapezoid area. For example, the water stress condition can be identified at the upper right corner where the largest LST occurs and the value of MSAVI is equal to 0.8. At the lower right corner, the saturated canopy indicative of the wet point appears to have lowest LST when the value of MSAVI is 0.8. On the other side, when the value of MSAVI is 0.1, it uniquely defines another extreme condition for dry point where the LST turns out to be the largest. This results in the situation that the lowest LST may appear when the value of MSAVI is 0.2 reflecting the saturated bare soil condition. The SEBTA may thus run through all pixels automatically to discern the wet and dry pixels estimate the ET according to such a VITT concept. In summary, wet point can be discerned when LST is lowest and MSAVI is 0.8 whereas dry point can be discerned when LST is largest and MSAVI is 0.1.
The Modified Soil Adjusted Vegetation Index (MSAVI) is a VI developed by Qi et al. (1994) to describe the effects of soil brightness in the background. It is: where ρ red is the red band (0.63 ∼ 0.69 µm) reflectance; ρ nir is the near red band (0.76 ∼ 0.90 µm) reflectance; ρ blue is the blue band (0.45 ∼ 0.52µm) reflectance; and ρ green is the green band (0.52 ∼ 0.60 µm) reflectance. Besides, L is adjustment factor that is set to address the minimum background effects (L = 0.5).

Dynamic determination of surface kinetic parameters in the SEBTA
In the process of thermal diffusion of surface water/heat flux, the surface kinetic parameters mainly refer to the surface roughness z oh and zero-plane displacement d (Farah et al., 2001). These surface kinetic parameters can be determined using a look-up table on a monthly basis according to the land data assimilation system developed by the Goddard Space Flight Center of National Aeronautics and Space Administration (NASA) in the US. The surface roughness and the zero-plane displacement can be estimated by surface effective height with the available empirical formulae as shown in Eqs. (27), (28) and (29). The roughness z oh can be finally obtained by Eq. (29) at the pixel (Wilson et al., 2003).
where H eff is the effective height of objects in the pixel (m), which can be estimated by considering object's height, coverage, and spatial distribution of vegetation and surface objects within the pixel. Equation (30) produces the effective vegetation height within each pixel.
where the H eff is the effective vegetation height within each pixel (m), VI is the vegetation index, H max and H min are the largest and the smallest effective vegetation height (m), H is the difference between the minimum effective height (H min ) and the maximum effective height (H max ), VI max and VI min are the largest and the smallest vegetation index, respectively. The effective height and vegetation index of different vegetation types are determined with respect to the surface parameter look-up table developed by NASA (see Table 3). The other surface objects which do not change with the seasons are taken into account as constant effective heights, based on a 1/100 000 land use database (Liu et al, 2005). This equation, however, is only used to compute the effective vegetation height within each pixel. As for other land covers (e.g., water body, urbanized land, rural residential land, constructed land and barren land) their roughness estimates (effective height) were obtained by using the lookup table (Wilson et al., 2003).
Even with dynamic change of effective height based on MSAVI parameter, it is not the full use of available data for determination of roughness parameters, as still effects of atmospheric turbulence on roughens parameters are not included. With the aid of remotely sensed MSAVI, however, integration between Eq. (30) and Table 3 may enable us to estimate the dynamic surface effective height, thereby producing the zero-plane displacement based on the empirical  Fig. 4). Topographically, elevation is generally high in the northwest and low in the southeast where the Yellow sea borders. Eastern region of this study area is the Taihang Mountains in Shanxi Plateau and northern region is the Mongolian Plateau and Yanshan Mountains, whereas east and southeast region is vast plains close to the Yellow Sea. The Haihe River Basin belongs to temperate climate zone from where the humid, semi-humid and semi-arid regions are affected by the continental monsoon. The multi-year average rainfall is 535 mm and the mean annual temperature is in a range of 1.5 ∼ 14 • C. The annual average relative humidity is in a range of 50 ∼ 70%, and the annual ET is in a range between 850 ∼ 1300 mm. With ample light and heat radiation for crop growth in the study area, the irrigated land is suitable for agricultural production. Thus, the Haihe River Basin is one of China's three major grain production bases. There is only one, the Yucheng Ecological Experimental Station (Yucheng station hereafter), with ground-based measurement of lysimeter and eddy co-variance measure-ment equipment in the study area. The temporal resolution of lysimeter observations is 1hr. The Yucheng station is run by the Chinese Academy of Sciences. It is located at the North China Plain (116 • 34 E, 36 • 57 N) in Yucheng County, Shandong Province. The elevation is 28 m. The site belongs to monsoon climate of warm temperate and semihumid environments. The annual average temperature is 13.1 • C, the precipitation is 582 mm, the total solar radiation is 5225 MJ m −2 , the sunshine hours are 2.640 h yr −1 . The area is covered mainly by farmland and grassland (NRC, 1992;Yu et al., 2006).

Data collection and analysis
The dataset collected and processed in this paper includes MODIS data, DEM data, land use data, climate data, and ground-based measurements at the Yucheng station. Application of the SEBTA to the Haihe River Basin requires using several MODIS products to estimate the regional water/heat fluxes; and they include: MODMG-GAD (1 km), MOD09GQK (250 m), MOD09Q1 (250 m), MOD09GHK (500 m), MOD09A1 (500 m), MOD11A1 (1 km), and MOD11A2 (1 km). Due to the geographical location of the study area, cloud contamination is inevitable. 8-day products of MODIS were used for LST measurements. To obtain sufficient amount of satellite images of the MODIS products for data analysis, one threshold was therefore set up to screen out any Terra MODIS product which has the cloud cover greater than 40% between 2005 and 2006. This ends up obtaining a total of 48 days of the MODIS products. Table 4 summarizes all the MODIS products that were used in this study.
For implementation, the data bands and quality control bands in various MODIS data products were extracted with MRT tools (MODIS Reprojection Tool) provided by NASA. The parameters in the SEBTA may be manipulated with mixed codes of ESRI Arc/Info ® and FORTRAN language. Year The selected    The projection of UTM/WGS84 whose central meridian is 117 • E/50 • N Zone was used in this computational process. The DEM data collected by Shuttle Radar Topography Mission (SRTM) was applied. The available DEM data in this study is actually SRTM3 whose resolution is 90 m. Besides, the meteorological data used in this study are dataset collected at 89 stations which are run by the National Meteorological Centre, China Weather Bureau. The temporal resolution of climate data used here is one day. This dataset includes average temperature, maximum temperature, minimum temperature, precipitation, average wind speed, amount of cloud and others. This dataset were vectorized and interpolated as grid datasets with UTM projection. All data were resampled to the uniform resolution of 1000 m for ET estimation.
The land use and land cover (LULC) data used in this paper was derived from the dataset with scale of 1:100 000 in a database developed by the Chinese Academy of Sciences. This LULC dataset was generated based on the proper interpretation of Landsat TM/ETM images and was validated with ground truth data. The climate data was collected from the National Meteorological Center of China Weather Bureau. The spectrum of climate data includes average temperature, maximum temperature, minimum temperature, precipitation, average wind speed, amount of cloud and others. All datasets were interpolated as grid datasets with UTM projection in advance to ease the application in geographical information system (GIS).  Comparison between the measured ET and the simulated ET at the Yucheng station.  Fig. 5 due to the wheat harvest activities, which led to a large difference between the measured and the simulated values. These values with big discrepancies were not used later on in producing Fig. 6.

Model validation
In this study, we used instantaneous estimations of ET i /(R ni − G i ) to compute the ratio, and ground measurements at Yucheng station as shown in Fig. 5a to validate the evaporative fraction. The Polar-orbiting satellite with MODIS terra generally overpasses China between local time 10:00 to 11:00 a.m. so that the instantaneous ET estimated by remote sensing can be extended to daily value or 24h average  Fig. 5b. In general, the simulated ET data by using the SEBTA shows a good agreement with the measured ET values. In order to avoid the pixel distortion caused by the satellite drift and other reasons, taking into account the impact of satellite resolution, the average values of simulation within about 1kmx1km area were extracted and compared with the observed values at Yucheng station for model validation. Table 5 lists both weekly and monthly ET which can be computed with crop coefficient, K c , retrieved with satellite data and the weekly or monthly reference ET. The comparison between the measured ET and the simulated ET between 2005 and 2006 at the Yucheng station shows that the correlation coefficient is 0.87. Randomness becomes phenomenal when the ET values are either relatively large or small. The offset of estimation errors of ET over different time scale might be apparent during crop growing period (March-November). In order to eliminate the impacts of rainfall and other factors on lysimeter ET measurements on weekly and monthly scales, Table 5 summarizes the comparisons in 2005 based on monthly and weekly records somparatively. The weekly average ET of 2005 is 3.19 mm day −1 as measured and 2.90 mm day −1 as simulated. When comparing the values between measured and simulated ET, the offset could be salient. In our case, the weekly measurements have rel-ative error of −6.3% whereas the monthly counterparts have 2.09%. From March to November in 2005, the total ET is 736.22 mm as simulated and 726.86 mm as measured. In summary, on a monthly basis, the relative error of ET estimation is 1.29% of 2005 and 4.05% of 2006, respectively. When different temporal scale is taken into account, the relative error is −11.34% based on daily averages, −6.30% based on weekly averages, and less than 4% based on monthly and quarterly averages.
The consistency index that is also called index of agreement based on pair wise comparison between simulated and measured values is actually a dimensionless index, which can be defined as Idx as below (Willmott, 1982): where O i is the measured value; O is the spatial average value; P i is the simulated value. In general, both simulated and measured values of ET show a good agreement, and the consistency index is up to 0.92. The discrepancies that are closely tied with energy closure issues inherently might partially arise from the use of energy balance approach for simulated values and observed values collected by a lysimeter using water balance approach. The other possibility is related to the linkage between the instantaneous ET from remote sensing and the hourly ET from lysimeter.

30
DEM data of the study area resulting in distinct terraces on each 500 m contour line. Figure   4 6(b) delineates the LULC patterns derived by the Landsat TM/ETM images (Liu et al., 2005).

Results and discussion
The inclusion of both DEM and LULC data into the coupled remote sensing and SEBTA aids in the estimation of regional ET, making the SEBTA no longer unnecessarily relies too much on the availability of local micro-meteorological measurements. It greatly improves the operational efficiencies when monitoring the water and heat fluxes under largescale heterogeneous terrain. Based on the SEBTA, this study demonstrates an application using two bands of the MODIS data with 250 m resolution and the MODIS 1km products of LST, vegetation index, albedo etc. Combined with LULC and DEM data we were able to derive the ET maps of 2005 in the Haihe River Basin. Hence, the selected MODIS data in 2005 support the DAET estimation leading to examine the DAET over four seasons (i.e., spring (March-May), summer (June-August), autumn (September-November) and winter (December-February)). Finally, the spatial pattern of the ET distribution can be summarized with the aid of this seasonal ET data.

The impacts of heterogeneous terrain on ET estimation
The bases for classification of "rough" and "moderate" terrains are based on the distinction in modern geography in China (Huang et al., 1999). The terrain within the study area is not gentle. It is considered the rough terrain in north-west with a large quantity of terrain components, and considered moderate terrain in southeast. Somewhere in between is the moderate rough terrain as part of the transitional zone. Overall, the Taihang Mountains and Shanxi Plateau areas are located at the western part, the Mongolian Plateau and the Yanshan Mountains are located at northern part, and the vast plains are located at the eastern and southeast parts of the study area. The plains and plateaus intermittently dominate the landscape changes from southeast to northwest in the study area accounting for 72% of the total area in which the mountainous area occupies about 28% of the total area. Since the DEM is in a grid format, we can analyze the spatial characteristics of the study area with GRID module in ARC/INFO software to generate the contour lines. Figure 6a shows the DEM data of the study area resulting in distinct terraces on each 500 m contour line. Figure 6b delineates the LULC patterns derived by the Landsat TM/ETM images (Liu et al., 2005). The gentle terrain with elevation below 500 m dominates the northern part of the coastal area of the North China Plain that is located at the eastern part of the study area where possesses good potential for crop growth. Geographically, this area is composed of the northern part of Shandong Province, the eastern part of Hebei Province, and the Cities of Beijing and Tianjin, which accounts for about 45% of total study area. The DAETs of 2005 being summarized in spring, summer, autumn and winter across the gentle terrain are 1.09 mm day −1 , 2.78 mm day −1 , 3.40 mm day −1 and 0.31 mm day −1 , respectively. With bigger terrain slopes in the northwestern of the study area with moderate rough terrain, the elevation changes sharply from 500 m to 1500 m in connection with the southern part of the Inner Mongolia, northwestern part of Hebei Province and the eastern part of Shanxi Province. Such an area with moderate rough terrain accounts for about 47% of the entire study area. The DAETs of 2005 being summarized in spring, summer, autumn and winter in this region are 0.98 mm day −1 , 1.61 mm day −1 , 2.70 mm day −1 , and 0.27 mm day −1 , respectively. The area with rough terrain where the elevation is greater than 1500 m is located at the northwest corner of the study area, which already reaches the central part of the Inner Mongolia Plateau, Loess Plateau, and the Taihang Mountains. This area accounts for about 8% of the entire study area. The DAETs of 2005 being summarized in spring, summer, autumn and winter in this region are 1.43 mm day −1 , 1.75 mm day −1 , 3.10 mm day −1 , and 0.29 mm day −1 , respectively. Although the elevation-based ET patterns can be further characterized seasonally later on, the present observational evidence confirms that the DAET turns out to be the biggest in the area where the elevation is below 500 m and becomes smaller in the area where the elevation is between 500 m and 1500 m. Yet the smallest ET occurs in the area where the elevation is greater than 1500 m.
Temporal analysis in accordance with the seasonal ET maps of 2005 indicates that the patterns of DAET are quite different over seasons. Whereas the DAET distributions in spring and winter are not terrain-dependent, they are indeed the case in summer and autumn. According to different ranges of elevation, the mean ET values falling in different ranges were computed in a GRID module of ARC/INFO software package, and this is the basic spatial analysis function in a GRID module to generate seasonal ET spatial data (spring, summer, autumn and winter). Figure 7 thus depicts such observations. For example, the summer DAET is smaller in the west and northwest of the Loess Plateau and the Inner Mongolia Plateau while it is much larger over the plains. The distribution of autumn DAET varies from the southeast with larger ET values to northwest with smaller ET values. Both ET distributions in summer and autumn are actually affected by the terrain and climate factors simultaneously.
The seasonal ET of 2005 may be further partitioned based on varying elevations (see Fig. 8). In Fig. 8, the ETSPR, ETSUM, ETAUM, and ETWIN represent the DAET within spring, summer, autumn, and winter, respectively. By looking into it in a greater detail, the autumn and summer DAETs have something to do with elevations. The DAET in areas below 50 m is larger and DAET in areas between 50 m and 1300 m is significantly smaller. Yet ET in areas above 1300 m becomes larger. The changes of DAET in areas above 800 m are obvious. Overall, the DAET in areas below 800 m is larger than 2.00 mm day −1 and in areas above 800 m is smaller than 2.00 mm day −1 . Such a descending trend in summer and winter DAET over elevation remains looming the same with a much tempered pace. Autumn DAETs exhibit relatively larger values over all elevations although they are relatively smaller in areas between 800 m and 1500 m. Yet it becomes larger than 3.00 mm day −1 in areas with lower elevation and the changes of DAET in areas below 800 m is not significant with regard to elevation. In summary, the impacts of the seasonal ET simulated by the SEBTA on terrain vary over time. Whereas the winter DAET has no obvious response to elevation changes, the summer and autumn DAET do have strong responses to elevation changes to some extent.

The impacts of different patterns of LULC on ET estimation
According to Fig. 6b, the predominant types of LULC in the study area include dry farmland, grassland, and woodland, which account for 51.38%, 25.1%, and 15.48% of total study area, respectively. The remaining types of LULC include paddy field, constructed land, water body, and unused land, which account for only 8% of the entire study area. The grassland is mainly distributed in the Inner Mongolia Plateau and the Taihang Mountains. The woodland is mainly distributed in the Taihang Mountains, the Yanshan Mountains and the Luzhong Mountains. The dry farmland is mainly distributed in the North China Plain, the Taihang Mountains and the Loess Plateau. The barren land is mainly distributed in the Inner Mongolia Plateau and the Bohai Bay. The constructed land, such as the mixed urban or builtup land, and water body are densely scattered around the plain areas.
The daily ET varied widely over different LULC patterns. With the simulated ET derived by the SEBTA, the DAET associated with different LULC can be computed as 1.73 mm day −1 of paddy field, 1.38 mm day −1 of dry farmland, 1.27 mm day −1 of woodland, 1.14 mm day −1 of grassland, 1.85 mm day −1 of water body, 0.85 mm day −1 of urbanized land, 1.35 mm day −1 of village and other constructed land, 1.26 mm day −1 of sparsely vegetated land, and 1.08 mm day −1 of barren land. The differences of ET over different seasons are salient. Tables 6a and b summarize such differences statistically with respect to the types of LULC and the height, respectively. Two extremes of DAET exist. The DAET of water body can be as big as 1.84 mm day −1 and the DAET of urban and rural settlements can be as small as 0.6 mm day −1 and 0.86 mm day −1 , respectively. The DAETs of barren land that is mainly distributed in the northern part of the study area and around the Bohai Sea are as big as 1.78 mm day −1 and 1.50 mm day −1 , respectively. Figures 7 and 9 collectively exhibit such intertwined phenomena of spatiotemporal variations of DAET associated with various vegetation and climatic regions. In Fig. 9, seasonal changes can be characterized over 10 types of LULC from which the yearly trends may be illuminated. It is observed that the temporal changes of spring DAET over different LULC was significant a year round in most seasons.  In spring and winter of 2005, the DAET decreased gradually from paddy land to dry farmland, and to woodland. Yet the DAET of grassland bounced back while the DAET of water body reached an extreme point. The DAETs of urban and rural constructed land are much smaller as compared to that of water body. The increasing trend of the DAET remains from constructed land to sparsely vegetated land. The DAET of barren land became smaller. The DAETs in autumn and summer of 2005 show slightly different trends related to the categories of other constructed land, sparsely vegetated land, and barren land. In Fig. 9, the trends of spring, summer and winter are quite similar. The maximum DAET of spring, summer and winter appears in water body always and the minimum DAET appears in urbanized land. In autumn, however, the maximum DAET appears in paddy fields and the minimum DAET appears in the cities, followed by rural settlements. The trend of ET curve for autumn DAET is markedly different with those in other three seasons. Overall, it shows that the impacts of LULC on DAET simulated by the SEBTA are significant in summer and autumn. With the smallest differences of ET over seasons, grassland thus plays an important role in economic development, water conservation, biodiversity and ecological safety.
ET variations also follow the climatologic regions in the study area changing from the humid Southeast, to semihumid Northeast, and to semi-arid areas Northwest areas. In temperate and humid regions generally four seasons are distinct. The reasons for having big ET variations over different seasons and different LULC in the study area is partially due to the Northern temperate climate across humid, semi-humid and semi-arid regions where the vegetation cover over four seasons is very different. During the growing seasons (i.e., spring and summer), the biome of grass, crop and deciduous forest can be quite productive such that the vegetation cover increase quickly. The growing period of grass and crops is much shorter than that of deciduous species as period of greatest biomass was attained in the summer. In winter, the vegetation canopy in northern China stops growing and it becomes very lean. Even the understory of shrubs and herbs in a mature deciduous forest can hardly survive. This renders such a different response to ET changes over seasons and to land cover as well.

The impacts of the inclusion of elevation, slope, and aspect on ET estimation
The effects of terrain characteristics on water vapor and heat fluxes were introduced and elucidated by the inclusion of DEM and LULC data in the SEBTA. It is worthwhile to investigate the impacts of the elevation, slope, and aspect on ET estimation in a greater detail. We picked up a sub-region where the range of elevation varies in between 8 m and 1532 m. With the SEBTA, the DAETs can be calculated based on each 100 m contour line in the study region. Figure 10 clearly reveals that the black curve that considers the elevation effect is quite different with that does not consider it when the fall data was considered. The elevation of study area ranges from 8 to 1532 m, and the ET values 33 to elevation changes, the summer and autumn DAET do have strong responses to changes to some extent. were computed based on scale of 100 m associated with and without the DEM effects. The land use on these subregions was derived based on the database developed by the Chinese Academy of Sciences. The 2006249 of MODIS images was used for this practice. Obviously, overestimation of ET could appear as the effect of elevation was not taken into account as evidenced by the red curve in Fig. 10. We also tested other seasons and the results are similar to Fig. 10. The spoon-type curves present in Fig. 10 also feature some inherent implications of the inclusion of terrain effects in the study area. Note that although the maximum height is 1532 m, Fig. 10 only addresses transect for simulation which does not cover this highest point of the study area. The ET in the plain areas where the elevation is smaller than 200 m is larger due to the seasonal crop growth and the ET in the area where the elevation is around 200 m becomes smaller due to the presence of barren land. It is known that the vegetation cover in the area where the elevation ranges from 300 m to 1000 m is denser. Combined with the larger wind effects and solar radiation, the ET could increase gradually.
The slope and aspect are deemed important parameters accounting for the heterogeneous terrain effects in the study area. Whereas the aspect mainly affects solar radiation and inclination angle of ground and prevailing wind direction, the increase of slope would aggravate these impacts further. Both of them make significant differences in terms of water and heat fluxes, thereby resulting in different ET estimation. The following sensitivity analysis picked up a subregion where the slope is above 5 • within the study area. In Fig. 11, it shows the net radiation flux can be significantly affected by the aspect. The net radiation flux estimated without considering terrain effects is uniform almost over all aspects. Within the simulation, we found that the maximum net radiation flux appeared on the aspect of 120 • , which is very consistent with the sun azimuth (118.7 • ) associated with the acquisition of an important role in economic development, water conservation, biodiversity and e 6 safety. the satellite images. This implies that the minimum net radiation appeared on the opposite aspect of 300 • at that time. Besides, the ET in northern slope of the southern hemisphere is smaller than that in southern slope of the northern hemisphere. Thus, the distribution of the ET with the changes of slope is also significant. The difference in terms of radiation energy due to the inclusion or exclusion of aspect and slope is truly sensitive to the ET estimation as shown in Figs. 10 and 11 collectively.

Final remarks
Retrospective reviews advise us that the testing areas of the SEBAL model were conducted in small plain regions such as Qattara Depression in an area of 20 000 km 2 in the Western Desert of Egypt and Nile Delta of Egypt in an area of 6950 km 2 (Bastiaanssen et al., 1998a, b). The validation of the SEBS model were carried out based on the datasets collected over a cotton field in Maricopa Farms in central Arizona, shrub dataset in the Lucky Hills, a shrubdominated ecosystem in southeastern Arizona, grasslands in the Kendall, and Efeda in the Barrax region in Spain. All of them covered areas no greater than 0.45 km 2 (i.e., the field is 1500 m east-west by 300 m north-south) (Su et al., 2002). Yet land use was considered by roughness parameters (Su, 2002). The validation of the S-SEBI model simply used a small fraction of a scene (LANDSAT5-TM, path/row for the 192/30) (Roerink et al., 2000). Although the validation accuracy of water-heat fluxes in the above three models in those selected regions is relatively high, none of them account for the challenge due to the presence of terrain complexity at a practical level.
The impacts of terrain factors (elevation, slope, aspect) and LULC patterns on the ET estimation have been a challenge in this field. Besides, the implementation of the ET algorithms to a bigger area where the heterogeneous terrain 37 es was used for this practice. Obviously, overestimation of ET could appear elevation was not taken into account as evidenced by the red curve in Figure   sted other seasons and the results are similar to Figure 10. appears to be influential is even more difficult. Either overor under-estimation may occur causing the credibility issues in ET simulation although some offsets might exist to cover the drawbacks due to the varying temporal scales applied for the ET estimation. In our paper, it was proven that in areas of higher elevation and larger shade slope, the ET can be over-estimated resulting in large discrepancies. In any circumstance, the inclusion of terrain factors (elevation, slope, aspect) and LULC patterns on the ET simulation in the context of the SEBTA was hard to be proven effectively in the present study because there is not any station located at different place with differing elevation, aspect, or slope. Yet the framework developed by this study can be transferable to future applications.

Conclusions
This present study develops a coupled remote sensing and SEBTA to estimate the actual ET under heterogeneous terrain, which is geared toward integrating the energy balance principle and aerodynamics turbulence theory with varying terrain and landscape effects. With the aid of DEM and LULC, the SEBTA enables us to account for the impacts of heterogeneous terrain and LULC on the ET estimation. DEM helps the calculation of net radiation and an adjustment on temperature difference. Yet the advection terms in relation to topography are hard to be included. The affects of kinetic parameters (roughness and zero-plane displacement) in relation to LULC and terrain factors (elevation, aspect, and slope) produced by the DEM database were synergistically woven together within the simulation steps. Besides, the dry and wet pixels embedded in MODIS images can be automatically discerned so that the DAET can be simulated by a more accurate way. Such a strategy not only expands the application potentials but also increases the maneuverability of the SEBTA implementation at a practical level. the area where the elevation ranges from 300 m to 1,000 m is denser. Combined w 5 larger wind effects and solar radiation, the ET could increase gradually. The slope and aspect are deemed important parameters accounting for the heteroge 10 terrain effects in the study area. Whereas the aspect mainly affects solar radiatio 11 inclination angle of ground and prevailing wind direction, the increase of slope 12 aggravate these impacts further. Both of them make significant differences in terms of 13 and heat fluxes, thereby resulting in different ET estimation. The following sens 14 analysis picked up a subregion where the slope is above 5° within the study area. In 15 11, it shows the net radiation flux can be significantly affected by the aspect. T 16 radiation flux estimated without considering terrain effects is uniform almost over all a 17 Within the simulation, we found that the maximum net radiation flux appeared on the 18 of 120°, which is very consistent with the sun azimuth (118.7°) associated wi 19 acquisition of the satellite images. This implies that the minimum net radiation appea 20 Fig. 11. The impacts of aspect and slope on R n .
In our case study, the 48 MODIS datasets collected between 2005 and 2006 were used for model validation to avoid cloudy and rainy days. 5 datasets were taken out due to the wheat harvesting events finally. With the 43 datasets, comparisons show a good agreement between the simulated and the measured DAET. Final validation of the SEBTA was confirmed by the consistency index (0.92) and the correlation coefficient (0.87) simultaneously. The selected 24 datasets of MODIS images in 2005 were used for ET analysis assonated with different elevation and land cover. Research findings clearly indicate that the responses of the seasonal ET simulated by the SEBTA to terrain factors are versatile over seasons. Whereas the winter DAET has no any response to elevation changes, the spring and autumn DAETs do have response to elevation changes to some extent and the summer ET has significant response to elevation changes. The DAETs simulated by the SEBTA in different seasons associated with differing elevation zones well reflect the seasonal characteristics in the context of the LULC patterns of the study area. The coupled remote sensing and SABTA therefore becomes a convenient tool for ecohydrological studies. Future work will be placed on the use of the SEBTA to monitor the drought impacts for larger area and conduct the holistic assessment of water resources availability under the impact of climate change. We will choose another site for further verification of this algorithm and validate ET change at different heights over different seasons in the future.