Evaluating the SSEBop and RSPMPT Models for Irrigated Fields Daily Evapotranspiration Mapping with MODIS and CMADS Data

: It is of great convenience to map daily evapotranspiration (ET) by remote sensing for agricultural water management without computing each surface energy component. This study used the operational simpliﬁed surface energy balance (SSEBop) and the remote sensing-based Penman– Monteith and Priestly–Taylor (RSPMPT) models to compute continuous daily ET over irrigated ﬁelds with the MODIS and CMADS data. The estimations were validated with eddy covariance (EC) measurements. Overall, the performance of RSPMPT with locally calibrated parameters was slightly better than that of SSEBop, with higher NSE (0.84 vs. 0.78) and R 2 (0.86 vs. 0.81), lower RMSE (0.78 mm · d − 1 vs. 0.90 mm · d − 1 ), although it had higher bias (0.03 mm · d − 1 vs. 0.01 mm · d − 1 ) and PBias (1.41% vs. 0.59%). Due to the consideration of land surface temperature, the SSEBop was more sensitive to ET’s change caused by irrigation before sowing in March and had a lower PBias (6.7% vs. 39.8%) than RSPMPT. On cloudy days, the SSEBop is more likely to overestimate ET than the RSPMPT. To conclude, driven by MODIS and CMADS data, the two simple models can be easily applied to map daily ET over cropland. The SSEBop is more practical in the absence of measured data to optimize the RSPMPT model parameters. respectively. It is more convenient to compute daily net radiation from CMADS data when there is no site of observed sunshine duration. The SSEBop Rn approach was chosen to drive the two models for actual evapotranspiration simulation for its higher accuracy.


Introduction
Accurately estimating regional evapotranspiration (ET) is crucial for the wide range of research in hydrology, climate, crop yield forecasting, and drought monitoring. Knowledge of continuous daily ET of cropland is a practical approach to monitor spatiotemporal variations in water use. It is the key to the rational allocation and management of agricultural water resources in semi-arid and arid regions, where water resources are scarce. Agricultural irrigation water is mainly from channels and pumping wells [1].
At the site scale, accurate quantification of ET can be achieved using in situ measurements, such as the observations from weighing lysimeters, Bowen ratio systems, large aperture scintillometer (LAS), and eddy covariance (EC) systems [1]. Among all of the in situ ways, the EC is a relatively popular method for both actual ET measurements and ET estimations' validation, such as the global FLUXNET network, which contains the worldwide EC observations and is widely applied in global water and energy fluxes research [2]. However, in situ methods are costly and only can obtain water and energy fluxes over small regions. For regional or large-scale ET estimations, remote sensing technology plays an irreplaceable role for its convenience of obtaining the land surface parameters spatially, such as surface emissivity, surface temperature (LST), surface albedo, and vegetation indexes, which are the necessary driving factors of ET simulations. For the past decades, many remote sensing ET models were developed and widely applied, and validated worldwide over various land-use types. The surface energy balance (SEB) and overestimated. Whether CMADS data applies to actual ET estimation over irrigated fields in northwest China's semi-arid and arid areas still needs careful investigation.
In this study, the CMADS daily variables were selected as the meteorological data for actual ET estimation. The RSPMPT and SSEBop models were adopted to compute continuous daily ET over irrigated fields in an oasis-desert region in the northwest of China. The modeled daily ET was also compared with the SSEBop global ET product provided by Senay et al. (2013) [15]. The results were validated using in situ observations from automatic meteorological and eddy covariance systems installed in the study area. Furthermore, the uncertainties resulting from CMADS data and the different schemes of models, the models' performances for cloudy days, and the error analysis are discussed.

Study Area
The Heihe River Basin (HRB) is the second-largest inland river basin in the arid region of northwest China. The study area is located in an oasis-desert region in the middle reaches of the HRB and is distributed in a range of about 150 km × 150 km. The elevation ranges from 2000 to 3000 m above sea level (a.s.l.), and the annual precipitation is 100 to 250 mm. Maize is the main crop of the irrigated fields that grows from April to October and consumes plenty of the limited water resources, covering approximately 12% of the entire HRB [22]. There is no crop in the fields from November to March. Figure 1a shows the standard false-color composite image from the GF1 WFV3 sensor in 2016 with a spatial resolution of 16 m. The red part is the central area of the oasis, and the blue point is the location of the meteorological station. The land-use type extracted from the MCD12Q1 International Geosphere-Biosphere Programme (IGBP) product is shown in Figure 1b; this region is mainly covered by croplands, grasslands, barren areas, etc. The grassland in Figure 1b is primarily alpine meadow. It is close to the upper reaches of HRB, where the maximum altitude is approximately 5000 m a.s.l. and the minimum air temperature approaches −40 • C. Over the past few decades, many scientific experiments have been carried out to better understand the ecological and hydrological process of the Heihe river basin [23].

EC and AWS Data
The observation system of meteorological elements gradient of superstation installed in the Daman irrigated area continuously collects data using eddy covariance and automatic weather systems through the HiWATER projects [24,25]. The ground-based measurements from 2013 to 2016 of this site were used to validate the estimations.
The raw EC observations were processed by EdiRe software, including spike detection, lag correction of H 2 O/CO 2 relative to the vertical wind component, sonic virtual temperature correction, coordinating rotation, correction for density fluctuation, and frequency response correction. The half-hourly sensible heat flux (H) and latent heat flux (LE) were provided in the final released EC measurements. The automatic weather system (AWS) observations were provided at every 10 min frequency, including the surface net radiation (Rn), the soil heat flux (G) measured using soil heat flow plate at a depth of 6 cm underground, mean soil temperature measured at the depths of 2 to 4 cm underground, soil moisture at different depths, wind velocity, air temperature and humidity at different heights, etc. The AWS measurements were averaged to 30 min frequency to keep consistent with the eddy data. The soil heat flow plate measured G was corrected by the Averaging Soil Thermocouple Probe (TCAV) [26] method based on the mean soil temperature and soil moisture measurements.
To solve the energy non-closure, the Bowen ratio correction method was adopted under the assumption that the ratio of latent heat to sensible heat is constant before and after closure correction [27]. Finally, the mean diurnal variation (MDV) method was used to fill the latent heat gaps to obtain continuous daily evapotranspiration measurements [28]. The average yearly ET from 2013 to 2016 is 759.5 mm, which is close to the value of 760.5 mm at the same site, according to Zhou et al. (2018) [29].
in 2016 with a spatial resolution of 16 m. The red part is the central area of the oasis, and the blue point is the location of the meteorological station. The land-use type extracted from the MCD12Q1 International Geosphere-Biosphere Programme (IGBP) product is shown in Figure 1b; this region is mainly covered by croplands, grasslands, barren areas, etc. The grassland in Figure 1b is primarily alpine meadow. It is close to the upper reaches of HRB, where the maximum altitude is approximately 5000 m a.s.l. and the minimum air temperature approaches −40 °C. Over the past few decades, many scientific experiments have been carried out to better understand the ecological and hydrological process of the Heihe river basin [23].

MODIS and CMADS Data
The MODIS remote sensing products, including the 500 m resolution daily land surface reflectance (MOD09GA) and 1 km resolution land surface temperature (MOD11A1), were used for ET model calculation. The Normalized Difference Vegetation Index (NDVI) was computed based on the near-infrared and red bands reflectance. This study also used the version 4 SSEBop Monthly Evapotranspiration product for comparison (From https://earlywarning.usgs.gov in 1 October 2020). As reported, it improves ET's spatial accuracy by subdividing and modeling each MODIS tile into 25 units [15].
To eliminate the noise caused by cloud on the vegetation index, the Savitzky-Golay filtering was employed to process the daily NDVI. The upper envelope defined by the local maximum was detected to obtain a new NDVI time series. The leaf area index (LAI) and fractional vegetation cover (fv) were calculated using the empirical formulations calibrated based on NDVI and observed values in the same study area [1].
The version 1.1 CMADS datasets (From http://www.cmads.org in 1 October 2020) contain the daily averaged wind speed at 10 m height (m/s), the daily maximum and minimum air temperature (T max and T min , k), daily mean air temperature (T mean , k), and daily mean relative humidity (RH) and solar radiation (R s , MJ·m −2 ·d −1 ) with 25 km resolution. The MOD11A1 and CMADS data were interpolated to a spatial resolution of 500 m by the bilinear method [21]. As the Penman-Monteith reference evapotranspiration equation needs the wind speed at the height of 2 m (u 2 ), the CMADS wind velocity was transformed into this height using the wind profile equation [30].

Daily Net Radiation
The daily net radiation can be expressed as the difference between the incoming net shortwave radiation (R nsd ) and the outgoing net longwave radiation (R nl ): where α is the albedo and was set as 0.23 for the hypothetic grass reference crop, and R s is the solar radiation that can be derived from CMADS. The outgoing net longwave radiation (R nl ) was computed based on the daily maximum and minimum air temperature, daily mean vapor pressure, and the extraterrestrial solar radiation (R a ), according to Senay et al. (2018) [16]. Senay et al. [16] introduced a new modification to estimate the solar radiation by replacing the clear-sky with the "average-sky" condition; k Rs is the adjustment coefficient and ranges from 0.16~0. 19.
The average-sky condition assumes a given pixel in a given day will be subjected to an average cloud cover for the rest of the day despite the clear-sky condition during the instantaneous moment of satellite overpass, avoiding a general overestimation of the net radiation, which would lead to an overestimation of ET [16]. The two approaches for solar radiation computation were compared with the AWS measured ground truth to determine a better solution.

RSPMPT Model
The RSPMPT model computes the total latent heat flux as the sum of the transpiration from the vegetation canopy [31] and the soil surface evaporation [32]. The main formula will vary depending on whether the wet canopy evaporation term is distinguished [17,18]. In this study, the simple form was used and calculated as follows: where R nc, R ns , and G are the daily net radiation (W·m −2 ) over canopy and soil surface and the soil heat flux (W·m −2 ), respectively. R nc and R ns were estimated by the fractional vegetation cover (fv) along with the net radiation, G was estimated by the ratio of soil heat flux to net radiation according to Su (2002) [5], λ is the latent heat of evaporation of water, ∆ is the slope of the saturation vapor pressure versus temperature curve [30], γ is the psychrometric constant of~0.066 (Kpa·K −1 ), ρ is the air density (kg·m −3 ), C p is the specific heat of the air (J·kg −1 ·K −1 ), and e s and e a are the saturated and actual vapor pressures (Kpa) and calculated by the daily mean air temperature and relative humidity (RH) from CMADS data. The relative surface wetness (f wet ) is used to determine whether the surface is wet or not and is equal to RH 4 . The f sm is a soil moisture constraint used as an indicator of surface soil moisture. The coefficient α pt is an empirical multiplier and was initially set to 1.26 for well-watered surfaces. The canopy resistance r c (s/m) was found to vary with different environmental variables and was modeled as a product of the response functions to the environmental factors according to Jarvis's model [33]. As soil moisture is not available in CMADS data, the solar radiation (f sd ), the vapor pressure deficit (f vpd ), and the air temperature (f ta ) factors were considered to scale the minimum canopy resistance (r cmin ) as follows: The r cmin was set as 36.6 s/m, which was optimized by MODIS and meteorological data in the same research area [34]. The aerodynamic resistance r a (s/m) controls heat and water vapor transfer from the evaporating surface into the air above the canopy. It was calculated mainly based on the surface aerodynamic roughness and the daily wind velocity from CMADS. The RSPMPT model can thoroughly explain how plants' transpiration rate varies with the increase or decrease in the environmental constraints under different conditions [35].

SSEBop Model
The SSEBop model assumes that the available net radiation mostly drives the surface energy balance process. Differences in land surface temperature can quantify a decline in ET due to water stress and other factors. SSEBop directly estimates pixel-based actual evapotranspiration (ET a ) by multiplying the ET fraction (ET f ) and reference ET (ET r ) as follows [15]: where the k max is the coefficient that scales the grass reference ET (ET r ) into the level of a maximum ET experienced by an aerodynamically rougher crop, and a recommended value for k max is 1.2. The ET r is calculated using the FAO-56 PM equation with the support of CMADS daily meteorological parameters [30].
The ET fraction is a critical parameter for the actual ET computation in SSEBop, and is calculated as: where the T s is the MODIS land surface temperature; dT is a predefined temperature difference (K) between the "hot-dry" and "wet-cold" surfaces for each pixel and calculated as: where r ah is the aerodynamic resistance to heat flow from a hypothetical bare and dry surface (110 sm −1 ), and the minimum value of dT is recommended as 6 K, reducing the ET uncertainty in cold climates. The T c is the estimated T s at the idealized "cold-wet" limit for each pixel and computed as: where c is a scene-specific factor that relates air temperature to T s_cold , which is the average land surface temperature on a well-watered, fully transpiring vegetation surface, T s_cold is obtained on all pixels where NDVI is greater than or equal to 0.7 in the image of a particular day [36]. If no pixels in the images of some days meet the NDVI conditions, the median value of all the effective c values within that year is used to replace the c value of those images. The minimum ET f is set to 0. The maximum ET f is 1.05, and the ET f exceeding 1.3 is suspected of being contaminated by clouds and is filled through temporal interpolation by assuming linear relationships between two consecutive satellite dates. The SSEBop model eliminates the subjectivity and uncertainty that could be introduced during manual hot and cold reference selection, and more details can be found in [14][15][16]37,38].

Daily Net Radiation
Daily (24 h) net radiation (Rn) is the energy boundary condition of the two models. The two approaches, one computed from the CMADS solar radiation (CMADS Rn), the other from the SSEBop equations under the "averaged sky" assumption (SSEBop Rn), were compared with the AWS measured daily net radiation. Generally, both the two methods could reflect the annual change of daily net radiation. However, there was an evident overestimation of daily net radiation based on the CMADS solar radiation, which had a higher root mean square error (RMSE) of 7.23 MJ·m −2 ·d −1 . In contrast, the daily net radiation from SSEBop had a lower RMSE of 3.12 MJ·m −2 ·d −1 from 2013 to 2016. Figure 2 also shows the overall Nash-Sutcliffe efficiency (NSE), coefficient of determination (R 2 ), mean bias error (bias), and percent bias (PBias) statistics for daily net radiation estimations. The NSE value close to 1 represents the overall high simulation accuracy of the model. The negative bias and PBias values indicated that both the two approaches overestimated daily net radiation. The bias error can be removed by applying the bias correction. It was found that both the NSE and RMSE improved after bias correction for the two approaches. The NSE increased to 0.77 and 0.72, and the RMSE decreased to 2.51 and 2.66 MJ·m −2 ·d −1 for SSEBop and CMADS Rn, respectively. It can be concluded that although the SSEBop Rn slightly overestimated the daily net radiation, it had lower RMSE, better NSE than those of CMADS Rn before or after the bias correction. As the assessment of outgoing longwave radiation had the same formula, the overestimation of CMADS Rn primarily resulted from the uncertainties of the solar radiation provided by CMADS data. As reported by Santhi et al. (2001) [39], the model performance was considered to be acceptable if the NSE and R 2 were larger than 0.5 and 0.6, respectively. It is more convenient to compute daily net radiation from CMADS data when there is no site of observed sunshine duration. The SSEBop Rn approach was chosen to drive the two models for actual evapotranspiration simulation for its higher accuracy.

Daily Net Radiation
Daily (24 h) net radiation (Rn) is the energy boundary condition of the two models. The two approaches, one computed from the CMADS solar radiation (CMADS Rn), the other from the SSEBop equations under the "averaged sky" assumption (SSEBop Rn), were compared with the AWS measured daily net radiation. Generally, both the two methods could reflect the annual change of daily net radiation. However, there was an evident overestimation of daily net radiation based on the CMADS solar radiation, which had a higher root mean square error (RMSE) of 7.23 MJ·m −2 ·d −1 . In contrast, the daily net radiation from SSEBop had a lower RMSE of 3.12 MJ·m −2 ·d −1 from 2013 to 2016. Figure 2 also shows the overall Nash-Sutcliffe efficiency (NSE), coefficient of determination (R 2 ), mean bias error (bias), and percent bias (PBias) statistics for daily net radiation estimations. The NSE value close to 1 represents the overall high simulation accuracy of the model. The negative bias and PBias values indicated that both the two approaches overestimated daily net radiation. The bias error can be removed by applying the bias correction. It was found that both the NSE and RMSE improved after bias correction for the two approaches. The NSE increased to 0.77 and 0.72, and the RMSE decreased to 2.51 and 2.66 MJ·m −2 ·d −1 for SSEBop and CMADS Rn, respectively. It can be concluded that although the SSEBop Rn slightly overestimated the daily net radiation, it had lower RMSE, better NSE than those of CMADS Rn before or after the bias correction. As the assessment of outgoing longwave radiation had the same formula, the overestimation of CMADS Rn primarily resulted from the uncertainties of the solar radiation provided by CMADS data. As reported by Santhi et al. (2001) [39], the model performance was considered to be acceptable if the NSE and R 2 were larger than 0.5 and 0.6, respectively. It is more convenient to compute daily net radiation from CMADS data when there is no site of observed sunshine duration. The SSEBop Rn approach was chosen to drive the two models for actual evapotranspiration simulation for its higher accuracy.

Validation of the Estimated ET with EC Measurements
The daily ET computed from RSPMPT, SSEBop with CMADS data, and derived from the version 4 SSEBop monthly global evapotranspiration product (hereafter, RSPMPT-CMADS, SSEBop-CMADS, SSEBop-Global) were compared with the eddy covariance system measured ET from 2013 to 2016 ( Figure 3 and Table 1).

Validation of the Estimated ET with EC Measurements
The daily ET computed from RSPMPT, SSEBop with CMADS data, and derived fro the version 4 SSEBop monthly global evapotranspiration product (hereafter, RSPMP CMADS, SSEBop-CMADS, SSEBop-Global) were compared with the eddy covariance sy tem measured ET from 2013 to 2016 ( Figure 3 and Table 1).  As shown in Figure 3a and Table 1, both the RSPMPT-CMADS and SSEBop-CMAD showed a good performance of daily ET compared to the ground truth. Overall, the pe formance of RSPMPT with locally calibrated parameters (fv and LAI) was slightly bett than that of SSEBop, with higher NSE (0.84 vs. 0.78) and R 2 (0.86 vs. 0.81), lower RM (0.78 mm·d −1 vs. 0.90 mm·d −1 ), although it had higher bias (0.03 mm·d −1 vs. 0.01 mm·d and PBias (1.41% vs. 0.59%) from 2013 to 2016. Figure 3b shows the change process monthly ET derived from different approaches and the EC observations. All the metho  As shown in Figure 3a and Table 1, both the RSPMPT-CMADS and SSEBop-CMADS showed a good performance of daily ET compared to the ground truth. Overall, the performance of RSPMPT with locally calibrated parameters (fv and LAI) was slightly better than that of SSEBop, with higher NSE (0.  Figure 3b shows the change process of monthly ET derived from different approaches and the EC observations. All the methods have the potential to deliver consistent changes in ET at different crop phenological stages. On the monthly scale, the NSE and R 2 of the two methods became higher, and the RMSE was close; the RSPMPT still had higher bias and PBias than SSEBop. Unlike other years, both the two models tended to overestimate ET from June to August in 2015. Table 2 shows the statistics of averaged measured and estimated ET in different crop growth stages from 2013 to 2016. From November to February, there is no crop in the field during the cold and dry winter, and the croplands have little water or energy for evaporation. In this study area, spring agricultural irrigation was carried out in mid-March. Due to the increase of soil moisture content, farmland ET gradually began to increase. In March, the average ET approached 40.7 mm, and SSEBop-CMADS had the best estimates, with a PBias of 6.7%. In contrast, the SSEBop-Global showed an apparent underestimation (with an average ET of 5.9 mm and a PBias of 85.4%). As crop growth requires a lot of irrigation during the growing season, both vegetation transpiration and crop evaporation begin to rise. Much higher ET values were found in both estimations or observations during the crop peak growing season from June to August. The three approaches had similar estimation accuracy, with PBias from −6.3% to 2.8% during this period. At the early and late stages of crop growth, roughly during April to May and September to October, the cropland ET was lower than that in peak growing season, and RSPMPT-CMADS had the best performance.  Figure 4 shows the cropland's evapotranspiration maps from 2013 to 2016 in the study area. The evapotranspiration spatial distribution obtained by different methods was inconsistent. The RSPMPT-CMADS and SSEBop-CMADS had more similar ET patterns, while the SSEBop-Global showed lower annual ET throughout this region.

Spatial Comparisons of the Estimated ET
Although the underlying surface is homogeneous, the maps demonstrate the spatial heterogeneity of evapotranspiration, which may be caused by the differences in irrigation intensity, irrigation types, crop types and conditions, and other environmental factors. In general, the RSPMPT-CMADS had spatially higher ET than the other two approaches and mostly ranged from 680~810 mm/y in the years 2013, 2014, and 2016. In 2015, both the RSPMPT-CMADS and SSEBop-CMADS had spatially more extensive ET than those of other years in the Daman irrigated area around the EC and AWS station.
The density scatter plots of the RSPMPT-CMADS and SSEBop-CMADS estimated annual croplands evapotranspiration from 2013 to 2016 are displayed in Figure 5. There was a moderate correlation between the two models' estimations with the correlation coefficients between 0.66 to 0.77 based on the statistic from more than 6500 pixels. The RSPMPT-CMADS computed higher ET than that of the SSEBop-CMADS from 2013 to 2016 in the study area. Relative large RMSE values were found between the two models, ranging from 74.0 to 107.2 mm/year, indicating the two approaches' different performances at spatial scale. The ET density told that more pixels were concentrated in evapotranspiration between 600~800 mm/year for the two methods, and the highest density of annual ET was around 700 mm/y for the croplands.

The Models' Performances for Cloudy Days
The most critical part to continuously monitor ET by remote sensing is to reconstruct the ET time series of cloudy days. The LST quality control documents in the MOD11A1 product were used to determine whether there are clear pixels or not around the EC site. Figure 6 and Table 3 show the scatter plots and the statistics for the estimated ET by the two models on clear and cloudy days. From 2013 to 2016, there were 805 clear days and 656 cloudy days. The RSPMPT-CMADS and SSEBop-CMADS models had the coefficients of determination (R 2 ) of 0.89 and 0.84 on clear days. On cloudy days, both models' performance decreased slightly with R 2 of 0.79 and 0.73, respectively. However, the root mean square errors (RMSE) of the two models were low on both clear and cloudy days. The relatively high NSE and R 2 for both clear and cloudy days indicates the effectiveness of the two models for estimating evapotranspiration continuously over irrigated fields. From the bias and PBias, both the two models showed slight overestimation of ET on cloudy days. Table 4 shows the canopy resistance (r c ) and ET f during the crop growing and non-growing seasons, which are the critical parameters of RSPMPT and SSEBop, respectively. Lower canopy resistance and higher ET f were found in the growing season. In both growing and non-growing seasons, the daily canopy resistance on cloudy days was higher than on clear days. The solar radiation environmental stress increased due to the decrease of the downward solar shortwave radiation on cloudy days. However, the ET f in SSEBop models did not show a significant decline on cloudy days. Therefore, in contrast, SSEBop is more likely to overestimate ET on cloudy days, with a PBias of −6.5%, which is lower than that of the RSPMPT (−2.2%). The RSPMPT estimated ET on cloudy days using smoothed NDVI and the meteorological data, whereas the SSEBop model estimated ET by filling the ET f on cloudy days through temporal interpolation based on adjacent valid ET f values. The

The Models' Performances for Cloudy Days
The most critical part to continuously monitor ET by remote sensing is to reconstruc the ET time series of cloudy days. The LST quality control documents in the MOD11A product were used to determine whether there are clear pixels or not around the EC site Figure 6 and Table 3 show the scatter plots and the statistics for the estimated ET by th two models on clear and cloudy days. From 2013 to 2016, there were 805 clear days an 656 cloudy days. The RSPMPT-CMADS and SSEBop-CMADS models had the coefficient of determination (R 2 ) of 0.89 and 0.84 on clear days. On cloudy days, both models' perfor mance decreased slightly with R 2 of 0.79 and 0.73, respectively. However, the root mea square errors (RMSE) of the two models were low on both clear and cloudy days. Th relatively high NSE and R 2 for both clear and cloudy days indicates the effectiveness o the two models for estimating evapotranspiration continuously over irrigated fields. From the bias and PBias, both the two models showed slight overestimation of ET on cloud days. Table 4 shows the canopy resistance (rc) and ETf during the crop growing and non growing seasons, which are the critical parameters of RSPMPT and SSEBop, respectively Lower canopy resistance and higher ETf were found in the growing season. In both grow ing and non-growing seasons, the daily canopy resistance on cloudy days was higher tha on clear days. The solar radiation environmental stress increased due to the decrease o the downward solar shortwave radiation on cloudy days. However, the ETf in SSEBo models did not show a significant decline on cloudy days. Therefore, in contrast, SSEBo is more likely to overestimate ET on cloudy days, with a PBias of −6.5%, which is lowe than that of the RSPMPT (−2.2%). The RSPMPT estimated ET on cloudy days usin smoothed NDVI and the meteorological data, whereas the SSEBop model estimated E by filling the ETf on cloudy days through temporal interpolation based on adjacent vali ETf values. The ETf variations caused by the changes of environmental factors are not con sidered in the interpolation process.

The Models' Performances for Growing Season
It is essential to identify the evapotranspiration estimation accuracy in crop growing season when the model is used in crop water-consuming monitoring [40]. Table 5 shows the error statistics for the evapotranspiration from RSPMPT-CMADS and SSEBop-CMADS models during growing and non-growing seasons. The simulation accuracy of the two models in the crop growing season is significantly better than that in the non-growing season. The RSPMPT-CMADS had higher NSE and R 2 , and lower RMSE than the SSEBop-CMADS in both growing and non-growing seasons. However, the SSEBop-CMADS had better bias and PBias, especially in the non-growing season, and the RSPMPT-CMADS significantly underestimated cropland evapotranspiration. During the non-growing season, soil water evaporation is the primary source of evapotranspiration. The canopy resistance in the RSRMPT model is substantial (~5000 s/m) as the LAI is very low, so the ET from vegetation decreases. The soil transpiration is scaled from its potential level based on the soil moisture constraint and relative surface wetness. The soil moisture constraint estimated by vapor pressure has more uncertainties than that estimated by soil moisture [41]. It takes a long time and a more extensive range for surface soil moisture change to affect the water vapor pressure near the surface. The soil moisture constraint estimated by vapor pressure is more suitable for evapotranspiration estimation over large scales and long time intervals (e.g., eight days). In contrast, the pixel-based land surface temperature can be sensitive to soil moisture. The irrigated fields with higher soil moisture always experience cooler land surface temperature than the dryland fields [40]. The better bias proved that when soil moisture changes in the non-growing season (e.g., irrigation before sowing or after rainfall), the SSEBop model has more potential to reflect the evolution of evapotranspiration. It is consistent with the statistics shown in Table 2 that the SSEBop model can effectively reflect the increase of evapotranspiration caused by irrigation around the AWS station before sowing in March.

Discussion
Both RSPMPT and SSEBop models are easy to apply at a regional scale because they do not need to compute each intermediate variable in the SEB equation and can compute daily ET without instantaneous estimations. Since CMADS directly provides daily meteorological parameters necessary for the two models, it is of great convenience to simulate ET by combining the two models with the CMADS data. The validation of daily net radiation shows that the overestimation of CMADS solar radiation leads to a lower negative percent bias (PBias), indicating the apparent overestimation of net radiation. It is recommended to use the maximum and minimum air temperature to compute the daily net radiation by CMADS under "average-sky" conditions.
As can be seen in Figure 3, both the two models overestimated ET in 2015 significantly. Since the daily net radiation estimations do not show apparent annual variation, the CMADS daily meteorological variables were further explored and compared to the AWS averaged measurements from 2013 to 2016. As shown in Figure 7, from 2013 to 2016, the wind speed was apparently higher than the other years. In addition, the air temperature (minimum, maximum, and averaged) and air pressure were slightly higher in 2015 than those in other years. Significantly higher wind speed may compute lower aerodynamic resistance and lead to higher canopy transpiration for RSPMPT and higher reference evapotranspiration for SSEBop.

Discussion
Both RSPMPT and SSEBop models are easy to apply at a regional scale because they do not need to compute each intermediate variable in the SEB equation and can compute daily ET without instantaneous estimations. Since CMADS directly provides daily meteorological parameters necessary for the two models, it is of great convenience to simulate ET by combining the two models with the CMADS data. The validation of daily net radiation shows that the overestimation of CMADS solar radiation leads to a lower negative percent bias (PBias), indicating the apparent overestimation of net radiation. It is recommended to use the maximum and minimum air temperature to compute the daily net radiation by CMADS under "average-sky" conditions.
As can be seen in Figure 3, both the two models overestimated ET in 2015 significantly. Since the daily net radiation estimations do not show apparent annual variation, the CMADS daily meteorological variables were further explored and compared to the AWS averaged measurements from 2013 to 2016. As shown in Figure 7, from 2013 to 2016, the wind speed was apparently higher than the other years. In addition, the air temperature (minimum, maximum, and averaged) and air pressure were slightly higher in 2015 than those in other years. Significantly higher wind speed may compute lower aerodynamic resistance and lead to higher canopy transpiration for RSPMPT and higher reference evapotranspiration for SSEBop. In previous studies, the two models' uncertainties have been thoroughly discussed [18,42]. The sensitivity analysis of the RMPMPT model showed that the canopy resistance was responsible for 95.7% of the variability in the ET results [18]. In this study, the Savitzky-Golay filtering was used to process the daily NDVI smoothing. Daily NDVI calculated the daily LAI, and the LAI and meteorological data estimated the canopy resistance. RSPMPT is easier to use because it does not require land surface temperature input. However, the parameters, LAI, and the fractional vegetation cover (fv) were also critical and sensitive for the model calculation [1,43]. At present, statistical regression methods are In previous studies, the two models' uncertainties have been thoroughly discussed [18,42]. The sensitivity analysis of the RMPMPT model showed that the canopy resistance was responsible for 95.7% of the variability in the ET results [18]. In this study, the Savitzky-Golay filtering was used to process the daily NDVI smoothing. Daily NDVI calculated the daily LAI, and the LAI and meteorological data estimated the canopy resistance. RSPMPT is easier to use because it does not require land surface temperature input. However, the parameters, LAI, and the fractional vegetation cover (fv) were also critical and sensitive for the model calculation [1,43]. At present, statistical regression methods are used to optimize the estimation of these parameters. In this study, the regression formula with local calibration in the same study area was used. ET estimation accuracy was improved compared with other LAI and fv empirical formulas (not shown).
The SSEBop model is sensitive to input variables, land surface temperature (LST) and reference ET, and parameters, differential temperature (dT), and maximum ET scalar (k max ), particularly during the non-growing season and in dry areas [42]. The necessary process of the SSEBop model is to determine the cold pixel coefficient (C factor) within the image range and calculate the "hot/dry" and "cold/wet" surface temperatures at the two idealized reference conditions for each pixel in the same period. Ji et al. (2019) defined the cold pixel as NDVI > 0.7 for vegetation cover based on the knowledge and experience with the use of satellite-derived vegetation indices in ET modeling and mapping [36]. In this study, the cold pixels were determined over the cropland, mostly located in the middle reach of the Heihe river basin. Besides, cold pixels can also be determined over all vegetation types within the image. For comparison, Figure 8 shows that the C values calculated over croplands were slightly higher than those of C calculated over all vegetation types. By comparing the minimum, maximum, middle, and average values of the C factors calculated by the two methods, it was found that the maximum value was almost the same. In contrast, the C factor's minimum values had a large gap, and the median and average values had a slight gap, which showed that the C factor calculated by all vegetation types was lower.
, 11, x FOR PEER REVIEW 14 of 18 used to optimize the estimation of these parameters. In this study, the regression formula with local calibration in the same study area was used. ET estimation accuracy was improved compared with other LAI and fv empirical formulas (not shown). The SSEBop model is sensitive to input variables, land surface temperature (LST) and reference ET, and parameters, differential temperature (dT), and maximum ET scalar (kmax), particularly during the non-growing season and in dry areas [42]. The necessary process of the SSEBop model is to determine the cold pixel coefficient (C factor) within the image range and calculate the "hot/dry" and "cold/wet" surface temperatures at the two idealized reference conditions for each pixel in the same period. Ji et al. (2019) defined the cold pixel as NDVI > 0.7 for vegetation cover based on the knowledge and experience with the use of satellite-derived vegetation indices in ET modeling and mapping [36]. In this study, the cold pixels were determined over the cropland, mostly located in the middle reach of the Heihe river basin. Besides, cold pixels can also be determined over all vegetation types within the image. For comparison, Figure 8 shows that the C values calculated over croplands were slightly higher than those of C calculated over all vegetation types. By comparing the minimum, maximum, middle, and average values of the C factors calculated by the two methods, it was found that the maximum value was almost the same. In contrast, the C factor's minimum values had a large gap, and the median and average values had a slight gap, which showed that the C factor calculated by all vegetation types was lower.  [16], the median value of the C factor within one year is suggested to replace the image without "cold/wet" pixels. In farmland area, the period when "cold/wet" pixels cannot be found is mostly non-growing season. A lower median value of the C factor may lead to a lower ET fraction (ETf) value in this period. Figure 9  According to Senay et al. (2018) [16], the median value of the C factor within one year is suggested to replace the image without "cold/wet" pixels. In farmland area, the period when "cold/wet" pixels cannot be found is mostly non-growing season. A lower median value of the C factor may lead to a lower ET fraction (ET f ) value in this period. Figure 9 shows the daily ET f values from 2013 to 2016, and the lowest values calculated over all vegetation types concentrated on non-growing season. A lower ET f means a lower ET estimate. From 2013 to 2016, the pixel-averaged cropland ET estimated by C factor calculated over all vegetation types was about 40 mm lower than that estimated by C factor calculated over croplands only. The possible reason is that alpine meadow surface temperature in high-altitude areas is lower than that of farmland. Simultaneously, the spatial resolution of CMADS is coarse, and the difference of maximum temperature between the two types is not as evident as that of surface temperature. It is worth noting that the version 4 SSEBop ET product is also significantly underestimated in the nongrowing season, which eventually leads to the underestimation of annual crop ET. This product improves the spatial accuracy of ET by subdividing and modeling each MODIS tile into 25 units. However, it is also an extensive area, close to 450 km, due to the MODIS scan width of 2330 km. To a large extent, if the altitude, climate, and vegetation types vary significantly in different parts, the selection of vegetation types is also crucial for determining cold pixels. vegetation types concentrated on non-growing season. A lower ETf means a lower ET estimate. From 2013 to 2016, the pixel-averaged cropland ET estimated by C factor calculated over all vegetation types was about 40 mm lower than that estimated by C factor calculated over croplands only. The possible reason is that alpine meadow surface temperature in high-altitude areas is lower than that of farmland. Simultaneously, the spatial resolution of CMADS is coarse, and the difference of maximum temperature between the two types is not as evident as that of surface temperature. It is worth noting that the version 4 SSEBop ET product is also significantly underestimated in the non-growing season, which eventually leads to the underestimation of annual crop ET. This product improves the spatial accuracy of ET by subdividing and modeling each MODIS tile into 25 units. However, it is also an extensive area, close to 450 km, due to the MODIS scan width of 2330 km. To a large extent, if the altitude, climate, and vegetation types vary significantly in different parts, the selection of vegetation types is also crucial for determining cold pixels.

Conclusions
Monitoring cropland ET continuously is essential for agricultural water resources management and drought monitoring, especially for the arid areas with limited water resources. Remote sensing-based ET models can simulate the spatiotemporal features of ET.
This study's main objective was to evaluate the performance of the RSPMPT and SSEBop models driven by MODIS and CMADS data over irrigated fields in the Heihe river basin of China. The eddy covariance observed ET was used to make a quantitative evaluation of the models.
This study demonstrated that both the RSPMPT and SSEBop models produced robust daily ET results. Overall, the SSEBop model's performance was slightly inferior to that of the RSPMPT model with local calibration. The performance result also strengthens the useful application of the CMADS data for daily ET estimation. The solar radiation is overestimated, and the wind speed has an apparent interannual mutation from 2013 to 2016 in the study area. This study also found that the cold pixel coefficients (C factor) in SSEBop determined over all vegetation types or croplands caused different ET estimation. The ET estimated using C factors from all vegetation types was lower than using C factors from croplands. During the non-growing seasons, the C factor cannot be determined from

Conclusions
Monitoring cropland ET continuously is essential for agricultural water resources management and drought monitoring, especially for the arid areas with limited water resources. Remote sensing-based ET models can simulate the spatiotemporal features of ET.
This study's main objective was to evaluate the performance of the RSPMPT and SSEBop models driven by MODIS and CMADS data over irrigated fields in the Heihe river basin of China. The eddy covariance observed ET was used to make a quantitative evaluation of the models.
This study demonstrated that both the RSPMPT and SSEBop models produced robust daily ET results. Overall, the SSEBop model's performance was slightly inferior to that of the RSPMPT model with local calibration. The performance result also strengthens the useful application of the CMADS data for daily ET estimation. The solar radiation is overestimated, and the wind speed has an apparent interannual mutation from 2013 to 2016 in the study area. This study also found that the cold pixel coefficients (C factor) in SSEBop determined over all vegetation types or croplands caused different ET estimation. The ET estimated using C factors from all vegetation types was lower than using C factors from croplands. During the non-growing seasons, the C factor cannot be determined from croplands. The C factor calculated over alpine meadow was too low for the croplands during this period due to the increase of soil evaporation caused by the irrigation before sowing. It is more reasonable to replace the C factor in non-growing seasons with the medium C factor values over croplands of one year. It proved that the selection of vegetation types is also crucial for determining cold pixels in the SSEBop model if the altitude, climate, and vegetation types vary significantly in different parts of the research area. On cloudy days, the SSEBop model is more likely to overestimate evapotranspiration than the RSPMPT model because its ET f time expansion method does not consider environmental factors' impact. During the non-growing season, the SSEBop model was more sensitive to the change of evapotranspiration caused by irrigation before sowing, which had a lower bias, and the random error accounted for the most proportion of total errors.
To conclude, this is the first comprehensive validation and comparison of the two simple and popular ET models driven by MODIS and CMADS data (e.g., RSPMPT-CMADS and SSEBop-CMADS) over irrigated fields in the Heihe river basin. The two approaches both had the potential to model continuous ET over irrigated fields at spatiotemporal scales and reflect the heterogeneity of croplands ET caused by different crop species, crop growth, and irrigation amount. The SSEBop is more practical over a large area because it has more capacity to capture the ET heterogeneity of irrigated fields and needs less empirical coefficient, especially in the absence of measured data, to optimize the RSPMPT model parameters.