Developing an Integrated Complementary Relationship for Estimating Evapotranspiration

The complementary relationship for estimating evapotranspiration (ET) is a simple approach requiring only commonly available meteorological data; however, most complementary relationship models decrease in predictive power with increasing aridity. In this study, a previously developed Granger and Gray (GG) model by using Budyko framework is further improved to estimate ET under a variety of climatic conditions. This updated GG model, GG-NDVI, includes Normalized Difference Vegetation Index (NDVI), precipitation, and potential evapotranspiration based on the Budyko framework. The Budyko framework is consistent with the complementary relationship and performs well under dry conditions. We validated the GG-NDVI model under operational conditions with the commonly used remote sensing-based Operational Simplified Surface Energy Balance (SSEBop) model at 60 Eddy Covariance AmeriFlux sites located in the USA. Results showed that the Root Mean Square Error (RMSE) for GG-NDVI ranged between 15 and 20 mm/month, which is lower than for SSEBop every year. Although the magnitude of agreement seems to vary from site to site and from season to season, the occurrences of RMSE less than 20 mm/month with the proposed model are more frequent than with SSEBop in both dry and wet sites. Another finding is that the assumption of symmetric complementary relationship is a deficiency in GG-NDVI that may introduce an inherent limitation under certain conditions. We proposed a nonlinear correction function that was incorporated into GG-NDVI to overcome this limitation. As a result, the proposed model produced much lower RMSE values, along with lower RMSE across more sites, as compared to SSEBop.


Introduction
According to the U.S. Geological Survey (USGS) Famine Early Warning Systems Network [1], the rate and amount of evapotranspiration (ET) plays a considerable role in the monitoring of water loss from agricultural lands. As noted by Senay et al. [2], ET may be used to show the current vegetation condition compared to the historical records. This comparison has the potential to help identify vegetation stress in time and space. ET estimation methods can be divided into two types: (1) ground-based ET methods that use standard meteorological data; and (2) ET models that use remote sensing data that must be combined with retrieval algorithms to estimate ET.
McMahon [3] classified the ground-based ET methods into six classes on the basis of application: 1) potential evapotranspiration (ETP); 2) reference evapotranspiration; 3) actual evapotranspiration; 4) open water evaporation; 5) lake/ storage evaporation; and 6) pan evaporation. We have focused on actual ET in this study because it can be representative of actual conditions, whereas reference evapotranspiration would require a vegetation resistance parameter and deep lakes would require water temperature data. In addition, we use the term "evapotranspiration (ET)" in this paper to include actual evapotranspiration except in places where the term "reference (crop) evapotranspiration" is used by other authors.
One approach to estimating ET with ground-based methods is the complementary relationship proposed by Bouchet [4]. The primary advantage of the complementary relationship is that it generally requires only meteorological data. Bouchet [4] suggested that as a surface dries, the decrease in ET is matched with an increase in potential evapotranspiration (ETP) as shown in Figure 1.
Such a relationship offers a simple and attractive approach for estimating ET using ETP without the detailed knowledge of surface properties. Examples of widely known models using this concept are the Advection-Aridity (AA) model by Brutsaert [5], the Complementary Relationship Areal Evapotranspiration (CRAE) by Morton [6], and the GG model proposed by Granger [7]. These three models have been widely applied to a broad range of surface and atmospheric conditions [8] [9] [10] [11] [12].
Granger [13], however, argued that the symmetric relationship in Bouchet [4] [14] showed that the radiometric surface temperature measurements can be successfully incorporated into Granger [13] equation. Similar to Crago [14], Anayah [15] proposed a modified version of the GG model using Priestley [16] equation instead of Penman [17] equation. The model proposed by Anayah [15] is hereafter called the modified GG model. The results of the modified GG model showed a decrease in Root Mean Square Error (RMSE) from 20% to as much as 80% compared to the recent studies of Mu [18], Mu [19], Han [20], and Thompson et al. [21]. On the other hand, Kahler [10] proposed an empirical constant, b, in Bouchet [4] hypothesis and demonstrated that b is generally greater than 1, based on their theoretical and experimental evidence, while the symmetric condition of Bouchet [4] hypothesis requires b = 1. More recently, Aminzadeh [22] extended the asymmetric complementary relationship with an analytical prediction of b for Kahler [10].  [30]. Figure 2 presents the results obtained from Kim [26]. These results are in agreement with Anayah [15], which showed that the modified GG model needs further improvements in dry conditions, and showed the lowest mean RMSE in both dry and wet sites. Overall, these results indicate that, among the groundbased methods [26], model can be used as a powerful methodology to estimate ET. While these findings are good within the realm of complimentary methods (or ground-based methods), some of the more commonly used ET estimation methods now use remote sensing data. If the complementary relationship and the corresponding methods, such as the model proposed by Kim [26], are to be accepted as operational models in field conditions, then the results should be compared and validated with remote sensing-based ET estimation methods. Taking into consideration of the improvements made with complementary relationship-based methods, this study examines the work of Kim [26] in comparison with a commonly used remote sensing method and measured ET data from 60 EC flux tower sites located across the USA.
Biggs et al. [31] grouped the remote sensing-based methods into three classes: vegetation-based methods, radiometric land surface temperature-based methods, and triangle/trapezoid or scatterplot inversion methods. Among them, the radiometric land surface temperature-based methods have a number of attractive features compared to the other classes: minimal ground data, ease of implementation, and operational application over large areas.
Radiometric land surface temperature-based methods use the fact that ET is a change of state in water that uses energy in the environment for vaporization and reduces surface temperature [32]. A subset of these methods is often called energy balance methods since they solve the energy balance equation. Moreover, these methods do not directly measure ET but must be combined with retrieval algorithms since data and technical requirements to solve the full energy balance equation can be challenging, especially in large regions. For example, the Surface Energy Balance Algorithm for Land (SEBAL) model [33] [34] requires the measurements of wind speed, iterative calibration, and review by an expert operator. Mapping EvapoTranspiration at high Resolution with Internalized Calibration Figure 2. Comparison of RMSE (mm/month) between different complementary relationship models for 29 dry and 30 wet sites in the US. NGG and GG-NDVI refer to the models of Han [20] and Kim [26], respectively. (METRIC) [35] needs high-quality meteorological data such as net radiation, air temperature, wind speed, and humidity. According to Allen [35], METRIC has higher accuracy for hourly reference ET than SEBAL, but the processing cost of METRIC is high.
As an alternative, FWESNET (USGS) has produced ET measurements from MODIS using the operational Simplified Surface Energy Balance (SSEBop) model [2]. The SSEBop setup uses the Simplified Surface Energy Balance (SSEB) approach developed by Senay [36]. The SSEB approach estimates ET using ET fraction scaled from thermal imagery in combination with a spatially explicit maximum reference ET. SSEB has an advantage in that it does not require air temperature and the knowledge of land cover types. Instead, the method uses the "hot" and "cold" pixel approach of Bastiaanssen [33] to calculate the ET fraction.
Later, Senay [37] enhanced SSEB to accommodate diverse vegetation and topographic conditions using a lapse rate correction factor. They successfully evaluated the results by comparing with METRIC and ET values computed from the water balance approach. As a result of the work by Senay [37], the enhanced SSEB model increased the correlation with METRIC from 0.83 to 0.90. Furthermore, Senay et al. [38] proposed a revised SSEB to handle both elevation and latitude effects on surface temperature using the difference between Land Surface Temperature (LST) and air temperature. Recently, Senay et al. [2] proposed an operational SSEB, renamed as SSEBop, which uses predefined boundary conditions for hot and cold reference pixels so that ET can be calculated as a function of LST and reference ET. The SSEBop approach has been validated comprehensively by comparing with 45 EC flux tower observations [2] and then with both MOD16 and 60 EC flux tower observations [39]. Later, Bastiaanssenet al. [40] applied SSEBop to determine ET in the Nile Basin, Ethiopia, for mapping water production and consumption zones. SSEBop ET data is now freely available through the USGS Geo Data Portal.
Despite the general consensus of using SSEBop for estimating ET, a detailed study of SSEBop conducted by Senay et al. [2] showed that the use of reference ET can introduce a significant difference of up to 20% in the magnitude of ET.
They also showed that the use of constant pre-defined differential temperature between the hot and cold boundary conditions can also create an inherent inaccuracy. Thus, it is important that SSEBop ET be validated and calibrated with available data such as EC flux tower data before using it to model ET.
The facts provided in the previous discussion indicate a need to further validate both Kim [26] and SSEBop models in the operational application of the complementary relationship in estimating ET. Therefore, the objectives of this

Methodology
GG-NDVI is the most updated model using the original GG model. GG-NDVI uses historical annual Normalized Difference Vegetation Index (NDVI) data and precipitation to improve the ET estimates of the modified GG model proposed by Anayah [15]. We then used the SSEBop model (Senay et al. [2]) to further validate GG-NDVI in comparison to an operational remote sensing model.

GG-NDVI Model
The first complementary relationship was proposed by Bouchet [4], who postulated that, as a surface dries, the actual ET decrease is matched by an equivalent increase in ETP. In spite of the fact that ET is negatively correlated with ETP, Morton [6] showed that the relationship has no defined shape. Granger [13] showed that the symmetrical relationship between ET and ETP only occurs when the temperature is near 6˚C and suggested the following complementary relationship formulation: where ET, ETP, and ETW are in mm/day, γ is the psychrometric constant (kPa/˚C), and ∆ is the slope of saturation vapor pressure-temperature (kPa/˚C) relationship. Thereafter, [7] developed the GG model based on Equation (1) using the concept of relative evaporation. Recently, Anayah [15] developed the modified GG model using the work of Granger [7]. The performance of the modified GG model improves when the Priestley [16] equation shown in Equation (2) is used to calculate ETW instead of Penman [17].
where α is a coefficient equal to 1.28, n soil R G − is net radiation (mm/day), and soil G is soil heat flux density (mm/day). Note that soil heat flux density is negligible compared to net radiation when calculated at daily or monthly time-scale [9].
ET is then estimated as a fraction of ETW using Equation (3): where G is the relative evaporation parameter derived from [7]. They proposed a unique relationship with a parameter called relative drying power (D). The unique relationship between G and D are described in Equations (4) and (5), respectively. 8.045 where a E is drying power of air (mm/day) given in Equation (6).
where U is wind speed at 2 m above ground level (m/s), which is adjusted using the work of Allen [41]; s e is saturation vapor pressure (mmHg); and a e is vapor pressure of air (mmHg). The performance of the GG model, including the modified GG model proposed later, decreased with increasing aridity. A possible reason is G in Equation (4), which was empirically derived from 158 sites representing wet environments in Canada. To improve the parameter G, GG-NDVI model (Kim [26]) used the latest version of the Fu equation (Li [27]). In particular, the Fu [42] equation is one of the formulations of the Budyko curve [43] and it is consistent with the complementary relationship [28] [29]. The corresponding analytical formulation of the Fu equation is given in Equation (7).
where P is precipitation (mm) and ETP is estimated using Penman [17]. Parameter ϖ is a constant and represents the land surface conditions of the basin, especially the vegetation cover [27]. Furthermore, Li [27] showed that ϖ is linearly correlated with the long-term average annual vegetation cover that can help improve ET estimates. Reference Yang et al. [44] showed that vegetation cover defined by M is calculated using Equation (8).
where NDVI min and NDVI max are chosen to be 0.05 and 0.8, respectively. An optimal value for the basin can be derived through a curve fitting procedure that minimizes RMSE between the measured and predicted evaporation ratio [27].
Li [27] proposed parameterization that is simply a linear regression between optimal ϖ and the long-term average M given as where a and b are constants that are found for each site.
To incorporate Equation (7) into the modified GG model, Kim [26] used the work of Yang [28] and Zhang et al. [29]. According to Zhang et al. [29], the Fu equation showed that the rate of change of ET with precipitation increases with ETP but decreases with precipitation. This is similar to the complementary relationship proposed by Bouchet [4]. Later, Yang [28] derived the consistency between the Fu equation and the complementary relationship using 108 dry regions in China. With this theoretical background, Kim [26] Note new G is the updated definition of relative evaporation, G, which includes the Budyko hypothesis and the vegetation index. To estimate new G , ETP is required and can be estimated using Equation (11) [17].
Having found new G from Equation (11) and estimated ETW from Equation (2), we can estimate ET of the proposed model from Equation (12).

SSEBop Model
The SSEBop algorithm (Senay et al. [2]) does not solve the full energy balance equation. This approach assumes that for a given time and location, the temperature difference between the hot and cold reference values of each pixel remains nearly constant throughout the year under clear sky conditions. Furthermore, the major simplification of SSEBop is based on the knowledge that the surface energy balance process is mostly driven by net radiation. With this simplification, the ET fraction, ETf, is calculated using Equation (13).

ETf
Th Ts Th Ts dT Th Tc Here, ETf is between 0 and 1, with negative ETf values set to zero; Ts is surface temperature derived from MODIS LST; Th is hot reference value representing the temperature of hot conditions; and Tc is the cold reference value derived as a fraction of maximum air temperature [2]. The difference between Th and Tc is dT with temperature units in Kelvin.
ET is estimated using Equation (14)

Data
First, we used the SSEBop ET data set from the USGS Geo Data Portal  We defined the climate class of each site using the aridity index of the United Nations Environment Programme (UNEP) proposed by Barrow [47]. The aridity index divided climate conditions to six classes: hyper-arid, arid, semi-arid, dry sub-humid, wet sub-humid, and humid. However, this work simplified the climate class definition to two classes, similar to the work of Anayah [15]: dry and wet. Using this simplification, 24 sites were identified as dry, compared to 36 sites under the wet class.

Results and Discussion
This study was conducted in two phases. Phase 1 is the validation stage in which comparisons are made between the SSEBop model and measured ET to assess the accuracy of the remote sensing method to estimate ET. In Phase 2, a comparison of estimated ET from GG-NDVI with observed data will be performed to identify the weaknesses of the GG-NDVI model, especially relative to the complementary relationship, and appropriate corrections will be proposed.

Phase 1: Validation of GG-NDVI
Capturing inter-annual variations of ET estimates is important. Although such variations are not significant when water is unlimited, estimating these variations in water-limited conditions is essential for water resources management. In this phase, ET has been estimated from both SSEBop and GG-NDVI and compared against measured monthly ET data from 2000 to 2007.    to the predefined cold boundary (Tc), which brings ETf closer to 1.0, resulting in a corresponding ET that is close to the maximum ET.
According to Table 1 and Figure 6, the mean RMSE of GG-NDVI ranged between 15 and 20 mm/month, and GG-NDVI showed lower RMSE than SSE- Bop every year. Although the magnitude of agreement (overestimation or underestimation) seems to vary from site to site and from season to season, Figure 6 confirms that the occurrence of an RMSE less than 20 mm/month with GG-NDVI  the previous studies of Hobbins [9], Xu [12], and Anayah [15].
Based on these results, we could conclude that GG-NDVI is a reliable approach for estimating ET, the novelty of GG-NDVI being that the Fu equation can be used to define relative evaporation in the original GG model using NDVI.
This approach showed a reasonable match between GG-NDVI and the 60 Ame-riFlux sites. However, GG-NDVI may not predict ET accurately when the vegetated cover changes significantly or is dense. For example, at Brooking in South Dakota, the mean RMSE of GG-NDVI was 42 mm/month, compared to 18 mm/month with all sites, and NDVI has a large seasonal vegetation cover as shown in Figure 7. A possible reason is that the relationship between NDVI and vegetation can be biased in sparsely vegetated areas with a Leaf Area Index (LAI) of less than 3. According to [50], the Soil Adjusted Vegetation Index (SAVI) is recommended instead of NDVI when LAI is less than 3. It should be noted that the LAI of Brookings is about 2.5. Furthermore, prior studies of Mu [19] and Yuan et al. [51] have demonstrated that NDVI is insufficient to represent vegetation under dense vegetation conditions. Recently, Zhang [52]   to Yang et al. [44], the relative infiltration capacity and the average topographic slope need to be taken into consideration when using the Fu equation, especially in small catchments. Therefore, more work is needed to generalize the relationship for the use of NDVI with changing vegetation cover within the Budyko framework. The next section will discuss options to improve the GG-NDVI model.

Phase 2: Enhancement to GG-NDVI
As described earlier, GG-NDVI performed slightly better than SSEBop in both dry and wet climate conditions, and GG-NDVI increased the predictive power with increasing humidity. One interesting finding is that RMSE from GG-NDVI increases slightly with the relative evaporation parameter as shown in Figure 8.  ( ) where f(G) is the correction function. We expect the correction function to be nonlinear, similar to an exponential function, since the magnitude of the difference between ET and ETW decreases exponentially as shown in Figure 1. In this work, we fitted 2772 data points to an exponential function similar to Equation (18). Multiple regression analysis was conducted to compute the values of the α and β coefficients.
Regression analysis found that α is 0.7895 and β is 0.9655. Hereafter, the GG-NDVI model with the proposed correction function given as Equation (17) is called the Adjusted GG-NDVI model.
To determine the accuracy of Adjusted GG-NDVI, comparisons were made between the results from the Adjusted GG-NDVI and GG-NDVI and between measured ET data and ET values from SSEBop. These comparisons are shown in Figure 9 and Table 2 across 60 sites. While ET from GG-NDVI at Blodgett in California ( Figure 9) showed deviations from measured ET, we can see that the Adjusted GG-NDVI produced ET estimates close to measured ET and reduced mean RMSE from 33 to 22 mm/month for Mize and 17 to 10 mm/month for Blodgett. In Table 2, overall RMSE across 60 sites for GG-NDVI and Adjusted GG-NDVI were found to be 18 mm/month and 15 mm/month, respectively. Figure 10, which presents a histogram of RMSE from the different ET models, shows a significant improvement attributed to the Adjusted GG-NDVI model. With Adjusted GG-NDVI, 38 sites have less than 15 mm/month of RMSE, compared to 26 sites with GG-NDVI. These results suggest that the use of the correction function in GG-NDVI can significantly improve accuracy in estimating ET. In addition, Equation (17) can be updated with the new definition of G as ( ) where the value of 2f(G) can vary between 1.64 and 3.04 as G varies based on site-specific conditions. The new formulation of the Adjusted GG-NDVI model described in Equation (19) clearly shows that the relationship between ET and ETP is not symmetric with respect to ETW, further confirming the earlier conclusions that the hypothesis of Bouchet [4] needs to be extended and applied    Adjusted GG-NDVI 7 15 34 with appropriate corrections.

Summary and Conclusions
ET estimation models using the complementary relationship are able to estimate ET in most instances. In particular, the model proposed by Anayah [15]  relationship models but also showed the need for further refinements, especially under dense vegetation conditions. On the other hand, remote sensing methods are more common as operational models under field conditions. In order to determine whether complementary methods such as GG-NDVI can compete and deliver accuracy similar to remote sensing methods, it is important to make appropriate comparisons. The objectives of this work were therefore twofold: (1) evaluate the recently developed ET estimation method, GG-NDVI, to see if it could deliver similar accuracy to the commonly used operational remote sensing method, SSEBop and (2) identify the inherent weaknesses of the original complementary relationship and make appropriate refinements to further improve the GG-NDVI model, especially under dense vegetation conditions. For this purpose, we selected 60 AmeriFlux sites located across the US.
The first phase of the analysis showed that the GG-NDVI model with the Budyko framework and relative evaporation was found to work reasonably well. Validation with 60 AmeriFlux sites indicated similar levels of accuracy for both SSEBop and GG-NDVI. R-square between GG-NDVI and measured ET ranged from 0.40 to 0.79, overall RMSE of GG-NDVI ranged between 15 and 20 mm/month, and GG-NDVI showed lower RMSE than SSEBop every year. Furthermore, the occurrences of RMSE less than 20 mm/month with GG-NDVI were more frequent than SSEBop. Based on these results, we concluded that GG-NDVI is a reliable approach for estimating ET.
The second phase of the analysis showed that the predictive power of GG-NDVI decreased with relative evaporation possibly due to the use of the symmetric complementary relationship in estimating ET. In order to identify the true relationship between ET and ETP with respect to ETW, an exponential correction function was proposed. This phase demonstrated that the inclusion of relative evaporation with a correction function greatly improved the performance of the Adjusted GG-NDVI. For example, 68% of Adjusted GG-NDVI sites had RMSE less than 15 mm/month compared 43% with GG-NDVI.
In essence, this study strengthens the idea that the use of vegetation cover information in the complementary relationship has increased ET estimation power. More importantly, this work showed that the symmetric relationship typically assumed with the complementary relationship may not be valid. Instead, the results show that the symmetrical relationship needs to be updated with a nonlinear correction function as proposed here. A key strength of this study is that the latest proposed version of the GG model, Adjusted GG-NDVI, overcomes limitations of both relative evaporation as proposed by Granger [7] and the assumption of a symmetric complementary relationship from the work of Bouchet [4]. Consequently, Adjusted GG-NDVI can lead to significantly increased accuracy of ET estimates under diverse climate conditions while producing comparable or even better results than the SSEBop operational remote sensing model.