BESS-STAIR: a framework to estimate daily, 30-meter, and all-weather crop evapotranspiration using multi-source satellite data for the U.S. Corn Belt

. With increasing crop water demands and drought threats, mapping and monitoring of cropland evapotranspiration 15 (ET) at high spatial and temporal resolutions becomes increasingly critical for water management and sustainability. However, estimating ET from satellite for precise water resources management is still challenging due to the limitations in both existing ET models and satellite input data. Specifically, the process of ET is complex and difficult to model, and existing satellite remote sensing data could not fulfill high resolutions in both space and time. To address the above two issues, this study presented a new high spatiotemporal resolution ET mapping framework, i.e., BESS-STAIR, which 20 integrates a satellite-driven water-carbon-energy coupled biophysical model BESS (Breathing Earth System Simulator) with a generic and fully-automated fusion algorithm STAIR (SaTallite dAta IntegRation). In this framework, STAIR provides daily 30-meter multispectral surface reflectance by fusing Landsat and MODIS satellite data to derive fine-resolution leaf area index and visible/near-infrared albedo, all of which, along with coarse-resolution meteorological and CO 2 data, are used to drive BESS to estimate gap-free 30-m resolution daily ET. We applied BESS-STAIR from 2000 through 2017 in six areas 25 across the U.S. Corn Belt, and validated BESS-STAIR ET estimations using flux tower measurements over 12 sites (85 site-years). Results showed that BESS-STAIR daily ET achieved an overall R 2 = 0.75, with RMSE = 0.93 mm d -1 and relative error = 27.9% when benchmarked with the flux measurements. In addition, BESS-STAIR ET estimations well captured the spatial patterns, seasonal cycles, and interannual dynamics in different sub-regions. The high performance of the BESS-STAIR framework is primarily resulted from: (1) the implementation of coupled constraints on water, carbon, and energy in 30 BESS, (2) high-quality daily 30-m data from STAIR fusion algorithm, and (3) BESS’s applicability under all-sky conditions. BESS-STAIR is calibration-free and has great potentials to be a reliable tool for water resources management and precision agriculture applications for the U.S. Belt, and even for input data. as compared to flux tower measurements over the three Mead sites from 2001 through 2012. The linear equation


75
The other critical requirement for accurate estimations of ET at high spatiotemporal resolutions is satellite input data at high resolutions in both space and time. This is challenging because existing satellite missions cannot satisfy the two conditions simultaneously. Data fusion techniques, which take multi-sensor data to generate fusion data with high resolutions in both space and time, provide a possible and scalable solution. Several such algorithms have been developed over the past decade Houborg and Mccabe, 2018;Zhu et al., 2010), and they have been successful for localized applications 80 Gómez et al., 2016;Wu et al., 2015). Notably, energy balance and thermal-based ET models such as ALEXI/DisALEXI and SEBS have been combined with the fusion algorithm such as STARFM and ESTARFM to generate daily 30-m ET estimations with favourable performance at several sites (Anderson et al., 2018;Cammalleri et al., 2013;Ma et al., 2018).

85
Here we propose and present a new ET estimation framework that combines BESS with a novel fusion algorithm SaTallite dAta IntegRation (STAIR) (Luo et al., 2018), for accurate ET estimation at both high resolution in time and space. BESS has demonstrated its high performance in estimating ET at medium to coarse resolutions, but the major obstacle of moving BESS's ET estimation to finer resolutions is the lack of key vegetation status variables at higher spatial resolutions, including leaf area index (LAI), and visible and near-infrared albedo (αVIS and αNIR). In BESS, these surface information are 90 critical to resolving spatial heterogeneity, while environmental information such as radiation, temperature, humidity and CO2 concentration are relatively homogeneous. To cope with the absence of high spatiotemporal resolution vegetation data, we propose to couple STAIR with BESS. STAIR is a genetic and fully-automated fusion algorithm to generate cloud-/gap-free surface reflectance product in high spatiotemporal resolution (Luo et al., 2018). Instead of manually selecting image pairs adopted by most other data fusion algorithms, STAIR automatically takes full advantage of time-series of daily coarse-95 resolution images and fine-resolution but less frequent images. Moreover, STAIR's high efficiency in computation allows scalability for large scale productions, which enable this new framework to deliver daily 30-m ET at regional and decadal scales.
The objective of this study is to address a fundamental issue in agro-ecological science and applications: lack of high 100 spatiotemporal gap-free ET data for decision-making. We implemented a new ET estimation framework BESS-STAIR and tested it at six study areas across the U.S. Corn Belt from 2000 to 2017. This is the first attempt to couple a satellite-driven LSM with data fusion technique to provide daily 30m-resolution ET estimations at regional and decadal scales. While existing frameworks retrieve clear-sky ET from satellite-observed LST and fill ET gaps for cloudy-sky days, BESS-STAIR simulates all-sky ET and LST as a result of crop biophysical properties. This manner has more referential significance for 105 crop modeling studies and has potential of breaking a new path to agro-ecological science and applications. We conducted comprehensive evaluation on the BESS-STAIR ET estimations with regards to the overall performance, spatial patterns, seasonal cycles and interannual dynamics, benchmarked on the ET observations from 12 eddy-covariance flux towers across the U.S. Corn Belt. The paper also discussed on the performance, advantages, limitations and potential improvements of the BESS-STAIR ET framework. 110 2 Materials and methods BESS-STAIR estimates cropland ET at 30-m resolution at daily interval ( Figure 1). BESS is driven by environmental variables (radiation, temperature, humidity, and CO2 concentration), plant structural variables (LAI, αVIS and αNIR), and plant functional variables (peak maximum carboxylation rate at 25 °C (peak Vcmax25) and Ball-Barry coefficients, for C3 and C4 120 plants, respectively). Among these key inputs, LAI, αVIS and αNIR characterize crop canopy structure, which are usually very heterogeneous. In the global BESS ET product (Jiang and Ryu, 2016), these vegetation variables are derived from MODIS satellite data at 1-km resolution; while in BESS-STAIR, they are derived from 30-m resolution surface reflectance fused from high spatial resolution Landsat data and high temporal resolution MODIS data by STAIR.

The ET estimation model: BESS
BESS is a sophisticated satellite-driven water-carbon-energy coupled biophysical model designed to continuously monitor and map water and carbon fluxes (Jiang and Ryu, 2016;Ryu et al., 2011). It is a simplified land surface model, including an atmosphere radiative transfer module (Kobayashi and Iwabuchi, 2008;Ryu et al., 2018), a two-leaf canopy radiative transfer module (De Pury and Farquhar, 1997), and an integrated carbon assimilation -stomatal conductance -energy balance 130 module. Specifically, the Farquhar model for C3 and C4 plants (Collatz et al., 1992(Collatz et al., , 1991, the Ball-Berry model (Ball et al., 1987), and the quadratic form of the Penman-Monteith equation (Paw U and Gao, 1988) are used for the simulation of carbon assimilation, stomatal conductance and energy balance, respectively. This carbon-water integrated module employs an iterative procedure to solve intercellular CO2 concentration, stomatal conductance and leaf temperature for sunlit and shade canopy. Instantaneous sunlit/shade GPP and sunlit/shade/soil ET and net radiation at Terra and Aqua overpass times 135 are simultaneously estimated, followed by a temporal upscaling procedure to derive daily GPP and ET using semi-empirical cosine functions (Ryu et al., 2012). The Priestley-Taylor equation is used to compute daily potential ET (PET) based on estimated daily net radiation and meteorological inputs.
A unique feature of BESS is that BESS takes full advantages of atmospheric and land products derived from multi-source 140 satellite data. By using MOD/MYD 04 aerosol products (Sayer et al., 2014), MOD/MYD 06 cloud products (Baum et al., 2012), MOD/MYD 07 atmospheric profile products (Seemann et al., 2003), along with gap-free atmospheric data provided by MERRA-2 reanalysis products (Gelaro et al., 2017), BESS calculates direct/diffuse visible/near-infrared radiation components at 0.05° resolution. By coupling CO2 concentration derived from SCIAMACHY and GOSAT satellite data (Dils et al., 2014) with those from OCO-2 satellite data (Hammerling et al., 2012), as well as NOAA long-term field observations 145 (www.esrl.noaa.gov/gmd/ccgg/trends/), BESS derives long-term continuous monthly CO2 concentration maps. Finally, in this study BESS uses air temperature and dew point temperature provided by ERA5 reanalysis products at 0.1° resolution (Hersbach and H., 2016). In addition to these environmental variables, BESS also highly relies on vegetation structural and functional variables. By using satellite-derived LAI, αVIS and αNIR, BESS quantifies the absorption of ultraviolet/visible/near infrared radiation by sunlit/shaded canopy through a canopy radiative transfer model. This model also upscales leaf level 150 (Vcmax25) to sunlit/shade canopy, which is used in the Farquhar photosynthesis model. Vcmax25 is a parameter depending on the plant functional type (Bonan et al., 2011;Kattge et al., 2009), and its seasonal variation is empirically parameterized by LAI (Ryu et al., 2011).

The data fusion algorithm: STAIR 155
STAIR is a generic and fully-automated method for fusing multi-spectral satellite data to generate high spatiotemporal resolution and cloud-/gap-free data (Luo et al., 2018). It fully leverages the complementary strengths in the high temporal resolution MCD43A4 nadir reflectance (daily but 500 m resolution) (Schaaf et al., 2002) and the high spatial resolution Landsat L2 nadir reflectance (30-m resolution but 16-day revisiting frequency)  time series data. STAIR first imputes the missing pixels using an adaptive-average correction procedure, and then employs a local interpolation 160 model to capture finer spatial information provided by Landsat data, followed by a time-series refinement step that incorporates the temporal patterns provided by MODIS data. This strategy allows higher efficiency in missing-data interpolation as well as greater robustness against concurrently missed MODIS and Landsat observation, which is a common situation during continuous cloudy/snowy days.

165
The algorithm starts from the imputation of the missing pixels (due to cloud cover or Landsat 7 Scan Line Corrector failure) in satellite images. For MODIS images, a Savitzky-Golay filter is first applied to reconstruct continuous time series. For Landsat images, a two-step approach is employed using both temporal and spatial information from clear-sky observations.
First, a temporal interpolation through a linear regression is applied as the initial gap-filling, based on the whole time series of images throughout a year. Second, an adaptive-average correction procedure is applied to remove inharmonic spatial 170 patterns between gap-filled and original data. The target image is partitioned into multiple segments, each of which contains one type of homogeneous pixels. The relative difference between a gap pixel and neighbourhood pixels of it within the same segment is calculated using clear-sky observations acquired in several dates close to the target image acquisition date. Based on the assumption that the relative difference remains roughly the same across different dates in a short time period (e.g., < 2-3 weeks), such difference is used to correct the filled values of the gap pixel derived from temporal interpolation so that 175 the spatial relationship between the gap-filled pixel and its neighbourhood pixels within the same segment is consistent with those in clear-sky observations. The STAIR fusion algorithm fully exploits the spatial and temporal information in the time series of gap-filled MODIS and Landsat images throughout the growing season (April -October). A nearest neighbour sampling is conducted for all the 180 MODIS images to achieve the same image size, pixel resolution and projected coordinate system with Landsat images.
Difference image is calculated for each pair of Landsat and resampled MODIS images, and a linear interpolation is applied to reconstruct the difference image for any given date when no Landsat image is available. Such difference image is used to correct the resampled MODIS image on that date and to generate a fused Landsat image. In this manner, the fused image captures the most informative spatial information provided by the high spatial resolution Landsat data and incorporates the 185 temporal patterns provided by the high temporal resolution MODIS data without any user interference. The fusion algorithm is applied to the six Landsat bands: blue, green, red, near-infrared (nir), shortwave infrared 1 (swir1), and shortwave infrared-2 (swir2).

Derivation of BESS inputs from STAIR data 190
At global scale, LAI, αVIS and αNIR can be obtained from MODIS and other satellite data, but for field-scale agricultural applications high spatial resolution data are needed to account for the spatial heterogeneity between fields or within a field.
At this point, we employed two approaches to estimate 30-m resolution daily LAI from STAIR fused surface reflectance data: an empirical approach based on linear relationship with vegetation indices (VIs) and a mechanic approach based on inversion of a canopy radiative transfer model (RTM). 195 First, we estimated LAI using the empirical approach, because of availability of field LAI measurements in the study area.
We calculated four VIs calculated from STAIR-derived spectral reflectance: Wide Dynamic Range Vegetation Index where ρB, ρG, ρR, ρN, and ρSW1 refer to the surface reflectance at blue, green, red, near-infrared, the first shortwave-infrared band, respectively. Subsequently, we used field measured LAI data collected using destructive method at Mead, Nebraska from 2001 through 2007 to build VI-LAI relationships (Gitelson et al., 2007). For each of the four VIs we build a linear regression between time series of VI and LAI for corn, soybean, and the combination of corn and soybean, respectively 205 (Table 1). At this point, the equation derived from the combination of corn and soybean was used for vegetation cover other than corn and soybean. Although this might cause bias for forest LAI estimation, it is not a concern in this study as we focused on crop ET only. We applied linear regressions to four VIs separately and averaged the four derived LAI as the final LAI estimation, with the expectation that such average would reduce uncertainty caused by individual VI-LAI relationship.
210 Table 1. Linear equations for LAI (y) as a function of VI (x) for corn, soybean, and the combination of corn and soybean.

VI
Corn Soybean Combination WDRVI y = 6.288x + 4.631 y = 4.584x + 3.432 y = 5.745x + 4.288 GWDRVI y = 8.964x + 5.875 y = 6.384x + 4.275 y = 8.110x + 5.395 EVI y = 10.569x -2.165 y = 8.116x -1.936 y = 9.665x -1.993 LSWI y = 9.156x + 1.070 y = 7.553x + 0.888 y = 8.944x + 0.982 Second, we inversed PROSAIL RTM using a look-up table (LUT) method. PROSAIL is an efficient and widely-used model to simulate canopy reflectance given a set of sun-object-view geometry, canopy structure, leaf biochemical, and soil optical parameters (Jacquemoud et al., 2009). It is a combination of the PROSPECT leaf hyperspectral properties model 215 (Jacquemoud et al., 1996;Jacquemoud and Baret, 1990) and the SAIL canopy bidirectional reflectance model (Verhoef, 1984(Verhoef, , 1985. PROSAIL is particularly suitable for grasslands and croplands (Darvishzadeh et al., 2008;Xu et al., 2019), and therefore used in this study. LUT is a robust and easy method to retrieve model parameters from observed canopy reflectance (Verrelst et al., 2018). It is based on the generation of simulated canopy reflectance database for a number of plausible combinations of model parameter value ranges, and the identification of parameter values in the database leading to the best 220 agreement between simulated and observed canopy reflectance. LUT is particularly suitable for big data processing (Myneni et al., 2002), and therefore used in this study.
We established a database by running PROSAIL with sampled parameter values listed in Table 2. For computation efficiency, we only sampled varied values for four parameters while others were fixed. These four free parameters, including 225 LAI (10 values), fraction of vegetation cover (6 values), soil brightness (5 values) and chlorophyll content (4 values), were chosen because they have been identified as the most sensitive parameters in canopy radiative transfer models (Bacour et al., 2002;Mousivand et al., 2014). Leaf inclination distribution function is also sensitive but we set fixed types "spherical", "planophile" and "plagiophile" for corn, soybean, and other biomes, respectively (Nguy-Robertson et al., 2012;Pisek et al., 2013). The fixed values of other parameters were set according to literature (Baret et al., 2007;Feret et al., 2008;230 Jacquemoud et al., 2009). Solar zenith angle at satellite overpass time can be calculated so we did not set it as a free parameter. Instead, we built a set of databases with solar zenith angle values (°) of 20, 25, 30, 35, 40, 45 and 50, respectively, representing the range during growing season in the study area. In PROSAIL, specific absorption coefficients and refractive index of leaf material are pre-measured hyperspectral data from 400 to 2500 nm with 1 nm interval (Feret et al., 2008), we averaged them over wavelengths to match Landsat 7 bands and assumed differences of spectral ranges between Landsat 5, 235 Landsat 7 and Landsat 8 have marginal influence on LAI retrieval. We did not use default soil spectrum in PROSAIL, but spatiotemporally averaged all cropland pixels spectral reflectance in April when no crop is planted across the study area to derive representative soil spectral reflectance. To retrieve LAI, we compared STAIR-derived surface reflectance (RSTAIR) with records in the canopy reflectance database simulated by PROSAIL (RPROSAIL) pixel by pixel. We used root mean square error (RMSE) as the cost function which was defined as: where λ = 1,2, … l indicates band number and l = 6 for STAIR. Ideally, the simulated reflectance in the database yielding the 245 smallest RMSE can be considered as the best simulation, and the corresponding LAI value can be considered as the solution for the satellite pixel. However, in reality the solution might not be unique, because different parameter combinations could derive similar reflectance simulations and errors in both satellite and model could further amplify this problem (Verrelst et al., 2018). For this reason, we chose top 10% small RMSE simulations in the database and considered the average of corresponding LAI values as the final solution. The threshold 10% was decided by evaluating LUT-retrieved LAI against 250 field-measured LAI at three Mead sites, and it was within a reasonable range from top 50 records to top 20% records suggested by previous studies (Duan et al., 2014;Weiss et al., 2000).
We further employed semi-empirical equations to calculate αVIS and αNIR (Liang, 2001) from STAIR-derived spectral 255 reflectance in six Landsat bands: where ρSW2 is the surface reflectance at the second shortwave-infrared band. For C3 crops/grasses, forests, and C4 crops/grasses, peak Vcmax25 values were set 180, 60 and 45, respectively (Kattge et al., 2009;Zhang et al., 2014). Ball-Berry slope and intercept are another two important parameters used in the stomatal conductance model, and their values were set 13.3 and 0.02 for C3 crops/grasses, 9.5 and 0.005 for forests, and 5.8 and 0.04 for C4 crops, respectively (Miner et al., 2017). 260 Distributions of C3 and C4 crops were obtained from Crop Data Layer (CDL) data (Boryan et al., 2011).

Evaluation of BESS-STAIR ET
The A total of 12 cropland sites scattered in six areas across the U.S. Corn Belt are registered in the AmeriFlux or FLUXNET network with publicly-available ET data ( Figure 2 and Table 3). These sites include both corn only and corn/soybean 275 rotation sites and both rainfed and irrigated sites, covering typical cropping patterns in the U.S Corn Belt. All of them were used in this study to ensure the representativeness of the validation for the precision agriculture applications in this region.

Performance of STAIR LAI 300
LAI is the key input of BESS. The accuracy of high-resolution LAI estimations determine the validity of high-resolution ET estimations. We evaluated VIs-based LAI and RTM-based LAI estimations derived from 30-m resolution STAIR fused surface reflectance data against field measurements. We also compared them with 500-m resolution MODIS LAI. Overall, STAIR-derived LAI agree well with measured LAI, with R 2 > 0.85, RMSE < 0.8 and mean bias error (MBE) ≈ 0 (Figure 3).
The RTM-based method which is calibration free yields same performance with VIs-based method which requires 305 substantial field measurements to build empirical relationships. Misclassification of CDL data between corn and soybean is

Performance of BESS-STAIR ET
BESS-STAIR daily ET estimations are in a highly aligned agreement with ground truth from the 12 flux-tower measurements ( Figure 4). Across all of the 12 sites, BESS-STAIR ET with RTM-based LAI achieves an overall coefficient 320 of determination (R 2 ) of 0.75, root-mean-square error (RMSE) of 2.29 MJ m -2 d -1 , relative error (E(|Xestimation-Xmeasurement|)/E(Xmeasurement), RE) of 27.9%, and no overall bias.  335 Figure 5. Site-by-site R 2 between flux tower measured and BESS-STAIR estimated daily ET for corn and soybean, respectively. Crop type is from CDL data. Figure 6 shows the comparison between BESS-STAIR daily ET estimations and flux tower measurements over site years with least data gaps in measurements. Across all of the 12 sites, BESS-STAIR well captures the seasonal characteristics of 340 ET observation from flux towers, as they exhibit generally consistent variations over the growing season. During the peak growing season (June, July and August), the radiation displays a dominant impact on measured daily ET, and it is reasonably estimated by BESS-STAIR ET as well. In most cases, measured daily ET do not show strong and fast response to precipitation and/or irrigation, possibly due to the plentiful water storage in soil. Two exceptions are US-IB1 (2006) and US-Ne3 (2012). In case of US-IB1, no precipitation is available in August and little in July. As a result, daily ET measurements 345 drop slightly quicker in August than other cases. Such anomaly is also depicted by BESS-STAIR ET. In case of US-Ne3, the severe drought in 2012 summer causes slightly lower ET values than the two adjacent irrigation sites (US-Ne1 and US-Ne2).
BESS-STAIR ET also captures this considerable reduction, although a slight bias is observed in July.  Although 2012 is a severe drought year, soil water content (SWC) at the rainfed site US-Ne3 still shows a relatively high level (> 0.2). As a result, ET/PET from both BESS-STAIR and flux tower are at the same level with the adjacent two 360 irrigated sites (US-Ne2 and US-Ne3). their peaks in July. Different seasonal cycles for corn, soybean and grass are also captured. Grass has the highest ET among the three vegetation types from April through June. Corn has higher ET than soybean in June and July, and decreases quickly since August. Soybean has the lowest ET from April through June, but has the highest ET in August. and a negative linear relationship (r = -0.58, p < 0.05) is observed between BESS-STAIR ET/PET and VPD. The relative stronger relationship between ET/PET and VPD than that between ET/PET and precipitation indicates atmospheric water demand is likely to contribute more to drought than soil water supply in this area.

Performance of BESS-STAIR ET
In this study, we have presented BESS-STAIR, a new framework for estimating croplands ET at field and daily scale, and we have demonstrated its high performance in the U.S. Corn Belt. The process-based biophysical model BESS, driven by 30-m resolution vegetation-related variables derived from STAIR fused surface spectral reflectance data ( Figure 3) and 410 medium resolution environmental inputs derived from MODIS and other satellite data (Figure 1), is able to produce gap-free ET and PET estimations on field-scale and at daily interval across space and time (Figure 4-7). Over the 12 sites across the U.S. Corn Belt (Figure 2), BESS-STAIR explains 75% variations in flux tower measured daily ET (Figure 4), with an overall RMSE of 2.29 MJ m -2 d -1 (equivalent to 0.93 mm d -1 or 26 W m -2 ), a 27.9% relative error, and stable performance across sties (Figure 5), as well as consistent seasonal dynamics with respect to flux tower measurements (Figure 6-7). 415 The error statistics of BESS-STAIR are commensurate with previous high resolution croplands ET mapping studies. Typical RMSE values include 25 W m -2 by TSEB-DTD (Guzinski et al., 2014), 35 W m -2 by METRIC (Irmak et al., 2011), 62 W m -2 by SEBS (McCabe and Wood, 2006), 0.60 mm d -1 by SSEBop (Senay et al., 2016), and 1.04 mm d -1 by SEBAL (Singh et al., 2008). Nevertheless, it is worth mentioning that those studies used original Landsat data and therefore suffered from 420 considerably large data gaps. In contrast, BESS-STAIR uses daily Landsat-MODIS fusion data free from any gaps, which leads to temporally continuous ET estimation at the field level, thus can meet the requirements of precision agriculture. In addition, it is worth mentioning that BESS-STAIR is calibration-free and therefore is scalable. It also indicates that the accuracy of BESS-STAIR ET is likely to further improve by using locally optimized driving force or parameter values.

425
BESS-STAIR is also comparable to other croplands ET mapping studies utilizing data fusion techniques. For example, DisALEXI-STARFM daily ET estimates were validated against the flux tower measurements over the three Mead sites (Yang et al., 2018). They reported error statistics around 1.2 mm d -1 RMSE and 29% relative error. BESS-STAIR's performance at these three same sites shows an average of 0.89 mm d -1 RMSE and 25.3% relative error (Figure 11). At monthly scale, the average RMSE and relative errors are only 0.48 mm d -1 and 14.3% ( Figure A1). In addition, BESS-STAIR 430 has a potential to apply to any croplands around the world back to 1984 when both high spatial resolution data (e.g., Landsat/TM) and high temporal resolution data (e.g., NOAA/AVHRR) were available. Figure 11. Scatter plots between ET measurements and ET estimations at three sites US-Ne1, US-Ne2 and US-Ne3. 435

Scientific advantages of BESS-STAIR ET
The efficacy of BESS-STAIR lies in several aspects. First, BESS is a water-carbon-energy coupled biophysical model. BESS employs atmospheric and canopy radiative transfer modules, carbon assimilation module, stomatal conductance 440 module, and energy balance module (Jiang and Ryu, 2016;Ryu et al., 2011). BESS integrates the simulation of carbon cycle, water cycle and energy cycle in the same framework. Such carbon-water-energy coupling strategy realistically and coherently simulates plant physiology and their response to the environment, specifically the carbon uptake and water loss by plants have been simulated synchronously through environmental constraints on stomatal conductance, with further constraints by available energy (Baldocchi and Meyers, 1998;Leuning et al., 1995). Many land surface models have already 445 adopted such strategy and have successfully simulated the evolution of terrestrial ecosystems (Ju et al., 2006;Sellers, 1997;Tian et al., 2010). However, this is not the case in commonly-used remote sensing models. Empirical methods, water balance methods, and Priestley-Taylor methods only focus on the water cycle. Energy balance methods, triangular space methods, and Penman-Monteith methods couple water cycle and energy cycle and consider ET in the context of energy partitioning. BESS, unlike these remote sensing models, constrains ET with regards to both energy requirement and carbon requirement, 450 thanks to explicit modeling of radiative transfer and stomatal behavior processes. For above reasons, BESS-STAIR ET does not only achieve high accuracy ( Figure A1 -A3), but also accurately capture responses to GPP, radiation, temperature, and humidity at daily scale (Table 4). Thus, BESS-STAIR has the potential to advance the understanding of crop responses to climate change through bridging remote sensing data and land surface models, which was first suggested by Sellers et al. (1997) more than 20 years ago. 455 Table 4. BESS-STAIR captures the correct response of daily ET to GPP, radiation (Rg), temperature (Ta) and humidity (VPD) as compared to flux tower measurements over the three Mead sites from 2001 through 2012. The linear equation slopes and correlation coefficients between ET and other factors are similar in flux tower measurements and BESS-STAIR estimations, for both the whole growing season (April -October) or only peak growing season (June, July and August). For 460 "flux tower" columns, Rg, Ta and VPD are from site measurements, while for "BESS-STAIR" columns, they are from satellite-derived coarse resolution inputs. The second strength is that BESS-STAIR is designed most sensitive to the variables which can be well-quantified from remote sensing data. BESS ET is most sensitive to solar radiation, followed by LAI (Ryu et al., 2011), as BESS ET is mainly 465 constrained by net radiation and GPP. In most cases, solar radiation is the predominant component of net radiation, while LAI determines the capacity of radiation absorption and subsequently determines GPP. BESS explicitly computes radiation components in high accuracy by driving an atmosphere radiative transfer model FLiES using MODIS cloud, aerosol and atmospheric profile products. Globally, BESS-estimated solar radiation has its R 2 about 0.85 and 0.95 for MODIS snapshots and 4-day averages, respectively . On the other hand, BESS-STAIR calculates high spatiotemporal 470 resolution LAI and albedo from fused surface reflectance data. Since Landsat and MODIS surface reflectance products are publicly-available and highly-reliable (Claverie et al., 2015;Masek et al., 2006), spatial heterogeneity and temporal dynamics of crop growing conditions are well captured (Figure 8). This study only uses reflectance data fused from Landsat and MODIS, but STAIR can be easily extended to further incorporate other types of data, such as Sentinel-2 (10 m resolution) and Planet Lab CubeSats (3 m resolution) (McCabe et al., 2017). By incorporating more high resolution 475 observations, the relevance of reconstructed high resolution image series can be further improved.
The third strength is that BESS-STAIR is able to perform under all-weather conditions. BESS-STAIR fills data gaps in surface reflectance, which has a smooth day-to-day variation even with changes in sky conditions . Based on filtered surface reflectance, LAI and albedo time series are well-reconstructed, and subsequently BESS-STAIR could 480 directly work under all-weather condition. In this manner, BESS-STAIR has no need to fill cloudy-sky ET using clear-sky ET estimations, which is error prone because the empirically-filled ET estimations usually lack sophisticated process-level model constraints and thus can have large uncertainties. Figure 12 shows that the estimation errors of BESS-STAIR ET do not change significantly under different sky conditions, with low to high "sky clearness index" referring more cloudy to more clear sky conditions. 485

Limitations and future improvements of BESS-STAIR ET
In this study, several inputs used by BESS have some limitations in terms of generality and accessibility. First, three plant functional parameters, peak Vcmax25, Ball-Barry slope and intercept are obtained from literatures, assuming constant given C3 495 or C4 plant type. Other land surface models tend to use the similar strategy by assigning fixed values to a given plant functional type (PFT) (Bonan et al., 2011;Kattge et al., 2009;Miner et al., 2017). The drawback of this strategy is the overlook of within-PFT variations and the feedback mechanisms between vegetation and its environment (Van Bodegom et al., 2014). These limitations might be mitigated by incorporating innovative leaf trait estimation techniques emerged in recent years, such as imaging spectroscopy (Serbin et al., 2015), sun-induced fluorescence (Zhang et al., 2018), and plant 500 optimization theory (Walker et al., 2017;Wang et al., 2017). Second, BESS-STAIR in this study uses CDL data which is only available in United States. Fortunately, BESS does not require specific crop types but only C3/C4 distributions, and the separation of the major C4 crop maize from other crops is practical using time-series satellite data (Cai et al., 2018;Zhong et al., 2016). It is noted that misclassification of C4 and C4 crops are likely to cause large bias in GPP, but relatively small bias in ET (Fig. A1 -A3). 505 Though BESS-STAIR is able to capture water stress impact on ET in the U.S. Corn Belt where atmospheric demands play a major role, its applicability to regions where soil supply dominates needs further investigation. Some studies suggest that optical signal as an indicator of drought performs at a longer time scale than thermal signal does (Otkin et al., 2017).
Drought first decreases soil moisture content due to enhanced ET induced by high atmospheric demand, then decreases ET 510 due to low soil moisture content, and finally causes damage to plants which changes surface reflectance. Accordingly, LAI may not serve as a relevant early warning of droughts. Furthermore, severe soil moisture stress may cause physiological deterioration in addition to structural damage that has been reflected in LAI. To address this issues for dry regions, we acknowledge that LST observations may provide essential adding values. At this point, the capacity of BESS-STAIR in estimating LST leads to a possibility of optimizing BESS-STAIR using satellite-derived LST. Recent advances of innovative 515 thermal observation platforms such as ECOSTRESS (Hulley et al., 2017), GOES-R (Schmit et al., 2017), and Sentinel-3 (Zheng et al., 2019) have provide great opportunity to integrate satellite-derived LST with the BESS-STAIR.
The BESS model itself in essence estimates instantaneous ET. The ratio of snapshot potential solar radiation to daily potential solar radiation is adopted as a scaling factor for the temporal upscaling of ET (Ryu et al., 2012). In this study, 520 BESS runs two times per day, utilizing radiation components derived from Terra/MODIS (around 11:00 AM) and Aqua/MODIS (around 1:00 PM) data, respectively. The two instantaneous ET estimates are separately upscaled to daily estimates and averaged. In spite of robustness of the upscaling algorithm (Ryu et al., 2012), bias cannot be avoided if the sky conditions at two overpass times are not representative for that day, which is natural and common in the presence of moving cloud. Since BESS is a time-independent model and can perform at any time during daytime, adding more snapshots to 525 account for the diurnal variations of radiation can solve this problem. Unfortunately, fine-resolution polar-orbiting satellite usually have similar overpass times (10:00 AM -11:00 AM and 1:00 PM -2:00 PM), so even adding more satellites is likely to bring redundant information only. Reanalysis radiation data covering diurnal cycle have limited accuracy and coarse resolution (Babst et al., 2008;Zhang et al., 2016b), so they may be unable to provide much added values as well.
Next-generation geostationary satellites, acquiring data with both high spatial and high temporal resolutions such as GOES-530 R and GaoFen-4 (Goodman et al., 2012;Xu et al., 2017), are expected to enable BESS-STAIR ET in hourly or sub-hourly interval and subsequently generate more realistic daily ET estimates.

Conclusions
In this study we presented BESS-STAIR, a new framework to estimate high spatiotemporal resolution ET that can be used 535 for field-level precision water resources management. BESS-STAIR couples a satellite-driven water-energy-carbon coupled biophysical model BESS with a generic and fully-automated fusion algorithm STAIR to generate gap-free 30-m resolution daily ET estimations. Comprehensive evaluation of BESS-STAIR ET estimations revealed: 1) reliable performance over 12 flux tower sites across the U.S. Corn Belt, and 2) reasonable spatial patterns, seasonal cycles and interannual dynamics. The proposed BESS-STAIR framework has demonstrated its ability to provide significant advancements with regard to daily 540 field-level estimations of ET at regional and decadal scales. We expect BESS-STAIR to become a solid tool for precision water resources management and other precision agriculture applications for the U.S. Corn Belt as well as other agricultural areas around the world, thanks to the global coverage of input data.
Data availability. The data generated in this study are available upon request. 545 Competing interests. The authors declare that they have no conflict of interest.
Author contribution. C.J. and K.G. designed the study, C.J. conducted the modeling and analysis, M.P. and B.P. provided research inputs during the analysis, Y.R. and S.W. provided guidance in the BESS model and the STAIR algorithm, 550 respectively. All the authors contributed to the writing of the manuscript. Landsat data free of charge. We also thank NASA freely share the MODIS products.

Appendix
Appendix 1. Key equations for energy balance in BESS. Subsript i represents sunlit or shaded canopy. 570 Canopy conductance Gc,i for sunlit/shade canopy is calculated as (Ball, 1988): where m and b are Ball-Berry coefficients, RH is relative humidityf, Ca is ambient CO2 concentration, Ac,j is carbon assimilation by sunlit/shaded canopy.
Canopy net radiation Rn for sunlit/shade canopy is calculated as (Kowalczyk et al., 2006): 575 where Qs,i and QL,i are net shortwave and longwave radiation for sunlit/shaded canopy, Tf,i is leaf temperature for sunlit/shaded canopy, Ta is air temperature, cp is specific heat, Gr is canopy conductance for radiation: At this point, Ac,j is calculated by the Farquhar photosynthesis model for C3 (Collatz et al., 1991) and C4 (Collatz et al., 1992) crops, Qs,i and QL,i are calculated by radiative transfer models for shortwave (Ryu et al., 2011) and longwave (Kowalczyk et al., 2006), respectively. Tf,i is solved through an iterative procedure (Jiang and Ryu, 2016). Subsequently, latent heat flux ƛEi 580 for sunlit/shaded canopy is calculated by the quadratic form of the Penman-Monteith equation : where ƛ is latent heat of vaporization, and = 1 2 where es(Ta) is saturated vapour pressure, ra is aerodynamic resistance, ρ is air density, γ is psychrometric constant, D is vapour pressure defecit, rc,i is canopy resistance, the reciprocal of Gc,i.