Hybridized Adaptive Neuro-Fuzzy Inference System with Metaheuristic Algorithms for Modeling Monthly Pan Evaporation

: Precise estimation of pan evaporation is necessary to manage available water resources. In this study, the capability of three hybridized models for modeling monthly pan evaporation (Epan) at three stations in the Dongting lake basin, China, were investigated. Each model consisted of an adaptive neuro-fuzzy inference system (ANFIS) integrated with a metaheuristic optimization algorithm; i.e., particle swarm optimization (PSO), whale optimization algorithm (WOA), and Harris hawks optimization (HHO). The modeling data were acquired for the period between 1962 and 2001 (480 months) and were grouped into several combinations and incorporated into the hybridized models. The performance of the models was assessed using the root mean square error (RMSE), mean absolute error (MAE), Nash–Sutcliffe Efﬁciency (NSE), coefﬁcient of determination (R 2 ), Taylor diagram, and Violin plot. The results showed that maximum temperature was the most inﬂuential variable for evaporation estimation compared to the other input variables. The effect of periodicity input was investigated, demonstrating the efﬁcacy of this variable in improving the models’ predictive accuracy. Among the models developed, the ANFIS-HHO and ANFIS-WOA models outperformed the other models, predicting Epan in the study stations with different combinations of input variables. Between these two models, ANFIS-WOA performed better than ANFIS-HHO. The results also proved the capability of the models when they were used for the prediction of Epan when given a study station using the data obtained for another station. Our study can provide insights into the development of predictive hybrid models when the analysis is conducted in data-scare regions.


Introduction
Out of several components of the hydrological cycle that affect regional water resources and agricultural production, evaporation is known as one of the essential components [1,2]. The importance of evaporation gains even more attention when the objective is the management of water resources in arid and semi-arid regions [3,4]. Evaporation is estimated using various indirect and direct methods, such as water balance, energy balance, mass transfer, Penman method, and evaporation pan [5]. Among them, the evaporation pan is a widely used and cost-effective method that reflects the amount of evaporation in the balance between the water and energy cycles, and measures the combined effect of several climate elements, such as solar radiation (SR), air temperature (TA), rainfall, relative humidity (RH), and wind speed (SW). Therefore, many researchers have attempted to estimate pan evaporation (Epan) using direct and indirect approaches [6,7]. Although a direct estimation of Epan seems to be more accurate, practical limitations and restrictions in instrumental devices encourage meteorologists to implement indirect approaches that work based on empirical, mathematical, and soft computing, as well as machine learning (ML) techniques [8].
In the case of the ANFIS method, many works have acknowledged its accuracy for Epan estimation [21,22]. ANFIS is a powerful tool for iterative optimization of complex real-world systems because ANFIS has the advantages of coupling both network-based models for the training process and fuzzy inference logic concepts for converting data patterns through membership functions. ANFIS also benefits from a fast calculation speed that allows for the finding of an optimal solution.
Apart from method selection, the selection of appropriate parameters is a critical step for the development of accurate and reliable Epan prediction models. Since each ML method, such as ANFIS, has several parameters, it is a difficult and time-consuming task to manually tune all parameters. Therefore, automatic parameter tuning using metaheuristic optimization algorithms has received attention from many researchers studying real-world problems. Examples of the metaheuristic optimization algorithms used for tuning a base method, such as ANN, LSSVM, SVR and ANFIS [23,24], include an electrostatic discharge algorithm [25], a water cycle optimization algorithm (WCA) [26], atom search optimization (ASO) [27], particle swarm optimization (PSO) [28], a cultural algorithm (CA) [29], a bee colony algorithm (BCA) [30], a genetic algorithm (GA) [31], biogeography-based optimization [24], and a firefly algorithm (FFA) [32]. Table 1 compares the hybridized ANFIS models with other ML methods for modeling evaporation as used in previous studies.
Owing to the capabilities of hybridized ANFIS models in capturing nonlinearity features of complex systems, this research aims to evaluate the performance of standalone ANFIS, along with three hybridized types of ANFIS, including ANFIS-PSO, ANFIS-WOA (Whale optimization algorithm), and ANFIS-HHO (Harris hawks optimization) in predicting Epan at three stations located in China. In this framework, limited data (Tmin and Tmax) and extraterrestrial radiation (Ra), which can be easily calculated with the Julian date, are used to develop the models. The major novelty of this study is to assess the potential of three hybridized models for estimating Epan. Furthermore, the contribution of the study relies on the evaluation of the models' ability in extrapolating Epan in one of the stations (Nanxian Station) in proportion to the other two stations (Jingzhou and Yueyang Stations).

Study Developed Model(s) Performance Comparison
Jasmin et al. [2] ANFIS and hybridized ANFIS with FFA, GA, and PSO The ANFIS-PSO model with R 2 = 0.99 and RMSE = 9.73 performed the best.

Case Study
The Dongting Lake basin (DLB) in southeastern China (20 • -32 • N and 108 • -123 • E) was selected as the case study area to illustrate the methodology proposed for Epan estimation ( Figure 1). Covering an area of around 261,400 km 2 , the basin is located in the south of the Yangtze River. It covers 12% of the Yangtze River Basin and is ranked as the world's second-largest freshwater lake. The DLB is drained by four main rivers, i.e., Li River, Yuan River, Zi River, and Xiang River. The DLB has a humid climate with many variations (mean annual =1380 mm) in precipitation that result in drought periods and flood events. The mean annual temperature of the basin is 17 • C, with the lowest temperature of 4.2 • C in January and the highest temperature of 28.9 • C in July. This region is mainly known for rice production, although production has been reduced due to drought events in recent years. The basin has an altitude gradient ranging from 30 m in plain areas to 2500 m in mountainous areas. Monthly minimum and maximum temperature data, as well as evaporation data, for the period between 1962 and 2001 (480 months) at the three selected stations were collected from the China Meteorological Administration (CMA). A summary of the statistical characteristics of the data is listed in Table 2. For the application of the machine learning models, data were divided into two sets, with 75% (360 months) of the data for model training and the remaining 25% (120 months) for model validation [13,28,33].

Adaptive Neuro-Fuzzy Inference System (ANFIS)
The adaptive neuro-fuzzy inference system (ANFIS) is an artificial intelligence method that was inspired by a combination of ANN and fuzzy logic. For an approximation function with three input variables and one output variable, a fuzzy rules base can be expressed as follows: where αi, βi, γi, and δi are the linear parameters in the consequent part of the ANFIS model; x 1 , x 2 , and x 3 are the input variables; and A i , B i , and C i are the fuzzy sets [36]. ANFIS consists of five layers which can be summarized as follows. In the first layer, all nodes are adaptive and their parameters (i.e., the premise parameters) are updated during the training process. For each node in the first layer, the linguistic label is calculated using the corresponding membership functions (MFs), as follows: For example, if the Gaussian membership function is used, it can be applied as follows: where σ and c are the parameters of the membership function (i.e., the premise parameters). For the second layer, all nodes are fixed (not adaptive), performing a simple multiplier; each one only calculates the firing strength (w i ) of the fuzzy rule.
The output of the third layer is calculated as the normalization of the firing strengths from the previous layer: For the fourth layer, the nodes are all adaptive and they perform the following output: where α i , β i , γ i , and δ i are the linear parameters (i.e., the consequent parameters) of the fuzzy rules. For the fifth layer, only one node is available and it provides the final response by calculating the overall output as a summation: ANFIS uses a hybrid training algorithm in two steps: (i) a forward pass, for which the least square method (LSM) is used for updating the linear parameters (i.e., consequent parameters) in the fourth layer; and (ii) the backward pass for which the gradient descent (GD) algorithm is used for nonlinear parameters (i.e., premise parameters).

Particle Swarm Optimization (PSO)
Particle swarm optimization (PSO) was proposed by Kennedy and Eberhart [37]. PSO is a metaheuristic algorithm inspired by a swarm and is generally based on simulation of the migration, natural behaviors, and aggregation of birds, insects, herds, and fishes during foraging. An important point to note is that each individual in the colony represents one particle. PSO is an effective tool for solving global optimization problems based on two major components: iterations and search [38]. The major goal of PSO is that, within the swarm, the best global and optimal solution should be sought among the undetermined number of particles using cooperation and information sharing. According to Kennedy and Eberhart [37], each particle in the swarm is defined by a pair of numbers containing two co-ordinates: flying speed, i.e., velocity (V i ), and position, i.e., location (X i ): We denote by P best the optimal position of each particle and G best the cluster (i.e., entire population) historical optimal position (i.e., location) as follows: During the PSO process, and in each individual flight, the particles change and update their positions and velocities as follows: where i is a particle in the group; N is the total number of particles in the swarm; d is the dimension of every particle; D is the number of variables; v k id and v k+1 id are the velocity of each individual particle i at the iteration k and k+1, respectively; x k id and x k+1 id are the positions of each individual particle i at the iteration k and k+1, respectively; θ 1 and θ 2 are random numbers ranging from 0 to 1; δ 1 and δ 2 are the cognitive and social acceleration coefficients; ω is an inertia weight parameter; p k id is the best position of the particle; and p k id is the swarm best position [39].

Whale Optimization Algorithm (WOA)
The whale optimization algorithm (WOA) is a swarm intelligence metaheuristic algorithm proposed by Mirjalili and Lewis [40], mimicking the physical movement and behavior of humpback whales in the ocean when swimming and attacking their prey. The overall WOA is formulated based on three optimization modes, namely: (i) encircling prey, (ii) spiral bubble-net feeding maneuver, and (iii) searching for prey. The algorithm process is divided into two important stages: exploration (i.e., spiral bubble-net feeding maneuver) and exploitation stages (i.e., searching for prey). To use WOA, similar to all metaheuristic algorithms, it is necessary to see the individual as an information carrier, for which each whale is considered an agent or a possible candidate solution. The agent, i.e., the individual whale, starts doing so using a solid commitment to look (i.e., keep seeking) for the best solution to the identified problem, i.e., the prey, using an iteration process [40].
Humpback whales possess the ability to detect even small prey and surround them (i.e., encircle). What is important to note is that, at the beginning of the process, the whale does not know the localization of the best solution, WOA hypothesizes that the identified target prey corresponds to the best candidate solution, or is relatively and approximately close to optimum. Consequently, once the best solution is decided, it becomes a haven for all remaining agents and will be forced to leave toward the best candidate (i.e., solution), and therefore, update their positions with respect to the best one. This first process is formulated as follows [40]: where → Z represents the distance separating the search agent to the position of the best candidate; t refers to the current iteration; where → µ is a variable with an absolute value, and it should be decreased from 2 to 0 during the iteration process; → θ is a random vector ranging from 0 to 1. Furthermore, to surround their prey, whales perform a spiral movement downwards as far as possible from the bottom to the top, simultaneously emitting an ensemble of bubbles in different sizes. This movement of the humpback whales is designated as a helix-shaped movement, and is expressed as follows: is the distance separating the prey and the whale; l is a variable ranging from −1 to +1; and b is a constant related to the shape of the logarithmic spiral [41]. The previous equation can be formulated as follows: where p is a random variable ranging from 0 to 1. Finally, during the exploration phase, a randomly searching strategy is adopted by the whale, in respect to the position of each whale compared to the other and not according to the best agent. Hence, the process can be formulated as follows: where → X rand is a random position vector selected from the current population (Mirjalili and Lewis, 2016; Wong et al., 2019) [40].

Harris Hawk Optimization (HHO)
The Harris hawk optimization algorithm (HHO) was developed and introduced by Heidari et al. [41]. Since it was proposed, it has become a famous and attractive metaheuristic optimization algorithm largely used for solving high complex nonlinear problems. The HHO can be classified into the category of population-based algorithms. The algorithm has two major phases: exploration and exploitation. During the exploration phase, we should highlight two important connected waypoints: (i) the candidate solution, which is the Harris hawk, and (ii) the best candidate solution, which corresponds to the targeted prey (i.e., the rabbit). During wild animal hunting (i.e., attack of the rabbit), the Harris hawk uses two strategies (Equation (27)) having the same chance q: (i) based on the position of other Hawk members of the population (q < 0.5), and (ii) the position of the prey (q ≥ 0.5) [41].
In the equation, we can see that there are four different positions, namely the actual position of the Harris Hawk (X(t)), the position of the Harris Hawk at the iteration t(X(t + 1)), the mean position of the Harris Hawk at the actual iteration (X m (t)), and the position of the prey (the rabbit), which is the X rabbit (t). More precisely, α 1 , α 2 , α 3 , and α 4 are random numbers ranging from 0 to 1; β and δ are the upper and lower bounds of variables; and X rand is a randomly selected Harris Hawk among all available individuals. The average position of the population is obtained using Equation (28): where X i (t) is the position of the ith Harris hawk and N is the total number of Harris hawks in the population. Before moving to the exploitation phase, we should consider a transition phase within which the energy of the prey is evaluated in a decreasing trend, which provides an opportunity for the Harris hawk to select a strategy among a large number of possible strategic choices depending on the energy loss of the prey (i.e., the rabbit), which steadily decreases as the number of iterations rises. Hence, the energy loss of the prey can be formulated as follows: where E is the escaping energy of the prey; E 0 is the prey energy before starting the flight forward; and T is the total number of iterations. From a mathematical point of view, the absolute value of E corresponds to two different cases: (i) if |E|≥1, we are in the exploration phase, and (ii) if |E|< 1, we are in the exploitation phase. We logically assume that the value of E 0 ranges in the interval from −1 to +1, and that its value should be updated at each iteration in decreasing or increasing ways. If the value of E 0 decreases from 0 to −1, we assume that the prey is weak, and if E 0 is higher than 0, we assume that the prey is strong [41]. During the exploitation phase, the prey (i.e., the rabbit) starts by drafting and adopting several escaping strategies in response to the different Harris hawk attack scenarios. Within the combinations of the escaping and hunting strategies, four possible scenarios can be formulated for the HHO. For this mathematical formulation, a numerical number (r) is used to demonstrate the success or failure of the escape plan adopted by the prey. Hence, r < 0.5 indicates the success of the escape, and conversely, r ≥ 0.5 corresponds to a failed escape plan. Based on these two values, a hard or soft besiege is performed, and the mathematical formulation of the two besieges is based on the values of (E). An absolute value of E higher than 0.5 (|E| ≥ 0.5) indicates that the HHO adopts a soft besiege, and an absolute value less than 0.5 (|E| < 0.5) leads to performing a hard besiege. The function of the soft besiege can be formulated as follows: The distance separating the location point of the prey and the location point of the Harris hawk corresponds to ∆X, and J is the random jump strength of the prey during flight; calculated as follows: where r 5 is a random number ranging from 0 to 1. The function of the hard besiege can be formulated as follows (r ≥ 0.5 and |E| < 0.5): The function of soft besiege with progressive rapid dives (r < 0.5 and |E| ≥ 0.5): If the result of the attack procedure is not positive, an update is made as follows: where LF is a levy flight; D and S are the number of dimensions of the problem and a random vector, respectively. LF can be expressed as follows: where β is a constant equal to 1.5; µ and υ are random values ranging from 0 to 1; and Γ is the gamma function. Finally, final formulation can be expressed as follows [42]:

ANFIS Optimization Using PSO, WOA, and HHO
To train ANFIS toward the aim of achieving a high level of efficiency and accuracy, two parameters should be properly adjusted. They are (i) linear parameters (i.e., consequent parameters) available in the fourth layer {α, β, γ, δ}, and (ii) nonlinear parameters (i.e., premise parameters) available in the second layer {c, σ} for the membership function. The LSM can be adopted to better identify the linear parameters, while the GD algorithm is used for nonlinear parameters. In this study, we used three metaheuristic algorithms for optimizing the premise and the consequent parameters. The number of parameters to be optimized can be calculated depending on the number of fuzzy rules and the number of MFs for each input variable (i.e., the linguistic label). The total number of premise parameters is equal to the total number of fuzzy labels multiplied by two (i.e., c and σ), while the total number of consequent parameters should be equal to the number of fuzzy rules multiplied by the number of input variables and one constant. Consequently, in the metaheuristic algorithms, the initial population (i.e., the total number of agents) is equal to the total number of parameters; thus, the initial particles (i.e., agents) are randomly obtained and the fitness function is then calculated, and the process will continue until the optimal condition with an acceptable error is achieved. In this study, the mean squared error (MSE) was used as the fitness function. Details of the HHA, WOA, and PSO algorithms are depicted in Figure 2 as the research flowchart. Water 2022, 14, x FOR PEER REVIEW 11 of 24 Figure 2. Flowchart of the methodology adopted in this study.  When the hybridized models were successfully developed, their performances were compared using the following metrics [43,44]: where Y c , Y o , Y o , and N are the calculated, measured, and mean of the measured Epan, and data quantity. The evaluation metrics were applied to data from three stations in China using several control parameters. The information on these parameter values is shown in Table 3. Population and iteration were set at 30 and 150 for 10 times to reach the robust outcomes.

Results and Discussion
The hybridized ANFIS models were compared for estimating the monthly Epan of the Jinzhou Station using various input combinations ( Table 4). The effect of periodicity was also investigated by adding the number of months over the year (α varies from 1 to 12) into the optimum input combination. From the first three input combinations, it is evident that the Tmax was the most influential parameter for Epan estimations for all models. For example, ANFIS with Tmax input had lower RMSE (0.6508 mm) and MAE (0.5210 mm), and higher R 2 (0.8956) and NSE (0.8922) compared to ANFIS, with only Tmin or Ra as input. Among the double input combinations, Tmin and Tmax had the highest accuracy (see combinations 4-6 in Table 4). Of the all input combinations, three inputs (Tmin, Tmax, and Ra) provided the lowest RMSE and MAE, and the highest R 2 and NSE in estimating monthly Epan for all four models. The last input combination comprised optimum inputs (Tmin, Tmax, and Ra) and α (Opt inputs, α). Table 4 reveals that the periodicity positively affected the accuracy of the ANFIS-HHO and ANFIS-WOA models, while involving α decreased the performance of the ANFIS and ANFIS-PSO models in the estimation of monthly Epan. Among the models developed here, the ANFIS-WOA provided the best accuracy, and the model with Opt inputs, α (Tmin, Tmax, Ra, and α), had the lowest RMSE (0.3127 mm) and MAE (0.2561 mm), and the highest R 2 (0.9726) and NSE (0.9704), followed by the ANFIS-HHO and ANFIS-PSO models. Table 4 shows that metaheuristic algorithms improved the performance of ANFIS in both the training and testing phases; the modeling error (RMSE) decreased by about 1.1%, 24.7%, and 29.1% in the testing phases of ANFIS-PSO, ANFIS-HHO, and ANFIS-WOA, respectively. The dominance of ANFIS-WOA can be clearly examined from the mean values of the statistics, which improved from 0.6661 to 0.5032, 0.5260 to 0.4030, 0.8622 to 0.9082, and 0.8592 to 9054 for the RMSE, MAE, R2, and NSE for the split scenarios from ANFIS to ANFIS-WOA, respectively. Table 5 shows the training and testing outcomes of the four hybridized ANFIS models in estimating the monthly Epan of the Nanxian Station. Similar to the Jinzhou Station, Tmax was identified as the most effective variable for Epan prediction in all models. Among the double input combinations, Tmin, Tmax, Tmax, and Ra generally produced the same level of accuracy and performed better than the input combinations of Tmin and Ra. Similar to the previous station, three input combinations with all three parameters (Tmin, Tmax, and Ra) had the lowest RMSE and MAE and the highest R 2 and NSE both in the training and testing phases of all models. Periodicity improved the estimation accuracy of the ANFIS-HHO and ANFIS-WOA models. Improvements in testing accuracy of ANFIS-WOA were 11.4%, 10.6%, 0.6, and 0.5% in terms of RMSE, MAE, R 2 , and NSE, respectively. It is clear from Table 5 that the metaheuristic algorithms improved the efficiency of ANFIS, and that the hybridized models performed much better in estimating monthly Epan using limited climatic input. The ANFIS-WOA model with three parameters and periodicity (Tmin, Tmax, Ra, and α) achieved the lowest RMSE (0.5886 mm) and MAE (0.4629 mm) and the highest R 2 (0.9526) and NSE (0.9507) in the testing phase. The WOA improved the accuracy of ANFIS with input Tmin, Tmax, and Ra by 23.9%, 20.7%, 4.5%, and 4.6% with respect to RMSE, MAE, R 2 , and NSE, respectively. The mean values of the comparison of RMSE, MAE, R2, and NSE statistics also endorsed the outperformed performance of ANFIS-WOA by improving these statistical values from 0.9029 to 0.7712, 0.6843 to 0.5893, 0.8484 to 0.9079, and 0.8461 to 0.9057 for the models ANFIS to ANFIS-WOA, respectively. Table 6 shows the modeling results for the Yueyang Station. Similar to the other two stations, Tmax was found to be the most influential variable on the performance of the models in estimating Epan. The combination of Tmin and Tmax had the best accuracy among the input combinations and produced the best accuracy in all models developed. Periodicity had a positive effect on the performance of all models. In this station, improvements in the accuracy of ANFIS were further achieved using the metaheuristic algorithms, with 11%, 17%, and 20.5% decreases in RMSE for the PSO, HHO, and WOA algorithms, respectively. The ANFIS-WOA model with three inputs and periodicity performed better than ANFIS-HHO, ANFIS-PSO, and ANFIS in estimating monthly Epan. The averages of the RMSE, MAE, R2, and NSE statistics also confirmed the dominancy of ANFIS-WOA over other ANFIS-based models in predicting monthly Epan. Figures 3-5 illustrate the scatterplots which compare the observed and predicted Epan by the best hybridized ANFIS models for the three stations. ANFIS-WOA had the least scattered predictions, followed by ANFIS-HHO, ANFIS-PSO, and ANFIS. Figure 6 displays the Taylor diagrams of the models of the Nanxian Station. This diagram is a useful approach in assessing three statistics (e.g., RMSE, correlation, and standard deviation) at the same time. The diagram revealed that the ANFIS-WOA model had the least RMSE and the greatest correlation and the nearest standard deviation to the observations, followed by the ANFIS-HHO models. Figure 7 shows the violin plots of the models for the Nanxian Station. The plots reveal that the most similar distribution to the observations belongs to the ANFIS-WOA model compared to other models. While known as a simple structure algorithm, WOA has fast convergence speed, high convergence accuracy, strong robustness, and stability [45]. WOA has the ability to balance exploration and exploitation to avoid local optima and reach a global optimal solution. Many studies have demonstrated that WOA not only minimizes the cost of solving engineering problems, but also provides better optimization efficiency compared to other population-based metaheuristic algorithms [45,46].   Nanxian Station. The plots reveal that the most similar distribution to the observations belongs to the ANFIS-WOA model compared to other models. While known as a simple structure algorithm, WOA has fast convergence speed, high convergence accuracy, strong robustness, and stability [45]. WOA has the ability to balance exploration and exploitation to avoid local optima and reach a global optimal solution. Many studies have demonstrated that WOA not only minimizes the cost of solving engineering problems, but also provides better optimization efficiency compared to other population-based metaheuristic algorithms [45,46].     . Taylor diagram of the predicted Epan using the hybridized ANFIS models in the testing phase using the best input combination-Nanxian Station. Figure 6. Taylor diagram of the predicted Epan using the hybridized ANFIS models in the testing phase using the best input combination-Nanxian Station. Figure 6. Taylor diagram of the predicted Epan using the hybridized ANFIS models in the tes phase using the best input combination-Nanxian Station.

Figure 7.
Violin plots of the predicted Epan using the hybridized ANFIS models in the testing ph using the best input combination-Nanxian Station. Overall, our results demonstrated that the hybridized models reduced the drawbacks of the standalone ANFIS and enabled a more accurate model Epan, in line with previous studies [14,21,23,31,32,34]. Table 7 shows the performance of the hybridized ANFIS models in estimating the monthly Epan of the Nanxian Station using climatic data of the Jinzhou Station. This type of comparison is essential, especially for the scare-data regions, and reveals that the models' outcomes with external input data are not as accurate as local input data. However, the models had acceptable accuracy in estimating Epan without local input data. For example, R 2 and NSE of the ANFIS-WOA model in the testing phase ranged from 0.8152 and 0.8126 to 0.9281 and 0.9286 for the worst and best input combinations. Again, the results demonstrated the efficiency of the metaheuristic algorithms for the improvement of the accuracy of ANFIS. The ANFIS-WOA performed the best in estimating monthly Epan using external input data. Table 8 compares the efficiency of the models in estimating the monthly Epan of the Nanxian Station using climatic data from the Yueyang Station. For this case, similar results were obtained; R 2 and NSE of the ANFIS-WOA model in the testing phase ranged from 0.8153 and 0.8125 to 0.9330 and 0.9283. The improvements in the accuracy of ANFIS when the metaheuristic algorithms were used were evident. A comparison between the results presented in Tables 7 and 8 reveals that the models performed better in estimating Epan at the Nanxian Station using the climatic input data of Jinzhou Station. The main reason for this might be related to the climatic characteristics of the stations. The Yueyang Station is located very near to Dongting Lake (Figure 1), and relative humidity may affect Epan. This study used the temperature input without considering this information.

Conclusions
This paper investigated the accuracy of four hybridized ANFIS models in modeling monthly Epan using limited climatic data from three stations in China. Various com-binations of minimum and maximum temperature and extraterrestrial radiation were considered as inputs to the models to investigate the effect of each variable on Epan. Among the input variables, Tmax was found to be the most effective variable on Epan in all stations and for all models. An examination of the effect of periodicity revealed that it could improve the models' accuracy in some cases. The hybridized models were also compared in estimating one station's evaporation using input data from other stations. It was found that the metaheuristic algorithms, specifically WOA and HHO, considerably improved the efficiency of ANFIS in modeling monthly Epan using limited climatic inputs in both applications. The ANFIS-WOA and ANFIS-HHO models with R 2 ≥ 0.81 and RMSE ≤ 1.04 mm performed the best using data from the all three stations, and successfully outperformed the other hybridized models in modeling Epan with local external temperature data. Our study can provide insights into the development of predictive models when analysis is based upon a narrow range of climatic variables.