Mapping reference evapotranspiration from meteorological satellite data and applications

Reference evapotranspiration (ETo) is an agrometeorological variable widely used in hydrology and agriculture. The FAO-56 Penman-Monteith combination method (PM method) is a standard for computing ETo for water management. However, this scheme is limited to areas where climatic data with good quality are available. Maps of 10-day averaged ETo at 5 km × 5 km grid spacing for the Taiwan region were produced by multiplying pan evaporation (Epan), derived from ground solar radiation (GSR) retrieved from satellite images using the Heliosat-3 method, by a fixed pan coefficient (Kp). Validation results indicated that the overall mean absolute percentage error (MAPE) and normalized root-mean-square deviation (NRMSD) were 6.2 and 7.7%, respectively, when compared with ETo computed by the PM method using spatially interpolated 10-day averaged daily maximum and minimum temperature datasets and GSR derived from satellite inputs. Land coefficient (KL) values based on the derived ETo estimates and long term latent heat flux measurements, were determined for the following landscapes: Paddy rice (Oryza sativa), subtropical cypress forest (Chamaecyparis obtusa var. formosana and Chamaecyparis formosensis), warm-to-temperate mixed rainforest (Cryptocarya chinensis, Engelhardtia roxburghiana, Tutcheria shinkoensis, and Helicia formosana), and grass marsh (Brachiaria mutica and Phragmites australis). The determined land coefficients are indispensable to scale ETo in estimating regional evapotranspiration. Article history: Received 22 June 2016 Revised 8 November 2016 Accepted 15 November 2016


IntroductIon
Fresh water is vital to terrestrial life forms. However, 96.5% of the water on Earth is salt water found in the oceans, and less than 0.5% of the world's water resources are available to provide the freshwater needs for human activities and terrestrial ecosystems (Shiklomanov and Rodda 2004). Evapotranspiration (ET), which includes water evaporation from land and water surfaces and transpiration from vegetations, play an essential role in energy and mass exchanges between the hydrosphere, atmosphere and biosphere (Brutsaert 1982). On global basis, the mean ET accounts for more than half of the total precipitation fallen on land surfaces (Chahine 1992;Oki and Kanae 2006). ET is also the major irrigation water and precipitation use on agricultural lands (Gowda et al. 2008). ET rates are affected by many factors, such as ground solar radiation (GSR), wind speed, air temperature and humidity, soil water content, vegetation type, growth stage, planting density, and management practices, etc. (Rosenberg et al. 1983;Monteith and Unsworth 1990), which make precise regional evaporative loss quantification rather difficult.
Pressures on water resources are increasing due to population growth, increased living standards, competition for water between agricultural and other sectors, and water pollution. Climate change further aggravates this rapidly decreasing freshwater availability worldwide (Jiménez Cisneros et al. 2014). Although increasing water stocks in natural and artificial reservoirs are helpful in increasing the available freshwater resources for human necessities, Oki and Kanae (2006) pointed out that the flow of water should be the main focus in water resources assessments because water is a naturally circulating resource that is constantly recharged. Reliable ET quantification is very important and valuable, particularly for freshwater resources planning and management.
A substantial number of remote sensing algorithms, varied greatly in main assumptions and input requirements, have been proposed over the past few decades to provide ET mapping at different time and spatial scales (Jiang and Islam 2003;Gowda et al. 2008;Kalma et al. 2008;Li et al. 2009;Jiménez et al. 2011;Liou and Kar 2014). Remote sensed land surface temperature (LST) is the most critical input for many of those proposed algorithms. However, accuracies in retrieving LST under cloudy conditions still needs to be improved (Jiang et al. 2004;Jiménez-Muñoz et al. 2009;Yu et al. 2009;Benmecheta et al. 2013;Li et al. 2013). In addition, many physically-based approaches require auxiliary ground and atmospheric observations that cannot readily be measured by remote sensing techniques (e.g., Anderson et al. 1997;Bastiaanssen et al. 1998;Mecikalski et al. 1999;Su 2002;Allen et al. 2007).
The FAO Irrigation and Drainage Paper No. 56 "Crop Evapotranspiration" (Allen et al. 1998) is the guidelines for most irrigation engineers to compute crop water requirements. It estimates ET from a specific crop field by first calculating reference evapotranspiration (ET o ) using local weather data to represent the weather effect on water consumption. To account for the influences of crop characteristics on ET, a crop specific coefficient (Kc) is used to scale ET o for crops grown under excellently managed and wellwatered standard conditions. Under non-standard conditions, a stress factor (Ks) is applied to account for the water and environmental stresses that may be imposed on the crop field. Therefore, ET = Ks × Kc × ET o .
The FAO Penman-Monteith combination method (PM method) for reference grass surface, Eq. (1) (Allen et al. 1998), is now the standard method most commonly selected for ET o computation (Pereira et al. 2015).
where ET o is reference evapotranspiration (mm day -1 ), R n is net radiation at the crop surface (MJ m -2 day -1 ), G is soil heat flux density (MJ m -2 day -1 ), T is mean daily air temperature at 2 m height (°C), U 2 is wind speed at 2 m height (m s -1 ), e s is saturation vapor pressure (kPa), e a is actual vapor pressure (kPa), Δ is the slope of vapor pressure curve (kPa °C -1 ), c is psychrometric constant (kPa °C -1 ). This method requires daily data of maximum and minimum air temperatures, va-por pressure, net radiation, soil heat flux, and wind speed at 2 m above ground level. However, the required inputs are available only at locations where fully equipped weather stations have been installed. For locations where part of the required inputs are missing, alternative procedures to estimate the missing climatic data are given in Allen et al. (1998), and daily maximum and minimum air temperatures are the minimum datasets required.
Since the evapotranspiration process is determined by the amount of energy available to vaporize water, ground solar radiation (GSR) is the primary driving force of evapotranspiration from ground surfaces, Jacobs et al. (2004), Bois et al. (2008), andde Bruin et al. (2010) applied daily solar radiation estimates derived from GOES and MeteoSat images as an input to the PM method, in conjunction with other ground measurements, to determine ET o . To avoid using remote sensed LST data and minimize the requirements on ground measurements, Syu et al. (2016) used GSR derived from geostationary satellite images by Heliosat-3 method, in conjunction with monthly averaged surface temperatures based on elevation and justified temperature lapse rates (TLR), to generate maps of 10-day averaged daily pan evaporation (E pan ) for Taiwan region at grid spacing of 5 km × 5 km. Instead of computing by the PM method using local weather data as inputs, Allen et al. (1998) pointed out that the reference evapotranspiration (ET o ) for periods of 10 days or longer can also be estimated by multiplying E pan with a properly determined pan coefficient (K p ).
Therefore, the work by Syu et al. (2016) may serve as the basis for an approach to provide ET o using geostationary meteorological satellite images. The ET o concept provides simple, reproducible means to estimate crop ET from weather-based ET o values and helped to enhance the transferability of the crop coefficients from one location to another. Therefore, the reference evapotranspiration concept has now been adopted to estimate ET from many different landscapes, such as forests (Hou et al. 2010;Pereira et al. 2015), wet lands (Drexler et al. 2008;Zhou and Zhou 2009), turf grasses (Jia et al. 2009;Wherley et al. 2015), and urban areas (Pannkuk et al. 2010;Sun et al. 2012;Nouri et al. 2013). Landscape coefficients (K L ) are the term commonly used to substitute crop coefficients (Kc) in estimating evaporation from landscapes other than field crops. It is also defined as the ratio of actual evaporation (ET) over reference evapotranspiration (ET o ) on corresponding landscape, i.e., In essence, K L equals Ks × Kc at crop fields.
The objective of this study is to develop and provide feasible remote sensing methods to estimate ET o , particularly for areas where the PM method is inapplicable due to lacking of good quality local weather data. Specific aims are (1) to derive the pan coefficient (K p ) suitable for generating reference evapotranspiration (ET o ) from pan evaporation (E pan ) derived from geo-stationary satellite imageries by Syu et al. (2016), (2) to validate the accuracies of ET o derived from E pan by comparing with ET o computed by the PM method, and (3) to derive landscape coefficients (K L ) at selected sites where long term latent heat flux measurements data are available.

study region and datasets used
Taiwan is an island located at the intersection of the Pacific Ocean and Asian continent and situated off the southeast coast of China by the Taiwan Strait. The north to south length of the island is about 400 km and about 150 km at its widest. The total area is close to 36000 km 2 . The central mountain range, which has over 100 peaks higher than 3000 m, runs from north to south and divides the island into eastern hills and western plains. About 70% of the island's land area is mountains or slope land while the remaining 30% is plains (Fig. 1). The climate is controlled mainly by orographic relief and by the alternation between the summer southwest monsoon and winter northeast monsoon. To serve as the basis for further analysis and comparisons, the entire study region was divided into 49 × 84 grids with each grid of an area 5 × 5 km 2 . Average elevation above sea level of each grid was derived from a 40-m resolution digital elevation model (DEM).
Images from geostationary satellites MTSAT-1R and MTSAT-2 from 2008 -2013 archived by the Meteorological Satellite Center of CWB were also retrieved and used for derivation of GSR and E pan at each grid of within the study region as described in section 2.2. Daily maximum and minimum temperatures collected from 2008 -2013 at 26 standard weather stations and up to 325 automated temperature/rainfall stations managed by Central Weather Bureau (CWB), Taiwan, were retrieved from databases archived and maintained by the Data Processing Branch of CWB. Before being archived, these temperature data were screened by the Data Processing Branch of CWB to ensure their quality and integrity. Averages of daily maximum (T max ) and minimum (T min ) temperatures for days 1 -10, 11 -20, and 21 -30 (28 or 29 for February) of every month from 2008 -2013 were computed from the retrieved CWB daily temperature datasets and used for computing ET o by the PM method as described in section 2.3.
Latent heat flux data measured at an agricultural site (SK), two forest sites (CLM and LHC), and a wetland site (GDP) were used in this study for derivation of landscape coefficients (K L ). Paddy rice (Oryza sativa) was the crop planted at SK site. The vegetation at the CLM site is mainly a subtropical cypress forest (Chamaecyparis obtusa var. formosana and Chamaecyparis formosensis) with average height about 10 m. Warm-to-temperate mountainous rainforest of mixed evergreens and hardwoods (including Cryptocarya chinensis, Engelhardtia roxburghiana, Tutcheria shinkoensis, and Helicia formosana) with averaged height about 17 m were vegetations at LHC site. The vegetations at GDP site, a wetland situated at the junction of Tamshui River and Jilong River, are mainly comprised by para grass (Brachiaria mutica) and reed (Phragmites australis).
Flux data at all four sites were measured using the Eddy Covariance (EC) method. Brief summaries about the sites and instrumentations used are listed in Table 1 wind vectors between 270 and 30 degrees (c.a. W to NE) were rejected because of inadequate fetch. The flux data quality used for analysis and comparison in this study was controlled by following the standard procedures adopted by FluxNet and AsiaFlux. Therefore, interferences from the surrounding landscapes were minimal.
Standard MODIS (MODerate-resolution Imaging Spectroradiometer) 8-day composite surface reflectance product at 500 m resolution (MOD09A1 collection 051), from 2008 -2013, were downloaded from the Earth Observing System Data Gateway. MODIS Reprojection Tool (MRT) software (Release 4.0) was used to re-project from native Sinusoidal (SIN) projection to the UTM-WGS84 reference system used in this study. Normalized Difference Vegetation Index (NDVI) values were then computed to monitor changes in vegetation at all four flux measurement sites. For each flux tower site, only one pixel with its center closest to the flux tower location was used because it covers the majority of the fetch in most cases.

derivation of remote sensed Pan evaporation (e pan )
Ten-day averaged daily pan evaporation estimates (E pan ) at each 5 × 5 km 2 grid within the studied region were derived from MTSAT images following procedures given by Syu et al. (2016). However, due to changes of satellites in operation, MTSAT-1R and MTSAT-2 images were used for the period from 1 st January 2008 to 30 th June 2010 and from 1 st July 2010 to 31 st December 2013, respectively. Since the stationary orbital positions of the two satellites were different (140.25°E for MTSAT-1R, 145°E for MTSAT-2), the third-order polynomial functions required for representative reference albedos computation under clear and thick cloud conditions at various sun-satellite configurations for each pixel were re-derived for MTSAT-1R images using all the available MTSAT-1R daytime images. For MTSAT-2 images, the required third-order polynomial functions followed those derived by Syu et al. (2016).
Ten-day averaged daily solar irradiance (R s ) maps, from 2008 -2013, were first produced following procedures for image preprocessing, retrieval of GSR, and integration for daily global irradiance as described in Syu et al. (2016). Corresponding E pan maps were then derived by applying the Hansen type empirical model, shown as Eq. (4), given by Syu et al. (2016).
Values of Δ and c were computed using Eqs. (7a) -(7e) using monthly mean temperatures estimated from height above sea level and corresponding regional monthly TLR functions provided by Chiu et al. (2014).

computation of reference evapotranspiration (et o )
Procedures recommended by Allen et al. (1998) were used to estimate the required inputs for Eq. (1) using only the 10-day averaged daily maximum (T max ) and minimum (T min ) air temperatures from CWB ground observations and solar irradiance (R s ) retrieved from MTSAT images as follows.
T dew was estimated using T min and the quantities relate to vapor pressure differences, (e se a ), were calculated as The quantities relate to net radiation (R n ) were calculated as where a is the albedo and assumed to be 0.23, and where v is Stefan-Boltzmann constant (= 4.903 × 10 -9 MJ K -4 m -2 day -1 ), R so is the clear sky radiation estimated from elevation (z) and extraterrestrial radiation (R a ) of the same 10-day period and where z (m) is the elevation above sea level, and where J is Julian day of the year. The slope of saturation vapor pressure curve (Δ), psychrometric constant (c), and latent heat of vaporization (m) were calculated using formulas given by Allen et al. (1998).
As suggested by Allen et al. (1998), soil heat flux was ignored (G = 0), and U 2 was set at constant wind speed of 2 m s -1 for all the ET o computations.

spatial Interpolation of daily Maximum and Minimum temperatures
Temperature has large spatial and temporal variability but measurements at weather stations are only point estimates. Thus, except for computing ET o at measuring weather stations, temperatures have to be interpolated for each 5 × 5 km 2 grid from observations at nearby weather stations. Among numerous spatial interpolation methods, Residual Kriging was considered to be the best method for temperature estimation because of the complex relations between temperature and orography as well as the climate characteristics of Taiwan (Chiu et al. 2009).
The spatial temperature trend with elevation was first removed by subtracting monthly reference temperature calculated based on the site elevation and the regional monthly TLR functions given by Chiu et al. (2014). The residuals of T max and T min were then spatially interpolated by Ordinary Kriging method (Goovaerts 1997) on every 5 km × 5 km grid in the study region for each 10-day period. The semivariogram required for each interpolation was computed individually for a 150 km active lag distance using a uniform lag class interval of 10 km. An isotropic Gaussian model was the preferred model in fitting the semivariogram because it was believed that strong spatial continuities should exist near the origin. However, occasionally, an isotropic Spherical or Exponential model was required for modeling the semivariogram. A search radius of 20 km with a maximum of 16 neighbors and 2 × 2 local grid block Kriging scheme were then applied to conduct the interpolation. The final T max and T min estimates were then restored as the sum of the residual predicted and reference temperature subtracted. All of the semivariance analysis, model fitting, spatial interpolation, and cross validation tests were executed using GS+ (Ver. 5.1.1, Gamma Design Software). In order to provide representative samples for use in section 2.5 within limited time, the spatial interpolation on T max and T min within the study region were conducted only for the middle 10-day period of each month (i.e., day 11 to 20) in the studied six years (2008 -2013), i.e., a total of 144 sets of spatial interpolation were conducted.

derivation and Validation of Pan coefficient (K p )
As shown in Eq. (2), good quality data pairs of E pan (from section 2.2) and ET o (from sections 2.3 and 2.4) are required to derive the pan coefficient (K p ) for converting the pan evaporation rate (E pan ) into reference evaporation rate (ET o ). Spatially interpolated T max and T min with great uncertainty can produce erroneous estimates of ET o and thus severely affect the K p value derived. Therefore, as discussed in section 3.2, the standard deviation of T max and T min estimates must be less than 0.5°C was set as the criteria to ensure qualified pairs of ET o and E pan were selected for K p derivation.
Data pairs belonging to January, April, July, and October (i.e., month in the middle of each season) were used as the training dataset for K p derivation. The data pairs from the rest 8 months were used to test the suitability of derived K p against variations in seasons, geographical regions, and elevations.

derivation of landscape coefficients (K l )
Data pairs of actual evaporative flux (ET) measured by EC method described in section 2.1 and corresponding ET o computed by Eq. (2), using E pan derived from MTSAT images (section 2.2) and K p value to be derived in section 2.5, were assembled for derivation of landscape coefficient (K L ). The K L of each landscape was derived, individually, by fitting a linear model of ET = K L × ET o through the assembled data pairs of each site.
Apparent seasonal changes of NDVI, derived from the MODIS data described in section 2.1, were not observed at LHC, CLM, and GDP sites (to be shown in section 3.3). Significant variations of K L with seasons were also not detected in the preliminary study (data not shown). Therefore, the landscape coefficients for LHC, CLM, and GDP sites were derived using all the data pairs available to each site, individually.
Changes in NDVI at the SK site revealed that the rice plantation followed a two-season cropping system commonly adopted in Taiwan. Other than the changes in leaf area index and vegetation vigilance, the soil water regime also varied significantly at different growth stages due to irrigation practices. During the period from field preparation to active tillering stage (ca. 30 days after transplanting), the paddy fields are generally flooded with water to a depth of 10 -15 cm. In the middle of the growing season (ca. 30 -90 days after transplanting), fields are intermittently irrigated to a water depth of 1 -2 cm at intervals of 7 -10 days. Fields are not irrigated from the yellow-ripening stage (ca. 30 days before harvest) to harvest. During the fallowing period, the fields are usually left dry and covered with stubs, straw residue and grasses.
NDVI of 0.6 can be used to divide the total paddy rice growing period (~120 days) into three stages that match the water management scheme described above (Su and Yang 1999). Therefore, to reduce the influences of changes in vegetation conditions and soil water regime on K L , the available data pairs at the SK site were subdivided into three stages as follows.
(3) Late Stage: from NDVI < 0.6 to start of next flooding for field preparation.

error Analysis
In order to evaluate the derived pan coefficient (K p ) and landscape coefficient (K L ) performances, the Coefficient of Determination (R 2 ), Nash-Sutcliffe efficiency (E, Nash and where O i is the i th measured value, P i is i th simulated value, N is the total number of available data, and O is the average of observations. R 2 indicates how well the data fit a statistical model. An R 2 of 1 indicates that the regression line perfectly fits the data, while an R 2 of 0 indicates that the line does not fit the data at all. E measures the proportion of total variance in the observed data explained by the simulated data and should optimally be one. If E > 0, the modeling results are acceptable; if E ≤ 0, the model should be rejected. d r , with upper and lower bounds of 1 and -1, indicates the magnitude of differences between the model predication and observed data relative to the average error in the observed data. If d r = 1, it is a perfect fit; if d r = 0.5, the magnitude of predication error equals the average error in the observed data. At values of 0 and -0.5, the error magnitudes are twice or fourfold the average-error in the observed data, respectively. Values of dr near -1.0 can mean that the model is fitted very poorly; but, it also can mean that there simply is little observed variability. The accuracy of T max and T min estimated using the Residual Kriging method and ET o derived from multiplying E pan by K p were evaluated by statistical indices including: mean absolute error (MAE) and root-mean-square deviation (RMSD) as well as mean absolute percentage error (MAPE) and normalized root-mean-square deviation (NRMSD).
where N is the number of data compared, X i and O i are estimated and observed values, respectively. For temperatures, the observed T max and T min were those 10-day averages computed from CWB daily temperature datasets. For ET o , the observed data were ET o derived using the FAO PM method at each effective grid. In order to explore the seasonal and geographical effects, the estimated ET o were further evaluated for different seasons (spring: March to May, summer: June to August, autumn: September to November, winter: December to February) and regions of the TLR classification. Small MAPE values indicate the estimates with few errors, while small values of NRMSD indicate more accurate predictions on a point-by point basis.

spatial Interpolation of temperature for et o computation
Before deriving K p for converting E pan into ET o , spatial interpolation of T max and T min between the measuring stations must be conducted to provide the required temperature datasets for ET o computation. Figure 2 depicts the typical omnidirectional semivariogram and fitted model for the residual values after removing the temperature variation trend with elevation, while Fig. 3 shows the corresponding cross validation test results. A Gaussian model was selected to ensure the weights assigned to nearby samples decreased rather slowly with increasing separation distances. The non-zero nugget accounted for the possible short scale variability and uncertainties in temperature measurement. Cross validation tests also indicated good agreements between the estimated and observed T max and T min values (Fig. 3). The MAPE and NRMSD were generally less than 8 and 10%, respectively, for T max and T min estimates.
The number and proximity of nearby samples are major factors related to the uncertainty of the Kriginged estimates. Typical variations of standard deviation (SD) of the T max and T min estimates against the number of nearby samples (n) were shown in Fig. 4. The standard deviations of the estimates decreased with increasing number of nearby samples. At n = 3, the decreasing trend changed from rapid mode to more gradual mode. However, wide ranges of standard deviation, due to variations of samples' proximity, were observed particularly when the numbers of the nearby samples were small. Therefore, the number of nearby samples is not a good index for determining whether or not the T max and T min estimates are suitable for ET o computation. Instead, using SD ≤ 0.5°C as the criterion, only the estimates with small error would be selected. Sensitivity analysis indicated that a 0.5°C error in T max and T min estimates would cause less than 7% error in ET o estimation within the study region.
As shown in Fig. 5, the point measurements of T max and T min could be extended using the Residual Kriging method to cover nearly 90% of the studied region, except at remote areas where temperature data were not available within the search radius (20 km). However, due to sparsely distributed temperature measurement stations, the high SD of T max and T min estimates at coastal and plain areas in the southwest region have resulted in the exclusion of ET o estimation in these areas. Increasing search radius offered little help, if not worse, to expand the area to be covered for ET o computation because the SD were also increased as results of incorporating more far away and irrelevant data. The above discussion illustrates several disadvantages in computing ET o by PM method. For example, its application is restricted only to regions that have good quality data available, and the value of computed ET o may be different using different spatial interpolation techniques and alternative methods to account for the missing data. Considering this, the E pan maps derived from satellite images may be a better data source for ET o estimation because of the ability to provide more complete coverage and consistent computation scheme.

Pan coefficient (K p )
For K p derivation, the distribution of E pan derived from satellite images versus ET o calculated using the PM method is shown in Fig. 6. A 1:1 linear relation, Eq. (10), exists between ET o and E pan . This relation is statistically significant, with an overall coefficient of determination of 0.94.
By definition, the fitted parameter (0.9976) is K p . Performance tests indicated that the derived pan coefficient (K p ) Fig. 2. Semivariograms for residuals of daily maximum (T max ) and minimum (T min ) temperatures for period from 11 to 20 January 2013. The dashed lines represent sample variances. Fig. 3. Cross validation test results for daily maximum (T max ) and minimum (T min ) temperatures for the period from 11 to 20 January 2013 based on the semivarigrams given in Fig. 2.  is applicable in different geographic regions and seasons ( Table 2). The validation test indicated that the ET o derived from E pan by multiplying the derived K p was good surrogate of the ET o computed by FAO PM method (Table 3). The surrogates were effective not only at all four geographical regions and seasons, but also at various heights in mountainous areas. Overall MAPE and NRMSD were only 6.2 and 7.7%, respectively. This is a big breakthrough in providing ET o for areas where even the minimum weather datasets (T max and T min ) are not available, either due to lack of nearby temperature measurement stations or high uncertainties in the interpolated temperatures. ET o maps for the study region can now be easily generated from E pan maps that are derivable from meteorological satellite images.

et Fluxes at Four long term eddy covariance Flux Measurement sites
Temporal variations in vapor fluxes (ET) measured using the EC method and NDVI values derived from MODIS observations at the four long term ET flux observation sites are shown in Fig. 7. Corresponding reference evapotranspiration (ET o ), generated by multiplying E pan maps with the K p derived in section 3.3, for these four sites are also plotted in Fig. 7. As expected, the ET o was higher than actual ET at all four sites.
At the LHC, CLM, and GDP sites, the latent heat fluxes increased gradually from January, reached peak values in summer, and decreased gradually toward December. However, the NDVIs at LHC and CLM sites remained relatively stable (~0.8) year round, indicating the vegetations maintained fully covered and experienced no severe stresses during the studied 3 years. The NDVIs at GDP site were stably maintained at about 0.5 year round, which also indicated no major changes in percent cover and vigilance in the vegetation. The less dense wet land vegetation at the GDP site may have resulted in the relative low NDVI values compared with those at the LHC and CLM sites. Since no apparent changes in vegetation conditions were observed, the amount of energy available for evapotranspiration might be the major factor dominating the seasonal changes of the observed ET fluxes at LHC, CLM, and GDP sites.
At the SK site temporal variations in NDVI values revealed that the vegetation cover and growth conditions were under a generally adapted two-season rice cropping system in Taiwan (Fig. 7). The ET measured was modified by weather conditions at the time (as shown by the changes of ET o ). Although parts of the flux data were missing due to inadequate fetch as the result of changes in wind direction, the ET fluxes still showed an observable general trend in the first crop season (i.e., the first half of the year). The ET increased gradually after transplanting (early February), reached a plateau in between maximum tillering (late March) to heading (late April), and then decreased gradually toward harvesting (mid June). A similar trend was not observable in the second cropping season due to huge gaps in missing data. However, a decreasing trend in ET from late September was still observable.

landscape coefficients (K l ) at test sites
To find K L , a linear model of ET = K L × ET o was fitted through all the available data pairs of ET and ET o at LHC, CLM, and GDP sites, individually (Fig. 8). The statistical measures for performance evaluation of the fitted linear relations were shown in Table 4. All fitted models yielded an acceptable agreement between the simulated and observed ET fluxes.
At the SK site, data pairs of ET and ET o were first divided into three stages, as described in section 2.6, to account for differences in soil water regime and vegetation conditions and then fitted with linear relations, individually (Fig. 8). The statistical measures for performance evaluation of the fitted linear relations are also shown in Table 4. The fitted model for the middle stage was acceptable. However, the less satisfied statistical measures for the early and late stages at SK site indicated that more ET flux measurements are required to improve the accuracy of the corresponding K L .

comments on the derived landscape coefficients (K l )
Although different in vegetation type and canopy structure, the derived K L values (~0.6) were comparable at LHC and CLM sites. At CLM site, it is a subtropical montane cloud forest. Other than vegetation transpiration and soil evaporation, direct evaporation of fogs and rainfalls intercepted by canopy also contributes to the ET measured. Therefore, the K L at CLM is slightly higher than that at LHC. However, constant K L throughout year at LHC and CLM sites indicated that the forests at these two sites have the ability to extract water from subsoil even when seasonal rainfall is reduced. However, the less than 1 value indicated that the evapotranspiration from these forests was limited by the boundary layer resistance and resistances due to xylem flow and stomata opening.
Although soil is commonly saturated at the GDP site (wet land), a K L value around 0.66 indicated that the tall and dense vegetation of the para grass and reed also increased resistances of transpiration from vegetation and evaporation from the wet soil surface.
For paddy rice fields (SK site), a 5 -10 cm of water inundation at early stage (starting from field preparation stage to active tillering stage) provided an open water surface which yielded a K L around 1. In the middle stage (from maximum tillering stage to milking stage), soil water condition was controlled by the intermittent irrigation water management scheme commonly adopted in Taiwan. Therefore, the K L value (~0.75) was smaller than those commonly reported values (~1) for actively growing crops with ample water supply (Allen et al. 1998). At the late stage (yellowing stage to harvesting and fallowing before field preparation for the next crop season), the lower K L (~0.57) resulted from the limited soil water availability and increased resistance due to wilting vegetation and plant residues left on the ground. Assuming a crop coefficient (Kc) of 1 if under well watered condition, the combined stress factor (Ks) for the middle and late growth stages were about 0.75 and 0.57, respectively, by definition (K L = Ks × Kc). These Ks values were within the range of those commonly reported (Allen et al. 1998).  An objective of water resource management is to obtain accurate estimates on water loss from crop evapotranspiration and distribute the available renewable freshwater resources accordingly. Since the derived K L for paddy rice used two years of field measurements, it may reflect the actual field conditions under the intermittent irrigation practices currently adopted by Taiwan rice farmers and are thus helpful for water resource management and paddy field irrigation scheduling.

conclusIons
The above results confirmed that maps of 10-day averaged ET o can be derived from meteorological satellite images. The estimated ET o was in good agreement with ET o derived from the FAO PM method. The proposed scheme simplified the algorithms for ET o derivation and removed the need for in situ climatic data as inputs for ET o computation. This is a particularly useful approach for locations where even the minimum datasets required for ET o computation by PM equations are not available.
Derivation of landscape coefficients (K L ) at four different landscapes has been used to illustrate the potential usages of the ET o map generated. The proposed scheme can easily be applied to other locations where long term latent heat flux measurements data are available. More flux measurements on the same landscapes but at other locations, or on different landscapes will provide data for validation and developing a more comprehensive K L database. Once the database has been established, more accurate estimation of actual ET from the region can be accomplished.