Long time series of daily evapotranspiration in China based on the SEBAL model and multisource images and validation

Satellite observations of evapotranspiration (ET) have been widely used for water resources management in China. An accurate ET product with a high spatiotemporal resolution is required for research on drought stress and water resources management. However, such a product is currently lacking. Moreover, the performances of different ET estimation algorithms for China have not been clearly studied, especially under different environmental conditions. Therefore, the aims of this study were as follows: (1) to use multisource images to generate a long time series (2001-2018) daily ET product with a spatial resolution of 1 km × 1 km based 15 on the Surface Energy Balance Algorithm for Land (SEBAL); (2) to comprehensively evaluate the performance of the SEBAL ET in China using flux observational data and hydrological observational data; (3) to compare the performance of the SEBAL ET with the MOD16 ET product at the point-scale and basin-scale under different environmental conditions in China. At the point-scale, both the models performed best in the conditions of forest cover, subtropical zones, hilly terrain, or summer, respectively, and SEBAL performed better in most conditions. In general, the accuracy of the SEBAL ET (rRMSE = 44.91%) was slightly higher 20 than that of the MOD16 ET (rRMSE = 48.72%). In the basin-scale validation, both the models performed better than in the pointscale validation, with SEBAL obtaining superior results (rRMSE = 19.15%) to MOD16 (rRMSE = 33.62%). Additionally, both the models showed a negative bias, with the bias of the MOD16 ET being higher than that of the SEBAL ET. In the daily-scale validation, the SEBAL ET product showed an RMSE of 0.92 mm/d and an r-value of 0.79. In general, the SEBAL ET product can be used for the qualitative analysis and most quantitative analysis of regional ET. SEBAL ET product is freely available at https://doi.org/ 25 10.5281/zenodo.4218413 (Cheng, 2020)SEBAL ET product is freely available at https://doi.org/ 10.5281/zenodo.4218413 (Cheng, 2020). The results of this study can provide a reference for the application of remotely sensed ET products and the improvement of satellite ET observation algorithms. Keyword: evapotranspiration; SEBAL; MOD16; accuracy validation; multiscale 30


Introduction
Evapotranspiration (ET) is the process of transferring surface water to the atmosphere, including soil evaporation and vegetation transpiration (Wang and Dickinson, 2012). This process is a key node linking surface water and energy balance. In the process of water balance, ET represents the consumption of surface water resources, and in the process of energy balance, the energy consumed by ET is called the latent heat flux (λET, W/m 2 , where λ is the latent heat vaporization), which is an important energy component 35 (Helbig et al., 2020;Zhao et al., 2019). Approximately 60% of global precipitation ultimately returns to the atmosphere through evapotranspiration (Wang and Dickinson, 2012). Therefore, accurately quantifying the ET of different land cover types is necessary to better understand changes in regional water resources. However, the methods for the estimation of ET based on point-scale or small-area-scale analysis, such as lysimeter and eddy covariance, cannot meet the requirement of global climate change research and regional water resource management (Li et al., 2018). Since the United States successfully launched the first meteorological 40 satellite in the 1960s, hydrological remote sensing (RS) applications have developed rapidly and have led to huge breakthroughs (Karimi and Bastiaanssen, 2015). Remote sensing technology with a high spatiotemporal continuity provides an effective means for regional ET estimation.
Satellite remote sensing provides a reliable direct estimation of ground parameters; however, it cannot measure ET directly (Wang 45 and Dickinson, 2012). Therefore, several RS-based algorithms for the estimation of ET have been proposed and reviewed (Pôças et al., 2020; Senay et al., 2020; Wang and Dickinson, 2012). These models can be divided into three types according to their mechanism: those based on surface energy balance residual (SEBR), those based on semi-empirical formulas (SEFs) and statistic methods.
SEBR-based models can be further divided into one-source models and two-source models (Wang and Dickinson, 2012). Onesource models do not distinguish vegetation from bare soil and regard the land surface as a system that exchanges energy and water 50 with the atmosphere. Examples of one-source models include the Surface Energy Balance Index (S-SEBI) (Roerink et al., 2000), the Surface Energy Balance System (SEBS) (Su, 1999), and the Surface Energy Balance Algorithm for Land (SEBAL) (Bastiaanssen et al., 1998a;Bastiaanssen et al., 1998b). These models have a theoretical basis, a simple principle, strong portability, and have been widely used ( using an empirical formula and a vegetation index obtained from remote sensing data to divide the land surface into vegetation and bare soil in each single pixel. Examples of two-source models include the Two-source Energy Balance (TSEB) (Kustas et al., 2003), Two-source Trapezoid Model for Evapotranspiration (TTME) (Long and Singh, 2012), and Hybrid Dual-source Scheme and Trapezoid Framework-based Evapotranspiration Model (HTEM) (Yang and Shang, 2013). Compared to one-source models, twosource models have a superior theoretical mechanism. SEF-based models using traditional semi-empirical formulas calculate λET data or remote sensing data) and in situ ET are also be widely used (Mosre and Suárez, 2021; Yamaç and Todorovic, 2020).
Since ET plays a critical role in the study of hydrology and ecology, ET products with a high spatiotemporal resolution are required.
Therefore, a growing number of ET products have been generated to meet research needs. These include MOD16, which is generated by NASA based on the Penman-Monteith algorithm and has a spatial resolution of 500 m × 500 m and a temporal resolution of 70 eight days (Mu et al., 2007;Mu et al., 2011). The GLEAM daily ET product with a spatial resolution of 0.25° × 0.25° has been generated by the University of Bristol, UK, based on the Priestley-Taylor method (Miralles et al., 2010). Additionally, Chen generated long time series daily ET datasets with a spatial resolution of 0.1° × 0.1° based on the SEBS algorithm (Chen et al., 2014a; Chen, 2019). However, there are few ET products which simultaneously meet the current research needs in terms of temporal and spatial resolution. Therefore, generating a kilometer-level daily ET product which can minimize the influence of mixed pixels is 75 critical.
Water resources management is essential for China as it has an unbalanced spatial and temporal distribution of water resources. ET, as a crucial component of terrestrial water cycle, is critical for understanding water resources budget in China. Therefore, a spatiotemporal continuous ET data is needed. Several studies evaluated the performance of various remote sensing-based algorithm

95
In order to improve ET products in China and better understand the performance of RS-based ET estimation models in China, in this paper, we aim to (1) generate a long time series daily ET product with a spatial resolution of 1 km × 1 km based on the SEBAL model and multisource remote sensing images, (2) validate the accuracy of the generated ET product in China based on flux tower observational data and hydrological data, and (3) compare the performance of the generated ET product with MOD16 datasets in China under different environmental conditions. 100

Generation of long time series daily ET product
In this study, a long time series daily ET product was generated based on SEBAL, which is a widely used one-source model (Gobbo et where R n is the net radiation flux, H is the sensible heat flux, and G is the soil heat flux (the unit of all three parameters is W/m 2 ). 6 In this paper, MODIS data (MCD43 surface albedo, MOD11 surface temperature (daytime), MOD13 NDVI) and meteorological data (air temperature) from the Global Modeling and Assimilation Office (GMAO) were used as input for surface parameterization (R n , G and H). The details of generating SEBAL ET can be referred to Appendix.

125
The spatial and temporal resolutions of the MCD43 surface albedo and the MOD11 daytime surface temperature are 1 day and 1 km × 1 km, while those of MOD13 NDVI are 16 days and 500 m × 500 m. In this study, MOD13 was resampled to 1 km × 1 km and processed by smoothing and gap-filling from time series to daily data (Vuolo et al., 2017). It should be noted that there are several missing or unreliable pixels in MODIS images which may cause by cloud or other reasons, these pixels were marked in quality control (QC) files. In this study, these anomalous pixels of MODIS dataset (MOD11, MOD13 and MCD43) were filled 130 referred to previous studies, the rules as follows: (1) the value of anomalous pixel will be computed by liner interpolation of the nearest reliable value after it or prior it; (2) if the anomalous pixel was found in the first or last day, it will be replaced by the closest reliable date value. The more details can be referred to the study of Mu et al. (2011) and Zhao et al. (2005). Land surface temperature is the crucial parameter for SEBAL ET generating, the Appendix 2 shows the ratio of interpolated pixels of land surface temperature (MOD11) data (31% ± 11%). The spatial and temporal resolutions of GMAO air temperature are 1 day and 0.25° × 0.25°, 135 respectively. The coarse-resolution GMAO data were non-linearly interpolated to a spatial resolution of 1 km × 1 km based on the four GMAO pixels surrounding a given pixel (Zhao et al., 2005). The spatial and temporal resolutions of wind speed are 1 day and 1 km × 1 km (China Meteorological Data Network, http://data.cma.cn). The final generated daily ET product has a spatial resolution of 1 km × 1 km and covers the period 2001 to 2018.

Point-scale validation
The eddy covariance method measures λET from the covariance between moisture fluxes and vertical wind velocity using rapid response sensors at frequencies typically equal to or greater than 10 Hz (Wang and Dickinson, 2012); it is regarded as the most effective method for the estimation of ET and has been widely used (Wang and Dickinson, 2012). In this study, eddy covariance tower-measured daily flux data from eight stations in China (Table 1) obtained in 2003-2010 were used to validate the modeled ET 150 (ET SEBAL , ET MOD ). The latent heat flux (λET) observed at the flux towers was converted into the observed ET (ET flux ). It should be noted that the energy balance closure issue, which indicates the sum of sensible heat (H), latent heat (λET) and soil heat flux (G) is not equal to net radiation (Rn), was often found in eddy covariance system. Therefore, the eddy covariance system measured value should be filtered and corrected. First, the data with Energy Balance Closure Ratio (ECR, Eq. 2) less than 80% were not selected for validation (Wang et al., 2019), and then, the remaining data with ECR more than 80% were corrected by using Bowen Ratio 155 energy balance correction (Eq. 3) (Chen et al., 2014b).
where Rn, G, H and λET are all eddy covariance system measured value, and λET cor is corrected value. To ensure a reliable evaluation, the pixel value where the flux tower located (area of 1 km × 1 km) was extracted for comparison with measured value 160 (Velpuri et al., 2013). The water demand is different under different environmental conditions. Therefore, it is necessary to understand the accuracy performance of ET products for different vegetation types when a single ET product is not comprehensive (Velpuri et al., 2013). In order to better understand the influence of different environmental conditions on the accuracy of the model, the modeled ET were validated for different terrain, climate zones, land cover types, and seasons. Additionally, MOD16 data were resampled to a spatial resolution of 1 km × 1 km and daily ET SEBAL and daily ET flux data were accumulated to eight days to match 165 the MOD16 data. ET SEBAL was validated at the daily scale and 8-day scale, respectively.

Regional-scale validation
Furthermore, the regional (basin-scale) ET was calculated using the water balance method (Eq. 34) to validate the modeled ET at 170 the regional scale.
where P (unit: mm) is the annual precipitation in the basin; Q (unit: mm) is the annual runoff in the basin, which includes surface

Accuracy estimation
The modeled ET values were compared with the observed ET (ET flux , ET WB ) to evaluate the performance of ET SEBAL and ET MOD , 180 respectively. The correlation coefficient (r), root-mean-square error (RMSE), relative root-mean-square error (rRMSE), and mean bias error (MBE) were selected to quantify the accuracy of the modeled ET. The equations for these parameters are shown below: where ET M is the modeled ET (ET SEBAL and ET MOD ); ET Ob is the observed ET (ET flux and ET WB ); and n is the number of samples.
r was calculated to evaluate the linear relationship between the modeled and observed ET; higher r-values mean a higher correlation.
RMSE and rRMSE were used to evaluate the performance of the model: smaller RMSE and rRMSE mean a higher accuracy. rRMSE

MOD16 data
The MOD16 ET product is one of widely used evapotranspiration dataset for water resources management and global change study, 195 which also performs accurate to some extent (

Auxiliary data
In order to ensure the objectivity of the comparison between the SEBAL and P-M models, MODIS satellite data were selected as  Armonk, New York, USA) were used for numerical calculation and analysis.

Validation of daily SEBAL ET at the point-scale using flux tower observations
The validation results for the daily SEBAL ET (ET SEBAL ) obtained using flux tower observational data are shown in Fig. 3. 225 Compared to ET flux , ET SEBAL showed a good performance in China; the two data showed a high consistency, with an r-value of 0.79 with 9896 samples. However, the bias of SEBAL was relatively high; the RMSE and rRMSE were 0.92 mm/d and 42.04%, respectively. As shown in the scatter diagrams in Fig. 3

Performance of the RS-based model for different land cover types
The validation results for different land cover types are shown in Fig. 5 indicated that both the ET models underestimated ET over all land cover types. Overall, the accuracy of SEBAL was higher than that of MOD16.

Performance of the RS-based model for different climate zones
The validation results for different climate zones are shown in Fig. 6. The results show that the r-value varied from 0.68 to 0.90 for SEBAL and varied from 0.61 to 0.94 for MOD16. Climate zones were found to influence the accuracy of the RS-based models. In tropical zones, both the two models showed poor accuracy, with RMSEs of 10

Performance of the RS-based model over different terrain types 280
The validation results for different terrain types are shown in Fig. 7. The results indicate that both models showed a negative bias (negative MBE) for all terrain types except mountainous areas, for which both models overestimated, with MBEs of 1. 19

Performance of the RS-based model in different seasons
The validation results for different seasons are shown in Fig. 8. SEBAL showed a negative bias in summer, autumn, and winter,

Validation at the basin-scale using the water balance method
Additionally, validation using hydrological data was performed to investigate the performance of the RS-based models at the basinscale. The results (Fig. 10) showed that both the models had a negative bias, with an MBE of -24.

Comparison of the spatial distribution of ET between SEBAL and MOD16
Regarding the modeled spatial distribution of ET, both the SEBAL and MOD16 models showed that the annual average (2001-335 2018) ET in China increased from the northwest to the southeast (Fig. 11(a), (b)). Fig. 11(d). The annual ET of SEBAL varied from 0 to 1600 mm in space, with a mean value of 482.27±192.31 mm, while that of MOD16 varied from 0 to 1200 mm, with a mean value of 359.61±59.52 mm. In general, compared to the ET value estimated using MOD16 and SEBAL, the ET value estimated using SEBAL was higher and showed a greater spatial difference of ET in China. For 84.07% of the total area of China, the annual ET estimated by SEBAL was higher than that estimated by MOD16; for 14.07% of the total area of China, the difference was more 340 than two times-these areas are mainly distributed in Southern China, where ET is relatively high, and the difference reaches more than 600 mm in some places. Only in 15.93% of the total area of the country was the annual ET estimated by SEBAL lower than that estimated by MOD16; these areas are mainly distributed in Northwest China, where ET is relatively low (Fig. 11(c), (e)).
Regarding the distribution of ET SEBAL , a bi-modal curve with the boundary of ~500 mm was shown in the Fig. 11d, it was likely contributed by the misestimation of part of regions. The ET SEBAL map was divided to two parts with 500 mm as threshold value, the 345 part of ET SEBAL below 500 mm was distributed in the Northwest China (Fig. 11f), whereas the part of ET SEBAL over 500 mm was distributed in the southeast (Fig. 11g). It should be noted that the vegetation cover in northwest of China are mainly grassland and a small part of cropland (Fig. 11h), and the ET SEBAL of grassland and cropland was underestimated by SEBAL model (Section 3.1).
In contrast, the ET SEBAL showed slightly overestimation of forest which is the main land cover types in southeast of China. Therefore, the distributed ET SEBAL around ~500 mm was underestimated or overestimated, and thus formed the bi-modal curve. 350

Summary of validation results and comparison with other studies
The ET SEBAL showed a relatively good performance in China as a whole, with an average r-value of 0.79 and an average RMSE of 0.92 mm/d. These results are close to those obtained in other studies. Rahimzadegan and Janani (2019) used SEBAL to estimate the actual ET of pistachio in Semnan, Iran, and found that the modeled value had a high consistency with the in-situ measured value (r 360 = 0.80); this value was slightly lower than the cropland validation obtained in the present study (r = 0.88, daily-scale). This difference is mainly due to differences in the validation method between these two studies. Rahimzadegan

Inaccuracy of input data
Both SEBAL and MOD16 used MODIS data as the main input images (e.g., MCD43 surface albedo, MOD13 NDVI, MOD11 surface temperature). However, the accuracy of these data is uncertain to some extent (Ramoelo et al., 2014). For instance, surface albedo is a critical radiative parameter, however, the complex algorithm-led remote sensing-based albedo products can contain errors respectively. In general, original MODIS data performed errors to some extent.
Additionally, gap-filling of missing or unreliable MODIS data may causes the errors to some extent. For example, spring and summer have the relatively frequent precipitation, which causes more unreliable pixels due to the cloud, and these pixels value were finally replaced by gap-filling of nearest date pixel value, therefore, the modeled ET value of these pixels was close to that of nearest date without precipitation. In fact, due to the high air humidity in rainy day, the evaporation and transpiration are relatively less than 395 that of nearest date (Ferreira and Cunha, 2020; Li et al., 2016). Moreover, it should be noted, due to the decrease of surface available radiation energy which was caused by cloud cover, the ET (both actual and modeled value) is also less than that of nearest date (Cheng et al., 2020). This may explain the reason of obvious overestimation at lower ET rates in spring, summer and other pixels affected by cloud. Furthermore, a relatively high bias of SEBAL ET was found in winter, the rRMSE reached 66.92% (the highest value among all situations). Due to the ice and snow cover caused by the frequent snowfall and low temperature in winter, which 400 will affect the remotely sensed information to a great extent, e.g., reflectance (Casey et al., 2017), and further affect the ET estimation.
Moreover, the underestimation was found at higher ET rates in the most of situations as shown in Figs

Errors in flux tower measurements
The eddy covariance system (flux tower observations) is the most commonly used observation system to calculate and analyze the energy and mass exchange between the surface and atmosphere (Wang and Dickinson, 2012). However, the typical error of ET estimation based on the eddy covariance system is about 5-20% (Culf et al., 2008;Vickers et al., 2010). Also, the eddy covariance system generally has an energy balance non-closure issue that, the sum of the soil heat flux, sensible heat flux and latent heat flux was found to be less than net radiation in most cases (Mu et al., 2011;Wilson et al., 2002). Recently, it was found that the nonclosure issue of the energy balance was explained by the energy fluxes from secondary circulations and larger eddies that cannot be captured by EC measurement at a single station (Foken et al., 2011). In this study, the Bowen ratio method (Eq. 3), which assuming 415 that the residual of the energy balance is attributed to sensible and latent heat flux and assigning the missing energy flux to them (Song et al., 2016;Wang et al., 2019), was used to enforce energy closure. Actually, this assumption is not very correct, which generally led the sensible and latent heat flux overestimation (Song et al., 2016), which may could explain that the SEBAL ET was generally underestimated when compared to flux tower observed ET (Fig. 9). The same issue was found in regional-scale validation, due to the ignoring of ΔS in the water balance computation process (although it's small), which could lead the regional ET 420 overestimation and further caused SEBAL ET underestimation in regional validation (Fig. 10).
Additionally, the mismatch of flux tower footprint and spatial resolution of SEBAL ET will causes errors as well. Generally, the footprint of flux tower varied from hundreds square meters to several square kilometers which determined by the height of the observation instrument, the intensity of the turbulence, terrain, environment and vegetation status (Chen et al., 2012;Damm et al., 2020;Schmid, 1994). Moreover, a footprint probability distribution function (PDF) could characterize the footprint in a fine spatial 425 resolution (Wang et al., 2019), but it may not suit for the coarse resolution in this study (kilometer-scale). In this study, the 1 km × 1 km area of pixel was used for matching the footprint of flux tower which was referred to the study of Velpuri et al. (2013), however, the footprint is not stable but varied with environment changed, e.g., vegetation height. Chen et al. (2012) reported that forest footprint has clear difference with grassland, the footprint of forest is much larger which reach kilometer-scale. In fact, forest footprint may more matching with the spatial resolution in this study. Therefore, it may explain that the SEBAL ET has the greatest 430 performance in forest but worst performance in grassland. Compared to the study of Velpuri et al. (2013), the grassland also showed the worst remote sensing ET estimation in US when using flux tower data for validation at a kilometer-scale.
In general, more ET sites could improve the reliability of the validation process. However, limited number of EC towers in this study also caused uncertainties in SEBAL ET validation. Although the validation of SEBAL was conducted in different situations (e.g., different climate zones and land cover types), and it should be noted that several situations only have one sites can be used, 435 e.g., cropland. Therefore, only one climate zone of cropland ET can be validated, and other situations were not considered. Moreover, in this study, the different classes of land cover types, climate zones, elevation and seasons were considered for SEBAL validation  The P-M algorithm calculates λET using a semi-empirical formula and A, and therefore avoids the direct calculation of H.
However, several studies have indicated that MOST has an error of 10-20% for the estimation of the boundary layer thickness (Foken, 2006;Högström and Bergström, 1996). Therefore, MOST is also a source of error in SEBAL. Due to the complexity of the sensible heat flux, SEBAL makes several assumptions to estimate H, which may introduce error into the ET estimation (Zheng et 465 al., 2016).
Additionally, the selection of the hot/cold pixel depends on the domain size (the actual size of the modeling domain/satellite imagery being used). For instance, the basin-scale selection of the hot/cold pixel with diverse vegetation cover and single vegetation cover, respectively, will lead to different results for dT. Theoretically, with the domain size increasing, there is a possible tendency that T s_hot increasing and T s_cold decreasing. For example, if T s_cold remains invariant and T s_hot increases under the condition of domain 470 size increasing, the H estimates will decrease and λET estimates could thus increase. In the study of Long et al. (2011), the results showed that a 2 K increase in T s_hot will result in a 9.3% increase but a 9.1% decrease in a and b, respectively, and further caused an 11.8% mean decrease in H. Recently, the study of Saboori et al. (2021) reported that the cold pixel performed more stable than hot pixel in time series, especially in winter, the hot pixel was highly varied may due to the similarity of NDVI over space, it could further explain the poor performance of SEBAL ET in winter. Seguin et al. (1999) demonstrated the method of hot/cold pixel 475 selection for the estimation of H generally has an accuracy of ∼50 W/m 2 . In this study, the ET was computed in the domain size of 1200 km × 1200 km which is a relatively large area. The performance of different domain sizes in ET computation were compared ( Fig. 12) and resulted in better performance are generally found in smaller area, which may due to the relatively limited spatial heterogeneity (Long et al., 2011). However, a larger domain size may have faster computational efficiency of ET in regional scale.
The trade-off between efficiency and accuracy (i.e., most suitable domain size) need be further studied. Overall, the domain size 480 employed in this study (1200 km × 1200 km) performed an acceptable accuracy. Although several algorithms have been proposed that use other methods to avoid the error caused by the selection of the hot/cold pixel , such as the SEBS (Su, 1999), these replaced the selection of the hot/cold pixel with the fitting of dry and wet edges. However, no evidence has been found that the method of

Data availability
The dataset that was generated using SEBAL with a spatial resolution of 1 km and a temporal resolution of 1 day can be used for various types of geoscientific studies, especially for global change, water resources management, agricultural drought monitoring, etc. The evapotranspiration (ET) dataset for China is distributed under a Creative Commons Attribution 4.0 International license. 500 The dataset is named SEBAL evapotranspiration in China (SEBAL ET) and consists of 18 years of data. More information and data are freely available from the Zenodo repository at https://doi.org/ 10.5281/zenodo.4243988 (Cheng, 2020).

Conclusions
In this study, we generated a long time series (2001-2018) ET product based on SEBAL and multisource images. We further conducted a comprehensive validation of the product and compared its performance under different environmental conditions in 505 China with the performance of the ET estimated using MOD16 data. The conclusions are as follows: (1) The ET product generated using SEBAL showed a good performance in China. Compared to flux tower observational data, the r-value of the SEBAL ET reached 0.79 for 9896 samples; the RMSE was 0.92 mm/d and the rRMSE was 42.04%. SEBAL underestimated ET as whole, with an MBE of -0.15 mm/d. The SEBAL ET product can adequately represent the actual ET and can 510 be used in research on water resources management, drought monitoring, ecological change, etc.
(3) Based on flux tower observational data and hydrological observational data, the ET estimated by SEBAL and MOD16 were validated at the point-scale and basin-scale. The results showed that, at the point-scale, the accuracy of SEBAL was 7.77 mm/8 d 525 for the RMSE, 44.91% for the rRMSE, and 0.85 for the r-value, and the accuracy of MOD16 was 8.43 mm/8 d for the RMSE, 48.72% for the rRMSE, and 0.83 for the r-value. At the basin-scale, the accuracy of SEBAL was 48.99 mm/year for the RMSE, 13.57% for the rRMSE, and 0.98 for the r-value. In general, SEBAL performed slightly better than MOD16 at the point-scale, while SEABAL had a larger accuracy advantage at the basin-scale.

530
(4) Overall, the SEBAL ET is higher than the MOD16 ET: for 84.07% of the total area of China, the SEBAL ET showed higher values. Additionally, the SEBAL ET is closer to the in-situ measured ET in most conditions, while the MOD16 ET performed better only in temperate zones, mountain areas, or the spring season. In general, the two models both have a good performance and can be used in the qualitative analysis and most quantitative analysis of regional ET. Furthermore, the combination of the two models can improve the overall ET estimation accuracy for use in applications with higher accuracy requirements. 535 Compared to the widely used MOD16 ET data, the SEBAL ET product showed a higher accuracy and temporal resolution. However, it still has a daily error of 42.04% (0.92 mm/d) at the point-scale and a yearly error of 13.57% (48.99 mm/year) at the basin-scale.
Therefore, the improvement of the SEBAL algorithm will be the focus of follow-up research. Moreover, the 1 km spatial resolution of the SEBAL ET product cannot meet the requirements of more detailed research. Due to the difficulty of simultaneously satisfying 540 the requirements for the spatial and temporal resolutions of remote sensing data, the fusion of multiple sources of remote sensing data may be the most effective way to improve the spatiotemporal resolution of daily ET products.

Acknowledgements:
This study was financially supported by the National

Conflict of interest statement:
The authors declare no conflict of interest. 550

Appendix 1：Description of generating SEBAL ET in details
The SEBAL model calculates the instantaneous λET of the satellite transit time as a residual based on the surface energy balance equation (Eq. 1) as follows: where R n is the net radiation flux, H is the sensible heat flux, and G is the soil heat flux (the unit of all three parameters is W/m 2 ).
In this paper, MODIS data (MCD43 surface albedo, MOD11 surface temperature, MOD13 NDVI) and meteorological data (air temperature) from the Global Modeling and Assimilation Office (GMAO) were used as input for surface parameterization (R n , G and H). The equations for R n are shown in Eqs. 2-5 below: where α is the surface albedo obtained from the MCD43 data; R s_down , R l_up , and R l_down are the downwelling shortwave radiation, downwelling longwave radiation, and upwelling longwave radiation, respectively (the unit of all three parameters is W/m 2 ). R s_down can be calculated using the Julian day (used to estimate the astronomical distance between the sun and earth), elevation (used to estimate atmospheric emissivity), and solar zenith angle at the time of satellite transit. R l_up and R l_down can be calculated using the surface temperature (MOD11), NDVI (MOD13, used to estimate surface emissivity) and air temperature (GMAO data), and 565 atmospheric emissivity based on the Stefan-Boltzmann law. The equations for R s_down , R l_up , and R l_down are given in Eqs. 3-5. where G sc is the solar constant (1376 W/m 2 ); θ is the solar zenith angle; τ sw is the atmospheric transmittance (Eq. 6) (Tasumi, 2000); 570 d r is the astronomical distance between the sun and earth (Eq. 7) (Bastiaanssen et al., 1998a); ε a and ε are the atmospheric emissivity (Eq. 8) (Bastiaanssen et al., 1998a) and surface emissivity (obtained from MOD11), respectively; σ is the Stefan-Boltzmann constant (5.67 × 10 -8 W/m 2 K 4 ); and T a and T s are the air temperature (unit: K; obtained from GMAO data) and surface temperature (unit: K, obtained from MOD11), respectively. where ρ air (unit: kg/m 3 ) is the air density (Eq. 11) (Smith et al., 1991); C p (unit: J/(kg×K)) is the specific heat of air at constant 585 pressure; dT (unit: K) is the difference between the aerodynamic surface temperature (T z0h ; unit: K) and the reference height temperature (T a , unit: K); and r a is the aerodynamic resistance (unit: s/m) (Eq. 12).  (13) where U r is the wind speed at height Z r , which can be calculated from the wind speed monitored by weather stations (U w , Eq. 14); Z r is 200 m in this study (Zeng et al., 2008); and z 0m is the surface roughness (unit: m, Eq. 15) ( Moran and Jackson, 1991 However, since it is difficult to calculate dT directly, the model assumes that there is a linear relationship between surface temperature (T s , unit: K) and dT, as shown in Eq. 16: SEBAL solves the values of a and b by selecting the hot and cold pixels; it assumes that the hot pixel represent pixel of dry cropland 600 with low vegetation covers or bare surfaces or saline alkali land covered by vegetation with zero λET (H = R n -G), and the cold pixel represent pixel with sufficient water supply, lush vegetation, and low temperature, with an H of zero (λET = R n -G). In this study, the hot and cold pixels were selected automatedly by following the certain rules (Long et al., 2011): for hot pixel, the pixels with high Ts (top 10%) and low NDVI (top 10%) in the image were selected first, and further to select the pixels with the land covers of cropland or bare surfaces (according to MOD12 land use product) from the pixels selected in last step, finally, the pixel 605 with highest Ts in these pixels was selected as the hot pixel. In contrast, for cold pixel, the pixels with low Ts (top 10%) and high NDVI (top 10%) in the image were selected first, and further to select the pixels with the land covers of dense vegetation (generally forest) from the pixels selected in last step, finally, the pixel with lowest Ts in these pixels was selected as the cold pixel. It should be noted that China area is made up by 28 tiles of remote sensing image (MODIS data), and each tile was computed independently, as well as hot and cold pixels selection independent in the ET generating process (Long et al., 2011). After hot and cold pixels 610 Several iterations were carried out until the value of H was stable. Then, Eq. 1 was used to calculate λET. However, it should be noted that all the estimated energy component was an instantaneous value including latent heat; therefore, the concept of the evaporation fraction (Λ) was used to temporally scale-up from the instantaneous value to the daily ET. The evaporation fraction was defined as the ratio of latent heat to available energy (e.g., R n -G) (Eq. 30). Several studies have indicated that the evaporation fraction can be regarded as constant throughout the day (Crago, 1996); therefore, the daily ET can be calculated as follows: where ET daily , R daily , and G daily are the daily evapotranspiration, daily net radiation, and daily soil heat flux, respectively. Finally, the daily ET value was calculated. More details about SEBAL can be found in Bastiaanssen et al. (1998a).