Recalibration of four empirical reference crop evapotranspiration models using a hybrid Differential Evolution-Grey Wolf Optimizer algorithm

Abstract

Accurate estimation of reference crop evapotranspiration (ET0) is essential for water resource management and irrigation scheduling. A multitude of empirical models have been employed to estimate ET0, yielding satisfactory outcomes. However, the performance of each model is contingent upon the empirical parameters utilized. This study examines the applicability of four empirical ET0 models, namely the Makkink (Mak), Irmark-Allen (IA), improved Baier-Robertson (MBR), and Brutsaert-Stricker (BS) models. Meteorological data from 24 weather stations across various regions in China were procured and employed to assess the ET0 simulation results. The study employed the Differential Evolution (DE) optimization algorithm, Grey Wolf Optimizer (GWO) algorithm, and a hybrid algorithm that combines DE and GWO algorithms (DE-GWO algorithm) to optimize the parameters of the four empirical models. The findings revealed that the optimization algorithms significantly enhanced the regional adaptability of the four models, particularly the BS model. The DE-GWO algorithm demonstrated superior optimization performance (RMSE=0.055-0.372, R2=0.912-0.998, MAE=0.037-0.311, and FS=0.864-0.982) compared to the DE (RMSE=0.101-2.015, i?2=0.529-0.997, MAE=0.075-1.695, and FS=0.383-0.967) and GWO (RMSE=0.158-0.915, i?2=0.694-0.987, MAE=0.111-0.701, and FS=0.688-0.947) algorithms. The DE-GWO-optimized BS model was the most accurate and improved, followed by the MBR model. The IA and Mak models also showed slightly better performance after optimization with the DE-GWO algorithm. The DE-GWO-optimized BS model performed better in the southern agricultural region than in other regions. It is recommended to utilize the DE-GWO to enhance the accurate prediction of empirical ET0 models across the nine agricultural regions of China.

Headnote

Abstract: Accurate estimation of reference crop evapotranspiration (ET0) is essential for water resource management and irrigation scheduling. A multitude of empirical models have been employed to estimate ET0, yielding satisfactory outcomes. However, the performance of each model is contingent upon the empirical parameters utilized. This study examines the applicability of four empirical ET0 models, namely the Makkink (Mak), Irmark-Allen (IA), improved Baier-Robertson (MBR), and Brutsaert-Stricker (BS) models. Meteorological data from 24 weather stations across various regions in China were procured and employed to assess the ET0 simulation results. The study employed the Differential Evolution (DE) optimization algorithm, Grey Wolf Optimizer (GWO) algorithm, and a hybrid algorithm that combines DE and GWO algorithms (DE-GWO algorithm) to optimize the parameters of the four empirical models. The findings revealed that the optimization algorithms significantly enhanced the regional adaptability of the four models, particularly the BS model. The DE-GWO algorithm demonstrated superior optimization performance (RMSE=0.055-0.372, R2=0.912-0.998, MAE=0.037-0.311, and FS=0.864-0.982) compared to the DE (RMSE=0.101-2.015, i?2=0.529-0.997, MAE=0.075-1.695, and FS=0.383-0.967) and GWO (RMSE=0.158-0.915, i?2=0.694-0.987, MAE=0.111-0.701, and FS=0.688-0.947) algorithms. The DE-GWO-optimized BS model was the most accurate and improved, followed by the MBR model. The IA and Mak models also showed slightly better performance after optimization with the DE-GWO algorithm. The DE-GWO-optimized BS model performed better in the southern agricultural region than in other regions. It is recommended to utilize the DE-GWO to enhance the accurate prediction of empirical ET0 models across the nine agricultural regions of China.

Keywords: reference crop evapotranspiration, empirical model, nine agricultural regions of China, hybrid algorithm, Brutsaert-Stricker

1 Introduction

Water scarcity has emerged as a pressing global concern, particularly for arid and semi-arid regions. Amidst the backdrop of exponential population growth and escalating pollution, nations are increasingly striving to optimize water resource management. This endeavor is being pursued through the implementation of innovative strategies and technologies aimed at enhancing the efficiency and sustainability of water use'11. According to the survey, population growth by 2050 will consume about double the current demand for food, and to meet food security in the coming decades, a major reform of agricultural water management is needed'21. Therefore, effective water management irrigation and water resources allocation are the best ways to improve the rational use of regional water resources. ET is an essential parameter in the hydrological cycle. An accurate ET estimation is the basis for efficient agricultural water use, which is crucial for water resource management and agricultural water-saving irrigation'31. ET can be directly measured using field instruments; however, this is expensive and difficult. Mathematical models for estimating the ET have been proposed by several researchers, including Priestley-Taylor, Blaney-Criddle, Turc, and Hargreaves'461. Such models can facilitate ET measurement. The FAO-56 Penman-Monteith (FAO-56 PM) model is the most accurate for estimating ET'781. The ET estimated by the FAO-56 PM model can be used as a benchmark for validating the other models'9111. However, the FAO-56 PM model requires a considerable number of parameters for ET estimation, including sunshine hours, air temperature, radiation, humidity, and other meteorological data. In practical settings, meteorological data are often missing, which affects the efficiency of ET estimations'121.

Hence, it is imperative to develop a simpler approach to accurately calculate ET. The Mak model is precise to a point where it can almost represent the universally accepted Penman-Monteith model when the ET process is applied'13151. Gunston and Batchelor reported that the Priestley-Taylor model is a suitable alternative to the FAO-56 PM model in tropical climates with high humidity'161. More than 60 models have been proposed for estimating the ET, including empirical models based on temperature, radiation, humidity, water pressure, and hybrid models'17201. Multiple empirical models have been integrated to validate their reliability. Tabari et al. compared various ET models under humid conditions in the northern parts of Iran and found that the radiation-based models were the best among the models studied, followed by the temperature-based and pan evaporation-based models'211. Zhao et al. used FAO-56 PM as the benchmark model and evaluated different ET models (including HG) in a semi-arid region of China and concluded that the HG was the best ET model in the studied region'221. Majidi et al. tested and ranked the applicability of various models by applying the Bowen ratio energy balance method and found that the Mak model (which is based on radiation and temperature data) had good applicability for evaporation estimation'231. Feng et al. investigated the applicability of three ET empirical models in southwest China and observed that the Priestley-Taylor model was more accurate than the Mak and Ritchie models'241. The results of the empirical models depend on the empirical parameters, and regions differ in their characteristics. Moreover, it is difficult to obtain precise ET estimates on different spatial scales. To address these challenges, optimization algorithms have been extensively used to optimize model parameters and improve the prediction accuracy'25-261. Wu et al. evaluated 12 daily-scale ET models (four temperature models, five radiation models, and three hybrid models) by using the whale optimization algorithm in the four climate zones in China, revealing that the hybrid empirical model was the most accurate'271. Yang et al. used the least squares method to optimize the Rohwer, Romanenko, Turc, Hargreaves, and Penman models and reported that the models performed better after calibration with the least squares method'201. Zhang et al. estimated ET in Northwest China using the Mak model, for which the parameters were optimized using the least squares method, and noted that the model increased when optimized using the least squares method'281. The Firefly Algorithm (FA) has strong robustness and reliability. It can be used to predict ET in the West African deserts. Inspired by nature, FA has significantly improved the performance of the classical Adaptive Neuro-Fuzzy Inference System (ANFIS) model optimization in practical applications'29301. However, no single optimization algorithm has good optimization effects for all constrained optimization problems because each evolutionary algorithm has its own advantages in solving constrained optimization problems. Therefore, hybrid optimization algorithms are useful in situations in which single optimization algorithms are ineffective. Mohammadi and Mehdizadeh coupled Support Vector Regression (SVR) and the Whale Optimization Algorithm (WOA) to build a new ET model. The results indicated that the performance of the hybrid model was better than that of the pure SVR model'261. Long et al. used a hybrid Grey Wolf Optimization algorithm to solve high-dimensional complex functions, and noted that this algorithm exhibited significantly higher accuracy and a more favorable convergence speed than Particle Swarm Optimization (PSO) and universal gravitational search algorithms. Accordingly, hybrid algorithms exhibit superior performance and precision compared with single algorithms.

Consequently, this research has chosen four evapotranspiration

(ET) models, namely, the Makkink (Mak), Irmark-Allen (IA), improved Baier-Robertson (MBR), and Brutsaert-Strieker (BS) models, and optimized them using the Differential Evolution-Grey Wolf Optimizer (DE-GWO), Grey Wolf Optimizer (GWO), and Differential Evolution (DE) algorithms. In order to assess the adaptability of these ET models across various agricultural regions and to validate the enhancement in performance achieved by the different optimization algorithms on empirical formulas, this study collected data from 24 locations across nine distinct agricultural regions in China. The objectives of this study are threefold: 1) to propose a hybrid optimization algorithm for the optimization of ET estimation models, 2) to compare the optimization accuracy of the hybrid algorithm with that of the individual constituent algorithms, and 3) to evaluate the adaptability of each model in the nine agricultural regions based on their respective accuracies.

2 Materials and methods

2.1 Study area

Based on China's Comprehensive Agricultural Regional Planning, the country is divided into 10 first-level agricultural regions and 38 second-level regions. The 10th region is the marine aquatic region; the first nine regions are comprehensive agricultural regions and are generally referred to as the nine major comprehensive agricultural areas, which are outlined as follows: the Northern Arid and Semiarid Region (NASR), the Northeast China Plain (NCP) area, the Yunnan-Guizhou Plateau (YGP) area, the South China (SC) area, the Sichuan Basin and surrounding (SBSR) area, the middle and lower reaches of the Yangtze Plain (MLYP), the Qinghai-Tibet Plateau (QTP) area, the Loess Plateau (LP) area, and the Huang-Huai-Hai Plain (HHP) area. In this study, a total of 24 stations were selected in these nine major agricultural areas.

2.2 Data collection

For model calibration and performance evaluation, meteorological data measured at the same sites were collected from the China Meteorological Data Sharing Service System (http:// cdc.nmic.cn) for the period from January 1, 1960, to December 31, 2019. The meteorological dataset (latitude, longitude, altitude, average temperature, relative humidity, wind speed (U), and maximum and minimum temperatures) were divided into two parts. All the data were collected from each site. The continuous daily meteorological data from each station were divided into two parts, according to the time series, as the training and test samples of the model. The first part, comprising 80% (1960-2007) of the dataset, was used for model calibration, and the second part, comprising 20% (2008-2019) of the dataset, was used for model testing.

2.3 Models for estimating ET 2.3.1 FAO-56 PM model

The FAO-56 PM model is based on the principles of aerodynamics and energy balance, and it comprehensively considers various factors that influence ET. The FAO-56 PM model did not require parameter correction to account for regional climatic differences. The accuracy of the FAO-56 PM model was verified using global meteorological data. Therefore, this study used the ET estimates provided by the FAO-56 PM model as a benchmark for comparison with other models. The FAO-56 PM model can be expressed as follows'31"11: crop surface net radiation, MJ/m2-d; G is the soil heat flux density, MJ/m2-d; T is the daily mean air temperature at a height of 2 m, °C; u2 is the daily mean wind speed at a height of 2 m, m/s; es is the saturated vapor pressure, kPa; ea is the actual vapor pressure, kPa; A is the slope of the curve of pressure versus temperature, kPa/°C; and y is the hygrometric constant, kPa/°C. 2.3.2 Empirical model

This study used the FAO-56 Penman-Monteith formula (PM) recommended by the FAO to verify the optimization results of the empirical model. The four empirical models (Mak, IA, MBR, and BS) were compared with the FAO-56 PM model and optimized to evaluate their applicability. Finally, the four empirical models were optimized using three optimization algorithms, and their performance and adaptability in different regions were observed through evaluation indicators. Each of these models is presented in Table 1.

2.4 Optimization algorithms

2.4.1 DE algorithm

The DE algorithm'361 is a stochastic and population-based optimization method to determine the overall optimal solution in a multidimensional space'371. The operating procedures are described as follows'381:

Initialization: The space is regarded as comprising N-dimensional (D) vectors, which implies that the population size is NP. The initial population is generated randomly.

where, i represents the z'th individual; j represents the jth component; x' represents the lower bound of the jth component; and x" represents the upper bound of the jth component.

Mutations: A mutation operator is employed to generate a new target individual.

Crossover: A crossover operation is performed between each individual and the offspring mutation vector they produce. More specifically, for each component, the algorithm probabilistically chooses either the offspring mutation vector or the original vector to create the trial individual.

Selection: The selection process employs a greedy strategy, where individuals with lower fitness function values are chosen from both the experimental and current populations to advance to the next generation.

2.4.2 GWO algorithm

The Grey Wolf Optimizer (GWO) algorithm is a metaheuristic approach devised by Mirjalili et al. Within this framework, grey wolves are classified into hierarchical roles: alpha (a), beta (/?), delta (<5), and omega (a) wolves. The alpha wolves hold dominance

over the lower-ranking wolves, who, in turn, follow the lead of their superiors'39"411.

1) Tracking, chasing, and approaching prey: During the hunting process, grey wolves strategically encircle their prey. This envelopment behavior can be mathematically modeled as follows:

where, X is the position of the grey wolf; t is the current number of iterations; jl is the position of the prey; and D is the distance between the grey wolf and its prey.

2) Pursuing, encircling, and harassing the prey until it stops moving: In the GWO algorithm, a represents the optimal grey wolf, P represents the second-best grey wolf, S represents the third-best grey wolf, and w represents the rest of the grey wolves. The algorithm is established on the basis that a is the optimal solution, and that/3 and the other wolves know the location of the prey. In the iterative process, a, /?, and S were used to guide the movement of w to achieve global optimization. The positions Xa,Xfi, and Xs of a, /?, and S are used to update the positions of all grey wolves using the following equations (7=1, 2, 3):

3) Attack toward prey: Grey wolves conclude the hunt by striking the prey when it ceases to move. In Equation (8), t denotes the current iteration count, while T signifies the predefined maximum number of iterations. As a diminishes from 2 to 0, the corresponding value of a also varies within the range [-a, a]. With an increasing value of a, the grey wolves venture further afield from their pack members in search of more promising prey, thereby encouraging a global search (|a|>l). Conversely, as a diminishes, the grey wolves draw closer to the prey, steering the wolves towards a more localized search (|a|<l).>

2.4.3 Hybrid DE-GWO algorithm

Taking into account the limitations of the Differential Evolution (DE) and Grey Wolf Optimizer (GWO) algorithms, including premature convergence, inadequate stability, and susceptibility to getting stuck in local optima, this research introduces a more effective hybrid optimization algorithm: the DE-GWO algorithm. This algorithm aims to enhance the global search capability and circumvent the aforementioned drawbacks associated with DE and GWO. The operational steps of the DE-GWO algorithm are detailed below.

1) Population initialization: The algorithm commences by configuring the essential parameters for the hybrid optimization algorithm, which include the population size N, the maximum number of iterations MAX_iter, the crossover probability pCR, the search dimension D, and the search range bounds ub (upper bound) and lb (lower bound), as well as the scaling factor F. It then initializes the parameters a, A, and C, and carries out the DE mutation operation on the population individuals to produce intermediate solutions. Following this, a competitive selection process is conducted to refine the initialized population individuals.

The iteration counter is then set to t= 1.

2) Iteration: The algorithm calculates the objective function value of each grey wolf in the population, sorts the grey wolves according to the size of the objective function value, and selects the top three optimal individualsXa,Xp, andX^.

3) Selection: The algorithm evaluates the distance of each grey wolf in the population from the optimal position and updates their positions based on this assessment. This operation is carried out using a cross-update strategy. Additionally, the parameters a, A, and C are utilized to perform the crossover operation, leading to the generation of new individuals. Following this, the positions of the top three grey wolves are refined with the optimal values obtained.

4) Termination condition: According to the count value, if the maximum number of iterations MAX_iter is reached, the algorithm exits and outputs the objective function value of the global optimal Xa, otherwise, it sets f to t+\ and restarts the execution from step (2). 2.5 Model evaluation metrics

The root mean square error (RMSE), determination coefficient (R2), mean absolute error (MAE), predictive score (FS), and global evaluation index (GPI) were used as metrics to evaluate the performance of each ET model.

where, m, is the rth daily value simulated by the model; v, is the rth standard daily value calculated by the FAO-56 PM model; u is the average value of ",; v is the average value of v,; n is the number of samples in the test set; g, is the normalized value of RMSE, i?2, MAE, and NSE; and ys is the median of the parameter.

3 Results and discussion

3.1 Recalibration of the ET empirical model

The parameters of Mak, IA, MBR, and BS models were optimized using DE-GWO, GWO, and DE algorithms, and the results were analyzed. Initially, the parameter a for the Mak model, as optimized by the three algorithms, remained within a stable range (0.620-0.754) across all study areas. The values were notably higher in the NCP, NASR, and QTP compared to the SC, SBSR, and other regions. For the Mak model optimized by the DE-GWO and GWO algorithms, the parameter b followed similar trends as a, but this was not the case for the DE algorithm, with b values showing greater variability in the NASR and NCP areas (-0.644-0.391).

In the IA model, the parameter a optimized by all three algorithms exhibited considerable variation, mirroring the characteristics of the Mak model's a, yet with less fluctuation in the NCP, NASR, and QTP areas compared to the SC, SBSR, and MLYP areas (-0.734-0.753). The parameter b for the IA model optimized by the DE-GWO algorithm showed the most significant variation in the NASR area. The b values of the IA model optimized by the three algorithms were less variable in the SC, SBSR, and other areas than in the NCP, NASR, and other areas (0.159-0.403), while the c values were considerably smaller and showed little variation (-0.036-0.110).

For the MBR model, the parameters a, b, and d optimized by the three algorithms were generally smaller than c and e. Although a, b, and d showed some variation across the nine areas, the extent of variation was relatively minor (-0.087-0.195, -0.010-0.201, and -0.010-0.201, respectively). However, the parameter c maintained a more stable range (2.051-5.309) in the SC, MLYP, HHP, and NCP areas. The parameter e showed a stable variation range in the LP and HHP areas compared to other regions.

Lastly, for the BS model, the parameter b optimized by the three algorithms varied across all nine areas (0.657-0.896), with the amplitudes of b in the SC, SBSR, and MLYP areas being relatively stable. The parameter a for this model, optimized by the three algorithms, fluctuated within a narrow range (0.140-0.880). Among the parameters, c was the most sensitive across the nine areas. The performance of the three optimization algorithms was not uniform. The c parameter of the model optimized by the DE-GWO algorithm showed the greatest variation in the LP area. The amplitude of c for the model optimized by the GWO algorithm was the largest in the SC area, and the variation in c for the model optimized by the DE algorithm was the most pronounced in the SC and SBSR areas (-18.244-73.099)

In summary, the three optimization algorithms provided relatively similar optimization results for various empirical models in the nine areas. Therefore, the parameters of the models optimized by the three optimization algorithms can provide a basis for obtaining highly precise ET estimates in the nine areas. 3.2 Comparison of different optimization algorithms

In this study, four daily ET prediction models (Mak, IA, MBR, and BS) were constructed and then optimized using three optimization algorithms (DE, GWO, and DE-GWO). Table 2 presents the i?2, RMSE, MAE, and FS values derived for the optimized models for the training and test datasets obtained at each of the 24 stations. The results revealed that the Mak model had the highest accuracy (RMSE=0.496-0.940, i?2=0.767-0.872, MAE=0.343-0.641), followed by the IA model (RMSE=0.500-0.794, i?2=0.748-0.846, MAE=0.425-0.720), MBR model (RMSE=0.803-1.752, i?2=0.588-0.852, MAE=0.664-1.504), and BS model (RMSE=2.649-3.594,i?2=0.383-0.451, MAE=2.226-2.947).

The DE, GWO, and DE-GWO algorithms improved the accuracy of the models to a certain extent. The DE-GWO-optimized BS model exhibited the highest accuracy and greatest improvement after optimization (RMSE=0.055-0.372, i?2=0.912-0.998, MAE= 0.037-0.311, FS=0.864-0.982), followed by the GWO-optimized BS model (RMSE=0.158-0.915, i?2=0.694-0.987, MAE=0.111-0.701, FS=0.688-0.947) and DE-optimized BS model (RMSE=0.101-2.015, i?2=0.529-0.997, MAE=0.075-1.695, FS=0.383-0.967). The DE-GWO-optimized MBR model was slightly less accurate than the DE-GWO-optimized BS model (RMSE=0.399-0.730, i?2=0.855-0.957, MAE=0.303-0.504, FS=0.289-0.787). Therefore, the DE-GWO algorithm resulted in greater improvement in the MBR model than the GWO algorithm (RMSE=0.441-0.814, i?2=0.811-0.956,

MAE=0.297-0.615, FS=0.273-0.742). The DE algorithm showed the least improvement among the three algorithms (RMSE=0.434-0.917, i?2=0.745-0.937, MAE=0.321-0.735, and FS=0.030-0.695). The DE-GWO algorithm slightly improved the accuracy of the IA model, but resulted in lower accuracy improvement levels for the MBR and BS models (RMSE=0.168-0.858, i?=0.825-0.981, MAE= 0.125-0.592, FS=0.040-0.756). The GWO-optimized IA model had relatively low accuracy (RMSE=0.291-0.907, i?2=0.705-0.938, MAE=0.179-0.711, FS=0.012-0.639), as did the DE-optimized IA model (RMSE=0.337-0.915, i?2=0.715-0.917, MAE=0.228-0.745, FS=0.042-0.584). The Mak and DE-GWO-optimized Mak models were the most accurate, followed by the GWO-optimized Mak model (RMSE=0.226-0.764, i?2=0.806-0.975, MAE=0.165-0.542, FS=0.145-0.523) and the DE-optimized Mak model (RMSE=0.242-0.775, i?2=0.792-0.969, MAE=0.188-0.566, FS=0.013-0.489). Although the DE-GWO-optimized Mak model had favorable accuracy levels, its accuracy improvement levels were the lowest among the optimized models (RMSE=0.226-0.737, i?2=0.821-0.975, MAE=0.167-0.49, FS=0.125-0.523).

The three optimization algorithms performed well. The DE algorithm has a simple mathematical structure, is easy to implement, converges rapidly, and is highly robust. It reacts quickly and requires short execution time. Additionally, it has a one-to-one selection strategy and an optimal mutation strategy, which allows it to converge faster than other algorithms, reducing the complexity of the population more quickly and, in turn, improving its search capability. If it is affected by interference, premature phenomena appear in later stages, which worsen the results. The GWO algorithm is more complex, converges more rapidly, has a better search capability, and is better at optimization than the DE algorithm. During the search process, all three individuals unavoidably fell into the local optimum, which considerably reduced the accuracy of the overall search results of the model, which in turn reduced the accuracy and optimization ability. The DE-GWO algorithm combines the advantages of DE and GWO algorithms. This improves the global search ability by avoiding the introduction of defects caused by the premature stagnation of individuals that have fallen into the local optimum. The DE-GWO algorithm breaks through the local limit value easily, is precise, converges well, and is better at optimization than the GWO and DE algorithms. The three optimization algorithms were the most effective for the BS model. The RMSE, R2, and MAE values derived for the BS model increased by 92%-96%, 120%-152%, and 92%-97%, respectively, after optimization using the DE-GWO algorithm; by 80%-84%, 107%-118%, and 81%-86%, respectively, after optimization using the GWO algorithm; and by 57%-86%, 83%-104%, and 68%-87%, respectively, after optimization using the DE algorithm.

In summary, the three optimization algorithms exhibited high performance in improving the ET prediction models. The values of the four indicators (RMSE, R2, MAE, and FS) for the models after optimization using the different algorithms are shown in Figure 1. The optimization algorithms can be ordered according to their performance as follows: DE-GWO>GWO>DE. 3.3 Performance of ET models in different agricultural areas

The adaptability of the optimized models differed in the nine areas, and the accuracy of the optimized models was higher than that of the non-optimized models. The terrain and climate of each area influenced the performance of each model. The RMSE, R1, MAE, FS, and GPI metrics were used to evaluate the applicability of each model for each area. The accuracy of each model is shown in the Taylor diagram in Figure 2; in the accuracy evaluation, the FAO-56 PM model was used as a benchmark. The estimates provided by the optimized models closely matched the results provided by the FAO-56 PM model, demonstrating the accuracy of the optimized models.

The DE-GWO-optimized Mak model had the highest accuracy and best adaptability in the MLYP area (RMSE=0.318-0.362, R2=0.942-0.960, MAE=0.237-0.283, FS=0.383-0.459, GPI=1.051-1.431, GPI ranking 5-16). The DE-GWO-optimized Mak model was slightly less accurate in the NCP area, and the GWO-optimized Mak model had the same accuracy as that of the DE-GWO-optimized Mak model in this area. The DE-GWO-optimized Mak model was less accurate in the NASR and LP areas than in the MLYP and NCP areas. Moreover, the DE-GWO-optimized Mak model had the second-highest accuracy in the NASR and LP areas. The optimized models performed moderately in the LJHP, SC, and YGP areas. The DE-GWO-optimized Mak model performed worst in the QTP and SBSR areas. The DE-GWO-optimized IA model was the best in the NCP area (RMSE=0.463-0.719,"2=0.845-0.921,MAE=0.338-0.536, FS=0.058-0.330, GPI=1.150-1.297, GPI ranking 6-9), followed by the NASR, MLYP, and LP areas, YGP, SC, and HHP areas, SBSR area, and QTP area. The DE-GWO-optimized MBR model was the best in the MLYP area (RMSE=0.463-0.719, i?2=0.845-0.921, MAE=0.338-0.536, FS=0.058-0.330, GPI=1.150-1.297, GPI ranking 6-9), followed by the NCP area, HHP and NASR areas, LP and YGP areas, and QTP and SBSR areas. The DE-GWO-optimized BS model was the best in the HHP area (RMSE=0.463-0.719, R2=0.845-0.921, MAE=0.338-0.536, FS=0.058-0.330, GPI=1.150-1.297, GPI ranking 6-9), followed by the NASR and LP areas. The DE-GWO-optimized BS model's performance in the NCP and SC areas was the same as that of the DE-optimized BS model's performance in the MLYP area; however, the performance in the NCP, SC, and MLYP areas was not as high as that in the HHP, NASR, and LP areas. The GWO-optimized BS model adapted moderately in the YGP area and poorly in the SBSR and QTP areas. Although the models had low precision in some areas, they still ranked highly in these areas. For example, the GWO-optimized Mak model had low adaptability in the SBSR area, but the GPI value at the Leibo site was 1.503, ranking fourth. Furthermore, the DE-GWO-optimized IA model had a GPI value of 1.885 at site 59 293, ranking it second.

Overall, the adaptabilities of the models in the nine areas were different. Some models had good adaptability in the NASR and MLYP areas but poor adaptability m the QTP and SBSR areas; their adaptability in other areas was relatively moderate. The adaptability of the independent empirical models was improved using optimization algorithms. The models optimized using the hybrid optimization algorithm outperformed independent models.

4 Conclusions

To improve the accuracy and applicability of empirical ET estimation models, this study selected 24 sites in nine major agricultural regions in China as representative research sites and used three optimization algorithms (DE-GWO, GWO, and DE) to optimize four ET empirical models. Thus, ET estimation models with high accuracy and universality were developed.

The optimized models performed differently in various agricultural areas. The optimization significantly improved the accuracy and applicability of the four empirical models. Compared with the GWO- and DE-optimized models, the DE-GWO-optimized models exhibited significant improvements in accuracy and adaptability. The optimization algorithms can be ordered as follows in accordance with their performance: DE-GWO>GWO>DE. The DE-GWO-optimized BS model exhibited the highest accuracy (RMSE=0.087-0.279, R2=0.966-0.995, MAE=0.062-0.208, FS= 0.371-0.684), followed by the DE-GWO-optimized MBR, DE-GWO-optimized IA, and DE-GWO-optimized Mak models. All three optimization algorithms significantly improved the performance of the models. In general, the optimized models exhibited good adaptability in the NASR and MLYP areas but poor adaptability in the QTP area. The performance of the models in other areas was relatively low, although their accuracy and GPI rankings were high. This study is mainly applicable to nine major agricultural regions in China, and the adaptability to various regions of the world may be slightly biased and insufficient. This study makes a certain contribution to the optimization algorithm in optimizing the parameters of the empirical formula; meanwhile, it can facilitate evapotranspiration research in various regions of China.

Acknowledgements

We thank the National Climatic Centre of the China Meteorological Administration for providing the climate database used in this study. This work was also supported by the National Natural Science Foundation of China (Grant No. 52309050 and Grant No. 32372680), the Key R&D and Promotion Projects in Henan Province (Science and Technology Development) (Grant No. 232102110264), the Key Scientific Research Projects of Colleges and Universities in Henan Province (Grant No. 24B416001), the Innovative Research Team (Science and Technology) in University of Henan Province (Grant No. 23IRTSTHN024), and the Youth Backbone Teacher Project of Henan University of Science and Technology (Grants No. 13450013 and 13450010).

References

[References]

[1] Wang J Z, Raza A, Hu Y G, Ali Buttar N, Shoaib M, Saber K, et al.

Development of monthly reference evapotranspiration machine learning models and mapping of pakistan - A comparative study. Water, 2022; 14(10): 1666.

[2] Mueller N D, Gerber J S, Johnston M, Ray D K, Ramankutty N, Foley J A. Closing yield gaps through nutrient and water management. Nature, 2012; 490: 254-257.

[3] Raziei T, Pereira L S. Estimation of ETo with Hargreaves-Samani and FAO-PM temperature methods for a wide range of climates in Iran.

Agricultural Water Management, 2013; 121: 1-18. [4] Hargreaves G H, Samani Z A. Reference crop evapotranspiration from temperature. Applied Engineering in Agriculture, 1985; 1(2): 96-99.

[5] Heydari M M, Tajamoli A, Ghoreishi S H, Darbe-Esfahani M K, Gilasi H. Evaluation and calibration of Blaney-Criddle equation for estimating reference evapotranspiration in semiarid and arid regions. Environmental Earth Sciences, 2015; 74(5): 4053^063.

[6] Priestley C H B, Taylor R J. On the assessment of surface heat flux and evaporation using large scale parameters. Monthly Weather Review, 1972; 100(2): 81-92.

[7] Charytanowicz M, Kulczycki P, Kowalski P A, Lukasik S, Czabak-Garbacz R. An evaluation of utilizing geometric features for wheat grain classification using X-ray images. Computers and Electronics in Agriculture, 2018; 144: 260-268.

[8] Yong S L S, Ng J L, Huang Y F, Ang C K, Kamal N A, Mirzaei M, et al. Enhanced daily reference evapotranspiration estimation using optimized hybrid support vector regression models. Water Resources Management, 2024; 38(11): 4213-4241.

[9] Allen R G, Pruitt W O, Wright J L, Howell T A, Ventura F, Snyder R, et al. A recommendation on standardized surface resistance for hourly calculation of reference ETo by the FA056 Penman-Monteith method. Agricultural Water Management, 2006; 81(1-2): 1-22.

[10] Zhao L, Zhao X B, Pan X L, Shi Y, Qiu Z M, Li X Z, et al. Prediction of daily reference crop evapotranspiration in different Chinese climate zones: Combined application of key meteorological factors and Elman algorithm. Journal of Hydrology, 2022; 610: 127822.

[11] Allen R G. FA056 Irrigation and Drainage Paper. Crop Evapotranpiration. Irrigation and Drainage, 1998; 300(56): 300.

[12] Djaman K, Irmak S, ASCE M, Kabenge I, Futakuchi K. Evaluation of FAO-56 Penman-Monteith Model with limited data and the Valiantzas Models for estimating grass-reference evapotranspiration in sahelian conditions. Journal of Irrigation and Drainage Engineering, 2016; 142(11): 04016044.

[13] Rawat K S, Singh S K, Bala A, Szabo S. Estimation of crop evapotranspiration through spatial distributed crop coefficient in a semi-arid environment. Agricultural Water Management, 2019; 213: 922-933.

[14] Petkovic B, Petkovic D, Kuzman B, Milovancevic M, Wakil K, Ho L S, et al. Neuro-fuzzy estimation of reference crop evapotranspiration by neuro fuzzy logic based on weather conditions. Computers and Electronics in Agriculture, 2020; 173: 105358.

[15] Cammalleri C, Ciraolo G. A simple method to directly retrieve reference evapotranspiration from geostationary satellite images. International Journal of Applied Earth Observations and Geo information, 2013; 21: 149-158.

[16] Gunston H, Batchelor C H. A comparison of the Priestley-Taylor and Penman methods for estimating reference crop evapotranspiration in tropical countries. Agricultural Water Management, 1983; 6(1): 65-77.

[17] Bormann H. Sensitivity analysis of 18 different potential evapotranspiration models to observed climatic change at German climate stations. Climatic Change, 2011; 104(3): 729-753.

[18] Gharbia S S, Smullen T, Gill L, Johnston P, Pilla F. Spatially distributed potential evapotranspiration modeling and climate projections. Science of The Total Environment, 2018; 633: 571-592.

[19] Li S, Kang S Z, Zhang L, Zhang J H, Du T S, Tong L, et al. Evaluation of six potential evapotranspiration models for estimating crop potential and actual evapotranspiration in arid regions. Journal of Hydrology, 2016; 543: 450-461.

[20] Yang Y, Chen R S, Han C T, Liu Z W. Evaluation of 18 models for calculating potential evapotranspiration in different climatic zones of China. Agricultural Water Management, 2021; 244: 106545.

[21] Tabari H, Grismer M E, Trajkovic S. Comparative analysis of 31 reference evapotranspiration methods under humid conditions. Irrigation Science, 2013; 31(2): 107-117.

[22] Zhao C Y, Nan Z R, Feng Z D. GIS-assisted spatially distributed modeling of the potential evapotranspiration in semi-arid climate of the Chinese Loess Plateau. Journal of Arid Environments, 2004; 58(3): 387^03.

[23] Majidi M, Alizadeh A, Farid A, Vazifedoust M. Estimating evaporation from lakes and reservoirs under limited data condition in a Semi-Arid

Region. Water Resources Management, 2015; 29(10): 3711-3733.

[24] Feng Y, Cui N B, Zhao L, Hu X T, Gong D Z. Comparison of ELM, GANN, WNN and empirical models for estimating reference evapotranspiration in humid region of Southwest China. Journal of Hydrology, 2016; 536: 376-383.

[25] Fang R M, Song S J. Daily reference evapotranspiration prediction of Tieguanyin tea plants based on mathematical morphology clustering and improved generalized regression neural network. Agricultural Water Management, 2020; 236: 106177.

[26] Mohammadi B, Mehdizadeh S. Modeling daily reference evapotranspiration via a novel approach based on support vector regression coupled with whale optimization algorithm. Agricultural Water Management, 2020; 237: 106145.

[27] Wu Z J, Chen X, Cui N B, Zhu B, Gong D Z, Han L, et al. Optimized empirical model based on whale optimization algorithm for simulate daily reference crop evapotranspiration in different climatic regions of China. Journal of Hydrology, 2022; 612: 128084.

[28] Zhang Q W, Cui N B, Feng Y, Gong D Z, Hu X T. Improvement of Makkink model for reference evapotranspiration estimation using temperature data in Northwest China. Journal of Hydrology, 2018; 566: 264-273.

[29] Tao H, Diop L, Bodian A, Djaman K, Ndiaye P M, Yaseen Z M. Reference evapotranspiration prediction using hybridized fuzzy model with firefly algorithm: Regional case study in Burkina Faso. Agricultural Water Management, 2018; 208: 140-151.

[30] Yaseen Z M, Ebtehaj I, Bonakdari H, Deo R C, Danandeh Mehr A, Mohtar W H M W, et al. Novel approach for streamflow forecasting using a hybrid ANFIS-FFA model. Journal of Hydrology, 2017; 554: 263-276.

[31] Chia M Y, Huang Y F, Koo C H. Support vector machine enhanced empirical reference evapotranspiration estimation with limited meteorological parameters. Computers and Electronics in Agriculture, 2020; 175(1): 105577.

[32] Makkink G F. Testing the Penman formula by means of lysimeters. J inst of Water Eng, 1957.

[33] Shiri J, Nazemi A H, Sadraddini A A, Marti P, Fakheri Fard A, Kisi O, et al. Alternative heuristics equations to the Priestley-Taylor approach: assessing reference evapotranspiration estimation. Theoretical and Applied Climatology, 2019; 138(1): 831-848.

[34] Baier W, Robertson G W. Estimation of latent evaporation from simple weather observations. Canadian Journal of Plant Science, 1965; 45(3): 276-284.

[35] Brutsaert W, Strieker H. An advection-aridity approach to estimate actual regional evapotranspiration. Water Resources Research, 1979; 15(2): 443^50.

[36] Storn R, Price K. Differential evolution - A simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization, 1997; 11(4): 341-359.

[37] Roy D K, Lai A, Sarker K K, Saha K K, Datta B. Optimization algorithms as training approaches for prediction of reference evapotranspiration using adaptive neuro fuzzy inference system. Agricultural Water Management, 2021; 255: 107003.

[38] Wu L F, Huang G M, Fan J L, Ma X, Zhou H M, Zeng W Z. Hybrid extreme learning machine with meta-heuristic algorithms for monthly pan evaporation prediction. Computers and Electronics in Agriculture, 2020; 168: 105115.

[39] Mirjalili S, Mirjalili S M, Lewis A. Grey Wolf Optimizer. Advances in Engineering Software, 2014; 69: 46-61.

[40] Rodriguez L, Castillo O, Soria J, Melin P, Valdez F, Gonzalez C I, et al. A fuzzy hierarchical operator in the grey wolf optimizer algorithm. Applied Soft Computing, 2017; 57: 315-328.

[41] Dong J H, Liu X G, Huang G M, Fan J L, Wu L F, Wu J. Comparison of four bio-inspired algorithms to optimize KNEA for predicting monthly reference evapotranspiration in different climate zones of China. Computers and Electronics in Agriculture, 2021; 186: 106211.