Improved Hargreaves Model Based on Multiple Intelligent Optimization Algorithms to Estimate Reference Crop Evapotranspiration in Humid Areas of Southwest China

: Reference crop evapotranspiration (ET 0 ) is an important indicator for precise regulation of crop water content, irrigation forecast formulation, and regional water resources management. The Hargreaves model (HG) is currently recognized as the simplest and most effective ET 0 estimation model. To further improve the prediction accuracy of the HG model, this study is based on the data of 98 meteorological stations in southwest China (1961–2019), using artificial bee colony (ABC), differential evolution (DE) and particle swarm optimization (PSO) algorithms to calibrate the HG model globally. The standard ET 0 value was calculated by FAO-56 Penman–Monteith (PM) model. We compare the calculation accuracy of 3 calibrated HG models and 4 empirical models commonly used (Hargreaves, Priestley–Taylor, Imark–Allen and Jensen–Hais). The main outcomes demon-strated that on a daily scale, the calibrated HG models ( R 2 range 0.74–0.98) are more accurate than 4 empirical models ( R 2 range 0.55–0.84), and ET 0-PSO-HG has the best accuracy, followed by ET 0-ABC-HG and ET 0-DE-HG , with average R 2 of 0.83, 0.82 and 0.80, average RRMSE of 0.23 mm/d, 0.25 mm/d and 0.26 mm/d, average MAE of 0.52 mm/d, 0.53 mm/d and 0.57 mm/d, and average GPI of 0.17, 0.05, and 0.04, respectively; on a monthly scale, ET 0-PSO-HG also has the highest accuracy, followed by ET 0-ABC-HG and ET 0-DE-HG , with median R 2 of 0.96, 0.95 and 0.94, median RRMSE of 0.16 mm/d, 0.17 mm/d and 0.18 mm/d respectively, median MAE of 0.46 mm/d, 0.50 mm/d, and 0.55 mm/d, median GPI of 1.12, 0.44 and 0.34, respectively. The calibrated HG models (relative error of less than 10.31%) are also better than the four empirical models (relative error greater than 16.60%). Overall, the PSO-HG model has the most accurate ET 0 estimation on daily and monthly scales, and it can be suggested as the preferred model to predict ET 0 in humid regions in southwest China regions. estimation accuracy of 3 improved HG models and the Priestley–Taylor, Imark–Allen and Jensen–Haise models on different time scales; and (iii) propose recommended models for ET 0 estimation in Southwest China at different scales.


Introduction
Reference crop evapotranspiration (ET0) is a major parameter for climate researches and agricultural water resources management [1]. This essential indicator is not only important to measure water balance and conversion, but also plays a significant role in the process of energy exchange between the earth and the atmosphere [2,3]. Therefore, accurate estimation of ET0 is valuable for the estimation of agricultural water demand, irrigation forecast, the management of regional water resources and budget [4,5].
In recent years, many mathematical models can be used to estimate the ET0, including FAO-56 Penman-Monteith (PM) [6], Imark-Allen [7], Priestley-Taylor [8], Hargreaves [9], Jensen-Haise [10] and Makkink [11]. Among them, the Penman-Monteith model is the most common one applied in the study of evapotranspiration, and there are many variations about it [12,13]. The PM model is applied to various environmental and climatic conditions without local or global correction, and the estimation result is accurate. Therefore, it is often used as the standard to verify other models [13]. The main limitation of the PM model is that it requires many meteorological elements [14]. The most ideal process, then, is to estimate ET based on as few input elements as possible without affecting the accuracy [15,16].
The radiations model has also attracted the attention of many researchers because it requires less input data and the temperature is easy to measure [17]. The famous Hargreaves (HG) model estimates ET0 only by inputting the maximum and minimum temperature and solar radiation, so this method is considered to one of the simplest and accurate methods to estimate ET0 [18][19][20][21]. Hence, when the data required by the PM model cannot be fully obtained, the short-term prediction of ET0 can be done by using the HG model. Liu et al. [22] evaluated 16 methods for estimating ET0 based on radiation models and among them, the HG model has more accurate average performance indicators in the arid, semi-arid, temperate, also in frigid and polar climates. At the same time, some scholars pointed out that the HG model needs global calibration, which is suitable for the estimation of ET0 in long time series. Almorox et al. [23] and Feng et al. [24] found that under low humidity or strong wind conditions (u > 3 m/s) the HG model estimation result is relatively low, while the estimation result ET0 is relatively high when the levels of humidity are high. Therefore, the HG model might need to be calibrated to accurately estimate ET0.
In recent years, with the development of AI technology, machine learning optimization algorithms have also provided some useful information to the field. Feng et al. [25] discussed four different machine learning models: extreme learning machine (ELM), generalized regression neural network (GRNN), back-propagation neural network (BPNN) and random forest (RF) model and predict soil temperature (TS) on different depths of 2 cm, 5 cm, 10 cm and 20 cm on the Loess Plateau and the results showed that ELM, GRNN, BPNN and RF models can accurately estimate TS at different depths. Zhu et al. [26] used particle swarm optimization (PSO) to optimize the extreme value learning machine (ELM) model under the condition of limited meteorological data and proposed a hybrid model (PSO-ELM) to accurately predict daily ET0 in the arid area of Northwest China. In the past few decades, several scholars have done a lot of work on the calibration of the HG model [27,28], however, these calibrations are specific to some locations and cannot be applied to other places with totally different climatic conditions. Moreover, the structure of the improved HG model is more complex than the original model [29][30][31]. With that being said, it is possible to shape the aims of this research as (i) calibrate the 3 empirical parameters (C, m, a) in the HG model based on artificial honeycomb (ABC), differential optimization (DE) and particle swarm (PSO) optimisation algorithms; (ii) evaluate the ET0 estimation accuracy of 3 improved HG models and the Priestley-Taylor, Imark-Allen and Jensen-Haise models on different time scales; and (iii) propose recommended models for ET0 estimation in Southwest China at different scales.

Study Area and Weather Data
The southwest region of China mainly includes the Northwest Sichuan Plateau (NSP), the Sichuan Basin (SB), the Guangxi Basin (GB) and the Yunnan-Guizhou Plateau (YGP). The SB is one of the main grain-producing areas in the country, with an area of about 260,000 km 2 and a population of 90 million. The famous Dujiangyan project located in the middle of the SB, with an irrigation area of about 67,000 km 2 . The NSP is one of the five largest pastoral areas in China and the largest animal husbandry base in Sichuan, covering an area of approximately 166,000 km 2 . The terrain of the GB is higher in the northwest and lower in the southeast, and the elevation is generally 80 to 200 m, mostly limestone, covering an area of about 70%, with obvious karst landforms. The YGP located in southwestern China and is one of the four major plateaus with a total area of about 500,000 km 2 . This specific distribution is shown in Figure 1. The meteorological data used in this article comes from the National Climate Centre of the China Meteorological Administration (http://data.cma.cn). The southwest region mainly includes 98 meteorological stations-more information is shown in Table 1. The obtained data are daily, including the maximum and minimum temperature, relative humidity, sunshine hours, wind speed and other aspects measured since 1961. At present, there are many methods to calculate ET0, in which the physical basis of Penman-Monteith model is relatively strict and the calculation accuracy is high. Therefore, FAO 56 Penman-Monteith formula (PM) recommended by FAO was used to verify the calibration results of Hargreaves model [32]. The specific expression of the formula is: where ET0 is reference evapotranspiration (mm/d), Rn is net radiation (MJ/m 2 d), G is soil heat flux density(MJ/m 2 d), Tmean is mean air temperature (°C), es is saturation vapour pressure, kPa, ea represents actual vapour pressure, kPa, ∆ is the slope of the saturation vapour pressure function (kPa/°C), γ means psychometric constant (kPa/°C), and U2 stands for wind speed at 2 m height (m/s). The detailed calculation process of this model is given in Allen, R.G. [6].

Priestley-Taylor Model
Priestley and Taylor [33] derived Priestley-Taylor based on solar radiation and soil heat flux. The specific expression is where ET0 represents reference evapotranspiration (mm/d), α is a constant 1.26, λ is the latent heat of vaporization (MJ/kg), ∆ is the slope of the saturation vapour pressure function (kPa/°C), Rn is the net radiation (MJ/m 2 d), G is soil heat flux density (MJ/m 2 d) and γ represents the psychrometric constant (kPa/°C).

Imark-Allen (IK)
Irmak [34] proposed the Irmak-Allen model, also based on solar radiation data. The final equation for this model is = 0.489 + 0.289 + 0.023 where Rn is the net radiation (MJ/m 2 d) and T means air temperature (°C).

Jensen-Haise Model (JH)
Still in the same area, Jensen and Haise [35] evaluated many evapotranspiration observations through soil sampling and derived the Jensen-Haise model, proposing the following final equation: where Rs is solar radiation (MJ/m 2 d), T is air temperature (°C).

Hargreaves Model
Hargreaves and Samani were the firsts to introduce daily air temperature range (TR = Tmax − Tmin) and Ra to estimate Rs [36,37]: where Rs is solar radiation, which has the same units as ET0 (MJ/m 2 d), KRS is the empirical coefficient fitted to Rs/Ra versus ∆ data with the value of 0.16, Ra is the extraterrestrial radiation (MJ/m 2 d) and ∆ goes for daily air temperature range in degrees Celsius. Based on Equation (5), Hargreaves and Samani [37] further obtained the original HG model: where C, a and m are empirical coefficients, with recommended values of 0.0023, 17.8 and 0.5, respectively. Ra is the extraterrestrial radiation (MJ/m 2 d), Tmean is mean air temperature (°C), and is the daily temperature range (°C).

Artificial Bee Colony
Artificial Bee Colony (ABC) divides bee colony into three categories by simulating the honeybee's honey collecting mechanism: honeybee, observation bee and reconnaissance bee. Assuming that the solution space of the problem is D-dimensional, the numbers of collecting bees and observing bees are both SN and are equal to the number of nectar sources [37][38][39][40][41][42][43]. The location of each nectar source represents a possible solution, and the amount collected from the source corresponds to the fitness of the corresponding solution, and there is a one-to-one correspondence between the honeybee and the nectar source. The honeybee corresponding to the i-th nectar source searches for a new one according to the following formula: where i = 1, 2, ..., The ABC algorithm compares the newly generated possible solution ′ = ′ , ′ , . . . , ′ with the original solution = , , . . . , , and uses a greedy selection strategy to retain better solutions. In cases like this, each observation bee chooses a honey source according to the probability, and the formula for this process is The fitness value of fitj is possible solution Xi, the observation bee searches for a new doable solution according to the probability formula. If the fitness value of a honey source is not improved in a given step, the honeybee will be discarded and the corresponding one will become a scout bee. The scout bee searches for new possible solutions by the following formula: where the interval of γ is a random number on [0, 1], and x d max are the lower and upper bounds of the D dimension.

Differential Evolution Algorithm
Differential Evolution (DE) is a heuristic random search algorithm based on group differences [44,45]. The calculation process mainly goes through 4 steps of initialization, variation, crossover and selection.
(1) Initialization Differential evolution algorithm needs to initialize the population, by following this expression: where Xi(0) is the i-th individual, j represents the j-th dimension; L j i x , and U j i x , are the lower and upper bounds of the j-th dimension, respectively, and rand (0, 1) is the space for the random number in the interval [0, 1].

(2) Variation
The individual mutation is achieved through a different strategy. This common difference strategy is to randomly select two different individuals in the population, and then scale their vector difference to perform vector synthesis with the individual to be mutated: where r1, r2 and r3 are three random numbers, the interval is [1, N] and F stands for the scaling factor that is a certain constant.
(3) Crossover The purpose of crossover operation is to select individuals randomly because differential evolution is also a random algorithm. The method for crossover operation is where CR is called crossover probability, which generates new individuals randomly by means of probability.

(4) Selection
Lastly, in the process of optimization, the greedy selection strategy is adopted to select the better individual as the new individual:

Particle Swarm Optimization
Particle Swarm Optimization (PSO) was first proposed by Eberhart and Kennedy in 1995 and its basic concept originated from the study of bird foraging behaviour [44][45][46]. Particles are used to simulate practical problems and each particle can be regarded as a search individual in N-dimensional search space. Also, the particles only have two attributes: speed and position. Speed represents the speed of moving, and position represents the direction of movement. The optimal solution searched by each particle individually is called the individual extremum and the optimal individual extremum in the particle swarm and is used as the current global optimal solution. It also iterates continuously, updating the speed and position, and finally obtains the optimal solution with satisfies the termination condition [47]. The formula for updating speed and position is where ω is the inertia coefficient. When it is large, the global and local optimization ability is strong; when it is small, the global optimization ability is weak, and the local optimization ability is strong. Therefore, the global and local optimization capabilities can be adjusted by adjusting the size of ω , C1 and C2 are acceleration numbers, which are the learning factors of individual particles.

Parameter Optimization Process
Equation (6) is reasonable in form for different types of climate regions, but many studies have found that these coefficients are different in regions with completely opposite climate conditions [48]. Therefore, the three parameters in the Hargreaves model need to be regionally corrected, and the original results should be kept as brief as possible.
The ABC, DE and PSO optimization algorithms are used to correct the above three parameters (C, m and a), requiring some specific steps for the proper correction of the Hargreaves model such as (1) Divide the calibration period and inspection period. The continuous daily meteorological data of each station, in this case, are divided into two parts used as calibration samples (L1) and validation samples (L2). (2) Determine the feasible region of the variable to be optimized. After analysis and debugging, the three parameters for this article were respectively taken: ∈ 5 × 10 , 0.02 , m∈ 0.02,2.0 and ∈ 2.0,85.0 .
(3) Determine the optimization objective function. To ensure that the Hargreaves model has high simulation accuracy and generalization ability at the same time after correction, the minimization of the following function F is taken as the optimization objective of the three optimization algorithms.
where CD1 and CD2 represent the Nash-Sutcliffe coefficients of the periodic rate and the inspection period respectively, and the calculation formulas used in this case are: where ETH(t) and ET0(t) are the ET0 of calculated by Hargreaves and PM formulas, and ET0 represents the value of ET0(t).
(4) The ABC, DE and PSO optimize methods are used with the Hargreaves model to find the minimum value of F function. The undetermined parameters are C, m, a, and the search interval are determined in the second step. The FAO recommended values are 0.0023, 17.8 and 0.5, respectively. When the F function is the smallest, the convergence is the best and the optimal result is output. The specific optimization process is shown below Figure 2:

Performance Evaluation
The determination coefficient (R 2 ), relative root mean square error (RRMSE), mean absolute error (MAE), and Nash-Sutcliffe coefficient (NS) were used to evaluate performances of the models [38,48].
Following the above equation, and are ET0 values at the i-th step obtained by PM and HG models, respectively. n is the number of the time steps and ET0,mean is the value of the PM calculation. The RRMSE is dimensionless, with the value from 0 (the perfect fit) to ∞ (the worst fit). The MAE is measured in mm/d, with the value from 0 (the perfect fit) to ∞ (the worst fit). The NS is also dimensionless, and its value also ranges from 1 (the perfect fit) to −∞ (the worst fit).
We adopted the Global Performance Index (GPI) to combine the results of all the above indicators and avoid the difference of individual indicators [49]. Such an index can be calculated as where R 2 and NS are the main indicators with equal to −1. For other indicators, is equal to 1, yj is the median value of the scale value of index j, yij stands for the scale value of index j of model i. The GPI value gets higher according to the accuracy of the model.

Calibration of Hargreaves
The distribution of three empirical parameters (C, a, m) of the HG model after calibration based on ABC, DE and PSO algorithms is shown in Figure 3a Parameter C is highest in GB, and it is overall a little lower in NSP and YGP than that in SB and GB. The geographical difference of parameter a is easily spotted, and after the optimization of the three algorithms, the main distribution intervals are [2.16, 42.23], [2.00, 46.03] and [2.00, 44.50]. The changes are large in the NSP and the YGP, and the difference between the optimization results in the SB and the GB is small. The parameter m is also slightly affected by regional differences. Moreover, the calibrated value of the three optimization algorithms is slightly higher than the recommended value of 0.5. In this case, the main distribution intervals are [0.18, 0.82], [0.17, 0.80] and [0.18, 0.80], and the difference among optimization results is insignificant. Overall, the YGP and the NSP showed complex topographic structures with large undulations. The parameter a is greatly affected by regional differences, followed by the parameter c, and the parameter m that is basically not touched by the terrain undulations, indicating that the accuracy of the HG model is mainly influenced by parameter a.  Table 2. It can be concluded that the parameter c calibrated in the 4 regions of the SW of the three optimization algorithms is larger in the GB, and the difference between the other regions is small. The change of the parameter m is not much and is closer to the FAO recommended value of 0.50. The parameter a fluctuation range is also larger in the NSP, but smaller in the GB.

Performances of HG Models on a Daily Basis
The accuracy of daily ET0 estimated by seven models in Southwest China and four regions is represented in Table 3. Observing the results, it can be said that in the Southwest region, the optimization models (ABC-HG, DE-HG and PSO-HG) are more accurate than 4 empirical models (HG, PT, IA and JH) and the PSO-HG model is the one with the highest accuracy. In SW area, PSO-HG model had the highest accuracy, followed by ABC-HG and DE-HG models, and and −0.01, respectively. Among the three optimization algorithms, the PSO algorithm was the one with the highest accuracy after calibration, because the particle speed determines the direction and distance of the particle movement in it, and its speed is dynamically adjusted according to the movement experience and other particles, to achieve the individual that does the best in the solvable space. Therefore, it is suggested to use PSO-HG model to estimate the daily ET0 of Southwest China and four regions, since it was the one with more accurate events in all the examinations.  Based on the HG model calibrated by ABC, DE and PSO optimization algorithms and 4 empirical models, the daily ET0 accuracy indexes estimated in SW, NSP, YGP, GB and SB regions are represented in According to the GPI estimated by all models, the accuracy of PSO-HG model was the highest in the whole southwest region, followed by ABC-HG and DE-HG model, with median lines of 1.54, 0.94 and 0.91, respectively. However, the accuracy of the JH model was relatively low, followed by HG and I-A models, with median lines of −0.79, 0.13 and −0.12. In summary, the accuracy of the optimized HG model in the southwestern region was higher than the empirical model, and the accuracy of the PSO-HG model is higher than the other models.

Performances of HG Models on a Monthly Basis
Based on ABC, DE and PSO optimization algorithms, the monthly average ET0 estimated in Southwest China and the four regions are properly explained in Figure 5 and the monthly average ET0 estimated by different models showed a parabolic change. For SW, NSP, SB, GB and YGP, the error range of HG model calibrated by optimization algorithm a were 15.40 to 3.40%, −28.17 to 2.10%, −14.32 to 3.01%, −33.75 to 5.42% and −12.41 to 10.76%. Also, the error ranges of 4 empirical models were −21.51 to 11.81%, −51.95 to 17.33%, −22.04 to 9.07%, −48.09 to 20.68% and −20.09 to 20.70%. These results showed that the accuracy of the HG model improved by the optimization algorithm in estimating monthly scale ET0 and it was obviously better than that of the commonly used empirical model. In addition, the accuracy of the estimated monthly ET0 of the seven models during the period from May to September was relatively high and the relative errors in SW, NSP, SB, GB and YGP were lower than 11.80%, 11.01%, 7.05%, 4.92% and 5.06%, respectively, which may be due to the complex topographic structure of the southwest region and the small difference in surface temperature. The 4 empirical models had higher accuracy in summer and lower accuracy in winter, because the terrain structure there is complex, and winter is mostly cloudy and rainy. On the YGP, the terrain structure is mostly karst landforms, caused by the different degrees of soil evaporation, which leads to differences in relative humidity. Overall, the PSO-HG model had the highest accuracy of ET0 estimated by SW, NSP, YGB, SB and GB, with average relative errors of −2.69%, 3.48%, −1.64%, 6.10% and −0.92%. Therefore, the PSO-HG model might be the best fit to estimate monthly scale ET0 in Southwest China. According to the GPI estimated by all models, the PSO-HG model was the most assertive one, followed by the ABC-HG and DE-HG models, with median lines that marked 1.13, 0.47 and 0.35, respectively. However, the JH model showed lower accuracy, followed by IA, P-T and HG model, with median lines as −0.82, −0.64, −0.04 and 0.02, respectively. Hence, in the southwest region, the accuracy of the optimized HG model is higher than that of the empirical model, and the accuracy of the PSO-HG model surpasses the other models. Figure 7 shows Taylor diagrams that represent the relationship among the standard deviation, root mean square error, and determination coefficient of each model, and intuitively shows the estimated monthly ET0 accuracy of the 7 models at each weather station site. This also shows that the monthly ET0 estimated by the PSO-HG model is the most precise one.

Discussion
This research found that a model improved via an intelligent optimization algorithm has a higher estimation accuracy than that of the original empirical models. The optimized HG models overcome the errors caused by topographical and climate factors (temperature, air pressure, and radiation), which possess regional characteristics that it could make corresponding adjustments in different elevations, topographic structures and land types, making its accuracy significantly improved. Wu et al. [49] used four models optimized by intelligent optimization algorithms (GA, ACO, CSA and FPA) to estimate the daily ET0 in different climate zones in China, and found a significant improvement in estimating values of daily ET0 in the temperate monsoon and (sub) tropical monsoon climates with the models optimized by the four algorithms, while the RMSE decreased by 14.0%, 25.1%, 31.4% and 33.1%, respectively. Additionally, under temperate continental and mountainous plateau climates, the RMSE also declined by 5.2%, 9.8%, 12.9% and 14.1%, respectively. Still in this realm, Hai et al. [50] studied 3 stations in Burkina Faso to get a better algorithm to improve the fuzzy model to predict regional daily ET0 and found that the estimation accuracy of the optimized hybrid model was revamped significantly. Based on the original model, and under the complete input of meteorological data, the RMSE value was between 0.24 and 0.40 mm/d, and the R 2 value was between 0.95 and 0.97. Yin et al. [51] used a genetic algorithm to optimize support vector machine hybrid model to predict ET0 in China arid areas, and the results showed that the RMSE and R 2 values of the support vector machine model upgraded by this new algorithm were 0.138 mm/d and 0.995 respectively using complete meteorological data. Its accuracy was also improved by about 10%, compared with the original support vector machine model. Fan et al. [52] used support vector machine firefly calculation (SVM-FFA), Copula based nonlinear quantile regression (CNQR) and empirical model to estimate daily radiation. The findings from this process showed that, compared with CNQR, the average SVMFFA decreased by 0.67% (0.01 MJ·m 2 /d), and the average R 2 increased by 0.43% (0.004 MJ·m 2 /d). As to the training time of the three methods, the calculation time of SVM-FFA (1.68 s) was less than that of CNQR (6.68 s), but the parameter optimization time of SVM-FFA (4.91 s) marked 105 times that of CNQR. However, some scholars think that the accuracy of the optimized models has not improved compared to traditional models [53]. The reasons for this phenomenon might be that in some biologically inspired algorithms coupled with the model course or in the optimization process, it may fall into a local optimum, resulting in the lower accuracy than the original model. Besides that, it may be that when the model is boosted, the level of the researcher subjective judgment involved in the parameter adjustment process affects the final prediction progress of the model.
It was found that the PSO-HG model has the highest accuracy in estimating ET0, followed by ABC-HG and DE-HG models. The main difference between PSO, ABC and DE optimization algorithms lies in the global convergence process of parameter calibration. The PSO optimization algorithms use nonlinear changes to adjust the samples, while ABC optimization algorithms adopt the dynamic adjustment. DE optimization algorithms involve crossover, mutation and selection operations, which helps to improve the performance accuracy of the model. It revamps the model parameters by looking for the optimal parameter solution by the nonlinear fitting. Therefore, the principle of the PSO algorithm is consistent with the nature of finding the optimal frameworks.

Conclusions
In this paper, we used ABC, DE and PSO optimization algorithms based on meteorological data obtained from 98 stations from 1961 to 1991 improve the HG model for estimating ET0 in Southwest China. The meteorological data from 1992 to 2019 was checked to evaluate the applicability of the calibrated HG model, HG model, PT model, Imark-Allen model, and JH model on different time scales, and it helped to shape some conclusions.
(1) On daily scale, the calibrated HG models were more accurate compared with the 4 physical models for daily ET0 estimation at the 4 sub-zones of southwest China.

Data Availability Statement:
The meteorological data used in this study are available here: http://data.cma.cn.