Questioning of empirically derived and locally calibrated potential evapotranspiration equations for a lumped water balance model

In the hydrological cycle, the loss of transpiration over the vegetation and the evaporation through the water bodies on the surface is called evapotranspiration. The maximum amount of this loss or its occurrence at the time of sufficient soil moisture refers to the potential evapotranspiration (PET), while actual evapotranspiration is limited to the current soil moisture content. Since many physical geographical factors (such as latitude, altitude, and vegetation) affect the PET, the mechanism of its occurrence varies depending on the region and can be associated with meteorological observations such as relative humidity (RH) and wind speed, primarily radiation and temperature (Xu & Singh 2002). In this context, however, direct measurement of the PET is not feasible, as in pan evaporation. Thus, the precise prediction of the PET has great importance especially in the determination of plant water consumption and irrigation water requirement. Accordingly, observations by lysimeters (direct method) or empirical formulas (indirect methods) can be taken as a basis. The question of which empirical method is more appropriate was addressed by a limited number of researchers up to the early 1980s (e.g., Thornthwaite 1948; Blaney & Criddle 1950; Makkink 1957; Hamon 1961; Jensen & Haise 1963; McGuinness & Bordne 1972; Priestley & Taylor 1972). Later, the FAO (Food and Agriculture Organization) report (FAO-56), edited by Allen et al. (1998), stated that the Penman–Monteith (Pen-Mon) equation, which is a combination of the approaches of Penman (1948) and Monteith (1965), is a comprehensive indirect method. In some validation studies conducted in specific regions throughout the world with various climatic conditions, it has been emphasized that the Pen-Mon equation is highly compatible with lysimeter measurements (Itenfisu et al. 2003; Allen et al. 2005; Jain et al. 2008). It is therefore currently used as a reference formula (Pandey et al. 2016).

After the release of FAO-56, researchers tended to compare the Pen-Mon equation with other available methods. Moreover, re-calibration procedures were applied to update the coefficients of existing equations and the regression-based equations were also suggested (Xu & Singh 2002; Irmak et al. 2003; Xystrakis & Matzarakis 2011; Tabari et al. 2013; Bogawski & Bednorz 2014; Pandey et al. 2016). For example, in a study performed by Pandey et al. (2016), it was detected that Irmak (Irm), Makkink (Makk), Turc (Turc), and Blaney–Criddle (Bl-Cr) equations were more compatible with Pen-Mon for the north-eastern region of India. Xystrakis & Matzarakis (2011) applied 13 empirical equations on seven meteorological stations in southern parts of Greece and stated that the absolute biases obtained from the McGuinness–Bordne (McG-Bor) and Hamon (Ham) equations were relatively less. Xu & Singh (2002) applied five equations on the Changins station in Switzerland and suggested that if the Pen-Mon method was taken as the basis, the coefficients of Priestley–Taylor (Prs-Tyl) and Rohwer (Roh) equations should be recalibrated. Kellner (2001) also pointed out that the Prs-Tyl is a more appropriate equation for humid regions with high latitudes. Tabari et al. (2013) tested 31 equations at the Rasht station in the north of Iran, which has a humid climate. In their comparative study, it was argued that the Bl-Cr is more closely related to the Pen-Mon in terms of various statistical indices.

The main reason why the above-mentioned studies were performed is that the Pen-Mon equation of FAO-56 has more data requirements than other empirical ones. In particular, the solar radiation measurement network is not well-distributed as in temperature stations. Moreover, the regional studies will also shed light on which empirical equations can produce more consistent predictions. Even if these regional investigations are generally the subject of agricultural practices, a detailed analysis about PET is also needed for water balance modeling (Oudin et al. 2005a, 2005b). For instance, the lumped rainfall–runoff models also require PET input in addition to precipitation. In these models, the physical aspects of the transformation of rainfall into runoff are usually conceptualized by means of the parameters and the runoff series can be simulated by applying the continuity equation to the representative storage elements. Paturel et al. (1995) proved that runoff predictions derived from water balance models showed sensitivity to precipitation at a primary level and the PET also had a significant influence on model sensitivity. In the water balance modeling studies, the PET input can be defined as a regression-based function of mean temperature (T) and RH variables (e.g., PET = aT^b(100−RH); PET = aT^b; PET=ae^bT) or any empirical PET equation can be used directly (Xu & Vandewiele 1995; Fowler et al. 2007; Okkan & Kirdemir 2018). On the other hand, Xu & Vandewiele (1995) stated that these extra free parameters embedded into modeling would make the calibration procedure more difficult and that geographic regionalization studies with empirical PET input would be more proper. Vandewiele et al. (1991) have confirmed that when an empirical PET is used, the parameters may be better associated with the physical properties of the basin.

Various rainfall–runoff model practices using diverse empirical PET inputs are found in the literature. However, the reason why the related empirical equation was preferred in those studies remained ungrounded (McKillop et al. 1999; Bárdossy & Das 2008; Caldwell et al. 2015). In a few studies, the impacts of varied PET inputs on the rainfall–runoff performances were investigated in depth. The most comprehensive among them was the study prepared by Oudin et al. (2005b). In their work, the predictions obtained from 27 empirical PET equations were separately evaluated as inputs in the rainfall–runoff modeling and the outputs derived from application were compared with observed runoff series. The study was carried out across a large area including the major basins in Australia, France, and the USA. In particular, the performance of the McG-Bor equation in the hydrological model was reported to be more reasonable than those obtained from other equations. In another study, Kannan et al. (2007) asserted that the Hargreaves method (Harg) yielded more reasonable results than those of the Pen-Mon in the performance evaluation of the SWAT (Soil and Water Assessment Tool)-based modeling study for Bedfordshire, UK. In a study conducted by Wang et al. (2006) for the northwestern Minnesota basin, unlike other studies, it was found that the responses of the hydrological model produced by Pen-Mon, Prs-Tyl, and Harg equations were close to each other.

As can be seen from the literature reviews, there is no general judgment about which empirical PET equation has a better impression on hydrological model outputs. In Turkey, while the Pen-Mon and Bl-Cr methods are generally used in the calculation of plant water consumption (e.g., Beyazgül et al. 2000; Okkan & Kirdemir 2018), there are no equations examined or proposed from the point of view of rainfall–runoff modeling. In addition, significant changes in meteorological variables have been observed in many regions in Turkey due to climatic change. The increase in population and industrialization capacity, as well as greenhouse gas emissions, which have been continuously increasing, have triggered a meaningful increase in temperature and the PET. In this context, some hydro-meteorological predictions of basins have been made under different scenarios in Turkey. However, comparative analyses on PETs are unfortunately not available in these studies (e.g., Ozkul 2009; Okkan & Fistikoglu 2014; Okkan & Kirdemir 2018).

The study, prepared on the basis of the various reasons mentioned above, deals with the evaluation of the PET equations which may constitute an input to a monthly water balance model arranged on the basin scale and the determination of alternative parsimonious equations in terms of meteorological variables. In this study, Gediz Basin, which represents an important reserve of agricultural activities in the Aegean region of Turkey, has been selected as the study area. It is thought that the study performed has a unique value in the context of the derived results and can be adapted by researchers to other neighboring regions. The remainder of the presented study is arranged as follows. The next section consists of the details about the study region and data, followed by the background of the PET equations employed. Then, we explain the fundamental mechanism of the implemented water balance model and its performance criteria. The results and discussions obtained from the employed modeling strategy are presented next and the final section provides brief conclusions.

In the study, Gediz Basin, which is located in the western part of Turkey, was selected as the study area (Figure 1). The basin, which has a drainage area of nearly 17,000 km², is fed by the Gediz River as well as several streams. In the basin, dominated by the Mediterranean climate, the annual precipitation regime is around 550 mm for the reference climate period of 1981–2010. The annual mean temperature for the same period is around 15 °C throughout the basin. The prevalent economic activity in the basin is agriculture and the reservoirs in the basin are therefore generally operated for irrigational purposes.

There are 39 meteorological stations operated by the Turkish State Meteorological Service (MGM) and the General Directorate of State Hydraulic Works (DSI) in the basin. The locations of these stations are given in Figure 1. In all of these stations, precipitation observations are made, while only 20 of them have temperature data. Nine flow gauging stations (FGSs) representing the basin were also identified. The data of the FGSs for the 1981–2010 water year period were obtained from the DSI. The station named Muradiye Bridge has the largest drainage area on the main branch and is located in the western part of the basin, while the Borlu, Topuzdamlari, Derekoy, and Acisu stations feeding the Demirkopru reservoir are also major FGSs. In addition, Kayalioglu and Hacihidir in the northwestern part of the basin, Hacihaliller in the southwest, and Taytan Bridge in the southeastern part of the basin were included in the study. The locations of the FGSs are also given in Figure 1, and general information about them is summarized in Table 1. The weights of the precipitation stations representing the FGSs listed were obtained using Thiessen polygons, while the areal mean, maximum, and minimum temperature series were computed by arithmetic mean approach. Thus, the inputs to be used in the water balance modeling for the drainage area of each FGS were compiled.

In the study, while mean, maximum, and minimum temperatures (T, T_max, and T_min) were compiled from the meteorological stations which have a regular distribution over the basin, ERA-Interim data sets having 0.75° × 0.75° resolution were used for other variables needed for various PET equations, denoted in Table 2. The compliance of these re-analysis data sets with the observations has been verified by different researchers for various regions of the world (e.g., Dee et al. 2011; Bao & Zhang 2013). In Figure 1, eight ERA-Interim grids that almost uniformly comprise the study area and the center coordinates of the grids can be seen. These data sets are served by The European Center for Medium-Range Weather Forecasts for several categories. In addition, the grid numbers representing FGSs according to those of drainage areas are indicated in Table 1.

Before the application, the level of compatibility between ERA-Interim grid data and observations of meteorological stations was investigated. Since the temperature stations were regularly distributed in the basin, the average temperatures obtained from 20 stations and the arithmetic mean temperatures of the eight re-analysis grids were compared. Figure 2(a) demonstrates that the relationship is of good quality in terms of average temperatures. On the other hand, surface pressure (Press), mean wind speed (W), and solar radiation (R_s) were measured at a rare station over the basin. Thus, an exemplary comparison was performed between the data of Akhisar meteorological station and the Grid 7 data for RH and R_s variables. Figure 2(b) and 2(c) support the fact that an adequate relationship can also be obtained with the data of RH derived from monthly mean temperature (T) and monthly mean dew point temperature (T_dew), and R_s served with ERA-Interim.

As predominant equations used in the study are particularly radiation based, it is essential to compile the radiation data (solar or net) at this stage. In the wake of extracting R_s solar radiation data from the ERA-Interim database for the majority of radiation-based PET formulas, the difference between the R_ns net short-wave radiation and the R_nl net long-wave radiation, which gives the R_n net radiation on the ground surface, is also needed for other radiation-based equations such as Pen-Mon and Prs-Tyl. In this instance, it can be assumed that approximately 23% of R_s is reflected due to the albedo (Allen et al. 1998; Bogawski & Bednorz 2014), and the remaining part causes net shortwave radiation (R_ns ≈ 0.77R_s). Here, R_nl values can also be estimated based upon the Stefan–Boltzmann law (see Bogawski & Bednorz 2014).

Various methods attributed to the prediction of PET can be found in the hydrology literature (Xu & Singh 2002; Xystrakis & Matzarakis 2011; Bogawski & Bednorz 2014; Pandey et al. 2016). These methods are generally evaluated in sub-categories as mass transfer, temperature-based, and radiation-based (Oudin et al. 2005a, 2005b). In this study, the effect of the PET predictions, which were separately generated from 21 equations, on the water balance modeling performances were investigated over the different FGSs in Gediz Basin, Turkey. Although many empirical PET equations and their modifications have been found in the literature, the focus of this study is on the equations that have been frequently cited. These equations are in different categories and they need several meteorological input sets. The employed equations such as Thornthwaite (Thw), Romanenko (Rom), Blaney–Criddle (Bl-Cr), Kharrufa (Khr), Hamon-1 (Ham1), and Hamon-2 (Ham2) can be evaluated in the category of temperature-based methods. Based on the research performed by Singh & Xu (1997), mass transfer methods such as Meyer (Myr), Rohwer (Roh), and Penman (Pen), which can be used in the estimation of surface evaporation, were also included in the study. The other utilized equations, such as Turc (Turc), Hargreaves (Harg), Makkink (Mak), Penman–Monteith (Pen-Mon), Priestley–Taylor (Prs-Tyl), Caprio (Cpr), Jensen–Haise (J-H), Irmak1 (Irm1), Irmak2 (Irm2), Irmak3 (Irm3), McGuinness–Bordne (McG-Bor), and Baier–Robertson (Bai-Rob), which demand more intensive data sets including solar, net, or extraterrestrial radiation, can be considered as radiation-based methods. Among them, as Pen-Mon and Prs-Tyl methods incorporate energy balance and aerodynamic water vapor mass transfer principles, they can also be known as combination methods (Oudin et al. 2005b). The formulas of empirical PET equations used in the study are given in Table 2. They were compiled from the references mentioned in the last column of Table 2. The weather data notations and their units are also indicated in this table. Additionally, Table 3 summarizes the information on several supporter variables such as the temperature-based functions, wind speed scaling function, and some required constants.

Budyko (1958) argued that actual evapotranspiration (E_act) is a function of the precipitation (P) and the PET, and this relationship, which is derived for the annual time scale, is also referred to in the literature as the Budyko curve. Zhang et al. (2008) arranged this curve for monthly time scale data and integrated it into a water balance model termed dynamic water budget model (dynwbm), which is both conceptual and lumped. On the other hand, in the basins where there are dense agricultural activities, the groundwater table can be affected. In addition, the physical properties (land use or topographic patterns) may induce alteration in the groundwater level. Hence, an additional parameter, which is termed groundwater efficiency, was incorporated into the available dynwbm model by Okkan & Kirdemir (2018) in order to better conceptualize the groundwater storage process and to improve runoff prediction performances. Similarly, this five-parameter version was taken into consideration in the study.

In the model, the total amount of precipitation falling on the basin in any month is made up of two components, namely, direct runoff and catchment rainfall retention X, respectively. In this partition, the parameter α₁, which controls the first Budyko-type curve (see Figure 3, E1), plays an active role and a larger α₁ value results in more rainfall retention and less direct runoff. The model also has a sensitive parameter termed maximum soil moisture capacity (S_max), representing the soil and vegetation characteristics of the basin. On the other hand, the parameter α₂ controls evapotranspiration efficiency. In the case of an increase in this parameter value, an increment also occurs in the part of the water allocated to actual evapotranspiration E_act (Zhang et al. 2008; Li et al. 2016). The same parameter also controls the variable defined as evapotranspiration opportunity y, which is assumed to be composed of the sum of soil moisture content S and E_act (Sankarasubramanian & Vogel 2002). Thereby, y and E_act are both organized through the other Budyko-type curves (see E4 and E6 in Figure 3). During the month i, the available water content W_i may be expressed by the sum of the soil moisture remaining from the previous month (S_i-1) and X_i, as well as by the sum of the soil moisture content, actual evapotranspiration, and the amount of recharge (Rec) draining into groundwater storage. After Rec, S, and E_act are taken from the budget calculations, and baseflow prediction is made with the parameter d, the balance equation is then mounted for the groundwater storage G, which is postulated to represent linear reservoir behavior, using groundwater efficiency parameter e. The detailed description of the equations existing in the original model structure can be accessed from Zhang et al. (2008). The conceptual flowchart of the implemented version of the dynwbm model, its computation steps, and the related parameter definitions are given in Figure 3.

The long-term mean statistics collected from empirically derived PET predictions were examined prior to analyzing sensitivity of dynwbm to PETs. In this context, the dendrogram technique, which is a clustering analysis, was used to make an inference of long-term mean statistics of PETs throughout the basin. With this approach, the hierarchical relationship between the clusters represented in the mean statistics is summarized in Figure 4. The distance between the clusters was calculated with the Euclidean distance formula and Minitab package was used in the analysis. The lowest monthly average PET predictions were obtained from the Irm3, Bai-Rob, Ham1, Thw, and McG-Bor equations (67 mm/month, left-hand clustering in Figure 4), while Pen, Roh, and Rom were the highest predictions in the group (123 mm/month, right-hand clustering in Figure 4). In particular, it was considered that equations such as Pen and Roh, which were recommended for evaporations of water surfaces, were sensitive to the wind speed and therefore produced overestimation. While the long-term means range from 65 to 125 mm/month, the predictions of the Pen-Mon, the reference formula, are almost in the center of the two groups mentioned above (∼95 mm/month).

After examining the long-term mean statistics, the effects of PET estimations on runoff modeling were assessed. In practice, data covering the 1981–1995 water year period were used for calibration, while data covering the 1996–2010 water year period were used in the validation phase. The optimization of the parameters of dynwbm during the calibration phase was based on the LM algorithm to minimize the RMSE. In the calibration, the optimization of the model with different initial conditions was repeated 30 times. Thus, we have tried to reduce the risk of getting the local minimums.

The model parameters determined for all FGSs and PET formula variations were evaluated at validation stage and the statistical performance metrics such as RSR and NS were calculated. For the example of Muradiye Bridge, dynwbm parameters calibrated under different PET inputs and the related validation performances are given in Figure 5. According to the coefficients of variation (C_v) regarding the calibrated parameters, the upper limit value of the soil moisture storage has a variability of up to 13% in the basin. As the drainage area of FGS grows, the uncertainty of the different PET equations over this parameter increases (especially Muradiye Bridge and Acisu). On the other hand, parameters d and e of the groundwater storage system were also subjected to significant variability depending on the PET method and flow regime. The fact that d and e were interdependent parameters in the same groundwater storage function made it difficult to generalize the PET-based uncertainties. In addition, the α₁ parameter, which allows precipitation to be converted to direct runoff, has no significant sensitivity (C_v < 3%). The increased inter-method variability in α₂ also increased the variability of E_act/PET ratio (15% < C_v < 25%).

In addition to the parameter-based interpretation, the main theme of the study is to examine the contribution of PET equations with different characteristics to the runoff generation capability. In this study, the question of which PET equations produced a higher quality runoff prediction was performed with the validation outputs of the dynwbm. In this context, only the NS results were handled because of their high correlation with RSR values (not shared in the paper). In order to scrutinize the results over the entire basin rather than an FGS-specific assessment, the NS performances obtained from nine FGSs were examined by means of a box-plot (Figure 6(a)). Moreover, it was thought that ranking of NS values (with ties) could provide an idea in evaluation (Figure 6(b)). As the small alterations among the NS values can bring about biases in the ranking, the values are rounded to two decimals. In Figure 6(a), it can be seen that both the medians and the inter-quartile ranges (IQRs) have significant dissimilarity depending on the method. In other words, the influence of the predictions obtained from the different PET equations on the water budget elements and hence the model performances cannot be denied. The median statistics varied from 0.77 to 0.81, and based on the criteria denoted by Moriasi et al. (2007), these values are attributed to ‘very good’ modeling. When the distances between IQRs and whiskers are also viewed, the wide range of indices extracted from the methods such as Bl-Cr, Rom, Ham2, Pen, Roh, Irm1, and Irm2 make their usability ineligible in the whole basin. It is quite apparent that Irm3, McG-Bor, and Prs-Tyl methods, which exhibit scattering above the threshold NS = 0.75, make a greater contribution to dynwbm for the basin. From the rank representations presented in Figure 6(b), the performance of the McG-Bor was found to be very close to that of the Irm3. In several respects, the notion is that McG-Bor is the most reasonable radiation-based formula for the scope of the study in terms of parsimony.

As seen from Figure 7, the concordance of the predictions made by PET_adj2 and PET_adj3 with McG-Bor is more prominent compared to those of both Ham1 and PET_adj1. These equations were integrated into the hydrological model with the specified coefficients and dynwbm's parameter optimization process on all FGSs was repeated. The validation performances obtained by the integration of the proposed PET equations into the hydrological model are given in Table 4 in comparison with the McG-Bor outputs. According to the results, PET_adj3 showed almost identical validation performance with McG-Bor in the majority of stations (the small biases for Taytan Bridge, Kayalioglu, and Hacihaliller are negligible). Moreover, the usage of PET_adj2 compared to McG-Bor formula is also functional and, at the same time, reasonable in case the RH data are not provided or computed. In other words, Table 4 has proved that the proposed equations requiring only temperature input are capable of competing with radiation-based equations.

Another idea is to write an intensive PET function containing four free parameters into the hydrological model considering the general structure defined in Equation (6). In this variation, which has been tried out but not comprehensively presented to the readers within the scope of the study, the current five-parameter dynwbm has become a model having nine free parameters to be calibrated. The common diagnosis monitored in those practices is the increase in the degree of freedom and a boosted performance in the calibration period, while a significant decrease in the validation performance is observed. An exemplary application supporting this finding is given in Figure 8 for Muradiye Bridge, which has the largest drainage area in the basin.

The performance of the conceptual water balance models depends on the input of PET in addition to precipitation. Thus, which PET equation should be used in the flow simulation stage is an essential issue, especially in arid basins. However, this problem has been discussed in relatively few studies (e.g., Oudin et al. 2005a, 2005b; Kannan et al. 2007). To contribute to the related literature, this study aimed to test various PET equations which can be presented as input to the dynwbm for Gediz Basin in Turkey, and to determine alternative equations that use less meteorological data. The accordance of dynwbm with the PET equations was questioned by means of the several indices in the validation period. In contrast to the idea defended in the study by Wang et al. (2006), it appears at first glance that the effects of different PET equations on the water balance model cannot be neglected. According to performance criteria, the radiation-based McG-Bor equation, which causes the hydrological model over the basin to behave in the ‘very good’ class, appears to be superior. According to the NS ranking assessment applied throughout the basin, it was found that the McG-Bor equation resulted in similar responses to the hydrological model as the temperature-based Ham1 equation. Since the Ham1 is an economical method in terms of data requirement and because of the difficulty of obtaining solar radiation in the basin, it is accepted as an alternative method. Another viewpoint is to recalibrate some of the constants defined in the Ham1 equation to simulate the McG-Bor predictions. In accordance with this alternative, PET_adj2 and PET_adj3 were proposed as a result of the various arrangements performed. When the PET predictions derived from these parsimonious equations were presented as inputs to the water balance model, and the validation outputs were evaluated following the calibration stage, it was detected that PET_adj3 in particular provided a rather similar performance to McG-Bor throughout the whole basin. Thus, it has been determined that the proposed equations have more practical usage compared to radiation-based equations.

Due to the limited studies conducted in the literature connected with the topic, it is thought that the presented study is an inspirational one in terms of the content and suggested methodology. Since all the deductions have been made for a semi-arid region in Turkey, the derived equations embedded into water balance modeling can also be adapted to problems such as missing flow completion, discharge simulation under climate change scenarios, and defining supply–demand relationships of a reservoir for the neighboring basins and regions dominated by Mediterranean climate characteristics. For future development, the aim is to widen the scope of the study by taking into account the climate change scenarios (representative concentration pathways) attributed to the Fifth Report of the Intergovernmental Panel on Climate Change. In addition, there will also be a focus on both model-based (e.g., downscaling model, PET model, rainfall–runoff model) and scenario-based uncertainties.