Modeling monthly pan evaporation process over the Indian central Himalayas: application of multiple learning artificial intelligence model

The potential of several predictive models including multiple model-artificial neural network (MM-ANN), multivariate adaptive regression spline (MARS), support vector machine (SVM), multi-gene genetic programming (MGGP), and ‘M5Tree’ were assessed to simulate the pan evaporation in monthly scale (EP m ) at two stations (e.g. Ranichauri and Pantnagar) in India. Monthly climatological information were used for simulating the pan evaporation. The utmost effective input-variables for the MM-ANN, MGGP, MARS, SVM, and M5Tree were determined using the Gamma test (GT). The predictive models were compared to each other using several statistical criteria (e.g. mean absolute percentage error (MAPE), Willmott’s Index of agreement (WI), root mean squared error (RMSE), Nash-Sutcliffe efficiency (NSE), and Legate and McCabe’s Index (LM)) and visual inspection. The results showed that the MM-ANN-1 and MGGP-1 models (NSE, WI, LM, RMSE, MAPE are 0.954, 0.988, 0.801, 0.536 mm/month, 9.988% at Pantnagar station, and 0.911, 0.975, 0.724, and 0.364 mm/month, 12.297% at Ranichauri station, respectively) with input variables equal to six were more successful than the other techniques during testing period to simulate the monthly pan evaporation at both Ranichauri and Pantnagar stations. Thus, the results of proposed MM-ANN-1 and MGGP-1 models will help to the local stakeholders in terms of water resources management.


Introduction
Evaporation plays the main role in environmental studies and water resources management. Therefore, the precise simulation of EP m is necessary for hydrological simulating, river flow forecasting, forestry, lake-ecosystems, agronomy and irrigation sciences (Burt, Mutziger, Allen, & Howell, 2005). The most effective meteorological factors on evaporation rate are relative humidity, air temperature, vapor pressure deficit, atmospheric pressure, vapor pressure, wind speed, and sunshine hours . Evaporation is generally measured using two methods: (i) direct methods such as Class A panevaporimeter and (ii) indirect methods include empirical equations (Ghorbani, Deo, Yaseen, Kashani, & Mohammadi, 2017). Doorenbos and Pruitt (1977) said that the Class A pan-evaporimeter performance was affected by CONTACT Zaher Mundher Yaseen yaseen@tdtu.edu.vn the turbidity of the water, watering of birds or other animals, human and instrumentation errors, and maintenance problems. Numerous studies reported the determination of the pan evaporation applying some empirical and semi-empirical formulations based on various meteorological information (Griffiths, 1966;Penman, 1948;Priestley & Taylor, 1972). However, there are limitations to the applications of these methods because of using different climatic factors. Therefore, alternative methods, which require less meteorological data are needed to predict the evaporation (Kisi, 2015;Kisi, Genc, Dinc, & Zounemat-Kermani, 2016). AI models has shown promising progress on the evaporation estimation (Tao et al., 2018). Several versions of AI models are explored for pan evaporation simulation including; ANN based models, Fuzzy based models, support vector machine (SVM), evolutionary computing, data mining and complementary wavelet-AI models (Adnan, Malik, Kumar, Parmar, & Kisi, 2019;Guven & Kisi, 2013;Kisi & Heddam, 2019;Qasem et al., 2019;Qutbudin et al., 2019;Rezaie-Balf, Kisi, & Chua, 2019;Sebbar, Heddam, & Djemili, 2019;Yaseen, El-shafie, Jaafar, Afan, & Sayl, 2015). Tabari, Talaee, and Abghari (2012) predicted daily pan-evaporation using CANFIS and MLPNN techniques in Iran. The research findings demonstrated that the MLPNN model provided a better accuracy than the CANFIS model. Kisi (2015) used least square SVM (LSSVM), M5Tree, and MARS for estimating pan evaporation in Turkey. They found that the MARS provided a significant superiority compared to other models. Deo, Samui, and Kim (2016) developed the ELM, MARS, and relevance vector machine (RVM) to predict evaporative using maximum and minimum temperatures, atmospheric vapor pressure, precipitation and solar radiation in Australia. The results indicated that the RVM model has high ability compared to the other techniques. Malik and Kumar (2015) hired the CANFIS, MLR, and ANN models for predicting the pan evaporation over different stations in India and the ANN models are more successful compared to other ones. Keshtegar, Piri, and Kisi (2016) used M5Tree, conjugate gradient (CG), and ANFIS models to simulate of evaporation. Results produced that the CG model is more satisfactory than the M5Tree and ANFIS techniques.  estimated pan evaporation using the FG, ANFIS-GP, and M5Tree models in China. Results proved the high ability of the FG model in pan evaporation estimation. Wang, Kisi, Hu, et al. (2017) predicted pan evaporation using MLPNN, GRNN, ANFIS with grid partition (ANFIS-GP), MARS, MLR, FG, LSSVM, and SS models in China. The results showed that the AI models provided the most accuracy than the SS and MLR techniques. Wang, Niu, Kisi, Li, and Yu (2017) Remesan, Shamim, & Han, 2008;Singh, Malik, Kumar, & Kisi, 2018). Moghaddamnia, Ghafari Gousheh, Piri, Amin, and Han (2009) simulated evaporation using the ANN and ANFIS models while applying GT to select suitable input variables in different area of Iran. Results showed the ANFIS and ANN models have high capability in the estimation of EP m values. Goyal, Bharti, Quilty, Adamowski, and Pandey (2014) identified suitable input variables using GT for FL, ANN, ANFIS, and least square-SVM techniques to estimate EP m in daily scale in India, and the SS and Hargreaves-Samani (HGS) equations were compared with the results obtained by AI models. The study showed the superior accuracy of the SCT over SS and HGS techniques. Malik, Kumar, and Kisi (2018) simulated daily EP m using RBNN, MLR, Griffiths, Stephens-Stewart, Priestley-Taylor, Christiansen, Penman, SOMNN, and Jensen-Burman-Allen models in India. The input combinations of MLR, RBNN, and SOMNN were determined using the GT. The results of the research indicated that the most accuracy was observed in the RBNN model.
The current research is conducted to accomplish the following aims: (i) An appropriate input combination is identified using GT before the construction of the proposed and comparable predictive models. (ii) Potential of the proposed model is investigated using climate information of Ranichauri and Pantnagar stations located in Indian central Himalayas. (iii) A comprehensive evaluation is performed on the obtained performance of the proposed predictive models.

Description of case study
The sources of the weather information used in this study are two stations (i.e. Ranichauri and Pantnagar) in India. The coordinates of the Pantnagar and Ranichauri stations are (79°38 0 E, 29°0 0 N) and (78°24 35 E, 30°18 40 N), respectively. These coordinates are depicted in Figure 1. The Pantnagar and Ranichauri stations are located at altitude of 243.8 and 2000 m above sea level, respectively. The recorded meteorological information was collected from the Crop Research Centres in Uttarakhand, India are the monthly minimum (T min ) and maximum temperatures (T max ), the wind speed (S w ), the hours of sunshine (H ss ), the monthly pan evaporation (EP m ), and the morning and afternoon relative humidity (RH 1 and RH 2 ). The EP m in monthly scale were obtained from Crop Research Centre (CRC) observatory. Figure 2(a and b) depict the box-whisker plots of the meteorological information collected over 27 years at Pantnagar station (January 1990 to December 2016), and at Ranichauri station for over 13 years (January 2000 to December 2012), respectively. These plots provide information such as the minimum, quartile-wise, and maximum values of meteorological data. The collected

Gamma test (GT)
Gamma test which includes continuous nonlinear models calculates the minimum standard error (SE) for a dataset. A data sample would be characterized as following formulas (Tsui, Jones, & De Oliveira, 2002); where, the x = (x 1 , . . . , x m ) is the vector of inputs; y i is the target; M and m are the number patterns and inputs, respectively. The GT extracts the k th (k = 1 . . . p) nearest neighbor lists (p is a user defined parameters) x N[i,k] (i = 1 . . . M) of the input vector (x i ). Finally, the GT calculates the following four statistics (Remesan, Shamim, Han, & Mathew, 2009): where, | . . . | indicates the Euclidean distance. Equation (3) indicates the gamma function of the output as : where, y N[i,k] is the output of the k th nearest neighbor of x i . To calculate the gamma ( ), a linear function is obtained using the p points { M (k), γ M (k)} as follows : where, A = gradient, = intercept (δ = 0), and y = output vector. Small values of the (near to 0) shows that the input parameter has been selected more logically. The standard error (SE) of is computed to measure the gamma value reliability. Small values of the SE shows that the gamma value is more reliable. V ratio reveals the predictability of model, and defined as : where is the gamma function and σ 2 (y) shows the output variance. The smaller V ratio indicates the predictability of the model output.  It will be possible to create a high-quality mathematical model (e.g. MM-ANN, MARS, MGGP, SVM, and M5Tree) if the mentioned four factors i.e. A, SE, , and V ratio are small. In this research, the best predictive variables were picked using the obtained lowest values of and V ratio (Malik, Kumar, Ghorbani, et al., 2019;Piri et al., 2009) for both stations.

Multiple model-artificial neural networks (MM-ANN)
In this research, the multiple model-artificial neural networks (MM-ANN), which is a hybrid model was used as to simulate the evaporation process. This hybrid model includes two levels of modeling learning. Level 1 manages the primary learning process, and the main inputs (meteorological variables) and output (pan evaporation) variables are used for training the candidate AI models, this is in one land. On the other hand, the Level 2 modeling process is initiated based on the Level one result. In this manner, a kind of binary learning process of machine learning modeling strategy is initiated. The outputs of level 1 learning process are utilized as input attributes while the original output (i.e. pan evaporation) is used as outputs in Level 2 (See Figure 3). The established modeling strategy in the current research evidenced its potential in multiple hydrological engineering application such as soil cation exchange capacity, streamflow prediction, inflow detection (Ghorbani, Khatibi, Karimi, Yaseen, & Zounemat-Kermani, 2018;Kashani, Ghorbani, Shahabi, Naganna, & Diop, 2020;Khatibi, Ghorbani, Jani, & Servati, 2018;Khatibi, Ghorbani, & Pourhosseini, 2017). Yet, for the evaporation process to be investigated. This is due to the fact; the evaporation process is highly stochastic and non-linear characterized and thus such a binary 'multiple model learning' is needed here.

Multivariate adaptive regression spline (MARS)
The MARS model represents a linear model that can auto-simulate parametric interactions and nonlinearities (Friedman, 1991). This model has a backward and forward stepwise plan, where the plan that will ensure an over-trained and complex model after series of splits but with a lower accuracy is chosen at the forward stepwise plan (De Andrés, Lorca, de Cos Juez, & Sánchez-Lasheras, 2011). However, the previously chosen nonobligatory variables are removed at the backward stepwise plan. This function uses two functions (e.g. a knot -a value of a variable) that represents the meeting point within range of inputs to map X to Y (Adamowski, Chan, Prasher, & Sharda, 2012). Two functions [(x-t) + and (t-x) + ] are deployed in MARS. Here, '+' represents only the positive as follows (Deo, Kisi, & Singh, 2017;Kisi, 2015): The MARS formula is written as following equation: where c i is the constant coefficients, and B i (x) indicates a weighted sum of function. To run the MARS model, the Salford Predictive Modeler 8 software was utilized.

Multi-gene genetic programming (MGGP)
The MGGP is the advanced stage of classical genetic programming (GP), which improves the fitness of classical GP by linearly combining low depth GP trees (Searson, 2015). Several types of research about different genotype GP have been reported in various fields by (Danandeh Mehr & Nourani, 2017;Guven, 2009;Shoaib, Shamseldin, Melville, & Khan, 2015;Zorn & Shamseldin, 2015). Figure 4 represents the tree structure while the function is as follows (Danandeh Mehr, Kahya, & Olyaie, 2013): The major inputs of GP include (i) training/validation pattern; (ii) function for fitness (iii) leaves and inner nodes for identifying the structure; and (iv) the syntax tree formation parameters (Danandeh Mehr & Nourani, 2017). For simple problems, the arithmetic parameters including plus (+), minus (-), divide (÷), and multiply (x) have been used whereas for more complex problems, other operators including exp, cos, sin, and tan are applied (Searson, 2015). For example, a pseudo-linear MGGP chromosome predicts the predictandŷ using several input variables as shown in Figure 5. The mathematical expression can be expressed as (Danandeh Mehr & Kahya, 2017;Danandeh Mehr, Kahya, & Yerdelen, 2014): whereŷ is the estimated value; T i is the i th gene output; b 0 is the noise; and b 1 , b 2 , . . . , b G are the gene weights.

Support vector machine (SVM)
The SVM is mainly useful in regression and classification tasks. The model performs using the concepts of structural risk minimization (SRM) and traditional empirical risk minimization (ERM). Both concepts are normally applied by conventional neural networks (CNN). All the input space operations in the potentially low-dimensional feature space are performed by the kernel function in SVM. The more information about the SVM and its applications can be found in some studies such as reference evapotranspiration modeling (Kişi & Cimen, 2009), pan evaporation simulation (Goyal et al., 2014), and drought modeling . In this research, the SVM model was conducted by considering a gamma regularization parameter value of 0.0064 and polynomial kernel function with 3rd degree.

M5tree model
The M5Tree determines a relationship between predictive and output variables by obtaining the parameter subspaces based on the decision tree (Sharafati, Khosravi, et al., 2019). There are two stages in the structure of the M5Tree model (Salih et al., 2020): (i) dividing the data into the subsets and form decision trees by reducing standard deviation (SDR); and (ii) fitting a linear function to the actual dataset. The SDR values are obtained by Equation (11) (Rahimikhoob, 2014): where, T i is the subset of dataset having the i th outcome of the potential set; T is a dataset that reach the node; and sd shows the standard deviation. The more information about the M5Tree model can be found in different studies such as reference evapotranspiration modeling (Rahimikhoob, 2014), pan evaporation modeling , drought forecasting , and sediment modeling (Goyal, 2014). In this research, the M5Tree was run using the Weka 3.9.2 software.

Model performance evaluation indicators
The potential of the MM-ANN, MARS, MGGP, M5Tree, and SVM techniques for simulating the EP were assessed using Nash-Sutcliffe Efficiency ( (Malik, Kumar, & Singh, 2019;Singh, Pal, & Singh, 2010). The mentioned performance criteria can be expressed as:

Selecting the GT based effective input variables
Among several learning process components, input variables selection is an essential one for modeling AI models. This is owing to the fact AI models behave differently in different cases based on the attributes of the simulated matrix. Hence, in the current study different input combinations were constructed incorporating various climatological variables and their positive influence on the pan evaporation at the inspected Ranichauri and Pantnagar meteorological stations (See Table 2). The input combinations were configured using the gamma test where the highly influential variables are determined. The statistical results of GT were reported in Table 3 for both stations. According to the results tabulated in Table 3 and with a fixed mask example (111111), it appears that the minimum value of = 0.070 and V ratio = 0.062 at Ranichauri station. Whereas, the minimum value of = 0.515 and V ratio = 0.068 at Pantnagar station. The mask example demonstrated the incorporation of the six input variables to predict the monthly scale pan evaporation. Hence, the following input variables were utilized for the prediction process (i.e. T max , T min , RH 1 , RH 2 , H ss , and S w ).

Simulation of EP m at Pantnagar station
The monthly pan evaporation was estimated using MM-ANN-1, MARS-1, MGGP-1, SVM-1 and M5Tree-1 models based on NSE, WI, LM, RMSE and MAPE for both training and testing stages at Pantnagar station. The values of NSE, WI, LM, RMSE and MAPE criteria during the training and testing periods for MM-ANN-1, MARS-1, MGGP-1, SVM-1 and M5Tree-1 models are presented in Table 4. As evaluated for Pantnagar station from Table 4 Table 4 exposed the MM-ANN-1 model performed the best simulation during the testing periods. Therefore, MM-ANN-1 model followed the best statistical criteria (i.e. minimum RMSE and MAPE values, and maximum NSE, WI and LM values) in testing stage and selected best model among other models. MARS-1 model followed the MM-ANN-1 as the second rank closely. The temporal variation between simulated and observed monthly EP m data, and their scatter plots (right side) for the MM-ANN-1, MARS-1, MGGP-1, SVM-1, and M5Tree-1 models was plotted in Figure 6 (a through e) for the testing period. In scatter plots, the regression line provided the coefficient of determination (R 2 ) as 0.955 for the MM-ANN-1 model, 0.959 for the MARS-1 model, 0.956 for MGGP-1 model, 0.946 for SVM-1 model, and 0.962 for M5Tree-1 model during the testing stage, respectively. The fitted regression line (RL) and the perfect line of fit (1:1) were close to each other for all techniques. The RL was located above the best fit (1:1) for MARS-1, MGGP-1, SVM-1, and M5Tree-1 models. This means that at Pantnagar station, the five models overpredict the monthly pan evaporation values. Also, the regression line was located under the 1:1 line for MM-ANN-1 models, which means the model under-predict the pan evaporation values at Pantnagar station.
According to Figure 7, all model symbols were very close to each other. In detail, SVM-1 was located as the furthest from the observed point. This introduces SVM-1 as the worst model. MM-ANN-1 was the closest   Figure 8. However, the density plot suggested a slightly different distribution of the simulated values.

Simulation of EP m at Ranichauri station
The MM-ANN-1, MARS-1, MGGP-1, SVM-1 and M5Tree-1 models were trained using (January, 2000to December, 2009) and tested using (2010to December, 2012 datasets. The adequacy of these models was assessed using the NSE, WI, LM, RMSE and MAPE values in the both training and testing phases. Table 4 Table 4 revealed that the MARS-1 model produced the best performances and results in the training stage while failed in the testing stage. The MGGP-1 model tracked the best statistical criteria (i.e. minimum RMSE and MAPE values, and maximum NSE, WI and LM values) over the testing stage and indicated the best performance compared to other models. MM-ANN-1 model follows the MGGP-1 as the second-best model. The scatter plots and temporal variations graphs of the observed against the simulated EP m of the MM-ANN-1, MARS-1, MGGP-1, SVM-1, and M5Tree-1 models over the testing phase were presented in Figure 9 (a through e). The RL in scatter plots marked R 2 as 0.910 for the MM-ANN-1, 0.910 for the MARS-1, 0.916 for the MGGP-1 model, 0.850 for SVM-1 model, and 0.852 for M5Tree-1 model over the testing period. The different conditions have been yielded between RL and 1:1 line for all the models. For the SVM-1 model (Figure 9(d)), the RL was exactly over the 1:1 line, it over-predicted the EP m values. For the M5Tree-1 model (Figure 9(e)), the regression line was exactly below the best fit line (1:1), it under-predicted the EP m values. In the case of MM-ANN-1, MARS-1 and MGGP-1 models (Figure 9(a-c)), the regression line divided the best fit line (1:1). It can be explained that the models over-predicted the smaller values ( < 2.5 mm/month for MM-ANN-1, < 3.0 mm/month for MGGP-1, < 2.5 mm/month for MARS-1) and under-predict the higher EP m values at Ranichauri station.
According to Figure 10, the highest correlation coefficient (CC) was observed in the MARS-1. But the standard   deviation showed much less than the observed values. The SVM-1 and M5Tree-1 model is far from being competing on the best performance, MGGP-1 can be considered as a high performed model with lowest RMSE and high correlation although it has a standard deviation lower than the observed values. Finally, MGGP-1 was the closest model to the observed values based on the standard deviation, correlation, and RMSE. This was giving superiority to MGGP-1 over the others. Considering the density plot (Figure 11), the MGGP-1 model, which identified as the closest model in the Taylor diagram yielded a similar box from the observed values.
The findings of the current research demonstrated that the MM-ANN at Pantnagar station and MGGP model at Ranichauri station achieved better performances in simulating EP m values. The MM-ANN-1 and MGGP-1 models were selected as the best ones based on minimum RMSE and MAPE values, and maximum NSE, WI, and LM values at the testing period. Various researchers reported that the soft computing models provide better performances during the testing period, and the best performing procedure was selected (Abdullah, Malek, Abdullah, Kisi, & Yap, 2015;Guven & Kisi, 2013). Malik, Kumar, and Kisi (2017) estimated EP m using SOMNN, CANFIS, MLPNN, RBNN, Griffith's (G), and Stephens-Stewart (SS) models at Ranichauri  and Pantnagar stations. They published that the specific heuristic approaches (i.e. MLPNN and CANFIS models) outperformed the other models. It can be clear that the previous investigations confirmed the multiple model strategies for modeling nonlinear time series phenomena. The findings of this research were also in close agreement with the study done by (Malik, Kumar, & Kisi, 2017). Future research can be devoted to the incorporation of other casual related climate variables such as vapor pressure deficit, atmospheric pressure or others.

Conclusion
The evaporation process of the hydrological cycle is characterized by highly complex and stochastic phenomena in nature. In the current research, the ability of multiple model strategies (MM-ANN, MARS, MGGP, SVM, and M5Tree) was evaluated in simulating monthly pan evaporation at two meteorological stations, Ranichauri and Pantnagar in India. The best input variables combination was selected using GT. The results of applied models were examined using some performance criteria and visual presentations. In general, we derived five conclusions from the current research as follows: • Six appropriate variables (RH 1 , RH 2 , S w , H ss , and T min , T max ,) were selected by GT as the optimistic input combination for simulating the pan evaporation in monthly scale at both studied stations. The nature of the EP m process in this region is highly stochastic and more climatic information is needed for building the prediction matrix of the proposed AI models. • At the Pantnagar station, the MM-ANN-1 model indicated the better performance compared to the MARS-1, MGGP-1, SVM-1 and M5Tree-1 models. However, the feasibility of the MGGP-1 model demonstrated batter results at Ranichauri station. This is normal for the case where AI models behave differently in different cases based on the actual internal mechanism between the predictand and predictors. • MARS model performed as the second accurate model following the MM-ANN at Pantnagar station, while the MM-ANN was following the MGGP model capacity in mimicking the EP m trend at Ranichauri station. • Overall, the attained results of the proposed MM-ANN and MGGP models demonstrated an optimistic intelligence approaches for this particular region where they would contribute to the water resources engineering practice, and management.