Multivariate modeling of river water withdrawal using a hybrid evolutionary data-driven method

Climate change and restrictions on water resources' accessibility affect the lives of human beings and all living organisms; and they pose serious challenges for humanity in different aspects (Jeihouni et al. 2019; Vanderhoof et al. 2019; Bissenbayeva et al. 2021). Various sectors of agriculture, industry, environment, etc. depend on water resources to meet the needs of communities. The agricultural sector, in addition to having a significant share in satisfying food needs, is the main sector of water consumption due to irrigation demand to produce agricultural products. Therefore, optimal management and planning in this sector can have a significant impact on saving more water resources for future generations (Vondracek et al. 2005; Mehta et al. 2013, 2014; Böhme et al. 2016; Ghalehkhondabi et al. 2017; Patel et al. 2018; Langat et al. 2019; Kumar & Mehta 2020). Low flows and droughts are also important for agricultural water use (Eris et al. 2019; Eris et al. 2020).

The amount of River Water Withdrawal (RWW) for agricultural consumption can be evaluated from various perspectives, including morphological, quantitative and qualitative hydrological, and land use factors (Msigwa et al. 2019; Magritsky et al. 2020). Due to the dynamic nature of rivers, the morphological factor is constantly changing and these changes can negatively affect the quantitative and qualitative status of rivers (Gostner et al. 2013) and agricultural lands (Hohensinner et al. 2018). Population growth, accompanied by increasing land use changes, and improper uses of surface and groundwater resources have not only reduced water resource quantities, but have also destroyed the quality of these valuable resources. Investigating the impact of land use on river water withdrawal can enhance the management and planning of water resources, which requires both a holistic view and engineering precision (Shank & Stauffer 2015; Msigwa et al. 2019; Shirmohammadi et al. 2020). The arid and semi-arid climatic conditions in the majority of basins in Iran highlights the limitations of water resources and the importance of proper management (Abbaspour et al. 2009; Sharafati et al. 2020).

In previous decades, the importance of surface water resources has led to the evaluation of various quantitative and qualitative hydrological variables of river flow rate, Total Dissolved Solids (TDS), and water temperature, etc. by various numerical methods (Mehta et al. 2013, 2014; Montaseri et al. 2018; Mehta & Yadav 2020; Mehta et al. 2020). Recently, hybrid data-driven evolutionary methods such as Wavelet–AutoRegressive (WAR) model, Wavelet–AutoRegressive Moving Average (WARMA), Wavelet–Linear Regression (WLR), Wavelet–Artificial Neural Networks (WANNs), Wavelet–Particle Swarm Optimization (WPSO), Wavelet–Adaptive Neural Fuzzy Inference System (WANFIS), Wavelet–Support Vector Machine (WSVM), and Wavelet–Gene Expression Programming (WGEP) have been extended and emphasized by many researchers and politicians for modeling, optimization, and management of the world's water resources (Yarar 2014; Barzegar et al. 2016; Montaseri et al. 2018; Zhang et al. 2018). In the estimation and modeling field of quantitative and qualitative characteristics of river water, some summaries of research are addressed. Adamowski & Sun (2010) extended ANN and WANN methods for flow estimating of the Kargotis and Xeros rivers, located in Cyprus with semi-arid climate conditions. In both rivers, the WANN models provide more real flow estimates than the ANN models. Barzegar et al. (2016) applied single and Wavelet–hybrid (W-hybrid) data-driven methods, including ANN, ANFIS, WANN, and WANFIS, to estimate the Electrical Conductivity (EC) value for Aji-Chay River, located in a cold semi-arid region in Iran. Respectively, the results reflected the superiority of Daubechies-4 (Db4) mother wavelet decomposition compared with the other wavelets, W-hybrid models compared with the single models, and the WANFIS model compared with the WANN model. The ability of evolutionary data-driven methods (e.g. GEP, Support Vector Machine (SVM), WGEP, and WSVM) has been investigated by Solgi et al. (2017) for estimating the daily and monthly flow rates of Gamasiyab River, located within a cold semi-arid basin in Iran. The results reflected the higher efficiency of W-hybrid methods compared with the single ones. Shafaei & Kisi (2017) predicted the daily flow of the Aji-Chay River, located in a cold semi-arid basin in Iran, applying SVM, ANN, and WANN models. The results reflected the superior performance of the WANN methods. Zhang et al. (2018) estimated four-station streamflow of the East River basin, located in a subtropical basin in China, through Multiple Linear Regression (MLR), ANN, and WANN methods. Based on the results, the WANN model performed well in comparison with the MLR and ANN models. Yaseen et al. (2018) estimated monthly Tigris river flow in Iraq, employing simple and W-hybrid data-driven evolution, as a new structure of ANN, which is called Extreme Learning Machine (ELM and WELM) methods. The results revealed that WELM models can be introduced as a reliable application for estimating river flow in semi-arid climatic conditions. Montaseri et al. (2018) used single and W-hybrid data-driven methods, including ANN, ANFIS, GEP, WANN, WANFIS, and WGEP, to estimate the amount of TDS in rivers of four basins with various climatic conditions (e.g. Dsa, Bsk, Bwk, and Bsh Köppen–Geiger climate categories), located in Iran. The results indicated the superior efficiency of W-hybrid methods compared with single methods. It also highlighted the importance of providing mathematical relationships derived from GEP and WGEP. Sun et al. (2019) predicted short-term flow rates of the Heihe and Pearl rivers in China, classified in arid-semi-arid and humid subtropical climate categories, in turn, using evolutionary data-driven methods (e.g. AR, ARMA, ANN, LR, WAR, WARMA, WANN, and WLR). Based on the results, they declared that W-hybrid models have higher performance than single models. Kumar et al. (2020) extended data-driven evolutionary methods including WANN and SVM to estimate the rate of perennial river discharge, located in India. The results indicated the appropriate and acceptable performance of evolutionary data-driven methods to estimate daily river discharge. Wang et al. (2020) modeled three water quality indicators of the Grand Canal, located in China, including Chemical Oxygen Demand (COD), Ammonia Nitrogen (NH₃-N), and Dissolved Oxygen (DO), using a single and hybrid form of W-PSO-SVR methods. The results indicated better performance of coupled wavelet models for estimating water quality indicators in tropical monsoon climatic conditions. Chen et al. (2020) and Rajaee et al. (2020) investigated single and W-hybrid data-driven methods to estimate river water qualitative variables. Their evaluations emphasize the applicability of the data-driven methods in the estimation of water quality indicators for rivers and the superiority of W-hybrid methods compared with the single data-driven methods for various climatic conditions. Kambalimath & Deka (2021) estimated daily stream flow in Malaprabha sub-catchment located in India via SVM and W-SVM models. Their results indicated that combining wavelet models with single SVM could improve the efficiency of daily stream flow estimating. Since in recent years, due to climate change, the Namak Lake basin has been facing a shortage of water resources, it is especially important to pay attention to the amount of water withdrawal for agricultural uses (Sheikh et al. 2020).

Despite the great improvement in the data-driven methods' application in hydrological-variable estimation, most of the previous studies focused on modeling the quantitative and qualitative characteristics of river flows and evaluating the elementary level of the hybrid models without estimating RWW, which is covered by the current research. The Namak Lake basin has an arid climate condition. As agriculture is affected by climate change and is the main consumer of water resources, scientists and policymakers should pay special attention to satisfying this sector to achieve sustainability. Consequently, estimating the amount of river water withdrawal for agricultural uses requires extreme attention. This study is conducted to assess three single data-driven models (ANN, ANFIS, and GEP) and their W-hybrid models (WANN, WANFIS, and WGEP) for RWW estimation in two study areas classified in various climate categories, while data de-noising, and evaluation of various factors (e. g. morphological, quantitative and qualitative hydrological, land use, and combined factors) is carried out. The novelty of the developed W-hybrid model's structure and estimated RWW variables with various factors and related variables, by collecting data with field sampling in the Lavasanat and Qazvin study areas, make current research distinguished from the previous literature. Besides, considering various climate zones provides a better comparison of the model's performance.

The Namak Lake basin with an area of 92,563 km² is located in the central part of Iran, which is one of the most critical basins in Iran. In the current research, the Lavasanat and Qazvin study areas were selected to evaluate river water withdrawal for agricultural uses (Figure 1). According to the Köppen–Geiger climate classification, the Lavasanat study area is categorized as arid desert cold (Bsk), and the Qazvin study area as a warm temperate with dry, hot summer (Csa) climate type. The Bsk is a climate in which average annual temperature is less than 18 °C and is too dry to support a forest, but not dry enough to be a desert, usually consisting of grassland plains. The Csa is a climate in which the coldest month is warmer than −3 °C but colder than +18 °C and summers are dry and hot. In selected study areas, field sampling was performed to determine the amount of river water withdrawal for agricultural uses, and data was collected by the Iran Water Resources Management Company (http://wrbs.wrm.ir/). The number and spatial distribution of samples are acceptable based on the significant results of the F-test at the level of 0.01 (sig = 0.000), applying SPSS software (IBM Corp., released 2017). The number of samples in the Lavasanat and Qazvin study areas is 291 and 198, respectively. The sampled variables included River Width (RW), River Depth (RD), minimum flow rate (Q_Min), maximum flow rate (Q_Max), average flow rate (Q_Mean), Water Temperature (WT), Electrical Conductivity (EC), pH, Cultivated Area (CA), and Orchard Area (OA).

RW and RD represent the morphological features, Q_Min, Q_Max, Q_Mean, WT, EC, and pH represent quantitative and qualitative hydrological characteristics, CA and OA signify land use features, and RW, Q_Min, Q_Mean, EC, and CA are regarded as combined factors. The statistical characteristics of the data for the selected study areas are given in Table 1. The RWW variable has average values of 182.08 and 154.59 (×10³ m³) in the Lavasanat and Qazvin study areas, respectively.

For more details see Mallat (1989).

Artificial neural networks (ANNs) are computational tools with a similar function to biological brain processes. The ANNs can be defined as a network of simple processors (or neurons) with three main layers including input, hidden, and output layers (Ferreira 2001). The mathematical form of neural networks is much simpler than their biological appearance (Rahmanpanah et al. 2020). In ANN methodology, the dataset is often separated into two main sections: training and testing. The training and testing sections contain the data for the model's learning to determine the weight of neurons and evaluate generalization performance, respectively. The training section should contain appropriate input data and ultimately the desired results. The selection of input variables is one of the main aspects of building a successful neural model in which keeping the cost of data collection down, eliminating the effects of duplicate data, and providing the model in the simplest possible way are prioritized. In modeling ANNs, determining input variables and collecting appropriate training data require more time and effort than network training. The special purpose of ANN modeling is to introduce the interaction between input and output variables (Montaseri et al. 2018; Chaplot & Birbal 2021). The Levenberg–Marquardt (LM) algorithm is one of the fastest ANN training algorithms, the structure of which is designed to achieve high training speed without the need to calculate the Hussein Matrix (HM) (Ghavidel & Montaseri 2014).

In our research, an ANN, containing a three-layer structure (including: input, hidden, and output layers), LM Back-Propagation Training Algorithm (BPTA), various transfer functions in the hidden and output layers, and a varied number of neurons (1–10) in the hidden layer, has been coded (in MATLAB software, 2018) and developed to model the amount of RWW for agricultural uses.

In the current study, the ANFIS model based on the Subtractive Clustering (SC) method called ‘ANFIS-SC’ with a varied radii value (0–1) has been coded (in MATLAB software, 2018) and developed to model the amount of RWW for agricultural uses.

The GEP, similar to the genetic algorithm (GA), can combine linear and simple chromosomes of constant length. GEP models are generated based on the Darwinian theory of natural selection (Ferreira 2001; Birbal et al. 2021). The most effective indicators were extracted by the GEP model formulations. The GEP modeling of RWW was based on several indices and parameters, namely the Functions set (F), Terminal set (T), Mutation Rate (MR), Inversion Rate (IR), IS Transposition Rate (ISTR), RIS Transposition Rate (RISTR), One-Point Recombination Rate (OPRR), Two-Point Recombination Rate (TPRR), Gene Recombination Rate (GRR), Gene Transposition Rate (GTR), the Linking Function (LF), Fitness Function Error Type (FFET), and Penalizing Tool (PT). The GeneXpro Tools 4.0 default values which are used in this for MR, IR, ISTR, RISTR, OPRR, TPRR, GRR, and GTR are 0.044, 0.1, 0.1, 0.1, 0.3, 0.3, 0.1, and 0.1, respectively. The sub-trees were linked by an addition function. The FFET was RRSE. The parsimony pressure was selected as PT. The number of chromosomes, the number of genes per chromosome, and the head size of the chromosome are 30, 3, and 7, respectively. More details about the technical formulation of the GEP approach can be found in Ferreira (2006).

The unique ability of the wavelet method for estimating variables has led to the development of wavelet hybrids with other methods. In this study, the wavelet was hybridized with all three methods, ANN, GEP, and ANFIS, to estimate RWW for agricultural uses, and then, the estimated results were compared. Figure 2 shows the used single and W-hybrid approach (ANN, ANFIS, GEP, WANN, WANFIS, and WGEP) structures, and Figure 3 displays the stages of the current research.

The data was decomposed and divided into main (A_l) and detail (D_l) sub series by applying one-dimensional Daubechies-4 (db4) mother wavelet analysis, which has been provided by many researchers, recently Shafaei & Kisi (2017) and Sun et al. (2019). The wavelet decomposed the variables of morphological, quantitative and qualitative hydrological, and land use factors including RW, RD, Q_Min, Q_Max, Q_Mean, WT, EC, pH, CA, and OA, in the Lavasanat and Qazvin study areas at level 2 and 1, respectively. The obtained results of wavelet analysis are illustrated in Table 2. The A_l (A₁ and A₂) series denote the main decomposed data and the D_l (D₁ and D₂) series denote the detailed decomposed or data noise. For instance, A₂ and A₁ with low-frequency approximation section at level 2 and 1 of the RWW variables for Lavasanat and Qazvin varied from −32.98 to +597.63 and −81.57 to +903.82, respectively.

In the current study, the three single data-driven methods ANN_ij, ANFIS_ij, and GEP_ij (i = 1 and 2 show the number of selected study areas, and j = 1,.., 5 show the number of selected effective factors on the RWW variable) are applied for estimating the RWW variable of Lavasanat and Qazvin. The ANN–Levenberg–Marquardt (LM-BP) algorithm with one hidden layer is applied for RWW modeling. The hidden nodes' numbers of ANN (1–10) and radii values of the ANFIS (0.10–0.80) models are examined by the trial-and-error method. The results of optimal single models based on R, RMSE, and NSE values and characteristics for ANN (the number of neurons and activation functions of hidden and output layers), ANFIS (radii values), and GEP models are presented in Tables 3 and 4 for the Lavasanat and Qazvin study areas, respectively. The transfer functions of the output layer are evaluated purelin or tansig or logsig for the ANN developed models of various factors in both study areas with Bsk and Csa climates. For example, the transfer functions of the hidden layers of the ANN₁₁ to ANN₁₅ models were obtained tansig-purelin, tansig-tansig, tansig-purelin, logsig-purelin, and tansig-tansig for the Lavasanat study area.

The values of R, RMSE, and NSE for morphological, quantitative and qualitative hydrological, land use, and combined factors in the Lavasanat study area with Bsk climate classification corresponding to the ANN₁₁ to ANN₁₅ models varied from 0.702 to 0.940, 53.549 to 125.508, and 0.335 to 0.879 respectively. These figures for the Qazvin study area with Csa climate classification corresponding to the ANFIS₂₁ to ANFIS₂₅ models ranged from 0.691 to 0.955, 80.644 to 200.817, and 0.436 to 0.909, respectively. In the current study, the Root Relative Squared Error (RRSE) with a pressure tree is selected as a fitness function for the GEP models. The R and RMSE for the GEP₁₁ to GEP₁₅ and GEP₂₁ to GEP₂₅ models of RWW estimations in the Lavasanat and Qazvin study areas have the values of [R (0.793, 0.928, 0.723, 0.899, 0.983), and RMSE (×10³ m³) (105.024, 64.525, 123.352, 69.824, 33.812)] and [R (0.843, 0.905, 0.813, 0.706, 0.966), and RMSE (×10³ m³) (153.283, 116.553, 162.991, 189.911, 78.017)], respectively. The results of R, RMSE, and NSE values for ANN, ANFIS, and GEP models to estimate the RWW variable by all aforementioned factors for the Lavasanat and Qazvin study areas are summarized in Tables 3 and 4 for the testing section.

The W-hybrid (WANN_ij, WANFIS_ij, and WGEP_ij) data-driven models are developed via wavelet tools for improving the performance of the single (ANN_ij, ANFIS_ij, and GEP_ij) data-driven models for estimating RWW in the Lavasanat and Qazvin study areas. In the current method, the first stage is to divide and decompose the used variables of morphological, quantitative and qualitative hydrological, and land use factors into the subseries of main and detail (A_l and D_l series) by applying a Wavelet-db4 tool. To make and develop the W-hybrid RWW estimated models, D_{l=1, 2} decomposed subseries are introduced as noise and removed from the models. After de-noising, decomposed subseries values are estimated separately with single models. The results of optimal W-hybrid models, WANN, WANFIS, and WGEP, by R, RMSE, and NSE values and specifications of models for various factors including morphological, quantitative and qualitative hydrological, land use, and combined, for the Lavasanat and Qazvin study areas with Bsk and Csa climate type are listed in Tables 3 and 4 for the testing section, respectively. The values of R and RMSE for estimating RWW corresponding to the WANN₁₁ to WANN₁₅ and WANN₂₁ to WANN₂₅ models varied from 0.773 to 0.986 and 20.615 to 93.470 in the Lavasanat, 0.670 to 0.963 and 62.734 to 176.643 in the Qazvin study areas, respectively. On the other hand, the values of R and RMSE corresponding to the WANFIS₁₁ to WANFIS₁₅ and WANFIS₂₁ to WANFIS₂₅ models for RWW ranged from 0.799 to 0.991 and 16.497 to 88.150 in the Lavasanat, 0.719 to 0.979 and 46.565 to 175.510 in the Qazvin study areas. Finally, the values of R and RMSE for the WGEP₁₁ to WGEP₁₅ and WGEP₂₁ to WGEP₂₅ models of RWW estimates in the Lavasanat and Qazvin study areas figure at [R (0.821, 0.932, 0.811, 0.919, 0.996) and RMSE (×10³ m³) (67.268, 54.659, 80.871, 50.796, 15.676)] and [R (0.898, 0.914, 0.885, 0.723, 0.990) and RMSE (×10³ m³) (105.532, 96.615, 105.018, 160.961, 44.332)], respectively.

Figures 4–6 show the observed and estimated RWW values for the single models and their W-hybrid models based on various factors for the study areas in the testing section. The R values of the developed data-driven models are close to 1, while the relation among R values for all applied factors of RWW estimation models in the study areas is as follows: R_GEP> R_ANFIS > R_ANN and R_WGEP> R_WANFIS > R_WANN. For instance, the percentage of WGEP performance improvement compared with GEP for morphological, quantitative and qualitative hydrological, land use, and combined factors are 35.95%, 15.29%, 34.44%, 27.25%, and 53.64% in the Lavasanat and 31.15%, 17.11%, 35.57%, 15.24%, and 43.18% in the Qazvin study areas.

The results comparison of the applied models for various factors in the Lavasanat and Qazvin study areas are shown in Figures 7 and 8, respectively. The results establish that the ANFIS model exceeds the ANN model's performance; the WANFIS model surpasses the WANNs; and also, the GEP and WGEP models have a better performance for RWW estimating than the two other single and W-hybrid applied models for both study areas with various factors and climates for the testing section.

The Taylor diagram reflects the degree quantification between the observations and estimations in terms of the R, RMSE, and the standard deviation (Taylor 2005). The comparison of the WGEP model's performance for estimating RWW with various factors in the Lavasanat and Qazvin study areas is shown in Figure 9. Based on the reflected results of Figure 9, the combined factor, including the RW, Q_Min, Q_Mean, EC, and CA variables, has the best performance in estimating RWW values compared with the other developed factors in both study areas.

One of the most popular charts that show many descriptive statistics of estimated RWW data is the box-plot, which is based on the five values of ‘minimum’, ‘first quartile (0.25%)’, ‘median (0.50%)’, ‘third quartile (0.75%)’ and ‘maximum’ values. Furthermore, depicting symmetry in the data is one of the tasks of this chart. It is worth noting that the degree of focus and even skewness of the data can be seen in this chart. Figure 10 displays the boxplots of estimated RWW values based on the combined factor for Bsk and Csa climate type, respectively. On account of the chart's results, the estimated RWW values of the Bsk climate category have the maximum values of the median and the lower variation coefficient of estimated data compared with the Csa climate category.

Figure 11 shows the confidence and prediction bands at 95% of the WGEP models (combined factor) in the Lavasanat and Qazvin study areas. The finding of the aforementioned intervals showed the reliability of the WGEP model's prediction with the combined factor for estimating the RWW variable. The prediction band indicates uncertainty in the true position of the curve (enclosed by the confidence bands), and also accounts for the scatter of data around the curve, while the confidence interval indicates how well the research data defines the best-fit curve. Estimating uncertainty enables water managers to analyze and discover a wide range of sustainable management practices and identify the ones that are most robust for all factors (Adib et al. 2019). On the other hand, for more details about the input variables' effect on RWW estimation, Tornado sensitivity analysis of the WGEP models (combined factor) at 5% and 95% is plotted in Figure 11. The highest and lowest sensitivity of the RWW estimation models at 95% are related to the (Q_Mean, RW) and (Q_Mean, EC) input variables in the Lavasanat and Qazvin study areas, respectively.

The unique structure of each data-driven method is the main reason for the models' efficiency difference. Gene and chromosome production in the GEP method makes the efficiency of this method considerably better than the other models. Also, the ANFIS method, by combining fuzzy rules and neural neuron structure, has relatively high efficiency compared with the ANN method. The efficiency of single models is significantly increased in all data-driven methods and climatic categories by combining wavelet theory and de-noising of data and creating complex nonlinear relationships in its structure. On the other hand, climatic categories only affect the models' performance values, not their priority. However, the results of single and W-hybrid methods could be acceptable for estimating the RWW variable in the Lavasanat and Qazvin study areas. Table 5 summarized the performance of applied factors to estimate RWW by scoring the studied areas in the Bsk and Csa climate categories for the testing section. The results indicated that the WGEP and ANN models respectively are the best and poorest models in both study areas without the effect of climate conditions. Also, a combined factor which includes RW, Q_Min, Q_Mean, EC, and CA variables was introduced as the best model to estimate RWW variables compared with other factors in both the Bsk and Csa climate categories. On the other hand, qualitative hydrological and land use factors were the weakest factors for estimating RWW variables in the study areas. The performance scoring of the applied models to estimate the RWW in the Bsk and Csa climatic categories are as follow: WGEP, WANFIS, WANN, GEP, ANFIS, and ANN models.

The performance of the W-hybrid models based on R-values for Bsk and Csa climatic categories with the combined factor for three ranges (30%_Min, 40%_Mid, 30%_Max) are listed in Table 6. The three ranges of 30%_Min, 40%_Mid, and 30%_Max are related to RWW 121.61, 121.61 < RWW < 220.50, and RWW220.50 for the Bsk climatic condition and RWW381.18, 381.18 < RWW < 1,568.19, and RWW1,568.19 for the Csa climatic condition in the testing section. The performances of the WGEP models in the three ranges are highly acceptable (R > 0.900) for estimating RWW in both climatic conditions.

On the other hand, the most important advantage of the GEP and WGEP compared with the other applied data-driven methods is in developing and extracting predictive equations. The driving mathematical equations from the GEP and WGEP models for estimating RWW in the Lavasant and Qazvin study areas are listed in Table 7. The mathematical equations can be used at various spatial–temporal scales. The R values of the validation section of the extracted mathematical equations were obtained as 0.865 and 0.897, respectively in Lavasanat and Qazvin for the combined WGEP models. The results indicated the capacity of applying the extracted equations in various basins, which has significant impact on river basin management.

Vafakhah & Bozchaloei (2020) analyzed the flow duration curves through ANN and SVR methods in the Namak Lake basin using Height (H), Area (A), Rangeland Area (RA), Drainage Density (DD), Permeable Formation (PF), and average Stream Slope (SS) variables. Their finding indicated that the R² values of ANN and SVR for the best models were obtained as 0.94 and 0.96, respectively. In our research the best values of R were 0.996 and 0.990 for Lavasanat and Qazvin, respectively. So, the sampling data of various related variables has the ability to estimate the hydrological variables in the Namak Lake basin. The results of the current study are in line with the findings of Adamowski & Sun (2010), Montaseri et al. (2018), and Kumar et al. (2020), regarding the superior performance of W-hybrid models in estimating quantitative and qualitative hydrological variables.

In recent years, lack of proper water resources management, growing demand for natural resource consumption, population growth, and climate change have put a huge tension on water resources. As a result, the scarcity of water resources makes it the most important challenge for mankind in the last decade. The main purpose of this study is to model river water withdrawal for agricultural uses in the Lavasanat and Qazvin study areas with different climatic conditions, while considering morphological, quantitative and qualitative hydrological, and land use factors. Saving time and money has been the reason for applying data-driven methods of hybrid evolution instead of traditional methods by researchers. Therefore, in this study, ANN, WANN, ANFIS, WANFIS, GEP, and WGEP models were used to estimate the amount of the RWW variable according to the various factors. Climate change and consequently changes in surface water resources with the influence of nonlinear multivariate conditions on morphological, hydrological, and land use factors have led researchers to use nonlinear models. Employing nonlinear methods such as hybrid data-driven evolution methods is one of the ways to overcome this problem. The results of this study showed that the combination of data-driven models with wavelet theory can improve the performance of models. Also, this study confirms the high dependence of RWW for agricultural uses on morphological, quantitative and qualitative hydrological, and land use factors. The best way to estimate the amount of RWW is to combine the more effective variables of the four mentioned factors. Since the used variables in this study were collected via field sampling from the study area, the researchers faced limitations in measuring and collecting all the variables related to the selected factors for estimating RWW in the study areas. Therefore, it is recommended to estimate the amount of RWW by measuring more variables related to morphological, quantitative and qualitative hydrological, and land use factors. On the other hand, the RWW for agricultural uses could be affected by more characteristics such as weather variables, population, etc. which are suggested to be included for future research. The use of other data-driven models such as MODWT and MODWT-MRA (for decomposing the data) to estimate the RWW could also be very useful in water resources management at the basin scale.

The authors thank Sari Agricultural Science and Natural Resources University (SANRU) for its financing support [Grant No. 02-1399-26].

None.

All relevant data are included in the paper or its Supplementary Information.