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

With increasing crop water demands and drought threats, mapping and monitoring of cropland evapotranspiration (ET) at high spatial and temporal resolutions become increasingly critical for water management and sustainability. However, estimating ET from satellites for precise water resource 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 presents a new high spatiotemporal resolution ET mapping framework, i.e., BESS-STAIR, which 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 m multispectral surface reflectance by fusing Landsat and MODIS satellite data to derive a fine-resolution leaf area index and visible/near-infrared albedo, all of which, along with coarse-resolution meteorological and CO₂ 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 across the US 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²=0.75, with root mean square error RMSE =0.93 mm d⁻¹ and relative error RE =27.9 % when benchmarked with the flux measurements. In addition, BESS-STAIR ET estimations captured the spatial patterns, seasonal cycles, and interannual dynamics well in different sub-regions. The high performance of the BESS-STAIR framework primarily resulted from (1) the implementation of coupled constraints on water, carbon, and energy in BESS, (2) high-quality daily 30 m data from the 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 resource management and precision agriculture applications for the US Corn Belt and even worldwide given the global coverage of its input data.

1 Introduction

Accurate field-level management of water resources urgently demands reliable estimations of evapotranspiration (ET) at high spatial and temporal resolutions. ET is the sum of water loss from the soil surface through evaporation and that from plant components through leaf transpiration and evaporation, and ET at cropland is usually considered for crop water use (Allen et al., 1998). ET consumes up to 90 % of total water inputs (precipitation plus irrigation) in agro-ecosystems in the western and Midwestern United States (Irmak et al., 2012). In the US Corn Belt, where more than 85 % of corn and soybean is produced in the US (Grassini et al., 2015), increasing vapor pressure deficit (VPD) and drought sensitivity have been recognized as severe threats to future crop security (Lobell et al., 2014; Ort and Long, 2014). The vulnerability to drought in this region is further exacerbated by elevated rates of grass-to-crop conversion and expansion of irrigated areas (Brown and Pervez, 2014; Wright and Wimberly, 2013). Furthermore, precision water resource management requires the capacity to account for spatial heterogeneity and to guide real-time decision-making (GAO, 2019). Accordingly, reliable tools are urgently needed to estimate, map, and monitor the total amount and spatial and temporal variations of cropland ET.

One critical requirement for the accurate estimations of ET at high spatiotemporal resolutions is reliable and advanced satellite-based models. This is challenging because the process of ET is complex and difficult to model. ET results from balance between atmospheric water demand and soil water supply, and it is also regulated by plants through canopy development and stomatal behaviors in order to optimize their water, carbon, and energy use strategies (Katul et al., 2012; Wang and Dickinson, 2012). A large number of satellite-based ET estimation methods have been developed based on different theories and techniques. In general, they can be grouped into many categories: statistical or machine-learning methods (Jung et al., 2010; Lu and Zhuang, 2010), water balance methods (Pan et al., 2012; Wan et al., 2015), energy balance methods (Anderson et al., 1997; Su, 2002), triangular or trapezoid space methods (Jiang and Islam, 1999; Li et al., 2009), Priestley–Taylor methods (Fisher et al., 2008; Miralles et al., 2011), and Penman–Monteith methods (Mu et al., 2011; Yebra et al., 2013). Kalma et al. (2008), Li et al. (2009), and K. Zhang et al. (2016) have provided detailed reviews of the pros and cons of different remote-sensing approaches.

Given the complexity of the ET process, we argue that a reliable ET model should include both necessary biophysical processes and high-quality multi-source observations to constrain ET estimations (Loew et al., 2016). While remote-sensing-based approaches tend to focus on constraints from various satellite data, land-surface models (LSMs) are proficient at including processes that account for interactions between environment and plant structure and functions. Given the gaps between remote sensing and LSMs, a distinct ET model, the Breathing Earth System Simulator (BESS), was developed (Jiang and Ryu, 2016; Ryu et al., 2011). Different from the above-mentioned remote-sensing models, BESS is a biophysical model, which adopts modules commonly implemented in LSMs but uses various satellite remote-sensing data as direct inputs. Specifically, BESS is a two-leaf water–carbon–energy coupled model driven by environmental and vegetation variables derived from multi-source satellite data. As the energy cycle, carbon cycle, and water cycle are jointly modeled and mutually constrained in BESS, it has produced a series of high-quality global long-term (2000–2017) products, including the 5 km resolution global radiation (R_g), photosynthetically active radiation (PAR) and diffuse PAR products (Ryu et al., 2018), and 1 km resolution gross primary productivity (GPP) and ET products (Jiang and Ryu, 2016), which enables tracking of crop growth and yields too (Huang et al., 2018). In particular, the 1 km resolution BESS ET product is able to capture the total amount and spatial and temporal variations in arid/semi-arid areas like Australia (Whitley et al., 2016, 2017), California (Baldocchi et al., 2019), and northwestern China (Wei et al., 2019). The fidelity of the coarse-resolution BESS ET product suggests its potential at fine resolutions.

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 (Gao et al., 2006; Houborg and Mccabe, 2018; Zhu et al., 2010), and they have been successful for localized applications (Gao et al., 2017; 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 algorithms such as STARFM and ESTARFM to generate daily 30 m ET estimations with favorable performance at several sites (Anderson et al., 2018; Cammalleri et al., 2013; Li et al., 2017; Ma et al., 2018).

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 high resolution in both 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, this surface information is critical for resolving spatial heterogeneity, while environmental information such as radiation, temperature, humidity, and CO₂ concentration is relatively homogeneous. To cope with the absence of high spatiotemporal resolution vegetation data, we propose coupling STAIR with BESS. STAIR is a generic and fully automated fusion algorithm to generate a cloud-/gap-free surface reflectance product at 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-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 spatiotemporal gap-free ET data for decision-making. We implemented a new ET estimation framework, BESS-STAIR, and tested it in six study areas across the US Corn Belt from 2000 to 2017. This is the first attempt to couple a satellite-driven biophysical model with a data fusion technique to provide daily 30 m resolution ET estimations at regional and decadal scales. While existing frameworks retrieve clear-sky ET from satellite-observed land-surface temperature (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 way has more referential significance for crop modeling studies and has the potential to forge a new path in agro-ecological science and applications. We conducted a comprehensive evaluation of 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 US Corn Belt. The paper also discusses the performance, advantages, limitations, and potential improvements of the BESS-STAIR ET framework.

2 Materials and methods

BESS-STAIR estimates cropland ET at 30 m resolution at a daily interval (Fig. 1). BESS is driven by environmental variables (radiation, temperature, humidity, and CO₂ concentration), plant structural variables (LAI, α_VIS, and α_NIR), and plant functional variables (peak maximum carboxylation rate at 25 ^∘C (peak V_cmax25) and Ball–Barry coefficients, for C₃ and C₄ 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.

Figure 1 The BESS-STAIR framework. The BESS ET estimation model and the STAIR data fusion algorithm are highlighted in green boxes. Blue boxes are satellite data, yellow boxes are ancillary data, and red boxes are key inputs to BESS. The output of BESS-STAIR is the 30 m resolution daily ET highlighted in a white box.

2.1 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 module. Specifically, the Farquhar model for C₃ and C₄ plants (Collatz et al., 1991, 1992), 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 CO₂ concentration, stomatal conductance, and leaf temperature for a sunlit/shaded canopy. Instantaneous sunlit/shade GPP and sunlit/shade/soil ET and net radiation at Terra and Aqua overpass times 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 it takes full advantage of atmospheric and land products derived from multi-source 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 CO₂ 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 (https://www.esrl.noaa.gov/gmd/ccgg/trends/, last access: 11 March 2020), BESS derives long-term continuous monthly CO₂ concentration maps. Finally, in this study BESS uses air temperature and dew point temperature provided by ERA5 reanalysis products at 0.1^∘ resolution (Hersbach, 2016). In addition to these environmental variables, BESS also highly relies on vegetation structural and functional variables. By using a 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 (V_cmax25) to sunlit/shaded canopy, which is used in the Farquhar photosynthesis model. V_cmax25 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 the LAI (Ryu et al., 2011).

2.2 The data fusion algorithm: STAIR

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 d revisiting frequency) (Masek et al., 2006) time-series data. STAIR first imputes the missing pixels using an adaptive-average correction procedure and then employs a local interpolation 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.

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 patterns between gap-filled and original data. The target image is partitioned into multiple segments, each of which contains one type of homogeneous pixel. The relative difference between a gap pixel and neighborhood pixels of it within the same segment is calculated using clear-sky observations acquired at 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 a difference is used to correct the filled values of the gap pixel derived from temporal interpolation so that the spatial relationship between the gap-filled pixel and its neighborhood 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-neighbor sampling is conducted for all the MODIS images to achieve the same image size, pixel resolution, and projected coordinate system with Landsat images. A 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 a 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 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).

2.3 Derivation of BESS inputs from STAIR data

At a 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 a linear relationship with vegetation indices (VIs) and a mechanistic approach based on inversion of a canopy radiative transfer model (RTM).

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 (WDRVI), Green Wide Dynamic Range Vegetation Index (GWDRVI), Enhanced Vegetation Index (EVI), and Land Surface Water Index (LSWI) for corn and soybean, respectively (Eqs. 1–3). These four VIs were chosen because they utilized information from different band combinations.

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{WDRVI}={\displaystyle \frac{\mathrm{0.1}{\mathit{\rho }}_{\mathrm{N}}-{\mathit{\rho }}_{\mathrm{R}}}{\mathrm{0.1}{\mathit{\rho }}_{\mathrm{N}}+{\mathit{\rho }}_{\mathrm{R}}}},\text{(2)}& {\displaystyle }\mathrm{GWDRVI}={\displaystyle \frac{\mathrm{0.1}{\mathit{\rho }}_{\mathrm{N}}-{\mathit{\rho }}_{\mathrm{G}}}{\mathrm{0.1}{\mathit{\rho }}_{\mathrm{N}}+{\mathit{\rho }}_{\mathrm{G}}}},\text{(3)}& {\displaystyle }\mathrm{EVI}=\mathrm{2.5}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{N}}-{\mathit{\rho }}_{\mathrm{R}}}{{\mathit{\rho }}_{\mathrm{N}}+\mathrm{6}{\mathit{\rho }}_{\mathrm{R}}-\mathrm{7.5}{\mathit{\rho }}_{\mathrm{B}}+\mathrm{1}}},\text{(4)}& {\displaystyle }\mathrm{LSWI}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{N}}-{\mathit{\rho }}_{\mathrm{SW}\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{N}}+{\mathit{\rho }}_{\mathrm{SW}\mathrm{2}}}},\end{array}$$

where ρ_B, ρ_G, ρ_R, ρ_N, and ρ_SW1 refer to the surface reflectance at blue, green, red, near-infrared, and the first shortwave-infrared band, respectively. Subsequently, we used field-measured LAI data collected using the 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 (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 LAIs as the final LAI estimation, with the expectation that such an average would reduce uncertainty caused by an individual VI–LAI relationship.

Table 1 Linear equations for LAI (y) as a function of VI (x) for corn, soybean, and the combination of corn and soybean.

Second, we inversed the 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 (Jacquemoud et al., 1996; Jacquemoud and Baret, 1990) and the SAIL canopy bidirectional reflectance model (Verhoef, 1984, 1985). PROSAIL is particularly suitable for grasslands and croplands (Darvishzadeh et al., 2008; Xu et al., 2019) and is 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 a 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 agreement between simulated and observed canopy reflectance. LUT is particularly suitable for big data processing (Myneni et al., 2002) and is 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 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 the literature (Baret et al., 2007; Feret et al., 2008; 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 the growing season in the study area. In PROSAIL, specific absorption coefficients and the refractive index of leaf material are pre-measured hyperspectral data from 400 to 2500 nm with a 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, Landsat 7, and Landsat 8 have a marginal influence on LAI retrieval. We did not use the 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.

Table 2 Parameter values needed to establish the canopy reflectance database by PROSAIL.

To retrieve LAI, we compared STAIR-derived surface reflectance (R_STAIR) with records in the canopy reflectance database simulated by PROSAIL (R_PROSAIL) pixel by pixel. We used root mean square error (RMSE) as the cost function which was defined as

$$\begin{array}{}\text{(5)}& \mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{l}}{\sum }_{\mathit{\lambda }=\mathrm{1}}^l{\left[R_{\mathrm{STAIR}}\left(\mathit{\lambda }\right)-R_{\mathrm{PROSAIL}}\left(\mathit{\lambda }\right)\right]}^{\mathrm{2}}},\end{array}$$

where $\mathit{\lambda }=\mathrm{1},\mathrm{2},\mathrm{\dots }l$ indicates band number and l=6 for STAIR. Ideally, the simulated reflectance in the database yielding the smallest RMSE can be considered the best simulation, and the corresponding LAI value can be considered 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 the top 10 % small RMSE simulations in the database and considered the average of the corresponding LAI values to be the final solution. The threshold of 10 % was decided by evaluating LUT-retrieved LAI against field-measured LAI at three Mead sites, and it was within a reasonable range from the top 50 records to the top 20 % of 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 reflectance in six Landsat bands:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{VIS}}=\mathrm{0.443}{\mathit{\rho }}_{\mathrm{B}}+\mathrm{0.317}{\mathit{\rho }}_{\mathrm{G}}+\mathrm{0.240}{\mathit{\rho }}_{\mathrm{R}},\text{(7)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{NIR}}=\mathrm{0.693}{\mathit{\rho }}_{\mathrm{N}}+\mathrm{0.212}{\mathit{\rho }}_{\mathrm{SW}\mathrm{1}}+\mathrm{0.116}{\mathit{\rho }}_{\mathrm{SW}\mathrm{2}}-\mathrm{0.003},\end{array}$$

where ρ_SW2 is the surface reflectance at the second shortwave-infrared band. For C₃ crops/grasses, forests, and C₄ crops/grasses, peak V_cmax25 values were set at 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 to 13.3 and 0.02 for C₃ crops/grasses, 9.5 and 0.005 for forests, and 5.8 and 0.04 for C₄ crops, respectively (Miner et al., 2017). Distributions of C₃ and C₄ crops were obtained from Crop Data Layer (CDL) data (Boryan et al., 2011).

2.4 Evaluation of BESS-STAIR ET

The BESS-STAIR ET estimations were evaluated against flux-tower ET measurements in the US Corn Belt. The US Corn Belt (Fig. 2) generally refers to a region in the Midwestern United States that has dominated corn and soybean production in the United States (Green et al., 2018), which currently produces about 45 % and 30 % of global corn and soybean, respectively (USDA, 2014). The region is characterized by relatively flat land, deep fertile soils, and high soil organic matter (Green et al., 2018). Most parts of the US Corn Belt have favorable growing conditions of temperature and rainfall. The majority of the croplands in the US Corn Belt are rainfed, with a small portion in the western part relying on irrigation.

Figure 2 Study areas. Red dots indicate 12 flux-tower sites scattered in six areas across the US Corn Belt. The background map indicates the percent each state contributes to the total national corn and soybean plantation area (USDA, 2018). The background image is a © NASA Blue Marble image.

A total of 12 cropland sites scattered in six areas across the US Corn Belt are registered in the AmeriFlux or FLUXNET networks with publicly available ET data (Fig. 2 and Table 3). These sites include both corn only and corn–soybean rotation sites and both rainfed and irrigated sites, covering typical cropping patterns in the US 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. For six sites, US-Bo1 (Meyers and Hollinger, 2004), US-Bo2 (Bernacchi et al., 2005), US-Br1 (Prueger et al., 2003), US-Br3 (Prueger et al., 2003), US-Ro2 (Turner et al., 2016), and US-SFP (Wilson and Meyers, 2007), level 2 half-hourly data were downloaded from the AmeriFlux website (http://ameriflux.lbl.gov/, last access: 11 March 2020). For three sites, US-IB1 (Matamala et al., 2008), US-Ro1 (Griffis et al., 2010), and US-Ro3 (Griffis et al., 2010), standardized gap-filled level 4 daily mean data were downloaded from the Carbon Dioxide Information Analysis Center data archive website (https://mirrors.asun.co/climate-mirror/cdiac.ornl.gov/pub/ameriflux/, last access: 11 March 2020). For the other three sites (Suyker et al., 2004), US-Ne1, US-Ne2, and US-Ne3, standardized high-quality gap-filled daily mean data were downloaded from the FLUXNET2015 website (http://fluxnet.fluxdata.org/data/fluxnet2015-dataset/, last access: 11 March 2020).

Table 3 Information of 12 flux-tower sites used for validation.

By comparing with eddy-covariance ET, we evaluated three ET estimations: BESS-STAIR with VI-based LAI, BESS-STAIR with RTM-based LAI, and BESS-STAIR with MODIS LAI. MODIS LAI refers to the MCD15A3H 500 m resolution 4 d composite LAI product downloaded from https://lpdaac.usgs.gov/tools/data-pool/ (last access: 11 March 2020). Flux-tower measurements usually have an irregular and dynamic footprint at scales from 100 m to 1 km (Fu et al., 2014), but for simplicity, only 30 m resolution BESS-STAIR pixels containing the flux tower were used for the direct comparison. With regard to flux towers, measurement data were directly used without energy closure adjustment. For AmeriFlux level 2 data, half-hourly data were averaged to daily ET only if no gaps exist during the day to avoid sampling bias caused by missing data. For AmeriFlux level 4 data and FLUXNET2015 data, gap-filled daily ET was used directly.

3 Results

3.1 Performance of STAIR LAI

LAI is the key input of BESS. The accuracy of high-resolution LAI estimations determines the validity of high-resolution ET estimations. We evaluated VI-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 the 500 m resolution MODIS LAI. Overall, the STAIR-derived LAI agrees well with the measured LAI, with R² > 0.85, RMSE < 0.8, and mean bias error (MBE) ≈0 (Fig. 3). The RTM-based method, which is calibration-free, yields the same performance as the VI-based method, which requires substantial field measurements to build empirical relationships. Misclassification of CDL data between corn and soybean is an important uncertainty source since both methods rely on crop types. During 2001–2007, 4 out of 21 site years (19 %) over the three Mead sites were misclassified. By using the correct classification (not shown), the accuracy of LAI estimations reaches R²=0.90 and RMSE =0.62 for the VI-based method and R²=0.89 and RMSE =0.68 for the RTM-based method. By comparison, the coarse-resolution MODIS LAI has relatively large errors, especially a negative bias (R²=0.55, RMSE =1.68, and MBE $=-\mathrm{0.97}$).

Figure 3 Scatter plots between LAI measurements and LAI estimations. LAI measurements are destructively collected at the three Mead sites. (a–b) STAIR-derived daily 30 m resolution LAI using the VI-based method and the RTM-based method, respectively. (c) 500 m resolution MODIS LAI.

3.2 Performance of BESS-STAIR ET

BESS-STAIR daily ET estimations are in highly aligned agreement with ground truth from the 12 flux-tower measurements (Fig. 4). Across all of the 12 sites, BESS-STAIR ET with RTM-based LAI achieves an overall coefficient of determination (R²) of 0.75, a RMSE of 2.29 MJ m⁻² d⁻¹, a relative error $\left(E\right(|X_{\mathrm{estimation}}-X_{\mathrm{measurement}}|)/E(X_{\mathrm{measurement}})$, RE) of 27.9 %, and no overall bias. Figure 5 further exhibits its performance over all of the 12 flux-tower sites. R² values range from 0.68 to 0.94 for corn and from 0.65 to 0.81 for soybean, highlighting the robustness of BESS-STAIR ET in the US Corn Belt. BESS-STAIR ET with VI-based LAI has similar performance (R²=0.75, RMSE =2.24 MJ m⁻² d⁻¹, and RE =27.4 %). Considering the relatively small difference between BESS-STAIR using RTM-based LAI and that using VI-based LAI, only the former one, which is calibration-free, is demonstrated in the following parts of this paper. By comparison, BESS ET with MODIS LAI shows larger errors (R²=0.65, RMSE =2.50 MJ m⁻² d⁻¹, and RE =30.2 %) compared to BESS-STAIR.

Figure 4 Density scatter plots between ET measurements and ET estimations. ET measurements are from eddy-covariance data collected at 12 flux towers. (a, b) BESS-STAIR ET with VI-based LAI and RTM-based LAI, respectively. (c) 500 m BESS ET with MODIS LAI.

Figure 5 Site-by-site R² 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 the fewest data gaps in measurements. Across all of the 12 sites, BESS-STAIR captures the seasonal characteristics of ET observation from flux towers well, 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 does not show a 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 the case of US-IB1, no precipitation is available in August and little in July. As a result, daily ET measurements drop slightly more quickly in August than in the other cases. Such an anomaly is also depicted by BESS-STAIR ET. In the case of US-Ne3, the severe drought in the 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.

Figure 6 Seasonal time series of flux tower measured and BESS-STAIR estimated daily ET for 12 selected site years. Daily radiation and precipitation/irrigation are overlaid except for US-Br3.

Figure 7 Seasonal time series of daily ET/PET derived from BESS-STAIR and the flux tower for US-Ne1, US-Ne2, and NS-Ne3 in 2012, along with measured daily mean soil water content (SWC).

Figure 7 shows the comparison between BESS-STAIR ET/PET and flux-tower ET/PET at three sites (US-Ne1, US-Ne2, and NS-Ne3) at Mead, Nebraska. Overall, BESS-STAIR agrees well with the flux tower in both magnitude and seasonal cycle. Although 2012 is a severe drought year, soil water content (SWC) at the US-Ne3 rainfed site still shows a relatively high level (> 0.2). As a result, ET/PET from both BESS-STAIR and the flux tower is at the same level as the adjacent two irrigated sites (US-Ne2 and US-Ne3).

3.3 Spatiotemporal variations of BESS-STAIR ET

BESS-STAIR daily ET demonstrates prominent spatial variations within the $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ area near the Mead site in Nebraska (Fig. 8). Because of the impact of drought, central pivot irrigated fields characterized by round-shaped plots generally display higher values than surrounding croplands, and croplands have much higher values than grasslands. Variabilities of ET between different crop fields and within individual crop fields are also observable. Such variabilities might be attributed to different irrigation strategies, varieties, and/or other management. By comparison, though 500 m resolution BESS-MODIS ET is able to capture the general spatial pattern, it has many mixed pixels and is unable to demonstrate gradients across field boundaries.

Figure 8 Daily ET (MJ m⁻² d⁻¹) derived from (a) BESS-STAIR and (b) BESS-MODIS at Mead (41.1–41.2^∘ N, 96.4–96.5^∘ W) on 1 August 2012. Circles indicate flux towers in this region.

Reasonable seasonal cycles for different land-cover types are revealed by BESS-STAIR monthly ET averaged from gap-free daily estimations. An example time series of monthly ET maps at Brooks Field during the growing season of 2000 is shown in Fig. 9. BESS-STAIR ET clearly captures the temporal dynamics throughout the growing season. All vegetation shows low values (e.g., < 2 mm d⁻¹) in April, May, September, and October but high values in June, July, and August (JJA), with 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 a higher ET than soybean in June and July and decreases quickly from August. Soybean has the lowest ET from April through June, but has the highest ET in August.

Figure 9 Monthly mean BESS-STAIR ET at Brooks Field (41.9–42.0^∘ N, 96.65–96.75^∘ W) during the growing season of 2000, along with a CDL land-cover map. The last subplot shows the average time series of corn, soybean, and grass. Circles indicate flux towers in this region.

BESS-STAIR is also able to produce long-term ET/PET estimation as an indicator of drought. Figure 10 shows an example time series of peak growing season ET/PET at Bondville from 2001 through 2017. Overall substantial interannual variability is shown, with regional average ET/PET values ranging from 0.76 in an extremely dry year (2012) to 0.91 in an extremely wet year (2015). A positive linear relationship (r=0.42, p < 0.1) is observed between BESS-STAIR ET/PET and precipitation, and a negative linear relationship ($r=-\mathrm{0.58}$, p < 0.05) is observed between BESS-STAIR ET/PET and VPD. The relatively 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.

Figure 10 Peak growing season (June, July, and August) BESS-STAIR ET/PET at Bondville (39.95–40.05^∘ N, 88.25–88.35^∘ W) from 2001 through 2017 along with two scatter plots between peak growing season precipitation and ET/PET and peak growing season VPD and ET/PET over the 17 years. Precipitation and VPD data are from Daily Surface Weather Data (Daymet) at Bondville, Illinois (https://modis.ornl.gov/cgi-bin/sites/site/?id=us_illinois_bondville&product=Daymet, last access: 11 March 2020), where VPD is derived using maximum air temperature and water vapor pressure. Circles indicate flux towers in this region.

4 Discussions

4.1 Performance of BESS-STAIR ET

In this study, we have presented BESS-STAIR, a new framework for estimating cropland ET at field and daily scale, and we have demonstrated its high performance in the US Corn Belt. The BESS process-based biophysical model, driven by 30 m resolution vegetation-related variables derived from STAIR fused surface spectral reflectance data (Fig. 3) and medium-resolution environmental inputs derived from MODIS and other satellite data (Fig. 1), is able to produce gap-free ET and PET estimations at field scale and at daily intervals across space and time (Figs. 4–7). Over the 12 sites across the US Corn Belt (Fig. 2), BESS-STAIR explains 75 % of variations in flux-tower measured daily ET (Fig. 4), with an overall RMSE of 2.29 MJ m⁻² d⁻¹ (equivalent to 0.93 mm d⁻¹ or 26 W m⁻²), a 27.9 % relative error, and stable performance across sites (Fig. 5), as well as consistent seasonal dynamics with respect to flux-tower measurements (Figs. 6–7).

The error statistics of BESS-STAIR are commensurate with previous high-resolution cropland ET mapping studies. Typical RMSE values include 25 W m⁻² by TSEB-DTD (Guzinski et al., 2014), 35 W m⁻² by METRIC (Irmak et al., 2011), 62 W m⁻² by SEBS (McCabe and Wood, 2006), 0.60 mm d⁻¹ by SSEBop (Senay et al., 2016), and 1.04 mm d⁻¹ by SEBAL (Singh et al., 2008). Nevertheless, it is worth mentioning that those studies used original Landsat data and therefore suffered from 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 and 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.

BESS-STAIR is also comparable to other cropland 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 of around 1.2 mm d⁻¹ RMSE and 29 % relative error. BESS-STAIR's performance at these three same sites shows an average of 0.89 mm d⁻¹ RMSE and 25.3 % relative error (Fig. 11). At a monthly scale, the average RMSE and relative errors are only 0.48 mm d⁻¹ and 14.3 % (Appendix Fig. A1). In addition, BESS-STAIR has the 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.

4.2 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, a carbon assimilation module, a stomatal conductance module, and an energy balance module (Jiang and Ryu, 2016; Ryu et al., 2011). BESS integrates the simulation of the carbon cycle, water cycle, and energy cycle in the same framework. Such a carbon–water–energy coupling strategy realistically and coherently simulates plant physiology and its 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 adopted such a 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 the 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, thanks to explicit modeling of radiative transfer and stomatal behavior processes. For the above reasons, BESS-STAIR ET not only achieves high accuracy (Figs. A1–A3), but also accurately captures 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 by bridging remote-sensing data and land-surface models, which was first suggested by Sellers et al. (1997) more than 20 years ago.

Table 4 BESS-STAIR captures the correct response of daily ET to GPP, radiation (R_g), temperature (T_a), 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 the peak growing season (June, July, and August). For “flux-tower” columns, R_g, 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 to be 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 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 at 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² at about 0.85 and 0.95 for MODIS snapshots and 4 d averages, respectively (Ryu et al., 2018). On the other hand, BESS-STAIR calculates high spatiotemporal 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 captured well (Fig. 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 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 (Liu et al., 2017). Based on filtered surface reflectance, LAI and albedo time series are well-reconstructed, and subsequently BESS-STAIR could directly work under all-weather conditions. 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 a low to high “sky clearness index” ranging from more cloudy to more clear-sky conditions.

Figure 12 BESS-STAIR estimated daily ET has a similar performance with a varying sky clearness index (the ratio of incoming radiation on the surface to that on the top-of-atmosphere). The lower and upper boundaries of boxes refer to the first and third quartiles of error statistics. The bars inside the boxes refer to median values. The whiskers indicate 1.5 times the distance between the first and third quartiles.

4.3 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 V_cmax25 and Ball–Barry slope and intercept, are obtained from the literature, assuming a constant given C₃ or C₄ plant type. Other land-surface models tend to use a 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 that it overlooks 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 that have emerged in recent years, such as imaging spectroscopy (Serbin et al., 2015), sun-induced fluorescence (Zhang et al., 2018), and plant optimization theory (Walker et al., 2017; Wang et al., 2017). Second, BESS-STAIR in this study uses CDL data which are only available in the United States. Fortunately, BESS does not require specific crop types, but only C₃/C₄ distributions, and the separation of the major C₄ 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 C₃ and C₄ crops is likely to cause a large bias in GPP but a relatively small bias in ET (Figs. A1–A3).

Though BESS-STAIR is able to capture water stress impact on ET in the US Corn Belt where atmospheric demands play a major role, its applicability to regions where soil supply dominates needs further investigation. Some studies suggest that an optical signal as an indicator of drought performs at a longer timescale than a 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 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 issue for dry regions, we acknowledge that LST observations may provide essential adding values. At this point, the capacity of BESS-STAIR to estimate LST leads to the possibility of optimizing BESS-STAIR using satellite-derived LST. Recent advances in innovative thermal observation platforms such as ECOSTRESS (Hulley et al., 2017), GOES-R (Schmit et al., 2017), and Sentinel-3 (Zheng et al., 2019) have provided a great opportunity to integrate satellite-derived LST with 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, BESS runs two times per day, utilizing radiation components derived from Terra/MODIS (around 11:00 solar time) and Aqua/MODIS (around 13:00 solar time) data, respectively. The two instantaneous ET estimates are separately upscaled to daily estimates and averaged. In spite of the 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 account for the diurnal variations of radiation can solve this problem. Unfortunately, fine-resolution polar-orbiting satellites usually have similar overpass times (10:00–11:00 and 13:00–14:00), so adding even more satellites is likely to bring redundant information only. Reanalysis radiation data covering the diurnal cycle have limited accuracy and coarse resolution (Babst et al., 2008; X. Zhang et al., 2016), so they may be unable to provide much added value as well. Next-generation geostationary satellites acquiring data with both high spatial and temporal resolutions, such as GOES-R and GaoFen-4 (Goodman et al., 2012; Xu et al., 2017), are expected to enable BESS-STAIR ET in an hourly or sub-hourly interval and subsequently generate more realistic daily ET estimates.

5 Conclusions

In this study we presented BESS-STAIR, a new framework to estimate high spatiotemporal resolution ET that can be used for field-level precision water resource 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 US 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 field-level estimations of ET at regional and decadal scales. We expect BESS-STAIR to become a solid tool for precision water resource management and other precision agriculture applications for the US Corn Belt as well as other agricultural areas around the world, thanks to the global coverage of input data.

Appendix A: Key equations for energy balance in BESS. Subscript i represents sunlit or shaded canopy

Canopy conductance G_c,i for sunlit/shaded canopy is calculated as (Ball, 1988)

$$\begin{array}{}\text{(A1)}& G_{\mathrm{c},i}=m{\displaystyle \frac{\mathrm{RH}}{C_{\mathrm{a}}}}{A}_{\mathrm{c},i}+b,\end{array}$$

where m and b are Ball–Berry coefficients, RH is relative humidity, C_a is ambient CO₂ concentration, and A_c,j is carbon assimilation by sunlit/shaded canopy.

Canopy net radiation R_n for sunlit/shaded canopy is calculated as (Kowalczyk et al., 2006)

$$\begin{array}{}\text{(A2)}& R_{\mathrm{n},i}=Q_{\mathrm{S},i}+Q_{\mathrm{L},i}-c_{\mathrm{p}}{G}_{\mathrm{r}}\left(T_{\mathrm{f},i}-T_{\mathrm{a}}\right),\end{array}$$

where Q_s,i and Q_L,i are net shortwave and longwave radiation for sunlit/shaded canopy, T_f,i is leaf temperature for sunlit/shaded canopy, T_a is air temperature, c_p is specific heat, and G_r is canopy conductance for radiation:

$$\begin{array}{}\text{(A3)}& G_{\mathrm{r}}=\mathrm{4}\mathit{\sigma }{\mathit{ϵ}}_f{T}_{\mathrm{a}}^{\mathrm{3}}/c_{\mathrm{p}}.\end{array}$$

At this point, A_c,j is calculated by the Farquhar photosynthesis model for C₃ (Collatz et al., 1991) and C₄ (Collatz et al., 1992) crops, and Q_s,i and Q_L,i are calculated by radiative transfer models for shortwave (Ryu et al., 2011) and longwave (Kowalczyk et al., 2006), respectively. T_f,i is solved through an iterative procedure (Jiang and Ryu, 2016). Subsequently, latent heat flux λE_i for sunlit/shaded canopy is calculated by the quadratic form of the Penman–Monteith equation:

$$\begin{array}{}\text{(A4)}& a\mathit{\lambda }{E}_i^{\mathrm{2}}+b\mathit{\lambda }{E}_i+c=\mathrm{0},\end{array}$$

where λ is latent heat of vaporization, and

$$\begin{array}{}\text{(A5)}& {\displaystyle }a={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{d^{\mathrm{2}}{e}_{\mathrm{s}}\left(T_{\mathrm{a}}\right)}{\mathrm{d}{T}_{\mathrm{a}}^{\mathrm{2}}}}{\displaystyle \frac{r_{\mathrm{a}}^{\mathrm{2}}}{\mathit{\rho }{C}_{\mathrm{p}}\mathit{\gamma }\left(r_{\mathrm{a}}+r_{\mathrm{c},i}\right)}},\text{(A6)}& {\displaystyle }\begin{array}{rl}b& =-\mathrm{1}-{\displaystyle \frac{\mathrm{d}{e}_{\mathrm{s}}\left(T_{\mathrm{a}}\right)}{\mathrm{d}{T}_{\mathrm{a}}}}{\displaystyle \frac{r_{\mathrm{a}}}{\mathit{\gamma }\left(r_{\mathrm{a}}+r_{\mathrm{c},i}\right)}}-{\displaystyle \frac{d^{\mathrm{2}}{e}_{\mathrm{s}}\left(T_{\mathrm{a}}\right)}{\mathrm{d}{T}_{\mathrm{a}}^{\mathrm{2}}}}\\ & {\displaystyle \frac{R_{\mathrm{n},i}{r}_{\mathrm{a}}^{\mathrm{2}}}{\mathit{\rho }{C}_{\mathrm{p}}\mathit{\gamma }\left(r_{\mathrm{a}}+r_{\mathrm{c},i}\right)}},\end{array}\text{(A7)}& {\displaystyle }c={\displaystyle \frac{\mathit{\rho }{C}_{\mathrm{p}}}{\mathit{\gamma }\left(r_{\mathrm{a}}+r_{\mathrm{c},i}\right)}}D+{\displaystyle \frac{\mathrm{d}{e}_{\mathrm{s}}\left(T_{\mathrm{a}}\right)}{\mathrm{d}{T}_{\mathrm{a}}}}{\displaystyle \frac{r_{\mathrm{a}}}{\mathit{\gamma }\left(r_{\mathrm{a}}+r_{\mathrm{c},i}\right)}}{R}_{\mathrm{n}}+a{R}_{\mathrm{n},i}^{\mathrm{2}},\end{array}$$

where e_s(T_a) is saturated vapor pressure, r_a is aerodynamic resistance, ρ is air density, γ is a psychrometric constant, D is vapor pressure deficit, and r_c,i is canopy resistance, the reciprocal of G_c,i.

Soil evaporation is calculated as (Fisher et al., 2008)

$$\begin{array}{}\text{(A8)}& \mathit{\lambda }{E}_i={\displaystyle \frac{\frac{\mathrm{d}{e}_{\mathrm{s}}\left(T_{\mathrm{a}}\right)}{\mathrm{d}{T}_{\mathrm{a}}}}{\frac{\mathrm{d}{e}_{\mathrm{s}}\left(T_{\mathrm{a}}\right)}{\mathrm{d}{T}_{\mathrm{a}}}+\mathit{\gamma }}}\left(R_{\mathrm{n},\mathrm{soil}}-G\right){\mathrm{RH}}^{D/\mathrm{1000}},\end{array}$$

where R_n,soil is net radiation for soil and G is ground heat storage calculated as $\mathrm{0.3}\times R_{\mathrm{n},\mathrm{soil}}$.

Figure A1 Time series of monthly mean ET from flux-tower measurements and BESS-STAIR estimations.

Figure A2 Time series of monthly mean GPP from flux-tower measurements and BESS-STAIR estimations. Significant underestimations in 2003 and 2011 for Ne1, in 2005 for Ne2, and in 2005 for Ne3 are due to misclassification of corn as soybean in CDL. Significant overestimations in 2006 for Ne2 are due to misclassification of soybean as corn in CDL.

Figure A3 Time series of monthly mean Rn from flux-tower measurements and BESS-STAIR estimations.

Figure A4 Examples of daily GPP (gC m⁻² d⁻¹) derived from BESS-STAIR at Rosemount (44.65–44.75^∘ N, 93.05–93.15^∘ W) in 2017. Circles indicate flux towers in this region.