A Hybrid Data-Driven Machine Learning Technique for Evapotranspiration Modeling in Various Climates

Abstract

In the current research, gene expression programming (GEP) was applied to model reference evapotranspiration (ETo) in 18 regions of Iran with limited meteorological data. Initially, a genetic algorithm (GA) was employed to detect the most important variables for estimating ETo among mean temperature (Tmean), maximum temperature (Tmax), minimum temperature (Tmin), relative humidity (RH), sunshine (n), and wind speed (WS). The results indicated that a coupled model containing the Tmean and WS can predict ETo accurately (RMSE = 0.3263 mm day⁻¹) for arid, semiarid, and Mediterranean climates. Therefore, this model was adjusted using the GEP for all 18 synoptic stations. Under very humid climates, it is recommended to use a temperature-based GEP model versus wind speed-based GEP model. The optimal and lowest performance of the GEP belonged to Shahrekord (SK), RMSE = 0.0650 mm day⁻¹, and Kerman (KE), RMSE = 0.4177 mm day⁻¹, respectively. This research shows that the GEP is a robust tool to model ETo in semiarid and Mediterranean climates (R² > 0.80). However, GEP is recommended to be used cautiously under very humid climates and some of arid regions (R² < 0.50) due to its poor performance under such extreme conditions.

1. Introduction

Evaluation of reference evapotranspiration (ETo) plays an important and undeniable role in irrigation scheduling, drought analysis, climate change studies, water level balance, agricultural and forest meteorology, long-term decision-making in food and water security policies, and optimum allocation of water resources [1,2,3]. Although several methods have been developed to predict ETo throughout the world, there is a limited number of models to estimate ETo where meteorological data is restricted or insufficient [1,2].

A solution to deal with this limitation is to use data-driven machine learning techniques, particularly genetic approaches, including genetic algorithm (GA) and gene expression programming (GEP). One advantage of the GEP to estimate ETo was that, unlike the artificial neural networks (ANN) method, the GEP generated an explicit model structure that can be easily comprehended and adopted [3].

The GA and GEP methods have been developed in various aspects of water resources such as streamflow forecasting [4], rainfall–runoff modeling [5,6,7,8], modeling transport streams with suspended sediment [9], predicting velocity in compound channels [10], characterizing risks in water supply systems [11], and modeling evaporation [12].

Although the GEP has been developed in water resources studies, the application of this technique for ETo modeling is limited. Some of the successful applications of the GEP to estimate ETo can be listed as follows.

Parasuraman et al. [3] modeled ETo using only ground temperature and net radiation. Irmak and Kamble [13] investigated evapotranspiration data assimilation with the GA and soil, water, atmosphere, and plant (SWAP) model for on-demand irrigation. The data assimilation methodology obtained from the present research can be considered a practical tool at the field scale for scheduling the irrigation estimated by remote sensing-based evapotranspiration. The results showed that the GA was effectively able to determine the terms included in the fitness function, and parameters were predicted reasonably, especially if only four variables were included. Traore and Guven [14] developed regional-specific numerical models of ETo using the GEP in Sahel. Statistically, the GEP was an effectual modeling tool for the successful computation of ETo under the study area. The results indicated that, using the GEP model, it would be possible to formulate an accurate and applicable numerical equation for each region by irrigation systems for the manual computation of ETo under the study area, where sufficient meteorological variable is often missing. Shiri [15] claimed that GEP outperforms empirical model to calculate ETo in hyper arid regions over Iran. Traore et al. [16] found that GEP is a robust tool to model ETo in in Jiangsu province, China. Mattar [17] demonstrated that GEP is a more accurate method than empirical equations to estimate ETo in Egypt.

Machine learning methods are needed in data sparse regions because Penman–Monteith equation requires a lot of data input, and many places do not have such data. Therefore, achieving similar results with fewer input data is desirable.

In addition, some studies reported greater accuracy of the GA/GEP rather than the ANN [18,19,20,21], Support Vector Regression (SVR) [20,22], Adaptive Neuro-Fuzzy Inference System (ANFIS) [1,2,23], and empirical models [2,24].

According to the literature review there are some investigations for which the superiority of the GA and GEP over other techniques such as the ANN, LP, ANFIS, and regression models have been reported. However, genetic approaches have not been simultaneously assessed under various climates. In addition, in most of previous studies, these methods were employed with many input data to estimate ETo. In different parts of the world, there is insufficient meteorological data due to the lack of synoptic stations or other limitations. Therefore, the current research seeks to evaluate ETo in four different climates in Iran to recognize the parameters with the most significant roles in ETo.

2. Materials and Methods

2.1. Designing Structure of Genetic Algorithm (GA) for the Current Research

In most of the previous investigations [1,2], a power function is presented as the optimal model to estimate ETo after modifying or calibrating the empirical models. Therefore, in the present research, various functions of the mean temperature (Tmean), maximum temperature (Tmax), minimum temperature (Tmin), relative humidity (RH), wind speed (WS), and sunshine (n) were defined and then combined by using summation and multiplying functions. The goal is to minimize the difference between the current functions with the FAO-Penman–Monteith (FPM) (

Table 1. Structure of genetic algorithm (GA) designed in this study.

Base	Formula
Mean Temperature	$ETo=a_1{\left(T_{mean}\right)}^{a_2}+a_3$
Differential Temperature	$ETo=b_1{\left(T_{\max }-T_{\min }\right)}^{b_2}+b_3$
Relative Humidity	$ETo=c_1{\left(RH\right)}^{c_2}+c_3$
Wind Speed	$ETo=e_1{\left(WS\right)}^{e_2}+e_3$
Solar Radiation	$ETo=f_1{\left(n\right)}^{f_2}+f_3$
Total (Summation)	${\left(ETo\right)}_{GA}={\displaystyle {\displaystyle \sum }_{i=1}^5}ET{o}_i$
Total (Multiplying)	${\left(ETo\right)}_{GA}={\displaystyle {\displaystyle \prod }_{i=1}^5}ET{o}_i$
Goal Function	$Minimize\left[{\left(ETo\right)}_{GA}-{\left(ETo\right)}_{FPM}\right]$

).

In the next step, a GA program was coded in MATLAB environment using the values shown in

Table 2. Parameters of genetic algorithm (GA) coded in this study by using MATLAB environment.

Value	Option	Value	Option
1.00E-06	TolFun	doubleVector	PopulationType
1.00E-06	TolCon	[2 × 1 double]	PopInitRange
10	InitialPenalty	20	PopulationSize
100	PenaltyFactor	2	EliteCount
1	PlotInterval	0.8	CrossoverFraction
@gacreationuniform	CreationFcn	Forward	MigrationDirection
@fitscalingrank	FitnessScalingFcn	20	MigrationInterval
@selectionstochunif	SelectionFcn	0.2	MigrationFraction
@crossoverscattered	CrossoverFcn	100	Generations
{[ @mutationgaussian ] [1] [1]}	MutationFcn	Inf	TimeLimit
Final	Display	-Inf	FitnessLimit
Off	Vectorized	50	StallGenLimit
Never	UseParallel	Inf	StallTimeLimit
1.00E-06	TolFun	doubleVector	PopulationType
1.00E-06	TolCon	[2 × 1 double]	PopInitRange

.

2.2. Parameters of Gene Expression Programming (GEP) Employed in the Current Research by Using GeneXpro Tools Version 5.0

Table 3. Parameters of gene expression programming (GEP) employed in this study by using GeneXpro Tools version 5.0.

Value	Option
+, −, ×, ÷, ^, √, ∛, ∜, sin, cos, ln, exp	Function set
500	Population size
95	Mutation frequency
50	Crossover frequency
217	Number of replication
30	Block mutation rate
30	Instruction mutation rate
40	Instruction data mutation rate
95	Homologous crossover
Ranking	Selection method
1	Recombination rate
0.2	Mutation rate
Generalize	Alternative method
Yes	Elitism
6	Maximum level of tree

shows the values of each option employed to design the GEP structures by using GeneXpro Tools version 5.0.

As shown in Table 3, the population size of 500 is the same as most previous research.

However, different functions are employed to enhance the results and to reduce the uncertainty in the research, and the number of replication is increased (almost 217). As a result, the elapsed time to run the GEP can be prolonged. In addition, the ranking method showed better performance compared to Roulet’s function (a fitness evaluation function).

2.3. Materials

The monthly averages of meteorological data were collocated from the Islamic Republic of Iran Meteorological Organization (IRIMO). These data contain mean, minimum, and maximum daily air temperature (°C), saturated vapor pressure deficit (kPa), mean and minimum relative humidity (%), wind speed (m s⁻¹) and direction, rainfall (mm month⁻¹), cloudy days, and sunshine (hr month⁻¹).

Table 4. Position and climate of the stations with length of collected data.

Station Name	ICAO Code	North Latitude	East Longitude	Altitude (masl)	Start Year	End Year	Climate
Ahvaz	40811	31°20′	48°40′	22.5	1961	2010	Arid
Arak	40769	34°6′	49°46′	1708.0	1961	2010	Semiarid
Bushehr	40858	28°58′	50°49′	9.0	1961	2010	Arid
Esfahan	40800	32°37′	51°40′	1550.4	1961	2010	Arid
Hamedan	40768	34°52′	48°32′	1741.5	1961	2010	Semiarid
Jiroft	40866	28°35′	57°48′	601.0	1989	2009	Arid
Kerman	40841	30°15′	56°58′	1753.8	1961	2010	Arid
Mashhad	40745	36°16′	59°38′	999.2	1961	2010	Semiarid
Moghan	40700	39°39′	47°55′	31.9	1984	2010	Semiarid
Qazvin	40731	36°15′	50°3′	1279.2	1961	2010	Semiarid
Rasht	40719	37°19′	49°37′	-8.6	1961	2010	Very humid
Sanandaj	40747	35°20′	47°0′	1373.4	1961	2010	Mediteranean
Shahrekord	40798	32°17′	50°51′	2048.9	1961	2010	Semiarid
Shiraz	40848	29°32′	52°36′	1484.0	1961	2010	Semiarid
Tabriz	40706	38°5′	46°17′	1361.0	1961	2010	Semiarid
Urmia	40712	37°40′	45°3′	1328.0	1961	2010	Semiarid
Yazd	40821	31°54′	54°17′	1237.2	1961	2010	Arid
Zabol	40829	31°2′	61°29′	489.2	1961	2010	Arid

shows the position of all 18 synoptic stations and their climates. It should be noted that Iran is in an arid and semiarid zone in Persian Gulf region, and the major part of Iran has the same climate mentioned in Table 4. Figure 1 shows location of the weather station in Iran.

Among all stations, there is 50-year period information for 16 regions. The accuracy of the Food and Agricultural Organization of the United Nations (FAO)−Penman−Monteith (FPM) is confirmed for these regions by using lysimeter measurements and with respect to the previous investigations [25,26,27,28,29,30,31,32,33,34]. In addition, there are 27-year and 21-year period datasets for Moghan and Jiroft, respectively, which due to the availability of lysimeter data, confirms the accuracy of the FPM for these two regions; they were also added to other 16 regions (with 50–year data).

The authors have used the last five years as the testing period and the rest of the data as training period.

Table 5. Selected models to estimate reference evapotranspiration including their references, formulae, and parameters.

Model	Reference(s)	Formula	Parameters
FAO Penman-Monteith (FPM)	[35]	$E{T}_o=\frac{0.408(R_n-G)+\gamma \frac{900}{T+273}u(e_s-e_a)}{\Delta +\gamma \left(1+0.34u\right)}$	H,φ,T,Tmin,Tmax,RH,u,n
Hargreaves-Samani (HS)	[36]	$E{T}_o=0.005508K_T{R}_a{\left(T_{\max }-T_{\min }\right)}^{0.5}\left(T+17.8\right)$	T,u,Tmin,Tmax,RH,n,φ
Jensen-Haise (JH)	[37]	$E{T}_o=0.408C_T\left(T-T_x\right)R_s$	T,Rs,H,Tmax,Tmin
WMO	[38]	$E{T}_o=\left(1.298+0.934u\right)\left(e_s-e_a\right)$	T,Tmin,Tmax,RH,u
Abtew (Ab)	[39]	$E{T}_o=0.01786\frac{R_s{T}_{\max }}{\lambda }$	T,Tmax,Rs
Makkink (Mk)	[40]	$E{T}_o=0.61\frac{\Delta }{\Delta +\gamma }\frac{R_s}{\lambda }-0.12$	T,Rs
Turc (Tu)	[41,42]	$E{T}_o=\left(0.3107R_s+0.65\right)\frac{T{a}_t}{T+15}$	T,RH,Rs
Blaney-Criddle (BC)	[43]	$E{T}_o=a+bP\left(0.46T+8.13\right)\left(1+0.0001H\right)$	H,T,n,RHmin,φ,u
Modified Hargreaves-Samani 1 (MHS1)	[44]	$E{T}_o=0.0005304R_a{\left(T_{\max }-T_{\min }-0.0123R\right)}^{0.76}\left(T+17\right)$	T,Tmin,Tmax,φ,R
Modified Hargreaves-Samani 2 (MHS2)	[45]	$E{T}_o=0.0009384R_a{\left(T_{\max }-T_{\min }\right)}^{0.424}\left(T+17.8\right)$	T,Tmin,Tmax,φ

shows the selected models with their references and parameters.

Among the numerous empirical methods to estimate ETo, nine models were selected that had the best performance under different climates on the basis of previous investigations [1,46,47,48,49,50], and the results were compared with the FPM.

ETo is the reference crop evapotranspiration (mm/day), R_n is the net radiation (MJ/m²/day), G is the soil heat flux (MJ/m²/day), γ is the psychrometric constant (kPa/°C), e_s is the saturation vapor pressure (kPa), e_a is the actual vapor pressure (kPa), Δ is the slope of the saturation vapor pressure–temperature curve (kPa/°C), T is the average daily air temperature (°C), u is the mean daily wind speed at 2 m (m/s), H is the elevation (m), φ is the latitude (rad), T_min is the minimum air temperature (°C), T_max is the maximum air temperature (°C), RH is the average relative humidity (%), n is the actual duration of sunshine (hr), R_s is the solar radiation (MJ/m²/day), R_a is the extraterrestrial radiation (MJ/m²/day), λ is the latent heat of vaporization (MJ/kg), and C_T, I, K_T, a, b, P, a_t, and T_x are empirical coefficients

Genetic algorithm (GA), genetic programming (GP), and gene expression programming (GEP) are different kinds of data-driven machine learning techniques to find an optimum solution for complex problems by artificial intelligence (AI). In a cell, the expression of the genetic information is an intricate procedure involving more than one hundred molecules. Two of the major players are DNA and proteins. DNA is the carrier of the genetic information and the proteins read and express the genetic information. GA employed biological evolution theory for computer applications (i.e., AI). In fact, GA is an oversimplification of biological evolution [51]. GP, introduced by [52,53], solves the problem of fixed length solutions (as stated for GA) by creating nonlinear entities. Each entity (parse tree) has a distinguished shape and size. GEP is the reliable modification of GP and GA, combining both the methodology of simple, linear chromosomes with fixed length (GA) and branched structures of various sizes and shapes (GP). GEP applies the same type of diagram representation of GP, however the entities (expression trees) evolved by GEP are the expression of a linear genome. The GEP is a constructive approach because the search operators of the GEP may constantly generate valid structures and are highly suited to genetic diversity. The GEP surpasses the old GP system for more than 100 times. In this study, the GA was employed to understand which parameters are more influential on the estimation accuracy of ETo. In addition, the GEP was used to model ETo considering the obtained results by the GA to achieve the maximum accuracy compared with the empirical methods. Therefore, a hybrid data-driven machine learning technique (considering both GA and GEP) was applied to model ETo in various climates.

To evaluate the accuracy of the models two indices were used as follows

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(X_i-Y_i\right)}^2}}{N}}$$

(1)

$$R^2=\frac{{\left[{\displaystyle \sum _{i=1}^n\left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}\right]}^2}{{\displaystyle \sum _{i=1}^n{\left(X_i-\overline{X}\right)}^2{\displaystyle \sum _{i=1}^n{\left(Y_i-\overline{Y}\right)}^2}}}$$

(2)

where, X_i and Y_i are the ith observed and estimated values, respectively; $\overline{X}$ and $\overline{Y}$ are the mean of X_i and Y_i; and N is the total numbers of data.

3. Results

3.1. Performance of Genetic Algorithm (GA) for Different Parameters

Table 6. Performance of genetic algorithm (GA) for different parameters.

Code	Parameters	Root Mean Square Error (RMSE) (mm/day)
GA1	Tmean	0.6823
GA2	Tmean,Tmin,Tmax	0.5049
GA3	Tmean,RH	0.5733
GA4	Tmean,n	0.6250
GA5	Tmean,WS	0.3263
GA6	Tmean,WS,RH	0.3087
GA7	Tmean,WS,n	0.3211
GA8	Tmean,WS,RH,n	0.3030

shows that although GA8 (combination of all parameters) was obtained with the highest accuracy (RMSE = 0.3030 mm day⁻¹) in Iran (average), it requires four parameters: Tmean, RH, WS, and n. On the other hand, GA5 is also obtained, which is approximately equal to GA8 (RMSE = 0.3263 mm day⁻¹) with less input data (only Tmean and WS). Meanwhile, GA5 estimated ETo throughout Iran with more accuracy than GA2, GA3, and GA4. This means that the WS was introduced as the most important factor (after Tmean) to control the dynamics of ETo process throughout Iran.

Table 6 reveals that a function of the Tmean and WS may predict ETo with acceptable accuracy. Therefore, this result is the basis for the development of the genetic models using the GEP. However, the Tmean is the first parameter to be measured in each station or region. Hence, the Tmean is considered as the input parameter in the entire GEP models in this study.

3.2. Performance of Gene Expression Programming (GEP) for Mean Temperature (Tmean) as Input Data

In the first step of the GEP—ETo—was estimated in all 18 stations using only the Tmean and then compared with the FPM (

Table 7. Performance of gene expression programming (GEP) for mean temperature (Tmean) as input data.

Region	The Best Structure	Training RMSE (mm/day)	Testing RMSE (mm/day)
Ahvaz	$ETo=\ln \left(58.598+Tmean\right)+\ln \left(\ln \left(Tmean+4\right)+Tmean+\sin \left(T+5+\sin \left(Tmean\right)\right)\right)$	0.5789	0.6839
Arak	$ETo=\sqrt{-\sin \left(Tmea{n}^2+11Tmean+28\right)-0.821-\sin \left(3Tmean\right)+{\left(1.386Tmean+8.318\right)}^{0.754}}$	0.5321	0.2542
Bushehr	$ETo=-0.011\cos \left(8Tmean\right)-\left(4.289+1.464Tmean\right)$	0.3230	0.2190
Esfahan	$ETo=\ln \left({\left(\ln \left(Tmean\right)\right)}^{2.197}+4^{\ln \left(Tmean\right)}+\sin \left(Tmea{n}^6+1.609\right)+\ln \left(2.398+{0.693}^{\ln \left(Tmean\right)}\right)\right)$	0.1200	0.3361
Hamedan	$ETo=\ln \left(\ln \left(\frac{\exp \left(4Tmean\right)}{\ln \left(\frac{3Tmean+3}{Tmean}\right)}\right)\right)$	0.2614	0.2844
Jiroft	$ETo=\exp \left(1.016\left(\ln \left(5.307-\cos \left(Tmean\right)+\cos \left(Tmean-1\right)\right)\right)\right)$	0.1980	0.2963
Kerman	$ETo=5-\cos \left(4-\cos \left(2-\cos \left(Tmean\right)\right)\right)$	0.2090	0.5285
Mashhad	$ETo={\left(Tmean+\sin {\left(Tmean+3\right)}^{1.609}+{\left(\sin \left(2.221\right)\right)}^{Tmean}\right)}^{\ln \sqrt{\ln \left(Tmean\right)}}$	0.1417	0.2707
Moghan	$ETo=\exp \left(3^{0.209{\cos }^{16}\left(Tmea{n}^3\right)}\right)$	0.2871	0.3207
Qazvin	$ETo=\sqrt[4]{Tmean-3}-0.945+\cos \left(\sin \left(\exp \left(\frac{Tmean}{3}\right)\right)\right)$	0.2508	0.3355
Rasht	$ETo=\ln \left(\frac{Tmean}{12-48Tmean}-10-1.115\sin \left(Tmean\right)-\cos \left(Tmean\right)\right)$	0.2391	0.1100
Sanandaj	$ETo=\sqrt{Tmean+0.041+\sin \left(0.414-Tmean\right)+\cos \left(\sqrt{Tmean+2}-6\right)}$	0.1478	0.2858
Shahrekord	$ETo=\sqrt{12.683-{\left(\ln \left(4Tmean-2.079\right)\right)}^{\cos \left(\sin \left(Tmea{n}^5\right)\right)}}$	0.2172	0.2667
Shiraz	$ETo=\cos \left(\sin \left(2-\sin \left(Tmean-3\right)\right)\right)+\frac{0.336}{Tmean-4.375}+4$	0.2251	0.2918
Tabriz	$ETo=\ln \left(\ln \left(\cos \left(Tmean+4\right)+4.255\right)\right)+\sqrt{1.302+Tmean}$	0.2873	0.3002
Urmia	$ETo=\sqrt{Tmean-\exp \left(\sin \sqrt{Tmean-8}\right)}$	0.2036	0.2396
Yazd	$ETo=\frac{5.071+\sqrt{6+Tmean}}{\sqrt{\frac{Tmean}{4}}\cos \left(\ln \left(\sqrt{Tmean-17}\right)\right)}$	0.4517	0.5519
Zabol	$ETo=\frac{4.256Tmea{n}^2}{\ln \left(Tmea{n}^3+2.545Tmean+2980.958\right)}$	0.7208	0.8229

).

According to Table 7, the best performance of the GEP belongs to Rasht (RMSE = 0.1100 mm day⁻¹), while, the worst accuracy was reported for Zabol (RMSE = 0.8229 mm day⁻¹). In 56% of regions the GEP was obtained with a RMSE < 0.3000 mm day⁻¹. Furthermore, the natural logarithm (ln) function was employed more than sinus, cosines, and particularly exponential functions in the GEP structures. It means that we can expect to have logarithmic relationship between ETo and Tmean more than other kinds of functions. Moreover, plus and minus functions were used more than multiplication and specially division.

3.3. Performance of Gene Expression Programming (GEP) for Mean Temperature (Tmean), Minimum Temperature (Tmin), and Maximum Temperature (Tmax) as Input Data

In the second step of the GEP, ETo was estimated in all 18 stations using the Tmean, Tmin, and Tmax was then compared with the FPM.

The results reveal that the best performance of the GEP belongs to Rasht (RMSE = 0.0884 mm day⁻¹), while the worst accuracy was reported for Zabol (RMSE = 0.8020 mm day⁻¹). In 61% of regions the GEP was obtained with a RMSE < 0.3000 mm day⁻¹. Furthermore, sine and cosine functions were employed more than the natural logarithm (ln) and particularly exponential functions in the GEP structures. It means that we can expect to have periodic relationship between ETo and temperature more than other kinds of functions. Moreover, plus and minus functions were used more frequent than multiplication sign and particularly division sign. In all regions, the accuracy was improved compared to the GEP models based on the Tmean only.

3.4. Performance of Gene Expression Programming (GEP) for Mean Temperature (Tmean) and Wind Speed (WS) as Input Data

With respect to Table 6, WS was introduced as the most important factor to control the variations of ETo. Thus, in the third step of GEP, ETo was estimated in all 18 stations using the Tmean, and WS was then compared with the FPM (

Table 8. Performance of gene expression programming (GEP) for mean temperature (Tmean) and wind speed (WS) as input data.

Region	The Best Structure	Training RMSE (mm/day)	Testing RMSE (mm/day)
Ahvaz	$ETo=3.189-\sqrt{WS-2+Tmean}+\sqrt{7+Tmean}$	0.2008	0.2128
Arak	$ETo=\ln \left(\sqrt{\sqrt{2.197\ln \left(2Tmean\right)}\left(\left(Tmean-4\right)WS\right)+1.386}\right)$	0.0963	0.1079
Bushehr	$ETo=\ln \left(17.167WS+137.339\right)-\frac{\exp \left(\cos \frac{1}{Tmean}\right)}{\exp \left(WS\right)+\cos \left(WS\times Tmean\right)}-\frac{\exp \left(\cos \left(4-Tmean-\frac{WS}{8}\right)\right)}{WS+\exp \left(\ln \left(WS\right)-0.8\right)}$	0.2018	0.1640
Esfahan	$ETo=\sqrt{Tmean}-\cos \left(WS\right)-0.006-\ln \left(\sqrt[4]{2-\cos \left(WS\right)}\right)$	0.0835	0.0925
Hamedan	$ETo=5-\exp \left(\exp \left(\frac{Tmean\sqrt{Tmean}}{12}\right)\right)-{0.518}^{\frac{W{S}^{3.5}}{11+\frac{8}{WS}}}$	0.1108	0.1457
Jiroft	$ETo=WS+3.001-\frac{\exp \left(\sqrt{5^{\sin \left(Tmean\right)}}\right)}{WS-13}$	0.1489	0.1401
Kerman	$ETo=\sqrt[4]{WS}\times \sqrt{\ln \left(WS+\frac{Tmean}{5}\right)}\times \ln \left(Tmean+\frac{\sqrt{Tmean}\times \ln \left(WS\right)}{5}\right)$	0.2489	0.4177
Mashhad	$ETo=\exp \left(\frac{Tmean}{28}\right)+\ln \left(\frac{0.571\left(Tmean+WS\right)}{-\ln \left(\frac{WS}{7}\right)}\right)$	0.1107	0.1725
Moghan	$ETo=\ln \left(\sqrt{0.303\ln \left(Tmean\right)\times Tmean\times WS}\times \ln \left(\ln \left(9WS\right)\times \left(Tmean-\cos \left(Tmean\right)\right)\right)\right)$	0.1082	0.0906
Qazvin	$ETo=\ln \left(WS\right)+\sqrt{\ln \left(WS\right)+Tmean-\ln \left(Tmean\right)+3}$	0.0827	0.1079
Rasht	$ETo=\sqrt{\sin \left(0.745Tmean\right)+\sqrt{5.064\left(\sin \left(WS\right)+4+WS\right)}}$	0.0486	0.0961
Sanandaj	$ETo=\sqrt{Tmean\times \left(\sqrt{WS}+0.745\right)}$	0.1740	0.1287
Shahrekord	$ETo=0.705\times \left(\sqrt{\left(Tmean+7\right)}+\ln \left(WS\right)-1-\cos \left(\ln \left(2+2Tmean\right)\right)\right)$	0.0583	0.0650
Shiraz	$ETo=\sqrt{\cos \left(2WS+Tmean+26\right)+\cos \left(Tmean+5\right)+2WS+\exp \left(\sqrt{WS+6}\right)}$	0.1474	0.1403
Tabriz	$ETo=\ln \left(\exp \left(\sin \left(Tmean\right)\right)+6+Tmean+\sin \left(\sin \left(Tmean\right)\right)\right)-\frac{\sin \left(Tmean\right)}{\ln \left(5+WS\right)}+{22025.465}^{\left(WS-7\right)\left(\sin \left(WS\right)+8\right)}$	0.1487	0.1620
Urmia	$ETo=\cos \left(\frac{WS-6.135}{Tmean}\right)+\cos \left(\frac{12.96}{0.75+WS+Tmean}\right)+\exp \left(\sqrt{\frac{WS}{5}}\right)-\frac{3}{Tmean}$	0.2049	0.1158
Yazd	$ETo=\sqrt[4]{\frac{Tmea{n}^2-8Tmean}{2}}\sqrt{\left(WS+0.008\right)}$	0.1790	0.1738
Zabol	$ETo=5.016-\frac{WS}{WS-1}-\frac{0.246-\frac{WS-2}{3}}{2.273-\cos \left(Tmean\right)}$	0.3289	0.3612

).

Table 8 represents that the best performance of the GEP belongs to Shahrekord (RMSE = 0.0650 mm day⁻¹), while the worst accuracy was reported for Kerman (RMSE = 0.4177 mm day⁻¹). In 83% of regions the GEP resulted with a RMSE < 0.2000 mm day⁻¹. Furthermore, the natural logarithm (ln) function was employed more than sinus, cosines, and particularly exponential functions in the GEP structures. Moreover, plus and minus functions were used more than multiplication sign and specially division. In all regions (with the exception of Rasht), the accuracy was improved compared to the GEP models based on the Tmean, Tmin, and Tmax. Therefore, for Rasht, wind speed is not required. It may correspond to minimum diurnal temperature rate (DTR) in very humid regions compared to other climates (role of relative humidity and saturated vapor pressure). However, wind speed-based GEP models predicted ETo in arid, semiarid, and Mediterranean climates more accurate than the other models.

It should be noted that the arctan has not acceptable accuracy compared to other functions. In addition, the results represent that ln and exp functions have a better performance than sine and cosine functions. Therefore, in future investigations, we can more focus on other functions than arctan. In addition, we can expect to have logarithmic relationship between ETo and temperature/wind speed more than other kinds of functions.

The obtained results are comparable and sometimes better than those reported by Wang et al. [54]. They resulted RMSE between 0.222 to 0.555 mm day⁻¹. In addition, Mehdizadeh [55] reported RMSE for GEP models from 0.46 to 2.08 mm day⁻¹.

3.5. Performance of Gene Expression Programming (GEP) for Mean Temperature (Tmean), Wind Speed (WS), and Relative Humidity (RH) as Input Data

The results of the GA revealed that use of the RH and n do not increase accuracy of the GEP significantly (Table 6). This is also confirmed by the GEP.

According to the results, the best performance of the GEP was for Esfahan (RMSE = 0.0730 mm day⁻¹), while the worst accuracy was reported for Kerman (RMSE = 0.4252 mm day⁻¹). In 89% of regions the GEP resulted with a RMSE < 0.2000 mm day⁻¹. Furthermore, sinus and cosines functions were employed more than the natural logarithm (ln) and particularly exponential functions in the GEP structures. Moreover, plus and minus functions were used more than multiplication sign and specially division.

A comparison of Table 8 indicates that adding the RH as input parameter not only did not increase the GEP’s accuracy, but also led to the reduction in the accuracy relevant to 56% of the regions. In addition, the improvements are inconsiderable. It should be noted that the best structures of the GEP are not a function of the RH in 44% of the regions (Arak, Bushehr, Mashhad, Moghan, Qazvin, Sanandaj, Shahrekord, Shiraz, Tabriz, and Yazd). It is an important result and confirms that the Tmean and WS have more value to be used as input variable for the GEP compared to the RH.

3.6. Accuracy of Empirical Models Against Gene Expression Programming (GEP)

Figure 2 and Figure 3 show that the best empirical models with respect to the smallest RMSE values for estimating ETo on a monthly scale lack a good performance at annual scale and this is alarming considering the importance of the agricultural water management to deal with water crisis in the world. However, the GEP may predict ETo in the most regions of Iran with considerable accuracy.

In arid regions (Figure 2), the best accuracy belongs to Esfahan (R² = 0.9217) and the worst accuracy belongs to Kerman (R² = 0.3636), due to poor performance of the GEP for peak events in this region. In the semiarid regions (Figure 3), the best accuracy belongs to Shahrekord (R² = 0.9403) and the worst accuracy belongs to Hamedan (R² = 0.8049).

The obtained results are comparable and sometimes better than those reported by Wang et al. [54]. They resulted R² between 0.639 to 0.944. In addition, Mehdizadeh [55] reported R² for GEP models from 0.084 to 0.969.

4. Discussion

The current research shows that the GEP is a robust tool to model ETo in semiarid and Mediterranean climates (R² > 0.80). However, the use of the GEP should be recommended cautiously in Rasht, Bushehr, and Kerman (R² < 0.50) due to its poor performance in extreme events of these areas. To enhance the accuracy of the GEP in extreme values, use of hybrid methods (coupling GEP with ANN, fuzzy logic, honey bee algorithm, etc.) can be recommended.

The results of this study has a potential to be compared with real data. However, reference evapotranspiration (ETo) refers to maximum evapotranspiration in the best growing conditions of a crop without any limitation (drought and salinity stresses). However, plant requirements evapotranspiration (actual evapotranspiration) refers to actual conditions of a crop in the field considering all limitations. Therefore, it is not acceptable to compare these to variables. To this end, a crop coefficient (Kc) must be considered to characterize actual evapotranspiration (ETa); Eta = Kc × ETo.

Considering complex processes of evapotranspiration, there are many difficulties to measure this parameter. Some of the methods for this end are lysimeter, remote sensing, eddy covariance, and Bowen ratio. Allen et al. [35] determined an accurate model (FPM) which has the highest adaptability with lysimeter measurements in all over the world. Therefore, in this study, such as most of previous investigations [13,18,19,20], the authors compared the outputs of the GEP with the FPM and the results (Figure 1 and Figure 2) appear that the accuracy of the GEP has been improved compared to empirical models.

It should be noted that we cannot see similar results for the stations. In fact, all the results are completely independent of each other. This is in line with previous investigations [1,2].

It is notable that, in future research, the findings obtained from the GEP models in comparison to what is noted in the biological literature as affecting this process, should be checked to achieve a reliable result to recommend the GEP for other regions and research.

The GEP has worse performance over extreme conditions. The reason is related to the structure of GEP models. In this study, the authors focused on the models in which we need one or two weather parameters (mean temperature and mean wind speed) to estimate ETo. However, in extreme conditions, sometimes other variables, like maximum and minimum temperatures, maximum wind speed, rainfall, minimum humidity, solar radiation, and vapor pressure deficit, are also important. Therefore, if we are looking to increase the accuracy of GEP model we should consider more variables. However, in many regions we have no access to all of weather parameters to apply them to support GEP models to estimate ETo.

The next step of this study is to train the GA/GEP using ETo calculated at one site with all the meteorological variables available. Then, use the resulting functions to estimate ETo at a site where there are insufficient meteorological variables. Indeed, if we do have enough meteorological information to calculate ETo, we can then train the GA/GEP on fewer input variables.

5. Conclusions

In this study, the capability of both GA and GEP was assessed using 50-year time series data under 18 regions in Iran with arid, semiarid, very humid, and Mediterranean climates. The GA suggested that use of a double-parameter basis including the Tmean and WS may predict ETo with good accuracy in arid, semiarid, and Mediterranean regions. However, in very humid regions, temperature-based models (Tmean, Tmax, and Tmin) are better alternatives to reduce uncertainty. Finally, the GEP is recommended to evaluate the annual dynamics of ETo process in arid, semiarid, and Mediterranean climates with reliable accuracy under conditions where meteorological information is limited. In this study, we represented the best structures of GEP models. In future research, when a researcher or expert is going to use an accurate model based on this study, s/he may employ the best structures of GEP whose have been considered for his/her climate region.