A Modeling Framework for Deriving Daily Time Series of Evapotranspiration Maps Using a Surface Energy Balance Model

Abstract

Surface energy balance models have been one of the most widely used approaches to estimate spatially distributed evapotranspiration (ET) at varying landscape scales. However, more research is required to develop and test an operational framework that can address all challenges related to processing and gap filling of non-continuous satellite data to generate time series of ET at regional scale. In this study, an automated modeling framework was developed to construct daily time series of ET maps using MODIS imagery and the Surface Energy Balance System model. The ET estimates generated from this modeling framework were validated against observations of three eddy-covariance towers in Oklahoma, United States during a two-year period at each site. The modeling framework overestimated ET but captured its spatial and temporal variability. The overall performance was good with mean bias errors less than 30 W m⁻² and root mean square errors less than 50 W m⁻². The model was then applied for a 14-year period (2001–2014) to study ET variations across Oklahoma. The statewide annual ET varied from 841 to 1100 mm yr⁻¹, with an average of 994 mm yr⁻¹. The results were also analyzed to estimate the ratio of estimated ET to reference ET, which is an indicator of water scarcity. The potential applications and challenges of the ET modeling framework are discussed and the future direction for the improvement and development of similar automated approaches are highlighted.

1. Introduction

Time series of remotely sensed evapotranspiration (ET) maps have extensive applications in agricultural, hydrological and environmental studies as they capture the spatiotemporal variability of vegetation consumptive use from field to continental scales. For example, spatial ET data have been used in agriculture sector for water right regulation, planning and monitoring [1], assessing irrigation and drainage performance [2,3,4], closing water balance at irrigation scheme levels [5] and managing agricultural water resources [6,7,8]. Recent studies have shown that remotely sensed ET can be used effectively for monitoring agricultural droughts [9,10,11] with the future potential of improving the performance of ET-integrated agricultural drought indices [12]. ET maps have been also used in assessing crop water productivity [13,14,15] and crop yield analysis [16,17].

Numerous studies have demonstrated the use of time series ET maps for ecological applications, such as capturing the progress of vegetation and wetland restoration [18], assessing the vulnerability of forest to fire and drought [19] and accounting water use from riparian vegetation and invasive species [20,21,22,23]. Remote sensing-based ET products have also been applied in improving the performances of hydrological models [24,25,26] and for climate studies to capture water feedbacks associated with seasonal cycles and soil moisture deficit at regional scales [27].

Among different approaches developed for mapping ET, the surface energy balance (SEB) approach has been widely used to acquire distributed ET at varying geographical scales [28,29,30,31]. Numerous SEB models have been proposed, including but not limited to Surface Energy Balance Index (SEBI) [32], Two-Source Energy Balance (TSEB) [33,34], Surface Energy Balance Algorithm for Land (SEBAL) [35], Simplified Surface Energy Balance Index (S-SEBI) [36], Surface Energy Balance System (SEBS) [37], Mapping Evapotranspiration at high Resolution with Internalized Calibration (METRIC) [38], Atmosphere-Land Exchange Inverse (ALEXI) [39], Regional ET Estimation Model (REEM) [40], Remote Sensing Evapotranspiration model (ReSET) [41], Operational Simplified Surface Energy Balance (SSEBop) [42] and Hybrid Dual-Source Scheme and Surface Energy Framework-Based Evapotranspiration Model (HTEM) [43]. Some of these models such as SEBAL and METRIC use manual selection of extreme pixels to compute sensible heat flux, which could result in variations in estimated ET [44] and may add uncertainty and errors based on the user’s experience [45]. Other models such as TSEB, SEBS and SSEBop do not require human intervention so that the associated uncertainties are minimized. The selection of the SEB model and the quality of input data are likely key factors to determine the accuracy of modeled ET [46].

Developing time series of ET maps requires complex, multi-step analyses to deal with issues associated with pre-processing of remote sensing data and post-processing of resulting ET products. The choice of the SEB model and satellite data could vary depending on intended applications of ET maps, availability and requirements of input data and availability of resources (time, money and expertise) to run the model. In general, the SEB-based ET estimation process can be divided into six steps: (i) collation of remotely sensed and ground-based input data, (ii) quality assessment of collected datasets and preparation of all necessary inputs for the selected SEB model, (iii) running the SEB model (including all modules and algorithms) to obtain the instantaneous ET at the time of satellite overpass, (iv) extrapolation of instantaneous ET to daily estimates, (v) filling the gaps due to cloud coverage over a portion of the map and (vi) interpolation of daily ET between image acquisition dates to obtain ET for longer time scales.

The first two steps are performed to ensure the quality of input data, a critical requirement for any remote sensing data analysis. A thorough QA/QC procedure for weather data as presented in [47,48] is necessary as the accuracy of final product depends on the quality of these datasets. The quality assurance of weather dataset is even more critical in case of SEB models as they are sensitive to weather parameters. For example, Webster et al. [49] found air temperature and wind speed as influential inputs for HTEM and SEBS models, whereas, S-SEBI was less sensitive to meteorological inputs.

For small-scale applications with similar climatic conditions, weather data from a single ground station are usually used as input in most SEB models. However, for regional applications with varying climatic conditions, distributed datasets are required. Several recent studies [50,51] have applied gridded weather datasets for mapping daily ET due to the ease of their application for regional studies. However, users need to confirm the integrity of the datasets before processing the SEB model. A study [52] found overestimation of reference ET due to biases in air temperature and wind speed in the widely used reanalysis data—North American Land Data Assimilation System when compared to reference ET estimates from the Texas High Plains ET Network [53]. The study recommended using weather station datasets within agricultural settings, whenever possible, for precise applications of time series ET information such as in irrigation scheduling. A few studies have explored the applicability of developing distributed weather data from the point measurements of a network of ground stations to account for the spatial variability of weather parameters [54].

The third step is to run the selected SEB model, which involves several sub-models to solve the SEB equation as shown in Equation (1).

$$\mathrm{LE}={\mathrm{R}}_{\mathrm{n}}-\mathrm{G}-\mathrm{H}$$

(1)

where LE is latent heat flux, R_n is net radiation, G is soil heat flux and H is sensible heat flux. All parameters are in units of W m⁻². Based on the sensible heat flux computation approach, SEB models can be categorized into single-source and two-source models. The sensible heat fluxes for soil and vegetation are computed separately in two-source models, while a single value for each pixel is computed in single-source models. Each approach has its own advantages and caveats. In theory, two-source models could provide more accurate ET over sparse vegetation as they close the energy balance separately for soil and vegetation. Timmermans et al. [44] found better accuracy from a TSEB model compared to SEBS across sparsely vegetated grasslands in the Southern Great Plains. Kustas et al. [55] reported that two-source performed better in sub-humid tallgrass prairie, whereas greater accuracy was found for a single-source model in semiarid rangeland.

As mentioned before, some single-source models require an additional step in running the model, which involves the manual selection of extreme hot and cold pixels by user. To remove the subjectivity in the selection of extreme pixels in SEBAL, Long et al. [56] introduced a trapezoidal approach to define boundary conditions for the selection of these pixels based on the relationship between vegetation fraction and surface temperature. Automated approaches have been proposed in [57,58,59] to replace human intervention. Alternative approaches are also applied by [60,61] to estimate ET from a cold pixel as a function of normalized difference vegetation index when an ideal cold pixel is difficult to find within a satellite image.

The fourth step is to extrapolate the instantaneous ET to daily values. Evaporative fraction (Λ) [35,37,62,63] and ETrF (fraction of reference ET) [38,64] are the common methods to obtain daily ET. Both of these methods assume the instantaneous Λ or ETrF is the same as the daily Λ or ETrF. However, a study [65] reported that this assumption was not satisfied when the fractional vegetation cover was close to a maximum. In the Texas Panhandle, Colaizzi et al. [66] found a better agreement of ETrF method for cropland and Λ method for bare soil when compared with lysimeter measurements. Chavez et al. [67] evaluated six extrapolation approaches on corn and soybean fields and found smaller error from Λ method when compared with eddy covariance measurements. Another study [68] found Λ method advantageous during several water stress events, whereas ETrF approach performed better under advective conditions [38,64], which could be significant in arid environments.

The fifth step is to fill the gaps caused by cloud coverage over a portion of the daily ET maps. One approach is to apply linear interpolation of nearest reliable values within an image [50]. This method is suitable when the nearest pixels are under the same land cover as that of missing pixels. However, it may not be appropriate when the area with data gap is large and encompasses heterogeneous terrain. Another approach includes the use of time-weighted interpolation of preceding and following images [69]. This method adjusts the vegetation development using normalized difference vegetation index (NDVI) across vegetated areas and residual soil moisture differences for the areas with bare soil surface. Anderson et al. [39] applied the available water for the root zone and soil surface layer to fill the gaps. The available water for the clear and cloudy days is used to estimate the daily water depletion due to ET from the root zone and soil surface layer and the fraction of available water is used to fill the gaps.

The final step is the interpolation of ET maps between consecutive satellite overpass dates to construct daily ET time series. Several interpolation and data-fusion approaches have been implemented for this purpose. A common approach is to apply linear interpolation of Λ or ETrF images between consecutive satellite overpass dates [70]. Another approach is to apply a curvilinear function using more than two Λ or ETrF images. For example, at least one cloud-free image for each month was used for spline interpolation within METRIC to obtain monthly and seasonal ET [38,71,72]. Singh et al. [70] evaluated the performance of several interpolation methods and found no significant difference in seasonal ET among cubic spline, fixed ETrF and linear interpolations. A backward-average iterative approach has been also proposed to estimate ET in between Landsat overpass dates [73].

While numerous studies have been conducted to address the issues related to specific steps involved in generating remotely sensed ET time series based on SEB models, only a few have focused on developing automated modeling frameworks, covering all hierarchical steps mentioned above. Such modeling frameworks, if validated, could have significant value in providing end-users with daily ET time series for practical applications in improving land and water management. Furthermore, a comprehensive and detailed documentation of the entire process of deriving daily ET maps at regional scales could be a useful resource to potential end-users who currently need to understand and select appropriate approaches for each of the six steps from many sources. Developing and documenting a comprehensive framework that generates complete ET time series from raw input data enables potential users outside the research community to utilize this framework for making more informed decisions and policies. The main goal of this study was to develop and document a modeling framework to construct daily time series of ET maps for the entire state of Oklahoma, USA. The performance of this framework was also evaluated by comparing its results with ET estimates of flux towers in Oklahoma. Finally, long-term variations in ET across Oklahoma were investigated.

2. Materials and Methods

2.1. Study Area

The study area covered the entire state of Oklahoma, USA, with an area of about 181,200 km² (Figure 1). Oklahoma Climate is classified as humid subtropical at most parts of the state and cold semi-arid at far west [74]. The state has nine climate divisions (CD) delineated based on precipitation and temperature gradients. The normal (1981–2010) annual precipitation is about 925 mm yr⁻¹, with significant spatial variation across CDs. While southeast (CD9) receives the largest amount of 1301 mm yr⁻¹ on average, the Panhandle (CD1) holds the smallest record of 520 mm yr⁻¹. The normal annual mean air temperature is 15.6 °C, with July and January being the hottest and coldest months, respectively. The southcentral (CD8) has the maximum mean annual temperature of 16.7 °C, whereas the Panhandle region has the minimum value at 13.6 °C. The top two land cover categories in Oklahoma are grassland (36.4%) and pastureland (11.3%) [75]. The elevation varies between 88 m above mean sea level at the southeast border with Arkansas and 1516 m at far-west border with New Mexico.

2.2. Modeling Framework

The modeling framework was designed to use daily images from the MODIS Terra satellite as input data. The single-source SEBS model [37] was selected as the SEB model for estimating energy fluxes. The main reason for the selection of SEBS over other SEB models was its applicability over large areas with heterogeneous surfaces [31]. In addition, this model does not require intermittent human intervention, which facilitates the automation process. A graphical illustration of the proposed framework is shown in Figure 2, followed by detailed explanation of specific approaches selected for each of the six computational steps mentioned before.

2.2.1. Step 1: Collation of Input Data

The daily surface reflectance (MOD09GA, [76]), daily land surface temperature (LST) and emissivity (MDO11A1, [77]) data were downloaded from the US Geological Survey Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/). Ground-based meteorological data included hourly air temperature, relative humidity, incoming shortwave solar radiation, wind speed and atmospheric pressure. These data were obtained from the Oklahoma Mesonet [78,79] weather stations installed across the state (Figure 1). The Oklahoma Mesonet is a world-class environmental monitoring network (https://www.mesonet.org/) consisting of 120 active stations with at least one station at each of the 77 counties in Oklahoma.

2.2.2. Step 2: Quality Assessment and Preparation of Inputs

The initial quality assessment of suitable MODIS images was based on cloud coverage. Images with less than 10% cloud cover were selected for further processing. Hence, any day when the cloud coverage was above 10% was assumed as a day with missing remotely sensed data. When the period of missing imagery was more than 10 consecutive days, images with less than 15% cloud cover were also included as acceptable quality. Then, cloud-covered pixels in each selected image were masked by applying a threshold of LST smaller than 250 K. These steps were repeated for all selected reflectance, LST and emissivity images. Since a single MODIS image tile was not sufficient to cover the entire state of Oklahoma, two image tiles (h09v05 and h10v05) were merged.

The quality assessment of each weather variables was performed as described in [47,48]. The solar radiation was checked against the upper limit under clear sky condition. Daily average temperature was compared against the average extreme temperatures to ensure the difference between them was within the acceptable range (2 °C) [48]. The quality of wind speed was maintained by considering gust factor threshold of more than 1. Relative humidity data were considered when the values were less than 100%. The missing weather data were filled by an average value of that parameter from four nearest Mesonet stations. Hourly alfalfa reference ET (ET_r) [48] was then computed at each station during the study period using the Bushland ET Calculator [80]. Daily ET_r estimates were obtained by summing 24-hour ET_r values. To incorporate the weather variability between the weather stations, spatial input data were generated by applying inverse distance weighted interpolation for all weather variables, including hourly and daily ET_r. As mentioned in the previous section, the Oklahoma Mesonet is a densely distributed weather station network, with about 1510 km² per station. This is a significantly finer spatial resolution than the 5000 km² per station value recommended by the World Meteorological Organization for evaporation stations on interior plains [81]. Hence, the adjustment of meteorological parameters with elevation was not considered during interpolation.

2.2.3. Step 3: The SEB Model

As mentioned before, the Surface Energy Balance System (SEBS) model of [37] was selected as the SEB model in the present study. However, other SEB models such as those reviewed in the Introduction section can be used in this step based on user resources, availability of input data and desired accuracy. Like other SEB models, SEBS estimates the latent heat flux (LE) as a residual of the land surface energy balance as shown in Equation (1). The R_n was calculated by applying the surface radiation balance equation:

$${\mathrm{R}}_{\mathrm{n}}=\left(1-\mathsf{\alpha }\right){\mathrm{R}}_{\mathrm{s}}+{\mathsf{\epsilon }}_{\mathrm{s}}{\mathsf{\epsilon }}_{\mathrm{a}}{\mathsf{\sigma }\mathrm{T}}_{\mathrm{A}}^4-{\mathsf{\epsilon }}_{\mathrm{s}}{\mathsf{\sigma }\mathrm{T}}_{\mathrm{S}}^4$$

(2)

where R_S is incoming shortwave solar radiation, α is surface albedo (dimensionless) estimated following [82], ε_a and ε_s are emissivities (dimensionless) of atmosphere and surface, estimated following [83,84], respectively. σ is the Stefan-Boltzmann constant (5.67 × 10⁻⁸ W m⁻² K⁻⁴), T_A is air temperature (K) and T_s is the surface temperature (K), estimated as a ratio of brightness temperature to ε_s^−0.25. The G was estimated by applying the relationship developed by [35]:

$$\frac{\mathrm{G}}{{\mathrm{R}}_{\mathrm{n}}}=\frac{\left({\mathrm{T}}_{\mathrm{s}}-273.15\right)}{100\mathsf{\alpha }}\left({\mathrm{c}}_1\mathsf{\alpha }+{\mathrm{c}}_2 {\mathsf{\alpha }}^2\right)\left(1-0.98{\mathrm{NDVI}}^4\right)$$

(3)

where c₁ and c₂ are calibration coefficients and were considered as 0.24 and 0.46, respectively.

SEBS uses similarity theories to estimate H: the bulk atmospheric similarity (BAS) theory for atmospheric boundary layer (ABL) scaling [85] and the Monin-Obukhov similarity (MOS) for atmospheric surface layer (ASL) scaling [86]. The ABL is a part of the atmosphere that is directly impacted by earth’s surface and responds to surface forcing with a timescale of an hour or less, whereas ASL is usually the bottom 10% of ABL [37]. During unstable conditions, an appropriate atmospheric (BAS or MOS) scaling is determined as presented in [83]. For stable conditions, functions given by [83,87] are used for ABL and ASL scaling, respectively. In the ASL, the similarity relationships for mean wind speed (u) and the difference between potential temperature profiles are derived using the MOS theory as:

$$\mathrm{u}=\frac{{\mathrm{u}}_{*}}{\mathrm{k}}\left[\ln \left(\frac{\mathrm{z}-{\mathrm{d}}_0}{{\mathrm{z}}_{0\mathrm{m}}}\right)-{\mathsf{\psi }}_{\mathrm{m}}\left(\frac{\mathrm{z}-{\mathrm{d}}_0}{\mathrm{L}}\right)+{\mathsf{\psi }}_{\mathrm{m}}\left(\frac{{\mathrm{z}}_{0\mathrm{m}}}{\mathrm{L}}\right)\right]$$

(4)

$${\mathsf{\theta }}_0-{\mathsf{\theta }}_{\mathrm{a}}=\frac{\mathrm{H}}{{\mathrm{k}\mathrm{u}}_{*}{\mathsf{\rho }}_{\mathrm{a}}{\mathrm{C}}_{\mathrm{p}}}\left[\ln \left(\frac{\mathrm{z}-{\mathrm{d}}_0}{{\mathrm{z}}_{0\mathrm{h}}}\right)-{\mathsf{\psi }}_{\mathrm{h}}\left(\frac{\mathrm{z}-{\mathrm{d}}_0}{\mathrm{L}}\right)+{\mathsf{\psi }}_{\mathrm{h}}\left(\frac{{\mathrm{z}}_{0\mathrm{h}}}{\mathrm{L}}\right)\right]$$

(5)

$$\mathrm{L}=-\frac{{\mathsf{\rho }}_{\mathrm{a}}{\mathrm{C}}_{\mathrm{p}}{\mathrm{u}}_{*}^3{\mathsf{\theta }}_{\mathrm{v}}}{\mathrm{kgH}}$$

(6)

where u_* is the friction velocity (m s⁻¹), k is the von Karman’s constant (0.41), z is the height above the surface (m), d₀ is the zero plane displacement height (m), z_0m is the roughness height for momentum transfer (m) estimated using an empirical relationship with NDVI [88], z_0h is roughness height for heat transfer (m), θ₀ is the potential air temperature at surface (K), θ_a is the potential air temperature at z (K), θ_v is the potential virtual temperature near the surface (K), ρ_a is the air density (kg m⁻³), C_p is the specific heat capacity of air (1013 J kg⁻¹ K⁻¹) and g is the gravitational acceleration (9.8 m s⁻²). ψ_m and ψ_h are the stability correction functions for momentum and sensible heat transfer, respectively and L is the Monin–Obukhov length (m).

The scalar roughness height for heat transfer, z_0h, is an important parameter to regulate the heat transfer between the land surface and the atmosphere and estimated as:

$${\mathrm{z}}_{0\mathrm{h}}=\frac{{\mathrm{z}}_{0\mathrm{m}}}{\exp \left({\mathrm{kB}}^{-1}\right)}$$

(7)

where kB⁻¹ is the Stanton number, a dimensionless heat transfer coefficient, estimated using a formulation from [89] as:

$${\mathrm{kB}}^{-1}=\frac{{\mathrm{kC}}_{\mathrm{d}}}{4{\mathrm{C}}_{\mathrm{t}}\frac{{\mathrm{u}}_{*}}{\mathrm{u}\left(\mathrm{h}\right)}\left(1-{\mathrm{e}}^{\frac{-{\mathrm{n}}_{\mathrm{ec}}}{2}}\right)}{\mathrm{f}}_{\mathrm{c}}^2+2{\mathrm{f}}_{\mathrm{c}}{\mathrm{f}}_{\mathrm{s}}\frac{\mathrm{k}\left(\frac{{\mathrm{u}}_{*}}{\mathrm{u}\left(\mathrm{h}\right)}\right)\left(\frac{{\mathrm{z}}_{0\mathrm{m}}}{\mathrm{h}}\right)}{{\mathrm{C}}_{\mathrm{t}}^{*}}+{\mathrm{kB}}_{\mathrm{s}}^{-1}{\mathrm{f}}_{\mathrm{s}}^2$$

(8)

The heat transfer coefficient in Equation (8) was formulated to account for three different land surface conditions. The first term follows the Choudhury and Monteith [90] model for full canopy, the second term accounts for the interaction between the vegetation and soil surface and the third term is for the bare soil surface given [83]. In this equation, f_c and f_s are canopy and soil fraction coverage, respectively, C_d is the drag coefficient for the foliage with a value of 0.2; C_t and C_t^* are the heat transfer coefficients of the leaf and soil, respectively. The value of C_t was taken as 0.03 and C_t^* was computed from Prandtl number and roughness Reynolds number (Re_*) [37]. The u(h) in Equation (8) is the horizontal wind speed at the canopy top (m s⁻¹) and h is canopy height (m) estimated as a ratio of z_0m to 0.136 [37]. The n_ec (within-canopy wind speed profile extinction coefficient) and Brutsaert term kB_S⁻¹ (for bare soil surface) were calculated as:

$${\mathrm{n}}_{\mathrm{ec}}=\frac{{\mathrm{C}}_{\mathrm{d}}\mathrm{LAI}}{\frac{2{\mathrm{u}}_{*}^2}{\mathrm{u}{\left(\mathrm{h}\right)}^2}}$$

(9)

$${\mathrm{kB}}_{\mathrm{s}}^{-1}=2.46{\left({\mathrm{Re}}_{*}\right)}^{0.25}-\ln \left(7.4\right)$$

(10)

where LAI is the leaf area index and estimated as a functional relation with NDVI [91].

SEBS requires estimation of H for dry (H_dry) and wet (H_wet) boundary conditions. Under dry conditions, the H_dry is equivalent to the available energy (R_n − G) as there is no evaporation due to the limitation of water availability and H_wet is calculated using the Penman-Monteith equation [92,93]. After computing H for boundary conditions, the relative evaporative fraction (Λ_r), the evaporative fraction (Λ) and ET are estimated. The steps and explanation are detailed in [37]

2.2.4. Step 4: Extrapolation of Instantaneous to Daily ET

The SEBS uses the Λ approach for scaling instantaneous ET to daily ET, assuming the Λ at the time of overpass is equal to the daily Λ. In this study, a modified approach was implemented where either Λ or ETrF is used for extrapolation of each pixel based on its NDVI value as shown in Equation (15).

$${\mathsf{\Lambda }}_{\mathrm{r}}=1-\frac{\mathrm{H}-{\mathrm{H}}_{\mathrm{wet}}}{{\mathrm{H}}_{\mathrm{dry}}-{\mathrm{H}}_{\mathrm{wet}}}$$

(11)

$$\mathsf{\Lambda }=\frac{{\mathsf{\Lambda }}_{\mathrm{r}}\left({\mathrm{R}}_{\mathrm{n}}-\mathrm{G}-{\mathrm{H}}_{\mathrm{wet}}\right)}{{\mathrm{R}}_{\mathrm{n}}-\mathrm{G}}$$

(12)

$${\mathrm{ET}}_{\mathrm{inst}}=\left(\frac{{\mathrm{R}}_{\mathrm{n}}-\mathrm{H}-\mathrm{G}}{\mathsf{\lambda }}\right)\times 3600$$

(13)

$$\mathrm{ETrF}=\frac{{\mathrm{ET}}_{\mathrm{inst}}}{{\mathrm{ET}}_{\mathrm{r}}}$$

(14)

$${\mathrm{ET}}_{24}=[\mathsf{\Lambda }\times {\mathrm{ET}}_{\mathrm{r}24} \mathrm{for}\mathrm{NDVI}<0.30\left] \mathrm{or} \right[\mathrm{ETrF}\times {\mathrm{ET}}_{\mathrm{r}24} \mathrm{for}\mathrm{NDVI}\ge 0.30]$$

(15)

where ET_inst and ET_r are the actual and reference ET at the hour of satellite overpass (mm h⁻¹), λ is the latent heat of vaporization (~2.45 MJ kg⁻¹). ET_r24 is the daily reference ET and ET₂₄ is the daily actual ET (mm d⁻¹). This modification was made to take the advantage of Λ and ETrF approach to better represent the water limited and energy limited conditions, respectively. The ETrF was estimated as a ratio of ET obtained from Step 3 to reference ET at the satellite overpass time (MODIS Terra satellite overpass local time around 10:30 AM).

2.2.5. Step 5: Filling the Gaps Due to Cloud Cover

Data-gaps due to cloud cover is a common issue in all space-borne satellites. In this study, crop coefficient (Kc) was used to fill the data-gaps. The Kc maps were created for all images as the ratio of ET₂₄ and respective daily ET_r. To fill the Kc of a cloud covered (missing) pixel for a specific image date, the Kc value of the same pixel from the preceding image date was first used. If the same pixel was missing in the preceding image, the Kc value was obtained from the next Kc map. The latter step was repeated if the next day was missing until a date was found with a Kc value estimated for the same pixel. This interpolation method was suitable to fill the data gaps as most of the selected images were less than 10 days apart during the crop growing season (April to October).

2.2.6. Step 6: ET for Longer Periods

After filling the data gaps in daily ET maps due to clouds, the ET maps needed to be created for days when the cloud coverage was more than 10% (or 15%) and thus no input imagery was available. To fill these gaps, the average Kc of the preceding and following images closest to the image date of interest was used. The Kc images were then multiplied with respective daily ET_r to obtain complete time series of daily ET maps. Construction of weekly, monthly, seasonal and annual ET maps was accomplished by summation of daily ET maps over corresponding periods. The processing of all steps was executed in Python language within ArcGIS environment.

2.3. Comparison with Flux Tower Data

Daily ET time series from the modeling framework explained above were compared against observed ET from three flux towers: US-ARc (35.5464 N, 98.0400 W), US-ARb (35.5497 N, 98.0402 W) [94] and US-AR2 (36.6358 N, 99.5975 W) [95]. The US-ARc and US-ARb were located close to each other over native grassland in central Oklahoma. The US-AR2 was located over planted switchgrass in northwest Oklahoma. The 30-minute flux data from the towers were downloaded from the AmeriFlux data archive (http://ameriflux.lbl.gov/) for years 2005 and 2006 for US-ARc and US-ARb sites and years 2010 and 2011 for the US-AR2 site. The flux tower data usually have the issue with energy balance closure, therefore, the closure error was corrected by maintaining constant Bowen-ratio following [96]. The corrected 30-min data were averaged to obtain daily data. The daily observed ET was then compared with the average values of 3 × 3 pixels (~1390 m × 1390 m) from the SEBS ET at the flux tower locations. It should be noted that the three flux towers used for validating the performance of the modeling framework in this study represent only two land covers (native and managed grassland). Hence, the performance of the framework may be different from what is documented here over different types of land covers not included in the present analysis.

For statistical analysis, correlation coefficient (r), the coefficient of determination (R²), mean absolute error (MAE), mean bias error (MBE) and root mean square error (RMSE) were used:

$$\mathrm{MAE}=\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|\mathrm{SEBS}\text{-}\mathrm{ET}-\mathrm{FT}\text{-}\mathrm{ET}\right|$$

(16)

$$\mathrm{MBE}=\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\left(\mathrm{SEBS}\text{-}\mathrm{ET}-\mathrm{FT}\text{-}\mathrm{ET}\right)$$

(17)

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left(\mathrm{SEBS}\text{-}\mathrm{ET}-\mathrm{FT}\text{-}\mathrm{ET}\right)}^2}$$

(18)

where FT-ET is the observed flux tower daily ET and SEBS-ET is the estimated daily ET from the SEBS model.

2.4. Application of the Modeling Framework

After evaluating the accuracy of the modeling framework, it was used to estimate annual ET maps over the entire state of Oklahoma, as well as its nine climate divisions (CD), during the 2001–2014 period. The annual ET were also compared with publicly available MOD16 ET dataset [97,98] over the same period, which covers the most recent drought episode of 2011–2014. The degree of water availability for each pixel and CD within Oklahoma was assessed by estimating the ratio of annual ET from the modeling framework and the reference ET. This ratio is an indication of the portion of the atmospheric demand that is supplied at each pixel and CD. Areas with smaller ratios represent water scarcity since the actual ET from the model is far from the potential limits of ET. The information on annual ET variations and water availability across Oklahoma can assist state water managers with making critical decisions based on long-term objective data from the implemented framework. As mentioned before, the validation dataset only represented native and managed grassland. About half (47%) of all lands in Oklahoma are under rangeland and grassland. With winter wheat being the most dominant crop, the majority of croplands have similar canopy characteristics. Nevertheless, the lack of representation of other land covers (e.g., 21% of forest in Oklahoma) should be considered in applications and interpretations of the results of the modeling framework.

3. Results and Discussion

3.1. Comparison with Flux Tower Data

The comparison with flux tower data showed good agreement between daily SEBS-ET and FT-ET. The modeling approach captured the spatial and temporal variations in ET. However, the model overestimated ET at all sites and years (Figure 3), with average MBE of 20.1 W m⁻². The range of MBE was between 1.7 W m⁻² at US-AR2 in 2011 and 29.3 W m⁻² at US-ARb in 2006 (

Table 1. Statistical indicators between SEBS and flux tower ET.

Site	Year	r	R2	MAE (W m−2)	MBE (W m−2)	RMSE (W m−2)
US-ARc	2005	0.78	0.61	39.6	19.1	40.1
	2006	0.77	0.59	36.7	27.5	49.2
US-ARb	2005	0.81	0.66	31.9	26.6	43.3
	2006	0.78	0.61	35.9	29.3	47.7
US-AR2	2010	0.61	0.37	29.4	16.4	41.7
	2011	0.62	0.39	24.7	1.7	34.1

). The mean MAE and RMSE were 33.0 W m⁻² and 42.7 W m⁻², respectively. The correlation coefficients varied from 0.61 to 0.81 and R² from 0.37 to 0.66.

The errors in the ET estimates of the modeling framework are due to errors generated in each of the six steps outlined in previous section. A major step for error introduction is step three, that is, the surface energy balance model. Previous studies have reported uncertain characterization of kB⁻¹ in water limited environments [99,100,101] and in low vegetation cover conditions [102]. Overestimating kB⁻¹ under these conditions would lead to overestimating z_0h, underestimating H and consequently overestimating ET [99]. The overestimation errors observed in this study were within the range of errors in previous studies when using MODIS as the input imagery to SEBS model. For example, [103] reported ET overestimation with MBE of 6.1 W m⁻²; [104] found MBE of 20.1 W m⁻² and RMSE of 34.7 W m⁻²; [105] reported MBE of 144.9 W m⁻² when comparing SEBS-ET from cropland and grassland with flux tower estimates; [106] found overall MBE of 31 W m⁻² and RMSE of 76 W m⁻²; and, [107] reported MBE of 95.1 W m⁻² and RMSE of 122.2 W m⁻² across several land covers and climatic conditions. While several studies have reported overestimation error from SEBS, the mean absolute error from the current study was smaller than the threshold of 50 W m⁻² suggested by [108].

Errors in other steps of the framework can contribute to biases in final ET estimates. A common source of error in estimating ET from satellite imagery is due to cloud contamination. A thin layer of cloud or a shaded area due to cloud presence over nearby pixels can result in underestimation of LST and consequently, overestimation of ET. In practical applications, it is impossible to remove all these contaminated pixels from the entire image even after applying the LST thresholds during quality control. In this study, there were days with underestimated LST due to cloud presence. For example, the LST at the flux tower pixel area dropped by 10.6 K from Day of Year (DOY) 113 to 114, while both DOYs were identified as cloud-free and no precipitation was recorded. The instantaneous T_A increased by 3.2 K over the same period. The smaller LST on DOY 114 affected ET estimation for this day and the following days until another cloud-free image was obtained for DOY 117 (Figure 3a).

A sensitivity analysis study [109] on SEBS model reported LST as the most sensitive parameter, with up to 70% error in H from irrigated fields expected with 0.5 K bias in LST. Another study [110] found that error in H varied between −41% and 152% when LST bias ranged from −4 K to 10 K. These studies show that a small bias in LST can significantly impact H and ultimately ET. The magnitude of error may depend upon the sensitivity of SEBS to LST, including other parameters such as T_A, u, ∆t [111] and could vary depending on whether the wet or dry limits have been reached [110,112].

In this study, the filtering criteria of less than 15% cloud cover limited the availability of cloud-free images. Applying this filter resulted in 125 and 154 cloud-free images for processing during 2005 and 2006, respectively. For the days with no cloud-free images, the ET estimate was dependent on the Kc approach explained before. However, the Kc approach may fail to account for the variability in pixel conditions, especially if land and weather conditions change dramatically during long periods of gaps in imagery. In this study, 10 to 15 cloud-free images each month were available for most months, which was assumed sufficient to capture general daily soil moisture and weather variations. In other periods, however, it was not possible to keep the length of gap periods short. For example, cloud-free images were not available for 17 consecutive days from DOY 270 to 286 in 2005, when larger differences between FT-ET and SEBS-ET were observed (Figure 3a,c).

The combined impact of LST bias due to cloud contamination and unavailability of cloud-free images significantly increase biases in ET estimation. The 15% cloud cover filter could be reduced to reduce cloud contamination issue but this would come at the cost of increasing the length of periods with no imagery at all. Increasing the filtering limit will have an opposite effect (more available imagery with larger cloud contamination within each image). Another solution is to manually inspect and select images. However, this increases the processing time and interrupts the automated nature of the ET modeling framework. Another factor that could play a significant role in increasing ET errors is the availability and quality of input weather data. Su et al. [113] reported about 40% increase in RMSE (from 73 W m⁻² to 102 W m⁻²) when using reanalysis dataset—Global Land Data Assimilation System within SEBS instead of ground-based weather data. In this study, the impact of this source of error is expected to be minimal since rigorous quality control was conducted on ground-based data and only less than 2% of data were missing during the study period.

As highlighted before, the daily ET results, uncertainty and potential biases of the proposed ET modeling framework were evaluated and discussed based on flux tower measurements over native and managed grassland at central and northwest Oklahoma. Flux tower data across other land covers were not available for comparison, thus the results from the framework may need further assessment to warrant the similar level of accuracy and uncertainty while applying the results to different land covers and climates across the state. In particular, the analysis and interpretation of results from current study may differ for vegetation with different canopy structure compared to grassland.

3.2. Application of the Modeling Framework

The automated operational ET modeling framework proposed in this study was used to create annual ET maps covering the entire state of Oklahoma for the period from 2001 to 2014. As expected, the annual ET followed the precipitation pattern and increased from southeast to Panhandle (Figure 4). When averaged over the entire 14 years, the southeast climate division (CD9) had the largest annual ET of 1272 mm yr⁻¹ and the Panhandle climate division (CD1) had the smallest annual ET of 588 mm yr⁻¹ (

Table 2. Average annual SEBS-ET, MOD16-ET, ET_r and the ratio of SEBS-ET to ET_r for all Oklahoma climate divisions (CD) during the 2001–2014 period.

Climate Division	SEBS-ET (mm yr−1)	MOD16-ET (mm yr−1)	ETr (mm yr−1)	SEBS-ET ETr−1
CD1 (Panhandle)	588	259	2140	0.27
CD2 (North Central)	918	364	1871	0.49
CD3 (Northeast)	1098	657	1521	0.72
CD4 (West Central)	790	338	2018	0.39
CD5 (Central)	1095	531	1700	0.64
CD6 (East Central) *	1175	736	1492	0.79
CD7 (Southwest)	845	363	2009	0.42
CD8 (South Central)	1163	599	1683	0.69
CD9 (Southeast) *	1272	798	1360	0.94
Oklahoma	994	516	1755	0.57

* These CDs have a forested area of more than 29%. The results presented in this table may not be accurate for these CDs since the flux towers used in validation did not include forest land cover.

). The reference ET (ET_r) had an opposite pattern, with CD1 having the largest amount at 2140 mm yr⁻¹ and CD9 the smallest (1360 mm yr⁻¹). This means that on average, about 94% of atmospheric demand was fulfilled at southeast, compared to only 27% in the Panhandle during the study period. In other words, water scarcity is a larger issue in CD1 compared to CD9 as available resources were not sufficient to keep up with atmospheric demand. The statewide average annual ET was 994 mm yr⁻¹, about 57% of the average annual ET_r of 1755 mm yr⁻¹.

The average annual ET comparison between MOD16 and SEBS indicated large differences across all Oklahoma CDs (Table 2). The differences between MOD16-ET and SEBS-ET varied between 37% at CD9 to 60% at CD2, with an average of 48% lower ET rates from MOD16. Three eastern humid CDs (CD3, CD6, CD9) had smaller differences between MOD16-ET and SEBS-ET compared to three western CDs (CD1, CD4, CD7). The difference between SEBS-ET and MOD16-ET is possibly due to a combination of overestimations from SEBS and underestimation from MOD16. The underestimation of ET from MOD16 has been reported in previous studies, particularly in semi-arid and arid climates [114,115].

The ratio of SEBS-ET to ET_r can be estimated on a pixel wise basis to provide information on water scarcity at a finer resolution for local water management and planning. This ratio is mapped in Figure 5. The general patterns are similar to those presented in Table 2, with western parts of the state under relatively larger water scarcity compared to the eastern parts. However, significant variability can be observed within some CDs. In CD1, for example, the western half of CD (Cimarron and Texas counties) had smaller ratios compared to the eastern half, suggesting a more severe water scarcity. CD2 was similar in terms of variations in the ET ratios across the CD. The surface water resources in western Oklahoma were visible in regions with a ratio value of more than 0.5. Examples include the riparian areas of Cimarron and North Canadian rivers in southwest of CD2, as well as Canton Lake and Foss reservoir in CD4 and the five reservoirs in CD7 (Lugert-Altus, Tom Steed, Lawtonka, Ellsworth and Fort Cobb). Maps similar to the one in Figure 5 can be developed at varying temporal and spatial scales to monitor changes in water availability more closely. The ratio of actual ET to reference or potential ET has been used in the past in monitoring water stress and drought, such as in the Evaporative Stress Index [17].

The inter-annual variations in ET were also examined for each CD and for the entire state. Figure 6 demonstrates deviations in SEBS-ET as percentage of the average annual ET during the 2001–2014 period. The impact of the 2011–2014 drought in western Oklahoma can be observed in this graph, with the maximum reduction in ET occurring in 2011 for the three western CDs of CD1, CD4 and CD7. The percent deviations from average was −22%, −21% and −33% for the same CDs, respectively. According to the U.S. Drought Monitor (USDM) [116], more than 80% of the three CDs was under extreme drought (D3 category) from June. The drought condition worsened in July and remained under D4 category until December 2011. The three eastern CDs of CD3, CD6 and CD9 were above average in 2011, with percent deviations of 9%, 8% and 10%, respectively. The USDM indicated almost no drought at CD3 in 2011, whereas CD6 and CD9 had less than 40% of their area under extreme drought from August to November 2011. The middle three CDs registered close to long-term average ET. The largest positive deviations for the three western CDs occurred in 2007, a year that was characterized by above normal precipitation.

The ET modeling framework proposed in this study can automatically generate time series of daily ET maps on a continuous basis, with several applications beyond those mentioned in previous sections. For example, ET maps over agricultural areas can be analyzed in conjunction with yield data to evaluate the water use efficiency. However, this modeling framework has some limitations that must be considered and improved in future applications. One limitation is the size of MODIS pixels, which practically hinders the possibility of using the ET data at field scale. This limitation can be overcome by modifying the framework to use satellite imagery at finer resolution (e.g., Landsat). Another challenge is identifying and removing cloud contaminated pixels. The filters used in this study were not always effective in identifying pixels that were covered by thin layers of cloud or were in the shadow of a cloud. Thus, further investigation and application of robust methods to examine cloud contamination are needed. Finally, there were periods when no images were available for several days due to clouds covering the entire scene. This negatively affects the ability to capture ET fluctuations during those periods. Data-fusion approaches can be implemented in the modeling framework as a potential solution to improving ET interpolation for days with missing images.

4. Conclusions

An ET modeling framework was proposed to automatically construct daily time series of ET maps across Oklahoma by integrating MODIS imagery, ground-based weather data and surface energy balance model. The comparison of the results with daily observations at three flux towers (two years of data at each site) showed good performance of the modeling framework with mean bias errors less than 30 W m⁻² and root mean squared errors less than 50 W m⁻². The results were then used to investigate spatial and temporal variations in ET across the state and its nine climate divisions (CD). The statewide annual ET varied between 841 and 1100 mm yr⁻¹ during the period from 2001 to 2014, with an average of 994 mm yr⁻¹. A large difference in ET was observed among CDs, with Oklahoma Panhandle (CD1) having the smallest and southeast (CD9) the largest average annual ET of 588 and 1272 mm yr⁻¹, respectively. The ratio of estimated ET to reference ET was used as an indicator of water scarcity at pixel and CD levels. The deviations in annual ET from the 2001–2014 average ET were also studied and found to be in good agreement with temporal and spatial variations in drought. The proposed ET modeling framework provided a pathway to construct daily time series of ET maps with potential for a range of applications. However, further improvements are necessary to resolve the issues highlighted in the current study.