Abstract

Reference evapotranspiration (ETo) is one of the most important factors in the hydrologic cycle and water balance studies. In this study, the performance of three simple and three wavelet hybrid models were compared to estimate ETo in three different climates in Iran, based on different combinations of input variables. It was found that the wavelet-artificial neural network was the best model, and multiple linear regression (MLR) was the worst model in most cases, although the performance of the models was related to the climate and the input variables used for modeling. Overall, it was found that all models had good accuracy in terms of estimating daily ETo. Also, it was found in this study that large numbers of decomposition levels via the wavelet transform had noticeable negative effects on the performance of the wavelet-based models, especially for the wavelet-adaptive network-based fuzzy inference system and wavelet-MLR, but in contrast, the type of db wavelet function did not have a detectable effect on the performance of the wavelet-based models.

INTRODUCTION

Reference evapotranspiration (ETo) plays an important role in hydrologic modeling and water balance studies. ETo is the evapotranspiration rate of a short green crop (grass) of uniform height completely shading the ground, and with adequate water status in the soil profile (Allen et al. 1998). ETo has a significant effect on irrigation water demands in arid and semi-arid regions, and many rainfall-runoff and ecosystem models use ETo as an important input variable. However, accurate estimation of ETo is challenging because ETo is a complex combination of evaporation from soil and transpiration from plants. The Food and Agriculture Organization (FAO) of the United Nations proposed a globally valid standard method (FAO56-PM) based on the Penman–Monteith equation (PM) for estimating ETo and which can be used in any climate and region (Allen et al. 1998). The major drawback of the FAO56-PM is that many input variables are required, including solar radiation, temperature, humidity and wind speed, and data records are often incomplete or unavailable for a given site. As a result, many researchers have examined or proposed different methods for estimating ETo based on limited data (Droogers & Allen 2002; Trajkovic 2005; Popova et al. 2006; Jabloun & Sahli 2008; Jain et al. 2008; Gocic & Trajkovic 2010; Martinez & Thepadia 2010; Shiri et al. 2014; Kisi et al. 2015).

Recently, there has been increased interest in using artificial intelligence (AI) methods for estimating ETo (Kisi 2007, 2008; Kumar et al. 2011; Gocic et al. 2015, 2016). Jain et al. (2008) examined an artificial neural network (ANN) for estimating ETo and showed that ANN was able to estimate ETo accurately, even with limited climate data. Kim & Kim (2008) found that generalized regression neural network models and a genetic algorithm could accurately estimate ETo. Landeras et al. (2009) used autoregressive integrated moving average (ARIMA) and ANN models for weekly prediction of ETo and found that both models worked well. Adeloye et al. (2012) used an unsupervised ANN, and a Kohonen self-organizing map (SOM) to develop models for ETo prediction and found that the SOM model estimated ETo with a good accuracy, even with limited climate data. Baba et al. (2013) showed that an adaptive network-based fuzzy inference system (ANFIS) and ANN were both able to accurately estimate ETo. El-Shafie et al. (2014) found that an ensemble neural network was a reliable tool in estimating ETo, and was better than traditional methods. Yassin et al. (2016) revealed that ANN models had a better performance than gene expression programming models for the estimation of daily ETo.

The wavelet transform (WT) is an advanced mathematical transform which has been used in recent years for meteorological and hydrological studies (Partal & Kisi 2007; Adamowski 2008; Adamowski & Sun 2010; Labat 2010; Shiri & Kisi 2010; Adamowski & Chan 2011; Sang 2013; Nourani et al. 2014; Araghi et al. 2015, 2017a, 2017b; Seo et al. 2015). Although there are many studies on ETo estimation using different kinds of AI or statistical methods, the number of studies using hybrid models based on WT and AI or statistical methods is relatively scarce. Partal (2009) used a combination of WT and ANN (called WANN) for modelling ETo and found that the combined model had a better performance than classical ANN approaches. Kisi (2011) showed that a hybrid model based on WT and linear regression (LR) had good accuracy for estimation of ETo. Falamarzi et al. (2014) revealed that the ANN and WANN could predict ETo with an acceptable accuracy, while the WANN had a better performance. Patil & Deka (2015) found that a combination of WT with some AI methods (i.e. ANN and ANFIS) could be employed as an accurate tool for estimation of daily ETo.

The main objective of this study was to compare ANN, ANFIS, multiple linear regression (MLR), and a combination of WT with each of these methods to estimate ETo in different climates and based on different input variables. There are some studies that estimate ETo based on wavelet-based models, but there is no comprehensive study comparing different wavelet-based models to estimate daily ETo in different climates; this is the primary contribution of the current study. Moreover, in this study, we examined the effect of wavelet function type and number of decomposition levels on the final results of the wavelet-based models, which has not yet been explored in the context of ETo studies.

THEORETICAL BACKGROUND

ETo estimation method

As briefly explained in the introduction above, we used the FAO56-PM for calculating ETo (Allen et al. 1998):  
formula
(1)
where ETo is reference evapotranspiration (mm day−1); Rn is net radiation at the crop surface (MJ m−2 day−1); G is soil heat flux density (MJ m−2 day−1) which is usually assumed to be zero; T is mean daily air temperature at a height of 2 m (°C); u2 is wind speed at a height of 2 m (m s−1); es is saturation vapor pressure (kPa); ea is actual vapor pressure (kPa); Δ is the slope vapor pressure curve (kPa °C−1) and γ is the psychrometric constant (kPa °C−1).
According to Equation (1), different input variables are required to calculate ETo. Saturation vapor pressure (es) at any temperature can be calculated using Equation (2). Also, the actual vapor pressure (ea) can be estimated using relative humidity (RH) (Allen et al. 1998):  
formula
(2)
 
formula
(3)
Using Equation (2), the slope of the saturation vapor pressure curve (Δ) can be calculated (Allen et al. 1998):  
formula
(4)
The psychrometric constant (γ) can be estimated using atmospheric pressure (P) (Allen et al. 1998):  
formula
(5)
 
formula
(6)
where P is atmospheric pressure (kPa); z is elevation above sea level (m); and γ is the psychrometric constant (kPa °C−1). If a record of P is available, it can be used instead of the values obtained from Equation (6).
Wind speed is usually measured at a height of 10 m at standard weather stations. Equation (5) can be used to convert wind speed at height z (uz) to u2 (Allen et al. 1998):  
formula
(7)
The net radiation (Rn) is the difference between the incoming net shortwave radiation (Rns) and the outgoing net longwave radiation (Rnl). Due to a lack of data for these types of radiation at many weather stations (usually because expensive measurement instruments are required), an estimation algorithm was suggested by Allen et al. (1998). It should be noted here that solar radiation (Rs) can be calculated with the Angstrom formula which relates Rs to extraterrestrial radiation (Ra) and relative sunshine hours (n) (Allen et al. 1998):  
formula
(8)
where Rs is solar shortwave radiation (MJ m−2 day−1); n is the actual sunshine hours; N is the maximum possible number of sunshine hours which is estimated using the geographical location of any site; Ra is extraterrestrial radiation (MJ m−2 day−1), and as and bs are Angstrom constants for any location which can assumed to be 0.25 and 0.5, respectively. More details about the FAO56-PM method can be found in Allen et al. (1998).

Artificial neural network

ANN, one of the well-known AI methods, can be expressed as a black-box mathematical model which transforms the input(s) into output(s) via an interconnected network of artificial neurons running in parallel. Multilayer perceptron (MLP) is one of the most frequently used ANN structures for modeling approaches in hydrometeorology. A typical MLP is made of one input, one output and one or more hidden layers. The proper number of hidden layers is chosen by the modeler based on the complexity of the problem. To produce an outcome in a neuron, each of the input values are multiplied by a weight, summing the product and then passing the sum through a transfer (or activation) function . Output of the jth neuron in a hidden layer can be written as (Yassin et al. 2016):  
formula
(9)
where , and are the input value(s), weight and bias parameters, respectively. The symmetric sigmoid (S-shaped) transfer function (Equation (8)) which was employed in this study is the most commonly used transfer function (Yassin et al. 2016):  
formula
(10)
where is the slope parameter. The outcome of the kth neuron in the output layer is typically calculated by a linear combination of results obtained from the neurons in the layer just before the output layer (Adamowski & Sun 2010; Falamarzi et al. 2014):  
formula
(11)
There are three main steps to develop an ANN model: training, calibration and validation. Training is the most important step in which the optimal values of the weights and biases are estimated. Comparing the final results of the model (outputs) and the observed values (targets) is the key to the training process. Mean square error (MSE) is the most frequently used index to evaluate the performance of the training step (Adamowski & Sun 2010):  
formula
(12)
where indicates the targets of modeling (i.e. the observation values). Different methods such as Levenberg–Marquardt (LM) and Bayesian regulation can be employed to perform the training step. In this study, the LM algorithm was used because it is rapid, accurate and reliable (Adamowski & Sun 2010). The LM method is a modification of the classic Newton algorithm for finding an optimum solution to a minimization problem. A schematic structure of an ANN is presented in Figure 1. Detailed information for MLP neural network models can be found in Hsieh & Tang (1998) and Adamowski & Sun (2010).
Figure 1

The schematic structure of a MLP ANN with one hidden layer (note: nI and nH are the number of neurons in the input and hidden layers, respectively).

Figure 1

The schematic structure of a MLP ANN with one hidden layer (note: nI and nH are the number of neurons in the input and hidden layers, respectively).

Adaptive network-based fuzzy inference system

ANFIS is the first fusion model based on ANN and the fuzzy inference system (FIS) in which the parameters of the FIS are estimated by ANN learning algorithms (Jang 1993). Generally, there are two main approaches for FIS, which are known as Mamdani (Mamdani & Assilian 1975) and Takagi-Sugeno (Takagi & Sugeno 1985). In Mamdani's approach, a fuzzy membership function (MF) is used in the consequent part (Q) of implication (i.e. if P, then Q), while Takagi-Sugeno's approach uses a linear or constant function, although both approaches can be interpreted in terms of fuzzy IF-THEN rules. The ANFIS structure is based on Takagi-Sugeno's approach and uses a hybrid learning rule combining back propagation gradient descent error digestion and a least-squares error method to determine a set of parameters (Nedjah & Mourelle 2005).

The ANFIS generally has five layers (Figure 2). The output of the ith node in layer 1 is the membership grade of inputs (x, y, etc.) using the MF (Jang 1993):  
formula
(13)
where is the MF for the Ai (Bi) fuzzy set. In this study, the Gaussian MF is employed, although there are various forms of MFs (Partal & Kisi 2007), such as triangular, trapezoidal, bell-shaped and Gaussian, which can be used to identify the membership grades. The Gaussian MF is expressed as (Jang 1993):  
formula
(14)
where is the parameter set for . Parameters in layer 1 are referred to as premise parameters. The output of the ith node in layer 2 is the product of incoming membership grades from layer 1 and it can be said that each node output in layer 2 denotes the firing strength of a rule (Jang 1993):  
formula
(15)
Figure 2

The schematic structure of the ANFIS.

Figure 2

The schematic structure of the ANFIS.

In layer 3, the ith node calculates the ratio of the ith rule's firing strength to the sum of all rules' firing strengths and therefore, the outputs of layer 3 are called normalized firing strengths. The output of the ith node in layer 3 can be expressed as (Jang 1993):  
formula
(16)
The output of the ith node in layer 4 , which represents the Takagi-Sugeno structure of the ANFIS, can be written as (Jang 1993):  
formula
(17)
where is the parameter set of this layer, which is composed of consequent parameters. In the last layer, the overall output of the ANFIS is acquired by summing all incoming signals (Jang 1993):  
formula
(18)
In ANFIS, the gradient descent method is used to determine the non-linear input parameters, while the linear output parameters are estimated using the least-squares method. More detailed information about ANFIS can be found in Jang (1993) and Nedjah & Mourelle (2005).

Multiple linear regression

MLR is one of the most widely used methods for modeling and prediction purposes. The MLR is defined as (Chatterjee & Hadi 2012):  
formula
(19)
where , , … , are regression coefficients (or parameters) which are unknown constants to be estimated from the predictor (, … , ) and response variables. Discrepancy in the approximation of Y is assumed to be a random error which is typically represented by . The MLR model is fitted such that the sum of squares of differences between observed and estimated values is minimized. In this study, nonlinear regression methods were not employed, as linear regression has been found to be a better approach to estimate ET (Kisi 2011). More details about MLR can be found in Chatterjee & Hadi (2012).

Wavelet transform

Fourier transform (FT) was one of the first methodologies to analyze frequencies in signals, but the technique has limitations, especially for signals with non-uniform periodic structures over time. Short-time Fourier transform (STFT) is a more appropriate method, in which the signal is broken down into several uniform sections using a window function. WT is an advanced revision of the STFT in which the window function is completely flexible and can be varied over time based on the compactness and shape of the signal (Hernández & Weiss 1996). Using WT, it is possible to analyze any signal with a uniform or non-uniform periodic structure (Mallat 2009). In general, there are two types of WT: continuous and discrete. The continuous WT (CWT) is defined as (Lau & Weng 1995):  
formula
(20)
 
formula
(21)
where is the WT of with a scale of s and a time shift of , represents the wavelet function and is the complex conjugate of . Flexibility of WT is founded on variations in s which is called the scale (Araghi et al. 2017a; Mallat 2009). In the discrete wavelet transform (DWT), the scale is dyadic (i.e. based on powers of two) and therefore, the calculation process can be reduced without any noticeable negative effect on the accuracy of the final results. DWT is written as follows (Daubechies 1990; Araghi et al. 2015):  
formula
(22)
where N is the length of . By applying the DWT to a signal, it decomposes the signal into two new ancillary signals called approximation (A) and detail (D) components, to express the low and high frequencies of the original signal, respectively. The decomposition process can be repeated at several levels, and at each level, component A is broken down into new A and D components. Following this, the summation of Ad and Dd components for decomposition level d is equal to the A component of decomposition level d−1 (i.e. Ad−1 = Ad + Dd). There are various types of wavelet functions (e.g. Haar, Meyer, Coiflets, Daubechies, and Symlets), but the Daubechies (db) is one of the most frequently used (Daubechies 1990; Araghi et al. 2015), and was employed in this study. Different types of db, called dbN are available, where N is the identification of type (e.g. db3). Increasing N adds complexity to the shape of the db wavelet function. More detailed information on DWT can be found in Daubechies (1990) and Mallat (2009).

STUDY SITES AND DATASETS

The second largest country in the Middle East, Iran, is located in southwest Asia between 25–40 °N latitude and 44–63 °E longitude, with a total area of 1,648,000 km2 and characterized mostly by arid and semi-arid climates (Araghi et al. 2018). The annual average precipitation in the country is around 250 mm, although local climatic regions may vary significantly, with annual precipitation levels ranging from 0 mm in central parts of Iran which is commonly desert (e.g. Lut and Central deserts), to 2,000 mm in coastal regions of the Caspian Sea in the north.

In this study, daily data from three weather stations in Iran with completely different climates were used. The variables including: minimum and maximum temperature (Tmin, Tmax), sunshine hours (S), wind speed (W), relative humidity (RH) and atmospheric pressure (P), were obtained during 2005–2009 from the three selected stations. The total number of data points for each variable was 1826 (i.e. 365 days/year for five years, plus 1 day for leap year 2008). Data from 2005 to 2008 were used for training and developing the models, and data from 2009 were used for testing the performance of models. Characteristics and locations of the selected weather stations are presented in Table 1 and Figure 3. The long-term monthly average of meteorological variables and ETo (calculated using FAO56-PM) at the three stations (Figure 4) indicate the distinctly different climatic status of each station.

Table 1

Identification and long-term climatic characteristics of the selected synoptic stations

Station name Longitude (°E) Latitude (°N) Elevation (m) Annual precipitation (mm) Annual temperature (°C) Climate type (De Martonne) 
Babolsar 52° 39′ 36° 43′ −21 892.6 16.8 Humid 
Mashhad 59° 38′ 36° 16′ 999.2 251.5 14.2 Semi-arid 
Yazd 54° 17′ 31° 54′ 1,237.2 57 19.5 Arid 
Station name Longitude (°E) Latitude (°N) Elevation (m) Annual precipitation (mm) Annual temperature (°C) Climate type (De Martonne) 
Babolsar 52° 39′ 36° 43′ −21 892.6 16.8 Humid 
Mashhad 59° 38′ 36° 16′ 999.2 251.5 14.2 Semi-arid 
Yazd 54° 17′ 31° 54′ 1,237.2 57 19.5 Arid 
Figure 3

Study region and the selected weather stations for this study (note: The climates of Babolsar, Mashhad and Yazd stations are humid, semi-arid and arid, respectively).

Figure 3

Study region and the selected weather stations for this study (note: The climates of Babolsar, Mashhad and Yazd stations are humid, semi-arid and arid, respectively).

Figure 4

Long-term monthly average of minimum temperature (a), maximum temperature (b), wind speed (c), relative humidity (d), sunshine hours (e), atmospheric pressure (f), precipitation (g), and reference evapotranspiration (h), at the studied stations during 1953–2014.

Figure 4

Long-term monthly average of minimum temperature (a), maximum temperature (b), wind speed (c), relative humidity (d), sunshine hours (e), atmospheric pressure (f), precipitation (g), and reference evapotranspiration (h), at the studied stations during 1953–2014.

METHODOLOGY

Model development

In this study we used three simple models, including ANN, ANFIS, and MLR, and three hybrid models, combining DWT and ANN, ANFIS and MLR, called WANN, WANFIS and WMLR, respectively. The hybrid models had the same structure as the simple models, but input variables to the hybrid models corresponded to the DWT of the original time series. For example, the original time series of Tmax in 2006 at the Mashhad station and its components decomposed via DWT using the db5 wavelet function and five decomposition levels are presented in Figure 5. According to our experience using WT in previous studies (Araghi et al. 2015, 2017a), and following the guidelines for selecting appropriate decomposition levels (de Artigas et al. 2006; Pradhan et al. 2006), we selected db5 as the wavelet function and five levels of decomposition were applied. Furthermore, db1, db5 and db10 wavelet functions and 1, 5 and 10 decomposition levels were examined for the DWT-based models with the highest accuracy at each of the stations, to see the effects of DWT settings on the final results of the models. It should be noted that an examination of DWT properties was not the main focus of this research; interested readers can refer to the following resources which have focused on optimal wavelet properties for hydrological and climatological applications in detail (Maheswaran & Khosa 2012; Sehgal et al. 2014a, 2014b; Sang et al. 2016; Yang et al. 2016; Peng & Wang 2017).

Figure 5

Original time series of daily maximum temperature in 2006 at the Mashhad weather station and its decomposed components via DWT using db5 wavelet function and five decomposition levels (note: The original time series is equal to the summation of D1 to D5 plus A5 components).

Figure 5

Original time series of daily maximum temperature in 2006 at the Mashhad weather station and its decomposed components via DWT using db5 wavelet function and five decomposition levels (note: The original time series is equal to the summation of D1 to D5 plus A5 components).

For each model, five different combinations of input variables were employed and each included Tmin and Tmax as it has been demonstrated that temperature has the highest correlation with ETo among all of the meteorological variables (Trajkovic 2005; Martinez & Thepadia 2010). Referring to daily temperature difference by using Tmin and Tmax serves as a surrogate for incoming solar radiation. The target values for each model were ETo values calculated via the FAO56-PM method.

Performance evaluation of models

To evaluate the performance of the models, four criteria were employed, including: the root mean square error (RMSE), the mean absolute error (MAE), the mean absolute percentage error (MAPE) and the coefficient of determination (R2). The mathematical expressions of these criteria are (Hyndman & Koehler 2006):  
formula
(23)
 
formula
(24)
 
formula
(25)
 
formula
(26)
 
formula
(27)
 
formula
(28)
where , and N are the model outputs, targets and length of data, respectively.

RESULTS AND DISCUSSION

Performance of the models

The performance of the models is presented in Tables 2 and 3, based on the different combinations of input variables. According to these tables, the RMSE values (mm day−1) at the Babolsar, Mashhad and Yazd stations range from 0.37 to 0.89, 0.55 to 0.95, and 0.39 to 1.15, respectively. The MAE values (mm day−1) for these stations vary from 0.31 to 0.73, 0.40 to 0.73, and 0.29 to 0.85, respectively, and the MAPE values fluctuate from 18.67 to 44.87%, 15.96 to 36.65%, and 8.11 to 23.11%, respectively.

Table 2

Model performance of ETO estimation based on non-wavelet-based models with different input variables (t: min and max temperature, w: wind speed, h: relative humidity, s: sunshine hours, p: atmospheric pressure) at stations with different climates

  RMSE (mm day−1MAE (mm day−1MAPE (%) R2 
Babolsar 
 ANNt 0.82 0.66 37.50 0.73 
 ANNt,w 0.67 0.49 21.54 0.84 
 ANNt,h 0.71 0.55 31.30 0.79 
 ANNt,s 0.45 0.35 21.87 0.92 
 ANNt,p 0.82 0.66 38.69 0.73 
 ANFISt 0.82 0.66 37.93 0.73 
 ANFISt,w 0.67 0.47 19.98 0.83 
 ANFISt,h 0.71 0.55 31.67 0.79 
 ANFISt,s 0.48 0.38 23.89 0.91 
 ANFISt,p 0.83 0.65 38.05 0.72 
 MLRt 0.89 0.73 43.98 0.68 
 MLRt,w 0.80 0.67 37.33 0.74 
 MLRt,h 0.75 0.62 37.27 0.76 
 MLRt,s 0.59 0.46 30.90 0.87 
 MLRt,p 0.89 0.73 44.87 0.69 
Mashhad 
 ANNt 0.79 0.60 22.16 0.91 
 ANNt,w 0.55 0.40 15.96 0.95 
 ANNt,h 0.74 0.54 19.19 0.91 
 ANNt,s 0.58 0.44 16.90 0.95 
 ANNt,p 0.76 0.57 21.72 0.90 
 ANFISt 0.77 0.57 20.88 0.90 
 ANFISt,w 0.57 0.42 16.40 0.95 
 ANFISt,h 0.74 0.54 19.14 0.91 
 ANFISt,s 0.59 0.44 17.08 0.94 
 ANFISt,p 0.77 0.57 20.25 0.90 
 MLRt 0.87 0.66 30.29 0.88 
 MLRt,w 0.86 0.70 34.63 0.89 
 MLRt,h 0.86 0.67 29.03 0.89 
 MLRt,s 0.76 0.59 29.89 0.90 
 MLRt,p 0.88 0.67 31.02 0.88 
Yazd 
 ANNt 1.01 0.77 19.24 0.85 
 ANNt,w 0.57 0.43 11.25 0.95 
 ANNt,h 1.02 0.77 18.75 0.85 
 ANNt,s 0.91 0.68 16.20 0.94 
 ANNt,p 0.94 0.73 17.97 0.87 
 ANFISt 1.02 0.77 18.83 0.85 
 ANFISt,w 0.52 0.40 10.72 0.96 
 ANFISt,h 1.01 0.77 18.38 0.85 
 ANFISt,s 0.93 0.71 17.05 0.88 
 ANFISt,p 0.92 0.71 17.63 0.88 
 MLRt 1.04 0.81 22.01 0.84 
 MLRt,w 0.68 0.56 16.59 0.94 
 MLRt,h 0.99 0.77 19.67 0.86 
 MLRt,s 0.97 0.76 20.09 0.86 
 MLRt,p 1.02 0.81 22.25 0.85 
  RMSE (mm day−1MAE (mm day−1MAPE (%) R2 
Babolsar 
 ANNt 0.82 0.66 37.50 0.73 
 ANNt,w 0.67 0.49 21.54 0.84 
 ANNt,h 0.71 0.55 31.30 0.79 
 ANNt,s 0.45 0.35 21.87 0.92 
 ANNt,p 0.82 0.66 38.69 0.73 
 ANFISt 0.82 0.66 37.93 0.73 
 ANFISt,w 0.67 0.47 19.98 0.83 
 ANFISt,h 0.71 0.55 31.67 0.79 
 ANFISt,s 0.48 0.38 23.89 0.91 
 ANFISt,p 0.83 0.65 38.05 0.72 
 MLRt 0.89 0.73 43.98 0.68 
 MLRt,w 0.80 0.67 37.33 0.74 
 MLRt,h 0.75 0.62 37.27 0.76 
 MLRt,s 0.59 0.46 30.90 0.87 
 MLRt,p 0.89 0.73 44.87 0.69 
Mashhad 
 ANNt 0.79 0.60 22.16 0.91 
 ANNt,w 0.55 0.40 15.96 0.95 
 ANNt,h 0.74 0.54 19.19 0.91 
 ANNt,s 0.58 0.44 16.90 0.95 
 ANNt,p 0.76 0.57 21.72 0.90 
 ANFISt 0.77 0.57 20.88 0.90 
 ANFISt,w 0.57 0.42 16.40 0.95 
 ANFISt,h 0.74 0.54 19.14 0.91 
 ANFISt,s 0.59 0.44 17.08 0.94 
 ANFISt,p 0.77 0.57 20.25 0.90 
 MLRt 0.87 0.66 30.29 0.88 
 MLRt,w 0.86 0.70 34.63 0.89 
 MLRt,h 0.86 0.67 29.03 0.89 
 MLRt,s 0.76 0.59 29.89 0.90 
 MLRt,p 0.88 0.67 31.02 0.88 
Yazd 
 ANNt 1.01 0.77 19.24 0.85 
 ANNt,w 0.57 0.43 11.25 0.95 
 ANNt,h 1.02 0.77 18.75 0.85 
 ANNt,s 0.91 0.68 16.20 0.94 
 ANNt,p 0.94 0.73 17.97 0.87 
 ANFISt 1.02 0.77 18.83 0.85 
 ANFISt,w 0.52 0.40 10.72 0.96 
 ANFISt,h 1.01 0.77 18.38 0.85 
 ANFISt,s 0.93 0.71 17.05 0.88 
 ANFISt,p 0.92 0.71 17.63 0.88 
 MLRt 1.04 0.81 22.01 0.84 
 MLRt,w 0.68 0.56 16.59 0.94 
 MLRt,h 0.99 0.77 19.67 0.86 
 MLRt,s 0.97 0.76 20.09 0.86 
 MLRt,p 1.02 0.81 22.25 0.85 
Table 3

Model performance of ETO estimation based on wavelet-based models with different input variables (t: min and max temperature, w: wind speed, h: relative humidity, s: sunshine hours, p: atmospheric pressure) at stations with different climates

Model RMSE (mm day−1MAE (mm day−1MAPE (%) R2 
Babolsar 
 WANNt 0.61 0.46 26.42 0.86 
 WANNt,w 0.67 0.51 30.01 0.83 
 WANNt,h 0.58 0.44 24.60 0.87 
 WANNt,s 0.37 0.31 18.67 0.96 
 WANNt,p 0.76 0.58 35.79 0.78 
 WANFISt 0.67 0.51 30.54 0.83 
 WANFISt,w 0.72 0.53 30.73 0.84 
 WANFISt,h 0.52 0.39 22.56 0.89 
 WANFISt,s 0.38 0.31 18.93 0.96 
 WANFISt,p 0.73 0.57 35.71 0.81 
 WMLRt 0.77 0.63 37.99 0.78 
 WMLRt,w 0.69 0.55 33.93 0.85 
 WMLRt,h 0.59 0.46 28.21 0.86 
 WMLRt,s 0.46 0.38 23.92 0.93 
 WMLRt,p 0.82 0.67 44.70 0.78 
Mashhad 
 WANNt 0.67 0.50 19.18 0.93 
 WANNt,w 0.56 0.45 18.76 0.97 
 WANNt,h 0.69 0.52 19.82 0.93 
 WANNt,s 0.62 0.48 20.32 0.94 
 WANNt,p 0.71 0.51 19.09 0.92 
 WANFISt 0.79 0.59 21.73 0.90 
 WANFISt,w 0.67 0.49 21.56 0.95 
 WANFISt,h 0.95 0.68 24.18 0.85 
 WANFISt,s 0.63 0.48 19.03 0.94 
 WANFISt,p 0.95 0.72 31.74 0.86 
 WMLRt 0.74 0.54 24.27 0.92 
 WMLRt,w 0.95 0.73 36.65 0.91 
 WMLRt,h 0.77 0.56 21.54 0.92 
 WMLRt,s 0.62 0.47 21.90 0.94 
 WMLRt,p 0.83 0.65 32.38 0.91 
Yazd 
 WANNt 0.95 0.72 18.37 0.87 
 WANNt,w 0.39 0.29 8.11 0.98 
 WANNt,h 1.06 0.82 21.42 0.85 
 WANNt,s 0.92 0.71 19.12 0.88 
 WANNt,p 0.92 0.71 18.11 0.88 
 WANFISt 0.99 0.76 18.55 0.86 
 WANFISt,w 0.44 0.32 8.71 0.98 
 WANFISt,h 1.06 0.82 21.20 0.84 
 WANFISt,s 1.15 0.85 23.11 0.82 
 WANFISt,p 0.89 0.68 16.57 0.89 
 WMLRt 1.01 0.78 21.91 0.85 
 WMLRt,w 0.56 0.42 13.25 0.96 
 WMLRt,h 0.95 0.73 19.25 0.87 
 WMLRt,s 0.91 0.70 17.84 0.88 
 WMLRt,p 0.98 0.78 21.30 0.86 
Model RMSE (mm day−1MAE (mm day−1MAPE (%) R2 
Babolsar 
 WANNt 0.61 0.46 26.42 0.86 
 WANNt,w 0.67 0.51 30.01 0.83 
 WANNt,h 0.58 0.44 24.60 0.87 
 WANNt,s 0.37 0.31 18.67 0.96 
 WANNt,p 0.76 0.58 35.79 0.78 
 WANFISt 0.67 0.51 30.54 0.83 
 WANFISt,w 0.72 0.53 30.73 0.84 
 WANFISt,h 0.52 0.39 22.56 0.89 
 WANFISt,s 0.38 0.31 18.93 0.96 
 WANFISt,p 0.73 0.57 35.71 0.81 
 WMLRt 0.77 0.63 37.99 0.78 
 WMLRt,w 0.69 0.55 33.93 0.85 
 WMLRt,h 0.59 0.46 28.21 0.86 
 WMLRt,s 0.46 0.38 23.92 0.93 
 WMLRt,p 0.82 0.67 44.70 0.78 
Mashhad 
 WANNt 0.67 0.50 19.18 0.93 
 WANNt,w 0.56 0.45 18.76 0.97 
 WANNt,h 0.69 0.52 19.82 0.93 
 WANNt,s 0.62 0.48 20.32 0.94 
 WANNt,p 0.71 0.51 19.09 0.92 
 WANFISt 0.79 0.59 21.73 0.90 
 WANFISt,w 0.67 0.49 21.56 0.95 
 WANFISt,h 0.95 0.68 24.18 0.85 
 WANFISt,s 0.63 0.48 19.03 0.94 
 WANFISt,p 0.95 0.72 31.74 0.86 
 WMLRt 0.74 0.54 24.27 0.92 
 WMLRt,w 0.95 0.73 36.65 0.91 
 WMLRt,h 0.77 0.56 21.54 0.92 
 WMLRt,s 0.62 0.47 21.90 0.94 
 WMLRt,p 0.83 0.65 32.38 0.91 
Yazd 
 WANNt 0.95 0.72 18.37 0.87 
 WANNt,w 0.39 0.29 8.11 0.98 
 WANNt,h 1.06 0.82 21.42 0.85 
 WANNt,s 0.92 0.71 19.12 0.88 
 WANNt,p 0.92 0.71 18.11 0.88 
 WANFISt 0.99 0.76 18.55 0.86 
 WANFISt,w 0.44 0.32 8.71 0.98 
 WANFISt,h 1.06 0.82 21.20 0.84 
 WANFISt,s 1.15 0.85 23.11 0.82 
 WANFISt,p 0.89 0.68 16.57 0.89 
 WMLRt 1.01 0.78 21.91 0.85 
 WMLRt,w 0.56 0.42 13.25 0.96 
 WMLRt,h 0.95 0.73 19.25 0.87 
 WMLRt,s 0.91 0.70 17.84 0.88 
 WMLRt,p 0.98 0.78 21.30 0.86 

According to the performance criteria, the most and least appropriate models for the Babolsar station were the WANN based on Tmin, Tmax and S, and the MLR based on Tmin, Tmax and P, respectively; for the Mashhad station these were the ANN based on Tmin, Tmax and W, and the WMLR based on Tmin, Tmax and W, respectively; and for the Yazd station these were the WANN based on Tmin, Tmax and W, and the WANFIS based on Tmin, Tmax and S, respectively. Applying the DWT approach did not enhance the performance of all models for all of the stations; model improvement depended on the input variables for each station. For instance, at the Babolsar station, applying DWT improved the model performance based on Tmin, Tmax and RH as input variables, while the performance of WANN and WANFIS models based on Tmin, Tmax and W was weak compared to that of the ANN and ANFIS models. These performances were completely reversed at the Yazd station. It can be concluded that the effect of applying DWT on a model's performance depends on the climate of the station and the input variables used. Overall, the performance of all models was acceptable with RMSE values of approximately 1 or less for all the studied stations.

In regions with limited climate data, the most available input variables were Tmin and Tmax. The performance of models based only on these input variables according to the RMSE, MAE and MAPE indices was from 0.61 to 1.02 mm day−1, 0.46 to 0.81 mm day−1 and 18.37 to 43.98%, respectively. These values represent a satisfactory performance for estimating ETo in the three climates studied, although in most cases, the WANN and MLR were the most accurate and least accurate, respectively. For example, a comparison of the performance of the MLR model based on Tmin and Tmax, and the WANN model based on Tmin, Tmax and W at the Yazd station in 2009 (Figure 6), indicates that the performance of the WANN model is much better than the MLR model. This figure is a good illustration of the effects of different input variables and the generally superior performance of wavelet-based models.

Figure 6

Observed and estimated ETO using the MLR model based on daily minimum and maximum temperatures (a) and the WANN model based on daily minimum and maximum temperatures and mean daily wind speed (b) in 2009 at the Yazd station.

Figure 6

Observed and estimated ETO using the MLR model based on daily minimum and maximum temperatures (a) and the WANN model based on daily minimum and maximum temperatures and mean daily wind speed (b) in 2009 at the Yazd station.

Effect of DWT settings on the final results

In the second part of the study, we examined the effects of DWT settings (type of db wavelet function and decomposition levels) on the final results of the wavelet-based models. To this end, the WANN, WANFIS and WMLR models with the best performance were selected for each station, based on the values represented in Table 3. All of these models were then run with db1, db5 and db10 wavelet functions and with 1, 5, and 10 decomposition levels. Results are presented in Figure 7. According to this figure, default settings of DWT, which were assumed for all of the wavelet-based models in this study, were appropriate settings. The most appropriate setting (i.e. the setting that yielded the lowest RMSE value) for each model and each station is identified with white text (Figure 7). The differences in RMSE between the ideal and the default settings are less than 0.5, although in four cases, db5 wavelet function and five decomposition levels represented the most accurate settings.

Figure 7

RMSE values for db1, db5 and db10 wavelet functions and 1, 5 and 10 decomposition levels in the study stations (notes: The best WANN, WANFIS and WMLR models were selected for each station, based on the values reported in Tables 2 and 3. The lowest RMSE value is identified with white text).

Figure 7

RMSE values for db1, db5 and db10 wavelet functions and 1, 5 and 10 decomposition levels in the study stations (notes: The best WANN, WANFIS and WMLR models were selected for each station, based on the values reported in Tables 2 and 3. The lowest RMSE value is identified with white text).

These results indicate that different types of db wavelet functions do not have a detectable effect on the final results of the models. However, decomposition levels equal to 10 can have a noticeable negative influence on the final results of the wavelet-based models, especially WANFIS and WMLR, and particularly when db10 was employed as the wavelet function. It can be concluded that medium decomposition levels could be the best choice for most cases.

CONCLUSIONS

The aim of this study was to compare the performance of three simple and three wavelet-based AI models for estimation of daily evapotranspiration in three different climates in Iran, based on different combinations of input variables. Results showed that, in most cases, the WANN and MLR models were more and less effective, respectively. Ultimately, the WANN model outperformed the other wavelet hybrid models (i.e. WANFIS and WMLR), which confirms the results of previous studies (Partal 2009, 2015; Falamarzi et al. 2014; Patil & Deka 2015). It was found that applying DWT generally had a positive effect on performance, but not in all three of the climates studied, nor based on all combinations of input variables. The other finding of this study was related to the effects of DWT settings on the final results of wavelet-based models. We found that the type of db wavelet function did not have a detectable effect on the performance of the models, but that applying a large number of decomposition levels had a noticeable negative influence on the performance, particularly for WANFIS and WMLR models. The combination of a medium decomposition level and a middle type of db wavelet function was found to be most appropriate in this study. This suggests the need for more studies with a similar approach using other meteorological and hydrological variables to find the most appropriate structure of wavelet-based models, and to analyze the effects of applying DWT to the models.

REFERENCES

REFERENCES
Adeloye
A. J.
,
Rustum
R.
&
Kariyama
I. D.
2012
Neural computing modeling of the reference crop evapotranspiration
.
Environ. Model. Softw.
29
,
61
73
.
Allen
R. G.
,
Pereira
L. S.
,
Raes
D.
&
Smith
M.
1998
Crop Evapotranspiration: Guidelines for Computing Crop Water Requirements
.
FAO
,
Rome
,
174
pp.
Araghi
A.
,
Mousavi-Baygi
M.
,
Adamowski
J.
,
Malard
J.
,
Nalley
D.
&
Hasheminia
S. M.
2015
Using wavelet transforms to estimate surface temperature trends and dominant periodicities in Iran based on gridded reanalysis data
.
Atmos. Res.
155
,
52
72
.
Araghi
A.
,
Mousavi-Baygi
M.
,
Adamowski
J.
&
Martinez
C.
2017a
Association between three prominent climatic teleconnections and precipitation in Iran using wavelet coherence
.
Int. J. Climatol.
37
,
2809
2830
.
Araghi
A.
,
Mousavi-Baygi
M.
,
Adamowski
J.
,
Martinez
C.
&
van der Ploeg
M.
2017b
Forecasting soil temperature based on surface air temperature using a wavelet artificial neural network
.
Meteorol. Appl.
24
,
603
611
.
Araghi
A.
,
Martinez
C. J.
,
Adamowski
J.
&
Olesen
J. E.
2018
Spatiotemporal variations of aridity in Iran using high - resolution gridded data
.
Int. J. Climatol.
38
,
2701
2717
.
Chatterjee
S.
&
Hadi
A. S.
2012
Regression Analysis by Example
,
5th edn
.
John Wiley & Sons Inc.
,
Hoboken, NJ
,
393
pp.
de Artigas
M. Z.
,
Elias
A. G.
&
de Campra
P. F.
2006
Discrete wavelet analysis to assess long-term trends in geomagnetic activity
.
Phys. Chem. Earth
31
,
77
80
.
El-Shafie
A.
,
Najah
A.
,
Alsulami
H. M.
&
Jahanbani
H.
2014
Optimized neural network prediction model for potential evapotranspiration utilizing ensemble procedure
.
Water Resour. Manage.
28
,
947
967
.
Gocic
M.
,
Motamedi
S.
,
Shamshirband
S.
,
Petković
D.
,
Ch
S.
,
Hashim
R.
&
Arif
M.
2015
Soft computing approaches for forecasting reference evapotranspiration
.
Comput. Electron. Agric.
113
,
164
173
.
Gocic
M.
,
Petković
D.
,
Shamshirband
S.
&
Kamsin
A.
2016
Comparative analysis of reference evapotranspiration equations modelling by extreme learning machine
.
Comput. Electron. Agric.
127
,
56
63
.
Hernández
E.
&
Weiss
G.
1996
A First Course on Wavelets
, 1st edn.
CRC Press
,
Boca Raton, FL
,
489
pp.
Hyndman
R. J.
&
Koehler
A. B.
2006
Another look at measures of forecast accuracy
.
Int. J. Forecasting.
22
,
678
688
.
Jang
J. S. R.
1993
ANFIS: adaptive-network-based fuzzy inference system
.
IEEE Trans. Syst. Man Cybern.
23
,
665
685
.
Kisi
O.
2011
Evapotranspiration modeling using a wavelet regression model
.
Irrig. Sci.
29
,
241
252
.
Kisi
O.
,
Sanikhani
H.
,
Zounemat-Kermani
M.
&
Niazi
F.
2015
Long-term monthly evapotranspiration modeling by several data-driven methods without climatic data
.
Comput. Electron. Agric.
115
,
66
77
.
Kumar
M.
,
Raghuwanshi
N. S.
&
Singh
R.
2011
Artificial neural networks approach in evapotranspiration modeling: a review
.
Irrig. Sci.
29
,
11
25
.
Landeras
G.
,
Ortiz-Barredo
A.
&
López
J. J.
2009
Forecasting weekly evapotranspiration with ARIMA and artificial neural network models
.
J. Irrig. Drain. Eng.
135
,
323
334
.
Mallat
S.
2009
A Wavelet Tour of Signal Processing
,
3rd edn
.
Academic Press
,
Burlington, MA
,
805
pp.
Mamdani
E. H.
&
Assilian
S.
1975
An experiment in linguistic synthesis with a fuzzy logic controller
.
Int. J. Man Mach. Stud.
7
,
1
13
.
Martinez
C. J.
&
Thepadia
M.
2010
Estimating reference evapotranspiration with minimum data in Florida
.
J. Irrig. Drain. Eng.
136
,
494
501
.
Nedjah
N.
&
Mourelle
L. d. M.
2005
Fuzzy Systems Engineering, Theory and Practice
.
Springer
,
Berlin
,
233
pp.
Nourani
V.
,
Hosseini Baghanam
A.
,
Adamowski
J.
&
Kisi
O.
2014
Applications of hybrid wavelet–artificial intelligence models in hydrology: a review
.
J. Hydrol.
514
,
358
377
.
Pradhan
P. S.
,
King
R. L.
,
Younan
N. H.
&
Holcomb
D. W.
2006
Estimation of the number of decomposition levels for a wavelet-based multiresolution multisensor image fusion
.
IEEE Trans. Geosci. Remote Sens.
44
,
3674
3686
.
Sang
Y.-F.
,
Singh
V. P.
,
Sun
F.
,
Chen
Y.
,
Liu
Y.
&
Yang
M.
2016
Wavelet-based hydrological time series forecasting
,
J. Hydrol. Eng.
21
,
06016001
.
Shiri
J.
,
Nazemi
A. H.
,
Sadraddini
A. A.
,
Landeras
G.
,
Kisi
O.
,
Fakheri Fard
A.
&
Marti
P.
2014
Comparison of heuristic and empirical approaches for estimating reference evapotranspiration from limited inputs in Iran
.
Comput. Electron. Agric.
108
,
230
241
.
Takagi
T.
&
Sugeno
M.
1985
Fuzzy identification of systems and its applications to modeling and control
.
IEEE Trans. Syst. Man Cybern. SMC
15
,
116
132
.