Comparative analysis of data-driven methods for daily reference evapotranspiration estimation of Southern Caspian Sea

Reference evapotranspiration (ET o ) as a component in the hydrological cycle is calculated using many methods. In this study, the capability of four data-driven methods including artificial neural network (ANN), adaptive neuro-fuzzy inference system (ANFIS), support vector machine (SVM), and M5 tree model has been evaluated for daily ET o estimation at the south of the Caspian Sea. Different combinations of climatic data, solar radiation ( R s ), mean air temperature ( T mean ), mean relative humidity (RH mean ), and wind speed ( U ) during 1991 – 2020 were used as input variables. The data were divided into training and test data. The values generated from the methods were compared with those of the FAO-56 Penman-Monteith as a standard method. The results indicated that the accuracy of ANFIS was increased for estimating ET o , especially in validation phase, when all climatic variables were used as inputs in the synoptic stations. Totally, based on the evaluation of the performances, it was created that the ANFIS with T mean , RH mean , R s , and U variables had the best accuracy, while the ANN, SVM, and M5 with only one input of U had the worst performance. The ANFIS with T mean and R s was recommended for modelling ET o if there are fewer climatic variables in these regions.


| INTRODUCTION
Reference evapotranspiration (ET o ) is one of the main variables for the control of reference crop growth. Accuracy estimation of ET o is necessary for many studies such as water resources management, drought management (Jiang et al., 2022), irrigation system design, water budget estimation, and meteorological and climatological studies (Raza et al., 2020). This water cycle component is a complex non-linear phenomenon that has a great deal of dependence on climatological factors. Therefore, research in this field is very important. Evapotranspiration is achieved by direct and indirect methods. It can be measured directly by some instrument such as lysimeters. However, construction of this equipment is difficult and expensive (Ma et al., 2015). In an indirect method, usually empirical, semi-empirical, or physical equations have been used to estimate ET o based on meteorological or climatological data (Lopez-Urrea et al., 2006). The Food and Agriculture Organization (FAO) recommends the use of the standardized FAO Penman-Monteith (PMF-56) equation for ET o estimation (Allen et al., 1998). Considering the definition of the reference crop (Doorenbos, 1977), the accuracy of this method has been confirmed in different climates and time periods (Allen et al., 2005). Therefore, in most studies, this method is known as a standard method for comparing other empirical or soft computing methods . The Penman-Monteith equation needs many meteorological factors, such as air temperature, relative humidity, solar radiation, vapour pressure deficit, and wind speed, to calculate ET o . The application and evaluation of this method in some areas is a challenge due to lack of dense weather stations and incomplete meteorological data (Valipour, 2017). Nonetheless, a large number of researchers compared the standard Penman-Monteith with other empirical methods for estimating ET o (Celestin et al., 2020).
Over the past few years, many soft computing approaches and data-driven methods like artificial neural networks (ANN), support vector machines (SVM), adaptive networks based on fuzzy inference system (ANFIS) and M5 model have been used in the field of various environmental sciences, such as water engineering, hydrological process, and meteorological and climatological studies (Pal & Deswal, 2009;Sekertekin et al., 2021).
Since the 1990s, ANN, based on the realization of the nervous systems, has been used extensively in hydrological prediction and simulation (Yang et al., 1996), optimization (Wen & Lee, 1998), water resources, and meteorological applications (Narvekar & Fargose, 2015). ASCE Task Committee went over the application of this method in hydrology problems (ASCE 2000a(ASCE , 2000b. Kisi (2008) investigated the accuracy of three different ANN methods, including the multi-layer perceptron (MLPs), radial basis neural networks (RBNNs), and generalized regression neural networks (GRNNs), in modelling ET o . The results demonstrated that the MLP and RBNN models could be used successfully in ET o estimation. In another study, Chauhan and Shrivastava (2009) evaluated the performance of four climate-based methods and ANN for ET o estimation in India. The results demonstrated that the ANN models performed better than the others. For the past 2 years, the use of ANN in ET o estimation has still been used by researchers (Majhi & Naidu, 2021;Walls et al., 2020).
The ANFIS method was first developed by Jang (1993) and has been used in many respects, for example, hydrology and water resources problems (Alquraish et al., 2021). Pulido-Calvo and Gutie'rrez-Estrada (2009) discussed the potential of ANFIS, ANN, and multivariate regression in ET o estimation in southern Spain. They found that the ANFIS performed better than the other methods. Ladlani et al. (2014) modelled the daily ET o in Mediterranean region of Algeria using ANFIS and multiple linear regression (MLR). They drew the conclusion that ANFIS simulated ET o better than the MLR. Keshtegar et al. (2017) studied the accuracy of subset ANFIS in estimating the daily ET o and compared it with M5, ANN, and ANFIS. They found that the accuracy of studied models was increased using subset ANFIS. The use of ANFIS is now recommended in modelling ET o in different climates (Roy et al., 2021).
The SVM is based on the concept of decision planes that define decision boundaries. This technique has a good distribution capability and has been proven to be a strong algorithm in many fields (Vapnik, 1998). Lately, the use of SVM has been considered in hydrology engineering, such as lake water level prediction (Wu et al., 2008) and stream flow simulations (Lin et al., 2006). Eslamian et al. (2009) investigated the ability of the SVM and ANN methods to model the ET o in a controlled greenhouse environment. They concluded that the SVM model was superior to ANN in modelling evapotranspiration. Rahimikhoob (2014) evaluated the performance of ANN and M5 methods with standardized PMF-56 method in an arid climate. The results showed that the ANN computed ET o better than the M5, but both methods performed well for this climate. The M5 tree model that was developed by Quinlan (1992) is a logical learning technique that deals with continuous numerical values, and its learning algorithm is called M5. This algorithm is the most popular categorization handled in the family of decision-making tree models. This model has been applied to study many problems in hydrology, such as flood forecasting (Solomatine & Xue, 2004) and meteorology (Kisi et al., 2021).
However, simultaneous use and comparison of these four data-driven methods for daily ET o estimation in different climates are limited in the literature (Rahimikhoob, 2014). Since ET o is a non-linear and stochastic component of the hydrologic cycle and has a great deal of dependence on meteorological factors; therefore, the study and comparison of many data-driven methods are necessary for estimating this component. From the literature review, it was considered that data-driven methods have been applied to modelling ET o in different regions and climates; however, the application of these methods has been studied finitely (Murat, 2011;Shiri et al., 2013). Because of the dynamic nature or even non-stationary of evapotranspiration, it is difficult to estimate ET o data accurately, especially in areas with limited meteorological data.
Therefore, the aim of this study is to peruse the usability of four data-driven methods (ANN, ANFIS, SVM, and M5 tree model) with different characteristics for estimating ET o in daily time step using different climatic input combinations of solar radiation (R s ), mean air temperature (T mean ), mean relative humidity (RH mean ), and wind speed at 2-m height (U) at the southern climates of the Caspian Sea. The gathered data of the three synoptic stations (from 1991 to 2020) were used for training and validation of the applied methods. The performances of these methods were compared with the standardized PMF-56 method. Some statistical criteria have been used for the model's performance evaluation. The best one with a smaller error and highest coefficient of determination was selected for future estimation of daily ET o . This inspection will help to find a correlation between meteorological variables and daily ET o and to recognize the most substantial climatic factors that effect on daily ET o .

| Study area and data
The research was carried out in the coastal zone of the southern part of the Caspian Sea (North of Iran). The Caspian Sea is the largest interior body of water on Earth, which has an area of 360,000 km 2 . Similar to many lakes, it does not feed into an ocean, but it looks like a sea in terms of its depth and size. That is why the lake is mentioned as 'Sea' in the literature. The Caspian Sea has about a third of the mean salinity of seawater, thus it is not actually a freshwater lake (Ghadiri et al., 2006). Because this sea has no outlets, the main factor in the loss of water in this area is the evaporation phenomenon. The coastline in the south of the Caspian Sea, expanding over 800 km, contains nearly 18% of the whole Caspian coastline. The Caspian Sea climate is mostly affected by its geographical status and the prevailing regional atmospheric circulation. Most of the precipitation occurs during December to March, and most of evapotranspiration occurs during June to September; therefore, there is a serious water scarcity during the growing season of crops due to the evapotranspiration phenomenon in this region. Evapotranspiration varies depending on the climatic conditions, the surface water availability, and the crop pattern in this area. Thus, accurate estimation of ET o is an important component in irrigation management in this coastal area. For this purpose, three synoptic weather stations, namely Gorgan, Babolsar, and Rasht, were selected in this study. These stations are located in the range of 36 54 0 to 37 19 0 N (latitude), from 49 37 0 to 54 24 0 E (longitude), and from À8.6 to 13.3 m above sea level elevation, covering three climate zones based on De Martonne Aridity index (Rahimi et al., 2013). These stations are on or near the coast of the Caspian Sea and are supposed to illustrate the general state of the southern Caspian Sea. The distances between Gorgan to Babolsar and Babolsar to Rasht synoptic weather stations are about 184 and 313 km, respectively. Table 1 shows the detailed geographical and climatological characteristics of these synoptic stations. The locations of these stations are shown in Figure 1. The long-term available daily meteorological records of solar radiation (R s ), mean air temperature (T mean ), mean relative humidity (RH mean ), and wind speed (U) were collected from IRIMO (Islamic Republic of Iran Meteorological Organization, 2021) during 1991-2020 were applied to calculate ET o . The quality control tests of raw data were accomplished according to Allen (1996). The statistical analysis was carried out for the data before using for the estimation of ET o . Moreover, the main parameters of the statistical analysis are illustrated in Table 1. The standard deviation value for RH mean is higher as compared with the other variables at the three stations. The daily ET o data for this study were divided into two parts: training and validation data sets. The data sets from January 1991 to December 2013 were used for training (80% of the data), meanwhile those from 1 January 2014 to December 2020 (20% of the data) were used for validation of the applied methods.

| Explanation of chosen methods
A physical base method (PMF-56) and four data-driven methods (ANN, ANFIS, SVM, and M5 tree model) were used in this study to estimate ET o . So, a summary of these methods is provided below.

| FAO-56 Penman-Monteith method
As there is no lysimeter for measuring ET o in the studied stations, daily ET o data were calculated by the FAO-56 Penman-Monteith (PMF-56) method. This method has been suggested as the standard equation for defining ET o by FAO (Doorenbos & Pruitt, 1977) and was applied as the benchmark model. According to Smith et al. (1992), this physical method estimates the ET o closer to the lysimetric values than other empirical methods. The PMF-56 can be used for daily or shorter time steps. For daily scales, Equation (1) is proposed:  where ET o is the reference evapotranspiration (mm day À1 ), R n is the net radiation at the crop surface (MJ m À2 day À1 ), G is the soil heat flux density (MJ m À2 day À1 ), T mean is the mean daily air temperature at a 2-m height ( C), u 2 is the wind speed at a 2-m height (m s À1 ), e s is the saturation vapour pressure (kPa), e a is the actual vapour pressure (kPa), e s À e a is the saturation vapour pressure deficit (kPa), Δ is the slope of the vapour pressure curve (kPa C À1 ), and γ is the psychrometric constant (kPa C À1 ). Several kinds of methods are suggested by Allen et al. (2005) to evaluate the parameters of Equation 1 from the observed climatic data. According to Allen et al. (1998), for daily estimates of the evapotranspiration rate of reference crop, the term G can be neglected, that is, we set G = 0. The net radiation at the surface (R n ) is given as the difference between incoming net short wave radiation and outgoing net long wave radiation. The computation of all data required for the calculation of the net radiation at the surface followed the procedure given in Chapter 3 of the FAO paper 56 (Allen et al., 1998).

| Artificial neural network
The ANN is distinguished as a set of similar structural components called neurons resembling the humans' neural system. In general, it is a biologically inspired network of artificial neurons put together to accomplish particular roles. An ANN has three fundamental elements including classes of connection weights, a threshold, and an activation function. In other words, the most popular employed ANN is the three-layer feed-forward back propagation (BP) training, which generally has, different layers such as input layer (where the data are offered to the network), hidden layer (with the number of neurons where the data are processed) and output layer (where the results are created), and some nodes in each layer. The structure of ANN depends on the nature of the problem and is distinguished by user. There are many ways to train the neural network, and the error BP algorithm was used in this research.

| Adaptive neuro-fuzzy inference system
A combination of ANNs and fuzzy inference systems (FISs), which was introduced by Jang (1993), creates the ANFIS system. This intelligence system uses a non-linear process from input to output, and it is a multi-layer form connected with a number of nodes that has three sections (Lȩski & Czogala, 1999). In the first part, a rule base including fuzzy rules is selected. The second part contains the database. In this case, a learning algorithm with the compatible membership functions (MFs) is used to create a set of fuzzy if-then laws from certain inputoutput subsets. For this purpose, the ANFIS receives primary FIS and modifies it with a propagation algorithm. Finally, in the last section, a logical inference system is used to remove it from the rules and input data to reach an acceptable output. The most famous of them, which are used in fuzzy systems, are Mamdani, Sugeno, and Takagi models. These systems recognize patterns and help to the adaptation of natural world. Figure 2 illustrates an ANFIS structure with two inputs (x 1 , x 2 ) and F I G U R E 1 Location of the three synoptic stations in the southern parts of the Caspian Sea.
one output (F) in which a circle represents a stable node, and a square shows an adaptive node. A description of the ANFIS layers is available in Lȩski and Czogala (1999).

| Support vector machine
A SVM theory is based on the classification technique, which is applied a linear model to divide the main data into a feature space through a non-linear mapping function (Vapnik, 1998). This technique maximizes the boundary between the two opposing classes and minimizes the risk of a structure. SVM can also be used in regression relationships. The regression method used in SVM is different from the usual regression methods (Mustafa et al., 2018). In general, in order to reduce the number of errors in the training phase, this model uses non-linear machines to find the most suitable hyper plane. In the next step, after constructing a Lagrange function and optimizing the target function, a non-linear regression function is obtained. The functions used for optimization are presented in the scientific resources (Tezel & Buyukyildiz, 2016). The non-linear part of the SVM is the recognition of Kernel functions. The simplest form of this function is the linear kernel (2005). Other types of kernels include the Polynomial, Gaussian radial basis function (RBF), and sigmoid. If the points could not be separated by a kernel, it would have to do so until it was properly segmented with other kernels. In this study, the kernel function of the RBF type was selected.

| M5 tree model
The tree model is a data-driven technique based on machine learning, presented by Quinlan (1992). The algorithm for learning this model is known as the M5. Decision trees are used to predict outputs from the input data. In these trees, each node (leaf) represents a class or internal decision that performs certain specified tests. The advantage of this model is that it is a white-box method in which the dependence between the variables is determined by using linear and non-linear functions. The M5 tree is the generalized regression tree and uses the linear regression functions in the leaves. In order to implement the model tree generation, three stages of splitting, pruning, and smoothing should be considered. Initially, an inductive algorithm is used to construct the decision tree. For this purpose, a splitting criterion is used. This criterion reduces the internal variation in the model. The criterion is based on the standard deviation (SD) behaviour of each class's values (Pal & Deswal, 2009). When some samples are left very small, the splitting process stops. As a matter of fact, in this stage, the M5 model calculates the error values based on the minimum error rate at each node. In the final step, the accuracy of the estimation of the tree model is improved using the smoothing of the sharp nodes. Further information about the M5 tree model can be found in the study by Quinlan (1992).

| Models' performance evaluation
Various techniques are recommended to assess the performance of models in the calibration and validation phases, which is generally considered as goodness of fit (Shcherbakov et al., 2013). Moreover, the performances of the data-driven methods were compared with the PMF-56 method. Five performance evaluation criteria were used in this study. The coefficient of determination (R 2 ) has been widely applied for model evaluation. It is used to describe how differences in a simulated variable can be explained by a difference in the measured variable. The range of this coefficient is between 0 and 1. If R 2 = 0, There is no linear relationship between observed and calculated values. If R 2 = 1, a very significant regression relationship exists. The root mean squared error (RMSE) evaluates the residual between observed (by PMF-56) and forecasted (using data-driven methods) ET o . The model and observational data are better suited if this statistic is smaller. The mean absolute percentage error (MAPE) has been widely used for relative error comparison in the estimation in respect of the actual rate of the variable. Consequently, The MAPE is an unbiased statistic for quantifying the predictive ability of the model. Smaller values of this index imply higher model efficiency (Lefebvre & Bensalma, 2015). The original index of agreement (IA) is a standardized measure of the degree of model estimation error and varies between F I G U R E 2 Adaptive neuro-fuzzy inference system architecture with two inputs and one output. 0 and 1 (Willmott, 1982). Perfect agreement would exist between observed and predicted data if IA = 1. This index is heavily sensitive to extreme values due to the use of squared error. Finally, the difference between the models was measured based on the Nash-Sutcliffe efficiency (NSE) coefficient (Nash & Sutcliffe, 1970). This coefficient can represent the accuracy of any modelling approach as very good (0.75 < NSE ≤ 1.0), good (0.65 < NSE ≤ 0.75), satisfactory (0.50 < NSE ≤ 0.65), acceptable (0.40 < NSE ≤ 0.50), and unsatisfactory (NSE ≤ 0.40). Altogether, the performance indices listed above are used to analyse the accuracy of the applied models in daily ET o estimation in this study. R 2 and IA with larger amounts, RMSE and MAPE with smaller amounts, and NSE point to higher model efficiency.

| Different modelling structures
In data-driven methods, the selection of input environmental variables is based on the recognition of the system physics (ASCE, 2000a(ASCE, , 2000b. Some variables have high linear correlation with ET o . Therefore, in this study, daily mean air temperature (T mean ), mean relative humidity (RH mean ), wind speed (U), and solar radiation (R s ) were used as inputs for estimation of ET o as output in all the proposed data-driven methods. From 10,950 daily data, 80% was used for training and 20% for validation. Because the four input variables and the output variable (ET o ) had different range, it was noticed to be essential to standardize the primary data to increase the training speed and the accuracy of the network. Therefore, the daily climatic inputs and output variables were normalized between 0 and 1 to improve the performance of the models (Equation 2) (Kumar et al., 2002).
where y is the normalized value of each variable, y o is real value of each variable, and y min , y max are the minimum and maximum rates of the actual values, respectively.
In this research, The ANN was implemented in the Neural Network Toolbox in MATLAB software (version 9.4.0.813654, R2018a). Similarly, the ANFIS was implemented by using FIS in this software during the whole process of the training and validation phases. Moreover, MATLAB software was used to write computational codes and develop learning and simulating algorithm for the SVM and M5 tree model.
Four different combinations of input variables were used to estimate the daily ET o using applied methods (ANN, ANFIS, SVM, and M5) at the three stations. The different combinations of input variables (as models No. 1 to 4) are presented in Table 2.
In this study, for the ANN method, after studying and testing various structures, the best one was obtained as 4-10-10-1. The most common network transfer function, which is known as the tan-sigmoid function, was selected as the network activation function. Moreover, for the calculation of the output, the feed-forward with learning algorithm back propagation error (BP) type with two hidden layers was used. In this algorithm, the output computation was done from the input to the output in the forward direction. The BP algorithm using the 10 nodes during the training phase gave the least mean square error and was used for network training.
Similarly, an optimal architecture was chosen after analysing different structures in the ANFIS method. In this case, the MF parameters were iteratively learned using a hybrid method. In accordance with a number of tests performed, the Gaussian type with two MFs was selected. A Sugeno's type FIS was selected for functioning, and the forward pass of hybrid learning algorithms was applied to identify the parameters. The resulting parameters are determined by the least squares estimate.
For the SVM method, the RBF kernel was used. The general performance of this method depends on an appropriate setting of the regulation parameters of kernel function. Therefore, two parameters including the penalty parameter (C = 10) and the regulation parameter (ϵ = 0.1) were obtained for the study stations.

| RESULTS AND DISCUSSION
Before presenting the main results of this paper, it should be noted that the sensitivity analysis of the estimation of daily ET o in the study area was previously performed by Bakhtiari and Baghizadeh (2012). Sensitivity analysis indicated that R s , RH mean , and T mean were effective parameters in the study area. However, the sensitivity of the variable U was less than other variables and was less important.
The results of the four different data-driven methods for estimating daily ET o are presented in Tables 3-5 in Bablosar Station were 0.10 and 6.02 mm day À1 , respectively, and the values related to the Rasht station were 0.05 and 6.85 mm day À1 , respectively. Since models No. 2-4 showed better performance at the three stations, the figures presented in Tables 3-5 were based on the results of these models. The values of the performance's indices for the cases of the best inputs and the best model were shown in bold. These tables present the statistical values of both training and validation phase, concerning various performance criteria. It can be observed that all the different data-driven methods have good performance during both training and validation phases.
For Gorgan synoptic station, during the training phase, the SVM method with four input factors (T mean , RH mean , R s, and U) obtained the best R 2 , RMSE, IA, and NSE performance criteria of 0.982, 0.312, 0.990, and 0.956 respectively, while the ANFIS method obtained the best MAPE of 12.284. As can be seen from Table 3, during the validation, the ANFIS model-4 performs the best on the whole, followed by the ANFIS model-3 (with three input factors including T mean , RH mean , and R s ) and ANFIS model-2 (with three input factors including T mean , R s , and U) models. ANFIS model-4 obtained the best values of RMSE = 0.381, MAPE = 13.313, and NSE = 0.946 (very good) during the validation phase. M5 model-3 is able to obtain the maximum ET o among all models created in this station.
Similarly, for Babolsar synoptic station, in the training phase, the SVM model-4 obtained the best RMSE, IA, and NSE statistics of 0.281, 0.997, and 0.979, respectively, while ANFIS model-4 got the best R 2 and MAPE statistics of 0.987 and 10.624, respectively. During the validation phase, the results indicated that both SVM and ANFIS models-4 obtained the best R 2 , IA, and NSE statistics of 0.982 and 0.995, and 0.97 (very good), respectively, whereas ANFIS model-4 obtained the best RMSE, and MAPE statistics of 0.202 and 9.826, respectively. For Rasht synoptic station, in the training phase, SVM model-4 obtained the best R 2 , MAPE, and NSE of 0.998, 14.879, and 0.936, respectively, while ANFIS model-4 obtained the best RMSE and IA statistics of 0.386 and 0.988, respectively. According to the figures in Table 3, we can conclude that ANFIS model-4 outperforms all other applied models in the validation phase.
Moreover, in the validation phase, as identified in Tables 3-5, the results of the ANN, ANFIS, and SVM method estimation were able to recognize a good, near estimate, as compared with those with M5 tree method; however, it can be achieved that the ANFIS method obtained the best RMSE between the observed and estimated ET o in all the study stations. In addition, it can be concluded that the ANFIS method obtained the best minimum absolute error between the observed and modelled minimum and maximum ET o in Babolsar and Rasht synoptic stations, and the SVM and ANFIS methods obtained the best minimum absolute error between the observed and modelled minimum and maximum ET o , respectively, in Gorgan weather station. From the other viewpoint, it can be achieved that the range of ET o (difference between the maximum and minimum values = 8.67 mm day À1 ) calculated from the ANN model-4 was near those of the observed values in Gorgan station (Table 3). Similarly, the nearest range of ET o to the observed values were 6.05 (in ANFIS model-4) and 6.24 (in ANFIS model-2) mm day À1 for Babolsar and Rasht stations, respectively (Tables 4 and 5).
According to the values in Tables 3-5, we can conclude that ANFIS model-4 had the best performance in the entire study region compared with the used models. The comparison between the observed and estimated ET o in the validation phase shows that the ANFIS method performs superior to the other data-driven methods in estimating ET o in all the study stations. The best scores were obtained using ANFIS model-4 in the validation phase for estimating Babolsar's ET o . The results indicate that the ANN model underestimated the daily ET o of PMF-56 in the Rasht station, especially in values greater than 3 mm day À1 (Figure 3).
Overall, the results of the data-driven methods indicated that during the training phase, the ANFIS or SVM method is suitable to get the best output, and during the validation phase, the ANFIS method can obtain the best output in terms of different performance statistics. Besides, these statistics illustrated that the degree of goodness or defect of estimation accuracy varies according to the different performance criteria during both the training and validation phases. ANFIS method with different input factors (model no. 1-4) can obtain a better output in terms of different evaluation criteria not only during the training but also during the validation phase. The statistics measures indicated that the results of the SVM method during the validation phase are lowerranking than the results of the training phase. ANN is in the middle level in both the training and validation phases compared with the other three methods, and is able to obtain the better RMSE and MAPE in validation phase in all the stations, especially ANN model-3 can obtain the maximum ET o among all models developed in Rasht synoptic station. Figure 4a-d shows the results of four method's regression error on the basis of RMSE during the validation phase for each station, which illustrates the effect of each meteorological variables on ET o estimation. It is clear that model-4 with input structure of T mean , RH mean , R s , and U has the lowest RMSE for all the study stations; therefore, they are the most relevant meteorological variables to predict daily ET o . Clearly, there is no significant difference between RMSE values for the first and even the second (with input variables including T mean , RH mean , R s ) most significant parameters. Moreover, it is found that U is a less relevant input variable since it has the highest RMSE. It is found that for all the stations, RH mean is the second unimportant parameter. This emphasizes the importance of selecting the appropriate input variables for predicting ET o at the study stations. It should be considered that a model with more simplicity is likely to be noted in terms of required input variables; therefore, considering the most relevant two input variables would be acceptable for many researchers in terms of optimum number of inputs to prepare a balance between the simplicity and high accuracy. In this case, the model two inputs of T mean and R s can be used in the study areas. It should be mentioned again that smaller rates of RMSE illustrate further exactness of the predicted ET o , and in a perfect case, they are zero. Thus, according to the results, it is found that the ANFIS method provides outperforming results to predict daily ET o than the ANN, SVM, and M5 tree methods.

| CONCLUSIONS
The permanent need to increase agricultural yield, together with the several more drought occurrences, requires a more careful evaluation of near-term irrigation demands and, accordingly, a more accurate estimation of evapotranspiration. Generally, climatic data in many regions is limited to a reduced number of variables measured; therefore, the current or future evapotranspiration estimation is restricted. Thus, in practical applications, data-driven methods may be a powerful tool for the estimation of precise evapotranspiration when a time series of a few years is accessible. The first question is what approach should be appropriately chosen from a variety of machine learning methods in practical applications. This is very important in usage, mainly in developing countries, where some weather data may be removed because of technical reasons.
In this study, four soft computing techniques were carried out based upon the data-driven methodology to identify the most relevant climatic variables for estimation of daily reference evapotranspiration (ET o ). The high number of climatic variables for estimating ET o has many disadvantageous. At first, data collection may be costly in some areas. Secondly, the various input climatic variables may reduce the capability of generalizability of a model. It is therefore helpful to apply methods, which allow simplifying the number of input variables. Therefore, the performances of several data-driven methods were investigated for modelling daily ET o time series in order to select the most important climate factors on ET o estimation in the Mediterranean and humid climates. The data-driven methods used in this study were the ANN, ANFIS, SVM, and M5 tree. The conventional standardized PMF-56 method, known for its accuracy, is also used as a benchmarking yardstick for comparison goals. The daily ET o data from actual weather observed data in the tree synoptic stations in the southern Caspian Sea were used to investigate various methods. Prior to data-driven modelling, the data set was divided into two subsets of training and validation. Four different combinations of climatic variables were delineated as models-1 to 4. Four standard performance evaluation criteria were used to assess the performances of each model. Based on the results obtained, it is noticeable that the most influential single variable for estimating daily ET o in Babolsar (with the humid climate) and Rasht (with the perhumid climate) is solar radiation, while for the Gorgan station (with the Mediterranean climate), the daily mean air temperature is the most effective factor for ET o estimation.
Furthermore, if the number of climatic variables was increased, the combination of the daily mean air temperature and solar radiation would indicate the optimal structure of the variables at all three stations. Finally, the combination of four variables, daily mean air temperature, mean relative humidity, solar radiation, and wind speed are the optimal combination for ET o estimation. The study indicated that the ANFIS method is able to compute better estimating accuracy in terms of different evaluation statistics during both the training and validation phases. ANFIS model-4 performs better than the others in all the study stations. The results of the SVM method during the validation phase were inferior to the results during the training phase. All-inclusive, the results indicated that the data-driven methods are strong tools for modelling daily evapotranspiration time series, and this may give valuable reference for researchers who apply data-driven methods for modelling ET o . The results of this study can provide a solution for accurate ET o estimation in the absence of complete climatological data. In further study, the artificial intelligence methods may be applied in estimating hourly evapotranspiration and comparing it with other empirical equations in different climates.
AUTHOR CONTRIBUTIONS BB suggested the topic and was the main supervisor, and helped in writing and editing the manuscript. AM-D acquired and analyzed the data, performed the computations, and prepared the initial draft. KQ supervised data driven methods analysis and edited the final manuscript. All authors provided critical feedback and helped shape the research, analysis and manuscript.