Improving the performance of rainfall-runoff models using the gene expression programming approach

Decision-making in water resources management, hydraulic structures, and flood control planning is of particular importance. Therefore, an acceptable estimate of runoff from rainfall is required. The results of rainfall-runoff models, divided into two categories of distributed and integrated, play an essential role in these decisions.

The effective management of the existing water supplies is a matter of utmost concern, especially in countries where freshwater resources are limited and scarce. In many parts of Iran, rivers serve as resources or potential resources of domestic, industrial, and agricultural water. In these regions, simulating the temporal variations in the flow of rivers is an essential requisite for the effective and sustainable management of basins (Adib et al. 2021a).

While process-based models use theoretical equations for the components of hydrological processes, conceptual models exploit simplified conceptualizations to describe each process considered (Carcano et al. 2008). Data-driven models, on the other hand, completely ignore the physics of the process under study and could be considered as alternatives to process-based or conceptual models due to their minimal information requirements. Moreover, based on the spatial distribution of the input data, models can be classified into lumped, semi-distributed, or distributed ones. Since the early 1960s, various process-based and conceptual hydrological models, such as TOPMODEL (Beven & Kirkby 1979), HEC-HMS (HMS 2000), SHE (Abbott et al. 1986), SWAT (Srinivasan et al. 1998), IHACRES (Young & Beven 1991), ARNO (Todini 1996), and SHETRAN (Ewen et al. 2000), have been developed to address a wide range of environmental and water resource problems. A review of previous studies (Clarke 1973; Singh & Woolhiser 2003; Adib 2008; Wu & Chau 2013; Paniconi & Putti 2015; Taormina & Chau 2015; Fatichi et al. 2016; Ashrafi et al. 2020; Farajpanah et al. 2020; Keum et al. 2020) provides an overview of different aspects of process-based and conceptual hydrological models. Despite considerable progress in developing such models, there are some factors, such as inadequate data availability, that may limit the use of process-based or conceptual hydrological models. Numerous studies have used rainfall-runoff models at the basin scale. Fernando et al. (2012) used GEP to combine the outputs of four lumped runoff models, i.e., LPM, LVGFM, SMA, and PDISC. This combination improved the runoff simulation performance. Salehpoor et al. (2018) used the SWAT and WEAP models for integrated water resource allocation in the Hablehroud Basin in the statistical period of 1998–2012. Lotfirad et al. (2018) exploited the IHACRES model to simulate the streamflow of the temperate and mountainous Navrood Watershed in the north of Iran. Their results indicated acceptable IHACRES performance in temperate areas. Lotfirad et al. (2019) employed the IHACRES model to estimate the runoff of the Hablehroud basin in the future and under climate change scenarios. Ahmadi et al. (2019) used SWAT, IHACRES, and ANN to estimate the runoff of the semi-arid watershed of Kan in Iran. The ANN, SWAT, and IHACRES showed the highest performance in flow simulation, respectively. Jaiswal et al. (2020) compared the performance of TANK, AWBM, and SWAT in the Chhattisgarh Basin, India; according to their results, the TANK model performed better than the other models. Onyutha et al. (2021) used seven lumped models to simulate the daily flow of the Kafu River in East Africa from 1952 to 1981. They calculated the flood return periods of these seven models and the average outcome of the models. The average of the seven models had good accuracy. Arsenault et al. (2015) studied four lumped runoff models in 429 sub-basins in the United States. Kunnath-Poovakka & Eldho (2019) used lumped GR4j, AWBM, and Sacramento models in five watersheds in India. They concluded that the GR4 J model performed better than the other two models.

The objective of the present study is to assess the simulation abilities of a semi-distributed and entirely conceptual model (SWAT) compared to four lumped and semi-conceptual models (i.e., HBV, IHACRES, AWBM, and SMA) to select an effective model to simulate flows of the Hablehroud River in north-central Iran. The prevailing climate of Iran is arid and semi-arid. The Hablehroud Basin represents the arid and semi-arid region and is discharged into the Garmsar Plain, which is one of the most important agricultural regions in Iran. Therefore, the hydrological simulation of this basin is essential. An ensemble of the outputs produced by the five above-mentioned models was also generated to improve the results.

The Hablehroud River in Tehran Province, north-central Iran, was selected in the present study. This river rises in the Alborz Mountain Range and generally flows south-westward through the eastern part of Tehran Province, discharging into the Garmsar Plain in the Kavir Desert (Figure 1). The Hablehroud River has a length of 119.5 km and a drainage basin of 3,261.2 km². The average minimum temperature, maximum temperature, precipitation, and discharge of the basin were 6.4 °C, 14.8 °C, 235 mm, and 7.5 m³/s, respectively, from 1985 to 2005. The regime of its principal tributaries, the Firuzkuh and the Namroud, is combined rain-fed/snow-fed. The highest flows occur in April and May on average. According to the Köppen-Geiger climate classification system (Kottek et al. 2006), the climate of the Hablehroud Basin is mid-latitude arid (BWk) in the south and mid-latitude semi-arid (BSk) in the north (B: arid; W: desert; S: steppe; k: cold). The daily-observed data of precipitation, minimum and maximum temperatures, sunshine hours, wind speed, and relative humidity recorded at eleven meteorological stations throughout the basin and the daily streamflow data recorded at the basin outlet were used. The FAO-56 Penman-Monteith (FAO-PM) equation (Allen et al. 1998) was used to estimate the daily reference evapotranspiration rates.

Land use is a critical factor affecting the runoff, evapotranspiration, infiltration, and surface erosion of a basin. Therefore, to prepare land use maps, Landsat 5 and the Enhanced Thematic Mapper (ETM+) satellite images were used, as shown in Table 1. Also, ENVI 5.3, Google Earth, and ArcGIS 10.3 were employed to process the satellite images and obtain land use information and the number of changes made. The satellite images were classified into six groups, as shown in Table 1.

ArcSWAT 2012, an ArcView extension and graphical user input interface for SWAT, was used to set up the model. Model calibration and uncertainty analysis were carried out using the sequential uncertainty fitting program of SUFI-2 (Abbaspour et al. 2004, 2007). In SUFI-2, SWAT output uncertainty was quantified using the 95PPU band ranging from the 2.5th percentile to the 97.5th percentile of the SWAT output. The percentage of the observed data bracketed by 95PPU is defined as the p-factor (Abbaspour et al. 2004; Yang et al. 2008). The average width of the 95PPU band, divided by the standard deviation of all observed data, is known as the d-factor (Abbaspour et al. 2007; Adib et al. 2021b). While the ideal value of the first factor is 100%, the second factor has an optimal value of zero. There is a trade-off between these two factors. Abbaspour et al. (2015) recommended the p-factor to be greater than 0.7 and the d-factor to be smaller than 1.5.

The region's cultivation pattern, planting and harvesting date, amount of water required, and irrigation resource were introduced as information and agricultural management to the SWAT model. A variety of plants are cultivated in the Hablehroud basin. However, it is impossible to simulate all of the plants. Thus, this study considered the dominant cultivation plants of the region, including wheat, barley, alfalfa, and potatoes as crops and apples, walnuts and peanuts as orchard products (Table 2).

The HBV hydrological model was developed by Bergstrom & Forsman (1973) at the Swedish Meteorological and Hydrological Institute. It is a widely used semi-distributed conceptual rainfall-runoff model. HBV has a simple structure and uses a few input variables to simulate the most important processes affecting the runoff produced in a basin. According to elevation and land cover, the model divides the basin into sub-basins. HBV uses precipitation, temperature, and potential evapotranspiration as input variables to simulate the streamflow and soil moisture content. HBV light is an updated variant of the original model with new capabilities such as simulating the groundwater contribution to the streamflow and considering a different basin reaction procedure through the delay parameter. Figure 4 represents a schematic of the HBV model.

The AWBM model was developed by the Australian CRCCH Institute to simulate the hydrological processes of catchments and has been one of the simplest and most practical continuous models in recent years in the world, especially in Australia (Boughton & Chiew 2007). Despite its simple structure, the AWBM hydrological model accurately calculates the excess rainfall. The data required in model calibration include precipitation, evapotranspiration, and observational discharge. Also, it is possible to obtain features such as automatic calibration through various optimization algorithms, e.g., the genetic algorithm, pattern search multi-start, uniform random sampling, Resenbrock multi-start optimizer, Resenbrock single start, SCE-UA, and pattern search.

If the storage moisture is negative, it is considered to be zero. If the storage moisture exceeds the reservoir capacity, the excess water is converted into runoff, and the storage moisture remains equal to the reservoir capacity. The model assumes that runoff results from two main resources: surface runoff and base water (Figure 5).

The model parameters are (1) the baseline flow index, (2) daily flow dryness constant, and (3) surface storage capacity (C1, C2, and C3) and the corresponding levels of the capacities (i.e., A1, A2, and A3). They are calculated by multi-regression method and automatic variables (Table 3).

The SMA model is a lumped rainfall-runoff model based on soil moisture to simulate the runoff produced by a basin. SMA divides the basin into two zones: a lower zone and an upper zone. It seeks to define the distribution of moisture and free water components in each zone using a set of parameters. Precipitation and temperature are the input variables. Parameters defining the states of soil moisture and parameters describing the relative permeability of the basin are used to estimate the input water, stored water, and output water of the basin (Sorooshian et al. 1993). The SMA model can estimate hydrological processes, such as evapotranspiration, percolation, interflow, and different forms of runoff from a basin. It can be used to forecast river flow or water supply or to assess the effects of climate change on hydrological variables. The model works best for large drainage basins and uses multiple years of records for calibration. Figure 6 illustrates a schematic of the SMA model.

Gene expression programming (GEP) is a method for developing computer programs and mathematical modeling based on evolutionary calculations. It is inspired by natural evolution. GEP was invented by Ferreira in 1999 and was officially introduced in 2001 (Ferreira 2002). It selectively abandons an unsuitable population and produces modified offspring. No functional relationship is considered in the beginning of the process. This method can optimize the model structure and components.

GEP works on the tree structure of formulas. Tree structures are created from a set of functions and terminals. The GEP steps include:

• (1)

  A set of terminals (problem variables and random fixed numbers),

• (2)

  A set of mathematical operators used in formulas,

• (3)

  Selecting the fitting function, and

• (4)

  Determining program control parameters (e.g., the population size and operator weights).

This study used GeneXproTools V.5.0 to apply the GEP model. The main idea was to construct a decomposition tree and mathematical functions for model interpretation (Table 4).

The Taylor diagram was employed to assess the accuracy of the models in the calibration and validation phases.

The models were calibrated and validated based on monthly data (daily mean averaged over the month). The observed discharge data during 1985–2001 at the basin outlet were used in calibration, and the discharge data during 2002–2005 were employed to validate the calibrated models. The SWAT model inputs, such as temperature and precipitation, were daily data. The SWAT was calibrated and validated based on monthly data (daily mean averaged over the month). The fitted, minimum, and maximum values of the model are shown in Table 5.

The genetic algorithm and Powell (GAP) method was used to calibrate the parameters of the HBV model (Pervin et al. 2021). The characteristics of the GAP method are listed in Table 6. The optimized values for the HBV parameters in the routines (i.e., the snow routine, soil moisture routine, response routine, and routing routine) are represented in Table 7.

The SMA model uses different methods, including the loss method (to model infiltration losses combined with the simple canopy and the simple surface method), SCS (to model the transformation of precipitation excess into direct surface runoff), Muskingum routing method, and linear reservoir method (to model base flow) in order to simulate the components of the rainfall-runoff process. The optimized values of the parameters in the aforementioned methods are presented in Table 8.

IHACRES was calibrated based on the daily mean temperature, rainfall, and discharge data. Finally, the model was calibrated and validated in certain periods. The fitted parameter values for the non-linear and linear modules of the model are shown in Table 9. The exponential loss parameter (p) and the threshold for the rain to give streamflow (l), were considered, respectively, one and zero.

To calibrate the AWBM model, the multi-start pattern search method was used. The Nash-Sutcliffe efficiency coefficient and baseflow method were considered as the first and second objective functions, respectively. The fitted parameter values are listed in Table 10.

The fitness statistics (i.e., RMSE, NS, R², and KGE) in the calibration and validation phases are given in Table 11. Table 11 shows satisfactory results for the SWAT, IHACRES, and HBV models in the calibration phase; the RMSE values are reasonably small, and the NS values are more significant than 0.50. In general, the river flow simulation can be said to be satisfactory if NS is more significant than 0.5 (Moriasi et al. 2007; Salehpoor et al. 2018). However, the negative values of NS in the validation phase indicate that HBV, IHARES, AWBM, and SMA failed to correctly simulate the process of rainfall-runoff in the basin under study. Also, only SWAT had good performance in the validation phase and outperformed the other models.

Figure 7 compares the monthly observed and simulated discharges in calibration and validation. While the overall fitness is satisfactory for calibration, the peak flows that occur in response to high rainfalls were generally underestimated by the models. According to Figure 7, the SWAT model had the best performance in simulating peak flows in the basin under study. The models are also compared in Figure 7 in the calibration and validation phases using the Taylor diagram (Adib et al. 2021c; Taylor 2001).

Figure 7 shows that the models developed for the basin under study had approximately the same standard deviation in the calibration phase. However, all the models underestimated the observed standard deviation in the calibration phase and had insignificant variability. For the validation phase, the HBV model and the observations have approximately the same standard deviation. It can be seen in Figure 7 that the SWAT model generally agreed best with the observations, with the highest calibration and validation correlation coefficients and lowest calibration and validation RMSEs.

The fitness statistics (i.e., RMSE, NS, R², and KGE) of the hybrid model in the calibration and validation phases are given in Table 11. According to the NS values, the hybrid model had satisfactory performance in both calibration and validation phases. It is seen from Table 11 and Figures 7 and 8 that the hybrid model outperformed all the hydrological models (the lowest RMES and the highest NS, R², and KGE in both calibration and validation phases and the most realistic model in the Taylor diagram).

As shown in Figure 8, the SWAT and hybrid models realistically calculated high flow and low flow close in both the calibration and validation phases. To more accurately evaluate the performance of the rainfall-runoff model and combined models, flows below the first quarter were selected as low flow, and flows above the third quarter were treated as high flow. According to Table 12, the performance of the models was much better in high flow simulation than in low flow simulation. In high flow simulation, the best performance was found for GEP, SWAT, and IHACRES models in calibration, respectively. For the validation phase, the best performance was shown by the SWAT, GEP, and HBV models, respectively. This is due to the higher variance of high flow than low flow; the simulation value is more likely to fall within this range at a higher variance.

The sensitivity analysis of parameters is an essential step in SWAT calibration. SUFI-2 uses multiple regression analysis to identify the relative significance of the parameter via the t-stats of the parameters (i.e., the parameter coefficient divided by the standard error) and the p-value of the t-stat. The present study started the calibration process with 29 parameters in SUFI-2. As can be seen in Figure 9, seven parameters, including available water capacity (SOL_AWC), initial SCS runoff curve number for moisture condition II (CN2), initial depth of water in the shallow aquifer (mm H2O) (SHALLST), base flow alpha factor (ALPHA_BF), soil evaporation compensation factor (ESCO), groundwater delay time (GW_DELAY), and snowfall temperature (SFTMP), were found to have the highest sensitivity to the monthly discharge of the basin as they had the largest absolute t-stat values and p < 0.1.

The local sensitivity analyses of the lumped models HBV, SMA, IHACRES, and AWBM were carried out to identify the most sensitive parameters with the greatest impacts on the model output (Figure 10).

The IHACRES model had fewer parameters than the other models. The parameter Tw had the highest sensitivity to the objective function (NS). Among the HBV parameters, UZL, K1, FC, TT, and MAXBAS were very sensitive with higher sensitivity to reduction than increase, based on the NS criterion. These results are consistent with Yaghoubi & MASSAH (2014).

In the AWBM model, the parameters BFI, C3, Kbase, and Ksurf were the most sensitive ones; Ksurf was more sensitive to reduction, and C3 was more sensitive to increase. This is consistent with Boughton & Chiew (2007) and Chouhan et al. (2016).

In the SMA model, the parameters GW1 storage, GW1 coefficient, impervious, and tension storage were more sensitive to reduction than increase. This is consistent with Bashar & Zaki (2005) and Ouédraogo et al. (2018).

The main purpose of this study was to provide an accurate method for estimating runoff in the arid basin of Hablehroud as a significant representation of other basins in Iran as most Iranian basins have arid and semi-arid climates. The novelty of the present study includes the use of lumped and distributed rainfall-runoff models in order to simulate the flow of the arid basin of Hablehroud and combine the outputs of the models with GEP and provide a mathematical relationship by exploiting the advantages of all five models used in runoff simulation. The ensemble model outperformed the other models in both the calibration and validation phases. This was also reported by Fernando et al. (2012); although Fernando et al. (2012) used only lumped hydrological models, the present study employed both lumped and distributed models.

Various models are available to simulate the complex rainfall-runoff process in watersheds. To find the most appropriate model for simulating the monthly flow (daily mean averaged over month) for the Hablehroud River in Tehran Province, north-central Iran, five hydrological models were calibrated and validated. The following conclusions could be drawn from the present study:

1. SWAT needs a variety of input data, given the limited availability of soil data in the Hablehroud Basin. However, the model fitness statistics (i.e., RMSE, RSR, NS) showed satisfactory results in both calibration and validation. The highest water consumption in the region was found in agriculture. In this research, only the SWAT could accurately model agricultural consumption. While the HBV, IHACRES, AWBM, and SMA models produced relatively acceptable results in the calibration phase, they failed to simulate the validation dataset.

2. In comparison to the SWAT model, IHACRES is a parsimonious model requiring only seven parameters. However, the performance of IHACRES in simulating the river flow in the study area in the calibration phase was quite comparable to that of the SWAT model. As shown in Figure 8, the ensemble, SWAT, and IHACRES models properly showed the fluctuations of river flow in the simulation period (wet and dry periods). The performance of the five hydrological models in high flow estimation was better than in low flow; SWAT and IHACRES were more accurate in estimating high flow than the other models. The ensemble model outperformed the five hydrological models in estimating high flow and low flow.

3. The GEP ensemble model outperformed the other models in both the calibration and validation phases since the GEP model seeks the best fit of the objective function by selecting the appropriate operators and constants. In other words, river streamflow is a nonlinear and complex phenomenon, and the GEP method has a high ability to model nonlinear processes. This study exploited GEP to combine the outputs of the hydrological models. In the review of the literature on combined lumped-distributed hydrological models, no study was found to have employed GEP. Therefore, the present study is a pioneer work in combining lumped and distributed hydrological models with GEP.

The authors wish to thank I.R. of Iran Meteorological Organization (IRIMO) and Iran Water Resources Management Company (IWRM) for providing the data used in this study. We should also like to thank the reviewers and the editors who have given up valuable time to review the manuscript and offer valuable comments.

The authors declare no conflict of interest.

This research received no external funding.

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