Assessing Accuracy of Vegetation Index Method to Estimate Actual Evapotranspiration

The estimation of actual crop evapotranspiration (ETa) maps using complex equations and remotely sensed shortwave and thermal infrared imagery can be challenging and may require input data that are not available. There is an opportunity, therefore to create a simpler and faster method to generate ETa maps using fewer input parameters for situations where limited input data is available or greater uncertainty in the resulting ET estimates are acceptable. We compared the estimates of ETa produced by a crop coefficient and NDVI-based (Kc-NDVI) method to a full energy balance (EB) method. Clear sky images from Landsat 7 and Landsat 8 were processed and used for the ETa estimations from maize during two growing seasons in eastern South Dakota, USA. The results showed that the ETa values from the Kc-NDVI method were lower than the ETa values from the EB method by 18% for 2015 and 11% for 2016 growing season. During study period the accuracy of ETa estimation decreased 17% with the Kc-NDVI method. For both years the mean bias error was 0.81 mm day -1 and the root mean square error (RMSE) was 0.37 mm day -1 . The average daily ETa of 5.3 mm day -1 . The Kc-NDVI method performed acceptable for ETa estimations, indicating that this method can be used to estimate ETa with minimum input parameters at focused regional and field scales for short time periods.


Introduction
The accurate estimation of crop evapotranspiration (ET) plays an essential role in irrigation water management such as in system planning and design, and irrigation scheduling [1]. ET varies relative to weather conditions including air temperature, solar radiation, wind speed, and air vapor pressure deficit and plant and soil conditions [2][3][4].
In irrigated agriculture a widely recommended method for estimating crop water needs or actual evapotranspiration (ET a ) is multiplying reference evapotranspiration (ET r ) with a crop coefficient (K c ) [3,5,6] (Eq. 1).
(1) ET r is estimated based on meteorological information from a local weather station using the Penman-Monteith equation [3,6]. Generalized values of K c can be taken from literature values [3,7] when appropriate. As an alternative to using K c values from the literature, there are several methods for measuring ET a directly to estimate K c values over homogeneous surfaces. Methods include weighing lysimeters, Bowen Ratio Energy Balance System (BREBS), Eddy Covariance (EC), scintillometers [8,9] or soil water balance methods. However, these methods provide point or near point measurements that may not fully represent the ET a from a larger population of fields other than where the measurement was conducted [10,11]. To overcome this problem of estimating ET a from a large number of fields, remote sensing-based ET estimation methods are increasingly used for estimating crop water use and K c values. A primary advantage of remote sensing-based methods is they provide ET a estimates at the same resolution of the satellite imagery, which enables the estimation of ET at a field-by-field or subfield basis at a regional scale [12][13][14].
Several models based on remote sensing techniques have been developed to estimate ET a at different scales [8]. Mapping EvapoTranspiration at High Resolution using Internalized Calibration (METRIC) Model [13,15,16] is one such model. METRIC utilizes shortwave and thermal wavebands along with ground-based weather information to solve the surface energy balance to estimate the K c through a series of steps, which includes estimates of the dominant atmospheric heat transport mechanisms. In the last decade the METRIC model has been used to estimate ET a at field and regional scales in different crops and vegetation types including cotton [17,18], wheat [10,19], banana orchard [20], soybean [21], maize [22], cover crops [23], alfalfa [24], pistacho [25], vineyard [26,27], olive orchard [28], sugarcane [29], and forest in the Amazon [30].
Another, simpler method using satellite imagery is using the Normalized Difference Vegetation Index (NDVI) to estimate K c [31,32] for ET a estimation using Eq. 1. NDVI indicates the density and robustness of surface vegetation [33] and reflects the actual crop conditions [32,34]. For well watered crops there is typically a linear, crop-specific correlation between NDVI and K c . For more than 30 years local regression functions for the NDVI and K c relationship have been established for agricultural crops (e.g. [16,[34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49]). We used two satellite-based approaches to estimate ET a for irrigation applications namely 1) the energy balance method using METRIC and 2) the K c vs NDVI method [9,[50][51][52]. The energy balance method (EB method) is complex, computational involved and data intensive and require trained personnel to complete. In contrast, the K c vs NDVI method, which will be referred to as K c-NDVI method henceforth, is simpler, less data intensive and can be completed within a shorter timeframe, and at the same spatial resolution as the energy balance [9,33,51,53]. The performance and a comparison between these methods for ET a estimation have not been clearly determined in eastern South Dakota. The objective of this study was to compare the accuracy of the K c-NDVI method to calculate ET a relative to the EB method calculated by the METRIC model over two growing seasons in eastern South Dakota. The study was carried out in eastern South Dakota during the 2015 and 2016 growing seasons (Figure 1 (a)). The study area had an average latitude of 44° 19' N and longitude of 96° 46' W and elevation of 500 m above sea level (Figure 1 (b)). Five maize fields located within 15 km of each other were studied each year. Different fields were used the two years due to the crop rotation (Figure 1 (c)). All fields were in a maize -soybean crop rotation system common to the region. Weather information was acquired from the Brookings Mesonet weather station operated by the South Dakota Climate Office in each growing season. The soils were silty clay loam with 0-2% slope (NRCS Web Soil Survey 2016). The actual maize plant population density was approximately 78,000 plants ha -1 and the fields were managed using common agricultural practices used in the region. The crop was not considered subjects to growthlimiting stress from pests, weed or nutrient deficiencies.

Study Area
The maize fields were around 64 hectares in size. Irrigation is uncommon in this area and none fields were irrigated. The normal average annual precipitation is 533 mm, of which ¾ typically falls during the growing season (April-October).

Landsat Images
Clear sky images were used for the ET a estimations ( Table  1). The images were downloaded from the United States Geological Survey (USGS) EROS Datacenter and processed using the METRIC model running in the ERDAS Imagine software environment [54]. The wedge-shaped gaps appearing within the Landsat 7 images as result of the SLCoff issue were removed using the Imagine built-in focal analysis tool [55].

Pixel Selection
Ten pixels in each field were randomly selected and their values for NDVI, K c and ET a were extracted. The same pixels were used throughout each growing season. The number of pixels (10) were assumed to be representative of each entire maize field.

METRIC Model and Input Parameters
METRIC model version 3.0 was used to estimate ET a . Please see [13,15,56] for a detailed discussion of the model calculations.
In the METRIC model four primary input parameters are used to estimate ET a namely the Landsat image (including shortwave and thermal bands), digital elevation map, land cover map, and weather data ( Figure 2). The elevation and land cover map were reprojected in meters to the same pixel size as the Landsat images (30 m x 30 m).

NDVI Calculations
The NDVI values range from -1.0 to +1.0, with water having negative values and dense vegetation having high positive values [57,58]. For Landsat 7 NDVI was calculated as: For Landsat 8 NDVI was calculated as: where and !" are the corrected spectral radiance in the near-infrared and red bands, respectively.

Crop Coefficient (K c ) Curves for NDVI Based Method
The alfalfa-based K c values from [7] for 2015 and 2016 crop growing seasons were used. For K c estimations this method divides the growing season into two periods, viz. percent of time from planting to effective cover and days after effective cover to harvest. The effective cover of maize for our study occurred in middle of July for 2015 and early July for 2016 based on field observations of the crop phenology.

Relationship Between NDVI and K c and Generation of ET a maps
A relationship between NDVI derived from NDVI maps and K c values from [7] at each overpass date was established. This relationship was used to develop a linear regression equation for both seasons. Those linear regression equations were used to generate K c maps using Model Maker tool of ERDAS Imagine. The K c values derived from the K c maps were multiplied by ET r to create ET a maps for both seasons using the K c-NDVI method. The ET r values were estimated based on weather data from the automatic Brookings weather station. In the final step, the ET a values from ET a maps were compared with ET a values obtained from the EB method for each overpass date and for each growing season.

Average Ratio of ET a K c-NDVI to ET a EB and Their Relationship
The average ratio of ET a K c-NDVI to ET a EB was calculated to quantify the accuracy and performance of the K c-NDVI method for ET a estimations. Figure 3 shows an example of ET a maps developed using the EB method and developed by K c-NDVI method on July 20, 2016. The ET a K c-NDVI method map generally shows higher ET a values compared to the ET a EB method. This is due to the calibration of the maps and to a lesser degree differences in resolution between the maps. Also there is a difference in how the colors are displayed between these two maps. The pixel resolution in the ET a K c-NDVI method is 30 by 30 m, while in ET a EB method the thermal pixel resolution for Landsat 7 is 60 by 60 m and for Landsat 8 is 100 by 100 m.

Figure 3. ETa maps generated using the EB method (left) and using the Kc-NDVI method (right) on July 20, 2016.
A similar comparison of ET a maps over agricultural areas generated by the METRIC model using energy balance and using vegetation index data were reported by Allen et al. [13] and Anderson et al. [59] in Twin Falls, Idaho. Mokhtari et al. [25], found that the METRIC-based ET is highly sensitive to surface temperature, but less sensitive to NDVI.
For the 2015 season, Figure 4 shows that the discrepancy between the ET a values were higher at the beginning and at the end of the growing season. The highest ET a values were showed in the mid-season (July 18) 7.9 and 7.7 mm day -1 for the EB method and the K c-NDVI method, respectively.
For the 2016 season, Figure 4 shows low ET a values at the beginning of the growing season at 2.8 and 1.7 mm day -1 for EB method and for K c-NDVI method, respectively. Moderate ET a presented at the end of the season for EB method was 4.2 mm day -1 and for K c-NDVI method was 3.0 mm day -1 . High ET a values were observed in the mid-season (July 12) with 8.9 mm day -1 for EB method and 8.7 mm day -1 for K c-NDVI method.
In general, the ET a values estimated with EB method were higher than the ET a values estimated with K c-NDVI method by 18 and 11% for 2015 and 2016 growing seasons, respectively. Because the K c-NDVI method overwhelmingly considers transpiration from green vegetation, and only to a small extent evaporation from bare soil, some underestimation during the shoulder periods of the growing season is common. These results coincide with those in previous studies reported by Anderson et al. [59] who reported that ET a calculated from vegetation index data always were found to underestimate seasonal ET a values in irrigated areas in Idaho.

Average Ratio of ET a K c-NDVI Method to ET a EB Method
The average ratio distribution of ET a K c-NDVI to ET a EB method for 2015 and 2016 crop growing seasons are shown in Figure 5. This figure shows that all average ratios are below 1, which is denoted by the thick blue line. This means that the ET a K c-NDVI values were lower than the ET a EB values during the two growing seasons. In early and late season the K c-NDVI method showed the far values from 1, while in the mid-season the values were close to 1. Indicating that K c-NDVI is more accurate for ET a estimations during the mid-season than early and late seasons, this reflects low vegetation cover, high soil evaporation, and leaf senescence [13,16,44,59]. Therefore, the K c-NDVI method gives less accurate estimation of ET a during early and late season periods. For irrigation scheduling purposes during periods with high crop water demand at the middle of the growing season, the K c-NDVI method may be acceptable. However, ET a values from K c-NDVI method need to be adjusted during early and during late season to get close or accurate estimates to ET a EB values. The adjustment factor (ET a K c-NDVI / 0.66 = ET a EB) for the 2015 growing season was 0.66 and (ET a K c-NDVI / 0.71 = ET a EB) for the 2016 growing season it was 0.71.
For the entire 2015 growing season the underestimation was 21% and for the mid-season only (July-August) (excluding early and late season) was 12%, while for entire 2016 growing season the percent of error was 13% and for the mid-season it was 7%. The total average error for the two growing seasons was 17%. This general percent of underestimation with the K c-NDVI method is may be acceptable in some applications and are within the 10-30% error for an experienced expert reported by Allen et al. [9]. The average error for both growing seasons during the midseason stage was less than 10%.

Relationship Between ET a EB Method and ET a K c-NDVI Method
An acceptable relationship was found between ET a EB method and ET a K c-NDVI method during the 2015 and 2016 seasons with coefficient of determination (r 2 ) of 0.97 ( Figure  6). The corresponding mean bias error was 0.81 mm day -1 and the root mean square error (RMSE) was 0.37 mm day -1 . The average daily ET a was 5.3 mm day -1 .
In this study, the K c-NDVI method performed acceptably for ET a estimations during the two growing seasons. This indicates the K c-NDVI method can be used to estimate crop water requirements at regional and field scale in regions where digital elevation, land cover map and thermal infrared data are not available and where higher uncertainty is acceptable.

Conclusions
ET a values calculated with K c-NDVI method were lower than the ET a values calculated with EB method by 18 and 11% for the 2015 and 2016 growing season, respectively. The ET a K c-NDVI values were less than the ET a EB values during the two seasons especially at the beginning and at the end of the seasons when the vegetation cover was incomplete. Soil evaporation is not fully captured by the K c-NDVI method. As a result, the ET a estimated using the K c-NDVI method underestimated the values by 17% compared to the EB method during the period of study. The K c-NDVI method is less accurate during the early and late portion of the growing season, however for irrigation scheduling purposes, this method may be acceptable.
The results showed a good relationship between EB method and the K c-NDVI method for ET a estimation throughout two growing seasons. The K c-NDVI method can be a reliable method to calculate ET a using minimum input parameters.