Application of A Simple Landsat-MODIS Fusion Model to Estimate Evapotranspiration over A Heterogeneous Sparse Vegetation Region

Abstract

A simple Landsat-MODIS (Moderate Resolution Imaging Spectroradiometer) fusion model was used to generate 30-m resolution evapotranspiration (ET) maps for the 2010 growing season over a heterogeneous sparse vegetation, agricultural region using the METRIC (mapping evapotranspiration with internalized calibration) algorithm. The fusion model performance was evaluated, and experiments were undertaken to investigate the frequency for updating Landsat-MODIS data into the fusion model during the growing season, to maintain model accuracy and reduce computation. Initial evaluation of the fusion model resulted in high bias stemming from the landscape heterogeneity and small landholdings. To reduce the bias, the fusion model was modified to be applicable pixel-wise (i.e., implementing specific pixels for generating outputs), and an NDVI-based (Normalized Difference Vegetation Index) coefficient was added to capture crop phenology. A good agreement that resulted from the comparison of the fused and non-fused maps with root mean square error (RMSE) of 0.15 mm day⁻¹ with coefficient of determination (R²⁾ of 0.83 indicated successful implementation of the modifications. Additionally, the fusion model performance was evaluated against in-situ observation at the pixel level as well as the watershed level to estimate seasonal ET for the growing season. The default METRIC model (Landsat only) yielded relative error (RE) of 31% and RMSE of 2.44 mm day⁻¹, while using the modified fusion model improved the accuracy resulting in RE of 3.5% with RMSE of 0.37 mm day⁻¹. Considering different data frequency update, the optimal fusion experiment (RMSE of 0.61 mm day⁻¹, and RE of 6.5%) required the consideration of the crop phenology and weekly updates in the early growing stage and harvest time, and bi-weekly for the rest of the season. The resulting fusion model for ET output is planned to be a part of ET mapping and irrigation scheduling systems.

1. Introduction

The motivation of this study stems from the requirement of high frequency and high-spatial resolution evapotranspiration (ET) monitoring in agricultural areas of semi-arid regions. Water scarcity in these regions poses challenges for maintaining environmental stability [1] as well as geopolitical security. Southern Iran is indicative of a growing example of regions where excessive pumping groundwater occurs, and which faces water scarcity. The combination of small agricultural land holdings, uncultivated areas and scattered pasture create a mosaic of a highly heterogeneous landscape. In addition to the dearth of information about crop water demand in these areas, in-situ data is also not available. Quantifying actual evapotranspiration in high resolution, and high frequency for assessing irrigation scheduling and prudent use of water resources will help the regional sustainability. Although heterogeneous landscape exacerbates the situation, using fusion ET model helps develop regional ET maps desired by the water managers.

Quantifying water budget components; in particular, rainfall and evapotranspiration, help improve water management and irrigation scheduling. While recent years have witnessed growth in rainfall monitoring using in-situ and satellite measurements [2], the availability of similar high spatio–temporal ET data at regional and local scale is still a challenge. Satellite data and remote sensing-based algorithms utilizing surface energy balance concepts facilitated generating spatial ET estimates [3,4,5,6]. The surface energy balance algorithm for land (SEBAL) [3] and its variant, mapping evapotranspiration with internalized calibration (METRIC) [5], have been extensively used and well established for assessing actual ET estimates [7,8,9,10]. While SEBAL underestimates ET under advective conditions [10,11,12], METRIC can account for advection and perform accurately in heterogeneous land use when proper extreme pixels (hot and cold pixels) are selected [13,14]. The METRIC algorithm can be used with Landsat data which are available at 30 m resolution and satisfies the need for high-spatial resolution. Despite the high-spatial resolution, Landsat data is insufficient due to the low-temporal resolution with a typical spanning interval of 16 days. This interval can be even longer under cloudy conditions [15]. Nevertheless, water use and biomass information for decision-making systems require a framework that generates daily ET estimates at field or subfield scales. Other satellites products with daily data coverage (e.g., MODIS (Moderate Resolution Imaging Spectroradiometer), geostationary satellites) have a coarse spatial resolution but can provide the required temporal continuity. Therefore, a data fusion approach that combines the spatial resolution of Landsat with the temporal frequency of MODIS helps provide daily field-scale ET estimates and can benefit local and regional decision-making framework.

The data fusion method was traditionally presented in the form of methods, such as the intensity-hue-saturation (IHS) [16], principal component analysis (PCA) [17], and high pass filter (HPF) [18], for integrating spatial characteristics of high-resolution panchromatic data with coarse resolution multispectral observation for generating new multispectral images [19,20]. These methods could only enhance the spatial resolution of the fused data, because using the panchromatic band was not effective in enhancing spatial and temporal resolution simultaneously [21]. To generate fused data with high spatio–temporal resolution Gao et al. [22] introduced the spatial and temporal adaptive reflectance fusion model (STARFM). The model blends MODIS and Landsat data and assumes a constant land cover type between observational periods. The model has been widely used in several studies to generate fused reflectance data [3,23,24]; however, the accuracy of the model depends on the availability of pure homogenous MODIS pixels, and the result can be misleading when applied in heterogeneous landscapes and small-scale farmlands [22,25]. The STARFM algorithm has been modified in a number of studies [25,26,27]. Weng et al. [28] used linear spectral mixing analysis to improve the original STARFM model and proposed the spatio–temporal adaptive data fusion algorithm for temperature mapping (SADFAT). This method required three consecutive MODIS and two Landsat images which made it unfeasible in many cases [29]. The spatial–temporal adaptive algorithm for mapping reflectance change (STAARCH) [25] and the enhanced STARFM method (ESTARFM) [27] were two other examples of trials to improve the original STARFM to detect land cover change and make it applicable in heterogeneous landscape. However, the difficulty in input data preparation (e.g., certain pairs of fine and coarse resolution images acquired on the same date) significantly limited their applicability [30].

While multi-sensor data fusion approaches have been widely used in reconstructing spectral reflectance for crop monitoring and phenology analysis [31,32], mapping forest-cover change and disturbance [25,33], environmental process detection and the estimation of biophysical parameters [23,34], the application of multi-sensor data in downscaling and generating fused ET mapping were limited and complex [35]. Hong et al. [36] suggested that MODIS to Landsat downscaling produces better results when regression-based approaches are employed. Singh et al. [37] developed the linear with zero intercept (LinZI) method for downscaling MODIS-based ET maps to the Landsat resolution scale. The regression model was not effective when the slope and intercept were calculated for a large area. Another data fusion method was proposed by Cammalleri et al. [35] to spatially disaggregate ET estimates, generated using the SEB (Surface Energy Balance) model and geostationary satellite data, with MODIS and Landsat imagery. The method was successfully applied over corn, soybean, and cotton fields [38,39], as well as irrigated vineyards [15]. Other approaches were also proposed for fused ET mapping using machine learning techniques [40] or modifying the STARFM model for ET applications [41,42]. Although these methods have a strong theoretical foundation and provide good estimates of fusion data in fine resolution, they are also computationally time-consuming and require the application of relatively complex algorithms. Therefore, such methods are often not used by local agencies and managers who look for simple techniques that can provide daily ET maps with acceptable or known uncertainty. A relatively simpler fusion model was proposed by Bharatti et al. [43] to generate daily and seasonal ET at a 30 m scale based on the high correlation reported between remotely-sensed vegetative indices and evapotranspiration [44,45]. Although a good performance of the model was reported in predicting crop water use [43], it needs to be analyzed and evaluated over other regions. More specifically, since there are a limited number of studies that evaluated fusion application under heterogeneous landscape of semi-arid areas with sparse vegetation and small landholdings [41], evaluation of simpler models in these regions is of interest.

In this study, because of chronic water shortage and the importance of water, weekly ET maps were desired by local water management agencies. The METRIC model was previously applied using Landsat imagery for estimating crop water use during the summer growing season of 2009–2010 and showed good accuracy with 4.7% error [10]. The model, however, showed high uncertainty (>30%) in the winter growing season due to the low-temporal resolution of available Landsat imagery. The low-temporal resolution occurred because of the lack of Landsat imagery particularly due to sensor malfunction and cloudy days. Therefore, a MODIS-Landsat fusion algorithm was utilized to generate ET maps with high spatio–temporal resolution. The feedback received from the users was that the data resolution was useful, but the fusion approach needed to be computationally efficient for routine/practical application. Building on this feedback, in this paper we report on the evaluation of the MODIS-Landsat fusion model [43]; the need for modifying it for the study region, and the generation of ET maps to estimate seasonal water use for local farmland over Gareh Bygone Plain (GBP) of southern Iran. The study period was the winter of 2010 for which in-situ measurements were available. The modifications to the fusion model were done with regards to the pixel selection to achieve computational efficiency and to reduce bias. Additionally, 24 experiments, considering the different frequencies required for updating Landsat-MODIS data into the fusion model, were undertaken to obtain the optimal strategy that could maintain model accuracy and reduce computations.

2. Materials and Methods

2.1. Study Site

The study area is Gareh Bygone Plain (GBP) which is located in southern Iran at 28°35′N and 53°53′E (Figure 1). The GBP covers an area of 18,000 ha and has mean annual precipitation of 219 mm that typically occurs between December and March. Air temperature reaches up to 46 °C in summer and drops to between −1 and −6 °C in winter. The land cover over GBP is sparse, with mixed landscapes, such as rangelands, shrublands, agricultural croplands, and bare soil. The landscape and climate of the study area is broadly representative of the southern/southeast Iran region. In the study area, corn and wheat are the dominant crop types during winter as well as summer growing seasons. The winter growing season starts at the end of June/early July and lasts through November/early December.

2.2. The METRIC Algorithm

The METRIC model is based on the surface energy balance equation [5]:

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

(1)

where LE is the latent energy consumed by ET (W m⁻²); R_n is the net radiation (W m⁻²); G is the soil heat flux (W m⁻²), and H is the sensible heat flux (W m⁻²). The model computes G as a ratio G/R_n using an empirical equation [5] as follows:

$$\frac{\mathrm{G}}{{\mathrm{R}}_{\mathrm{n}}}=({\mathrm{T}}_{\mathrm{s}}-273.15)\times (0.0038+0.0074\mathsf{\alpha })\times (1-{\mathrm{NDVI}}^4)$$

(2)

where T_s is surface temperature (K); $\mathsf{\alpha }$ is surface albedo and NDVI is the Normalized Difference Vegetation Index. Sensible heat flux is a function of air density ($\mathrm{P}$_air, kg m⁻³), the specific heat constant of air at constant pressure (C_p = 1004 J kg⁻¹ K⁻¹), the aerodynamic resistance (r_ah, s m⁻¹), and the temperature difference (dT, K), and computed as:

$$\mathrm{H}={\mathrm{P}}_{\mathrm{air}}\times {\mathrm{C}}_{\mathrm{p}}\times \frac{\mathrm{dT}}{{\mathrm{r}}_{\mathrm{ah}}}$$

(3)

The model computes LE based on other components (i.e., H, R_n, and G), and the instantaneous ET (ET_inst) is then determined as:

$${\mathrm{ET}}_{\mathrm{inst}}=3600\frac{\mathrm{LE}}{\mathsf{\lambda }}$$

(4)

The instantaneous ET fraction (ET_rF) is then estimated as ET_inst divided by the reference ET (ET_ref) at the image acquisition time.

$${\mathrm{ET}}_{\mathrm{rF}}=\frac{{\mathrm{ET}}_{\mathrm{inst}}}{{\mathrm{ET}}_{\mathrm{ref}}}$$

(5)

The ET_rF is typically considered to be constant throughout the day [5]. Combining daily reference ET (calculated from weather data) with ET_rF (calculated from METRIC) results in actual daily evapotranspiration for the days after the Landsat overpass day until the next available image [5].

2.3. The MODIS–Landsat Fusion Model

We derived ET maps considering a growing season using the MODIS–Landsat ET fusion scheme [43]. The scheme applies linear regression to integrate Landsat-derived ET fraction using the METRIC model, and the VTCI (Vegetation Temperature Condition Index) derived from the MODIS imagery. The VTCI is derived from the plot of surface temperature and vegetation index, and the trapezoidal shaped distribution (Figure 2), and can be used to describe the water stress condition of a pixel. The VTCI is calculated using the following equations:

$$\mathrm{VTCI}=\frac{{\mathrm{LST}}_{\max }-\mathrm{LST}}{{\mathrm{LST}}_{\max }-{\mathrm{LST}}_{\min }}$$

(6)

$${\mathrm{LST}}_{\max }={\mathrm{a}}_1+{\mathrm{b}}_1\times \mathrm{NDVI}$$

(7)

$${\mathrm{LST}}_{\min }={\mathrm{a}}_2+{\mathrm{b}}_2\times \mathrm{NDVI}$$

(8)

In the above, LST is land surface temperature corresponding to an NDVI value, and the LST_max and LST_min are maximum and minimum land surface temperatures, respectively, from the available satellite maps over the study domain. The coefficients (a₁, b₁, a₂, and b₂) are estimated from linear regression utilizing maximum and minimum values of NDVI and the associated LST values within the image domain.

The maps of VTCI corresponding to two days ‘a’ and ‘b’ were generated based on Equation (6). from the available MODIS imagery. Linear regression was performed between VTCI_a and VTCI_b for a pixel ‘i’ and a specific land cover as:

$$(\mathrm{VTCI}{)}_{\mathrm{b}.\mathrm{i}}=\mathrm{e}\times (\mathrm{VTCI}{)}_{\mathrm{a}.\mathrm{i}}+\mathrm{f}$$

(9)

Assuming that ET values are available at day of year ‘a’ (DOY_a), the model estimated ET_rF for DOY_b is developed assuming a linear relationship between ET_rF and the VTCI as:

$$({\mathrm{ET}}_{\mathrm{rF}}{)}_{\mathrm{b}.\mathrm{i}}=\mathrm{e}\times ({\mathrm{ET}}_{\mathrm{rF}}{)}_{\mathrm{a}.\mathrm{i}}+\mathrm{f}$$

(10)

where, (ET_rF)_b,i is the ET fraction values at 30 m resolution for a pixel ‘i’ at DOY_b using METRIC and Landsat imagery; ‘e’ and ‘f’ are coefficients from the linear regression of MODIS-based VTCI_a and VTCI_b, respectively.

For an available Landsat scene corresponding to DOY_a, a high-resolution ET_rF map is generated using the METRIC model. This map is then used in a linear regression equation to generate a high-resolution ET_rF map corresponding to DOY_b using the coefficients derived from the regression of MODIS VTCIs for DOY_a and DOY_b. This is the standard approach for using the fusion model and was applied over the study region. The model algorithm, as well as the modifications (described in the next section), are summarized in Figure 3.

2.4. Modifications to the Fusion Model

The fusion model has been successfully applied in prior published studies (e.g., [15,35,43]). For example, the fusion model ET was compared with eddy covariance measurements at a field site in Florida and RMSE (Root Mean Square Error) of 0.44 mm day–1, R² of 0.80, and 2% relative error for the cumulative growing season were reported [43]. Similar good performance was reported using MODIS–Landsat fusion model over corn, and soybean field sites in Iowa, with an average error of 0.58 mm day⁻¹ [35]. While, for a vineyard, the error has been about 1 mm day⁻¹ [15]. Indeed, this good performance reported in the literature was the motivation for pursuing the fusion modeling assessment for the study region.

However, initial evaluation of the model when applied over the study region resulted in significant bias. Reviewing the model spatial plots, it was apparent that the bias was likely a result of the landscape heterogeneity. As the spatial resolutions of the satellite sensors varied from 30-meter in Landsat to 1-kilometer in MODIS, and the fusion model needed input parameters from both resolutions, adequately dealing with the spatial scale issue was a challenge. The issue would be trivial over a homogeneous surface, but for a heterogeneous landscape, the difference in the satellite data products and measurements (e.g., NDVI data) can yield different values over the same surface with different resolution. Examples of the NDVI maps developed using Landsat and MODIS data over the study domain are presented in Figure 4.

The study region featured a heterogeneous landscape with small landholdings, and the sub-pixel variability considering MODIS and Landsat data introduced high uncertainties in calculations. Three different spatial resolutions were used in the fusion model’s inputs: MODIS LST at 1-km, MODIS NDVI at 250-m, and Landsat data at 30-m resolution. Incidentally, the different grid spacing (resolution) between MODIS NDVI and LST was another source of bias in calculating MODIS VTCIs in Equations (6)–(8). A green pixel with high NDVI value and neighboring fallow land pixels with low NDVI can be combined in one corresponding LST pixel representing the temperature of a mixed landscape. However, this can result in a low correlation between MODIS NDVI and LST used to derive the coefficients in Equations (6)–(8).

To reduce the bias stemming from different resolutions under a heterogeneous landscape, we applied a selective filter approach. First, the MODIS data were resampled to 30-m resolution using the nearest neighbor method. The filter scans a MODIS–Landsat pair for a specific DOY and compares the NDVI values from Landsat scene with the neighboring pixels from the 30-m resampled MODIS image. It only selects the areas and pixels that have the NDVI values as obtained from MODIS and Landsat within a pre-specified variation. The variation is a user-specified input and was taken as 20% based on average value during trial and error procedure, and evaluating the values presented in the literature as a difference threshold between MODIS and Landsat vegetation indices. A similar value was reported by Sesnie et al. [46] in the Sonoran Desert whose vegetation cover and its sparsity resembled some parts of our study area. Tong and He [47] reported an R² range from 0.33 to 0.38 between MODIS and Landsat-5 NDVI in grassland at a national park in Canada. Ke et al. [48] compared Landsat 8 with multiple satellite sensors and reported that NDVI values have a mean bias error of between 5% and 10% for four study sites in China and Korea that covered forest, crop, urban, and water body cover. Busetto et al. [49] compared Landsat-5 and MODIS NDVI at an agricultural landscape of Northern Italy and found the RMSE between the two ranged from 0.065 to 0.094 (6.2% to 21.4% of RE) for different dates. That is, we considered 20% to be the implicit bias when comparing the Landsat and MODIS fields further justifying the pre-specified variation threshold. The filter scans spatio–temporally through the growing season and applies the similarity requirement for all available MODIS–Landsat pairs to capture the effect of changes in crop phenology. These sub-areas and pixels were then used as the input to the fusion model to derive the required coefficients. A schematic of this process is shown in Figure 5.

Due to the finer NDVI resolution as compared to the LST in MODIS data, a relatively smaller bias is expected for the former over the heterogeneous landscapes. Recognizing this advantage, and considering the heterogeneity due to small landholdings in the study region, an NDVI-based coefficient was added to the MODIS VTCI calculations. As a result, Equation (6) was updated by considering the ratio of average NDVI of DOY_b to average NDVI of DOY_a for the calculation of VTCI_b as:

$$\mathrm{VTCI}=(\frac{{\mathrm{NDVI}}_{\mathrm{avg}-\mathrm{b}}}{{\mathrm{NDVI}}_{\mathrm{avg}-\mathrm{a}}})\frac{{\mathrm{LST}}_{\max }-\mathrm{LST}}{{\mathrm{LST}}_{\max }-{\mathrm{LST}}_{\min }}$$

(11)

The modified fusion model then was applied over the study domain. Experiments were conducted by considering different strategies in updating MODIS–Landsat input pairs. These aspects are described in the following section.

2.5. Frequency Update Experiments

As a default approach, we used four Landsat imagery coupled with bi-weekly (16-day) MODIS data updates for generating the 30 m ET maps with the fusion model. At the base level, the fusion model requires a pair of MODIS images (that can be selected across any time window during the growing season) and one Landsat images. The pair of MODIS images is used to derive coefficients in Equations (6)–(9), and the coefficients were used with Equation (10). This, paired with the ET_rF resulted from Landsat and METRIC model, was used to generate a fusion ET map at 30 m pixel size for the days in between the availability of MODIS pair. Accordingly, experiments were undertaken to select different dates to update MODIS LST, NDVI data from which the coefficients were calculated to be used with the available Landsat imagery. For example, one experiment could consider DOY_a as the beginning of the growing season and DOY_b at the harvest time to generate the coefficients and use it with the ET_rF map at the DOY_a to generate the fusion ET for the whole season between DOY_a and DOY_b. Alternatively, one can consider more frequent updates of MODIS and Landsat data for the fusion model. A total of 26 set of experiments were undertaken for updating satellite data for the fusion model. These included updates with Landsat only data (non-fusion approaches) or using one to four Landsat images and considering MODIS updates for every 8-, 16-, 24-, 32-, 40-, and 48-day interval. We also considered one hybrid frequency update by using weekly (8-day) updates at an early stage and harvest time and bi-weekly (16-day) during maturity. These experiments aimed to identify a practical pairing that yielded high accuracy and relatively low data processing effort.

2.6. Data and Model Evaluation

2.6.1. Data

Four cloud-free Landsat scenes between July and December 2010 corresponding to path/row 161/40 and DOY 189, 285, 301, and 325 were obtained from the USGS (United States Geological Survey) global visualization viewer (http://glovis.usgs.gov/). Landsat 5 TM imagery includes six bands from the visible to mid-infrared range with 30 m resolution (bands 1–5 and 7) and a thermal band with 120 m resolution (band 6). The USGS resamples the thermal band and distributes it at 30 m pixel size; therefore, the 30 m thermal band was used in this study. A subset of each image was extracted to focus on the study site for the analysis (Figure 1). The digital elevation model (DEM) was a required input to the METRIC model and was processed from ASTER (Advanced Spaceborne Thermal Emission and Reflection radiometer) data. The fusion model also required LST and NDVI as input, and were obtained from MODIS. The MODIS level 3 LST data at 1 km resolution (MOD11A1) and NDVI with 250 m resolution (MOD09GA) were retrieved from publicly available sources (https://modis.gsfc.nasa.gov). In addition, daily weather data were available from meteorological stations at the Gareh Bygone and the Fasa weather observatories (http://farsmet.ir/Amars.aspx). Pan evaporation data from the weather station were used with pan coefficient (0.75 for the station) to generate reference ET for the study region, which was consistent with the ET_ref resulting from ASCE Penman–Monteith standardized method [50] as reported by Jamshidi et al. [10] for the study region.

The in-situ ET measurements are reported by Pakparvar [51] for the 2009 and 2010 growing seasons. Availability of these measurements was one of the main reasons we chose the 2010 growing season as our test period in this observational data void region. The experiments were conducted over several farms to measure the water use by the cultivated crops (corn) throughout the growing seasons. For this measurement, the water balance approach was adopted by direct measurements of the amount of irrigation, precipitation, soil water content, and drainage water. Pakparvar et al. [52] used the measurements to calculate an average volume of water loss through evapotranspiration per hectare of agricultural fields. Knowing the total acreage of cultivated farmlands in the study domain, the cumulative seasonal crop water use was computed. We used the data provided by Pakparvar et al. [52] for comparing the result at the regional scale (i.e., cumulative water use) and for calculating the crop coefficients during the growing stage. The coefficients were compared against the available measurement-based coefficients for similar climatic region of the province [52,53,54,55] as well as against the remotely sensed results [10]. These coefficients were then used with the reference ET values to generate daily time series of actual ET for the growing season. A summary of the vegetation characteristics and the measurements is presented in

Table 1. Characteristic of the cultivated farmlands.

Crop Type	Planting Date	Harvesting	Average Irrigation Depth (mm)	Seasonal Crop Water Use (Mm3)	Initial Kc	Mid-Season Kc	End-Season Kc
Corn	Early July	Early December	1000	5.44	0.52	1.15	0.40

.

2.6.2. Model Evaluation

To evaluate the fusion model, two approaches were adopted. First, the ET estimates were evaluated at the pixel level for the study area. Applying the fusion model resulted in the generation of 30 × 30 m² ET maps which were then integrated to be used with in-situ ET measurements. The ET value from each pixel corresponding to agricultural land cover was extracted from the fusion maps. These values were averaged spatially by using arithmetic mean of all the extracted values to generate an average value for actual ET corresponding to each specific day during the growing season. These values were used to generate a time series of daily actual ET for the growing season. To compare and visualize the resulting fusion ET time series, these values were plotted against ET times series that were generated based on the in-situ measurements in the study area [51,52] described in the previous section. Additionally, the resulting ET maps from the fusion model were compared with those resulting from using only the METRIC model with no integration of MODIS products. To compare daily ET values resulting from the two maps (i.e., the fusion model and the non-fusion maps), two categories including cultivated areas (agricultural farmlands) and sparse vegetation (shrubland and pasture) were defined. The pixel values corresponding to each category were extracted, averaged, and compared.

Second, building on the availability of in-situ measurements including crop water use at the study area for the growing season [51,52], we evaluated the ET estimates integrated across the watershed in the study area. For this, the ET values corresponding to agricultural farmlands from the resulting fusion maps were extracted and summed over the entire season. The cumulative value was then converted to the volume of crop water use by multiplying it by the acreage of cultivated area. The calculated value was then compared against the cumulative crop water used from the in-situ measurements.

To quantify the performance of the fusion model in terms of ET estimates, we computed the coefficient of determination (R²). Additionally, the error components including RMSE (root mean square error), and RE (relative error) were also used to compare the results from the model.

$${\mathrm{R}}^2=1-\frac{(\sum ({\mathrm{O}}_{\mathrm{i}}-{\overline{\mathrm{O}}}_{\mathrm{i}})({\mathrm{P}}_{\mathrm{i}}-{\overline{\mathrm{P}}}_{\mathrm{i}}{)}^2}{\sum ({\mathrm{O}}_{\mathrm{i}}-{\overline{\mathrm{O}}}_{\mathrm{i}}{)}^2\sum ({\mathrm{P}}_{\mathrm{i}}-{\overline{\mathrm{P}}}_{\mathrm{i}}{)}^2}$$

(12)

$$\mathrm{RMSE}=\sqrt{\frac{\sum ({\mathrm{P}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}{)}^2}{\mathrm{n}}}$$

(13)

$$\mathrm{RE}=100\frac{\sum (|{\mathrm{P}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}|}{{\mathrm{O}}_{\mathrm{i}}}$$

(14)

where, ${\mathrm{O}}_{\mathrm{i}}$ and ${\mathrm{P}}_{\mathrm{i}}$ represent the observed and estimated ET of the i_th pixel, respectively, and ${\overline{\mathrm{O}}}_{\mathrm{i}}$ and ${\overline{\mathrm{P}}}_{\mathrm{i}}$ represent average values of the corresponding variable, n is the number of data samples considered.

3. Results and Discussion

3.1. Initial Results with the Original and Modified Models

Initial work with the fusion model highlighted an issue in estimating ET values over the study domain with heterogeneous landscapes. Figure 6a and Figure 7a show a sample correlation for LST–NDVI and the resulting VTCI using the default approach (i.e., before any modifications). Using such a framework in regions with mixed agriculture and small farmlands (<500 m², and 150 m² on avergae), it was obvious that the heterogeneity is an important factor introducing bias in MODIS pixels over the GBP landscapes. As a result, LST–NDVI trapezoidal (or triangular) distribution did not form, and a low correlation was found between LST and NDVI to derive the coefficients for use in Equations (7) and (8). Using the selective filter (Figure 6b, Figure 7b), reduced the bias caused by heterogeneous landscapes and yielded a better accuracy that made the fusion model applicable to such regions. Adding the NDVI coefficient (Equation (11)) increased the weight of the NDVI in calculating the VTCI, and further helped achieve improvements in the results (Figure 7c).

Table 2. Coefficient of determination (R²) of acquired vegetation temperature condition indexes (VTCIs) for different dates before and after modifications.

–	Original Method	Using Selective Filter	Using Selective Filter and Weighting Coeff.
213–227	0.27	0.63	0.70
228–242	0.31	0.67	0.76
243–257	0.33	0.69	0.72
258–271	0.38	0.70	0.72
272–285	0.35	0.71	0.75
286–299	0.30	0.69	0.73
300–315	0.32	0.64	0.72

shows the improvement in R² that was acquired for VTCI after using the selective filter during the growing season. Adding this simple step helped decreased the error and overcame a potential limitation in applying the algorithm. The improvements (Table 2) are more notable for the early and during harvest period. That likely occurred due to the maximum mid-season greenness which reduced heterogeneity and caused a more regionally homogeneous landscape.

3.2. Pixel Level Evaluation Results

After applying the modifications, the daily ET maps were first generated based on the METRIC model with no integration of MODIS products, and then by using the fusion model with more frequent MODIS and Landsat data updates. As the default, we generated the ET maps by using four Landsat imagery and 8-day MODIS updates. The METRIC-estimated ET maps using only Landsat data at DOYs 189, 285, 301, and 325 were compared to its corresponding values from the fusion model-estimated ET. These comparisons are shown in Figure 8. Considering cultivated areas, results showed that the METRIC model had closer estimates to the observations compared to the results from the modified fusion model. We cannot directly comment on the resulted values in sparse vegetation (as our observations were carried out from the cultivated areas), but it was observed that the fusion model underestimated the values, and the underestimation was more significant in sparse vegetation areas (

Table 3. Evapotranspiration (ET) values resulted from using mapping evapotranspiration with internalized calibration (METRIC) model and modified fusion model.

DOY	Observed ET (mm/day)	METRIC ET (mm/day)		Fused ET (mm/day)
		Cultivated Area	Sparse Vegetation	Cultivated Area	Sparse Vegetation
189	3.5	3.2	1.1	2.75	0.4
285	7.2	6.8	3.5	6.2	2.4
301	4.3	4.1	1.8	4.0	1.6
325	2.9	2.5	1.0	2.3	0.3

).

Considering the domain, the results showed good agreement between the two sets of the resulting products, both visually and the statistical sense (RMSE of 0.15 mm day⁻¹ and average R² of 0.83). The high-resulting correlation indicated that the modifications were successfully applied to the fusion scheme and eliminated/reduced the potential limitation in applying the scheme over the study region.

Comparing the computed ET time series at the pixel level (i.e., averaged values over all the cultivated farmland during the growing season) resulted in a significant difference when the fusion model is used (Figure 9). Analysis of the METRIC model when it was used standalone for generating ET time series against in-situ measurements yielded a poor result with RMSE ranging from 2.88 to 3.44 mm day⁻¹, when using available Landsat imagery. Using four Landsat images was inadequate to capture the phenology change due to crop growth that occurred between the gap dates of the available imagery. The inadequacy of Landsat data alone as an input to the METRIC model was specifically observed between DOY 205 and 280 when the planted crop (mostly corn) shifted from early tilling stages to leaf and stem development stage that contributed to the underestimation of vegetative indices (i.e., leaf area index and NDVI) in the model calculations. Adding MODIS data with 8-day updates helped to capture the phenology changes that occurred during the dates with missing Landsat data. Using the fusion model decreased the RMSE to 0.37 mm day⁻¹ when the default approach (i.e., four Landsat images with 8-day updates of MODIS data) was used. However, using such a scheme would not be optimal for the end users in the study region because of the computational requirements, and also because the plant phenology change might not be significant over 8-day interval. Therefore, we tested different update frequency, and the results from which are presented in the following section.

3.3. Effects of Data Update Frequency on the Fusion Model

The resulting fusion based ET time series with different data frequency update are compared in Figure 10 for ET time series at the pixel level, and in

Table 4. Crop water consumption and ET estimates for different experiments ¹.

Landsat Images	No MODIS	48-day	40-day	32-day	24-day	16-day	Hybrid	8-day
1 DOY = 189	---	(a) ** 2.53 53.5 4.21	(b) ** 3.29 39.5 3.11	(c) ** 3.47 36.3 2.86	(d) ** 3.6 33.8 2.66	(e) ** 4.03 25.9 2.04	(f) ** 4.2 22.8 1.5	(g) ** 4.26 21.7 1.4
2 DOY = 189,325	---	(h) ** 3.06 43.7 3.44	(i) ** 3.51 35.4 2.79	(j) ** 3.95 27.3 2.15	(k) ** 4.28 21.4 1.68	(l) ** 4.63 14.9 1.17	(m) ** 4.73 13.0 1.1	(n) ** 4.77 12.3 0.97
3 DOY = 189,285 325	(o) ** 2.83 48.0 3.78	---	---	(p) ** 4.30 21.0 1.65	(q) * 4.73 13.1 1.03	(r) ns 4.95 7.1 0.69	(s) ns 5.00 6.5 0.61	(t) ns 5.1 6.2 0.55
4 DOY = 189,285 301,325	(u) ** 3.75 31.0 2.44	---	---	(v) * 4.49 17.5 1.38	(w) * 4.94 9.1 0.72	(x) ns 5.02 6.0 0.53	(y) ns 5.19 4.2 0.41	(z) 5.25 3.5 0.37

¹ Each experiment is denoted with a letter (The first letter in each cell). The second number in each cell is the cumulative seasonal water use estimated corresponding to the experiment (Mm³). The third and the forth numbers are the relative error (%) and the RMSE value (mm/day) corresponding to each experiment. ** = significant at P ≤ 0.05, * = significant at P ≤ 0.01, ns = not significant.

for estimating the seasonal crop water use in the study watershed. Statistical analysis highlighted that lower update frequencies yielded poor accuracy with relative error values as high as 53% when using only one Landsat image and updating the MODIS data for every 48-day (Figure 10a). Higher frequency updates yielded better accuracy; however, no significant improvement (based on ANOVA at 95% confidence interval) was seen between the hybrid, the 8-day, and 16-day frequency updates with the implementation of three to four Landsat data (Figure 10r–t,x,y). This explicitly highlighted the consideration of crop condition and phenology in the calculations. Using the fusion model with three Landsat images and 16-day frequency updates of MODIS data helped reduce the bias, optimize the processing time, and maintained accuracy. Additionally, the hybrid experiment yielded accurate results as the 8-day experiment with R² of 0.84, while using fewer computations. In the early growth stage, when crop growth is rapid and the crop coefficient changes significantly, updates can be done on a weekly basis. During maturity, the crop condition does not change appreciably, and updates every two to three weeks (depending on the crop type) could be sufficient to maintain the model accuracy. In addition, consideration should be given to the harvest time when a change in land cover is expected.

3.4. Comparison at the Watershed Level

The estimated cumulative crop water use from the fusion ET values for the season led to a similar conclusion as the previous section. The results are presented in Table 4. It should be noted that in the table, each experiment is denoted with a letter (The first letter in each cell), the second number in each cell denotes the cumulative seasonal water use estimated correspond to the experiment, the third and the forth numbers are relative error (%) and the RMSE value (mm/day) and correspond to each experiment. For the growing season, the total volume of net applied water was measured as 5.44 Mm³ which accounts for the ET_a [52]. Using a non-fusion approach with three Landsat images (experiment ‘o’) predicted the seasonal crop water use as 2.83 Mm³. The use of the fusion scheme reduced the relative error (RE) from 48% in experiment ‘o’ to 3.5% in the high-frequency fusion approach (experiment ‘z’). The fusion approach with 8-day updates provided estimated crop water use as 5.25 Mm³, and yielded the highest accuracy with 3.5% of RE. According to the ANOVA, the resulting accuracy was not significantly changed at 95% confidence interval when the 16-day and hybrid experiments were used, while lower frequency updates required even lesser computation, they yielded less accurate estimates (RE ranged from 22% to 53%). Thus, considering the accuracy and reduced computational burden, from the user perspective a 16-day or hybrid approach is recommended.

4. Conclusions

The potential of using existing MODIS–Landsat fusion algorithms for generating high spatio–temporal ET data in sparse, highly heterogeneous vegetation regions, motivated this study. The default fusion algorithm had a relatively high bias and did not account for different special resolutions of inputs because of heterogeneous landscapes caused by small land withholdings and farmlands of the study area. The issue of small landholdings and heterogeneity is more generic especially in arid and semi-arid farming regions. To reduce the bias, simple modifications were done to apply the fusion model by adding a selective pixel-wise filter and NDVI-based coefficient. The modified approach successfully reduced the bias over the heterogonous landscape and was used to assess weekly ET over the study domain (GBP region of southern Iran). The fusion model did not require complex, time-consuming processing steps, and yielded good results. Thus, the modification provided an adaptable framework to generate ET maps for local-scale applications. Another advantage of the fusion model was that it did not require processing MODIS data in the METRIC model, which can be challenging especially over heterogeneous landscapes for locating the “pure” end member pixels.

The question of frequency of imagery update and corresponding computation time was considered for ET estimations. Different experiments using available data over the growing season were analyzed to assess the impact on balancing accuracy and the practicality of the computational efforts. It was concluded that adding more information beyond bi-weekly updates of MODIS data made the algorithm more complex and did not result in significantly higher accuracy. The knowledge of plant phenology (crop growing cycle) helped the selection of data update frequency. The optimal fusion approach used weekly updates during initial crop growth stage and bi-weekly updates for the rest of the season. This led to ET estimates which agreed reasonably well with available observation (RE=4.2% and RMSE of 0.41 mm day⁻¹). Results highlight the potential and value of satellite fusion models in generating high-resolution ET estimates in data-sparse regions for water and agriculture management. However, attention should be paid to the geographical characteristic, vegetation covers and land cover compositions of the study area, and as a result; additional pre-processing steps may be required before applying the fusion scheme. This simple methodology should be tested with additional case studies over different regions to further demonstrate its utility.