Towards a hybrid algorithm for the robust calibration of rainfall–runoff models

Conceptual rainfall–runoff (CRR) models are frequently exerted by hydrologists and water resources planners to interpret and manage both climate change and anthropogenic activities affecting the watershed system. Since the defined parameters in these models are often not directly measured on account of the physical constraints and scaling problems, the effective implementation of modeling depends primarily on how and to what extent the model is calibrated (Goswami & O'Connor 2007; Zhang et al. 2009; Qin et al. 2017). Early calibration studies were performed manually due to the lack of inadequate usage of optimization tools. In parallel with the developments in the field of optimization sciences, the model calibrations have been carried out with automatic search approaches for nearly 30 years (Piotrowski et al. 2017). The first notable applications of automatic calibration were performed through basic techniques such as simplex and pattern search algorithm in the category of local search algorithms. In the same category, gradient-type algorithms (such as Gauss–Newton) have been also preferred for the calibration of CRR models owing to their rapid convergence features (e.g., Hendrickson et al. 1988). Then, the non-convex nature of the calibration problem in the CRR models was demonstrated, and it was stated that the multiple local minima presence in the models could made the exploring of global optimal point complicated. It was also detected that the obtained results were hypersensitive to the initial parameter estimations (Duan et al. 1992). These constraints and increases in degrees of freedom based upon model structures have led hydrologists to global search algorithms with stochastic characteristics rather than local search algorithms since the 1990s. The metaheuristics such as genetic algorithm (GA), shuffled complex evolution (SCE) algorithm and simulated annealing (SA) algorithm have been the first global search techniques adapted to hydrological model calibration (Duan et al. 1992; Franchini & Galeati 1997; Wang 1997; Thyer et al. 1999).

Even though the GAs can overcome a great deal of trouble encountered in local search methods, the relatively high number of variables controlling the crossover, mutation and selection operators can confuse the user in the decision-making phase. In addition to this matter, the fact that it does not guarantee the same global solution for each run has triggered the development of more practical and stable metaheuristics (Li et al. 2016). The differential evolution algorithm (DEA) and the particle swarm optimization (PSO) algorithm, which are the elementary alternatives to GAs, have been used in the relevant literature (e.g., Goswami & O'Connor 2007; Zhang et al. 2009; Piotrowski et al. 2017). As the number of nature inspired metaheuristics is increasing day by day, it will be unavoidable to adapt them to the scope of hydrological model calibration. For instance, relatively new algorithms such as artificial bee colony (ABC), ant colony optimization (ACO), artificial immune system (AIS), cuckoo search (CS), harmony search (HS), invasive weed colonization (IWA) and successful history-based DEA with linear population size reduction (L-SHADE) have been found to be implemented to water resources issues in a limited number (e.g., Zhang et al. 2009; Arsenault et al. 2014; Chen et al. 2015; Asgari et al. 2016; Piotrowski et al. 2019).

In the above paragraphs, the effectuations of algorithms ranging from local search ones to population-based metaheuristics are mentioned within the scope of the calibration of CRR models. Apart from a few exceptions, it is widely accepted that the derived results from algorithms such as GA, SCE, DEA and PSO are more stable than those of local search algorithms. However, some contradictions have been experienced when questioning the superiority of them against each other. In other words, a study in the literature states that a certain algorithm is fairly prosperous, while another study argues that different methods may produce similar results in both statistical and hydrological terms. For instance, Wang et al. (2010) used two variants of GA and the SCE to calibrate a distributed rainfall–runoff model and demonstrated that the calibration performance obtained from three algorithms was quite similar for a catchment in Taiwan. In a study conducted by Goswami & O'Connor (2007), algorithms such as GA, SA, SCE and Nelder–Mead Simplex were tried to calibrate the soil moisture accounting and routing the model for various basins in Ireland and China. Though their results did not imply statistical significance, it was reported that SA exhibited slightly better performance. Piotrowski et al. (2017) tested both the diverse variants of PSO and DEA and much more sophisticated metaheuristics in the calibration of two CRR models but emphasized that the responses from employed algorithms did not differ significantly. In contrast, some algorithms have been claimed to have obvious advantages over other ones in various studies. In contrast to the findings of Goswami & O'Connor (2007), Cooper et al. (2007) underlined that SA was found to be weaker than other adopted algorithms. In a study performed by Zhang et al. (2009), the Soil and Water Assessment Tool (SWAT) was calibrated by GA, SCE, PSO, DEA and AIS algorithms. According to their results, if the control parameters are well adjusted and the number of fitness functions calls is increased, the appropriate results can be achieved through GA. However, it has also been proved in the same study that the PSO can give reasonable results with low population size and few runs. In the study prepared by Arsenault et al. (2014), 10 algorithms were utilized to calibrate parameter sets for three models on several basins. Among them, the weakest performances were offered with GA, PSO, DEA, CS and HS. The above-mentioned inferences will seem to continue to be made.

So, what could be the reasons of these dissimilarities? According to Kavetski & Clark (2011), incorrect tuning of the variables controlling the algorithm, the variability of the tested hydrological models, the study regions and thus the data attributes/quantities may lead to the aforementioned troubles. Besides, despite the robust structures of the metaheuristics, in a comprehensive work prepared by Tolson & Shoemaker (2007), they were shown to be markedly influenced by the dimension problem. This situation has been broadly discussed by various researchers, and it is stated that the calculation process with metaheuristics can move away from hydrology practice (Qin et al. 2018). According to Qin et al. (2017, 2018), the rapid convergence feature of gradient-based algorithms should be utilized to effectively reduce the computational costs, and they should be made more resolute to get rid of local minima traps. One of the strategies to be applied for this purpose is to run consecutively the gradient algorithm by assigning multiple random parameter values to the CRR model initially. Qin et al. (2017) indicated that Gauss–Newton techniques with the multi-start version produced more convergent solutions than SCE. In terms of the frequencies of convergence to the global optimum, it is expressed that there are model-based uncertainties and extra effort is required to obtain a robust case.

Another remedy is to integrate global search properties of metaheuristics with local search capabilities of gradient algorithms (Karahan et al. 2013). In this context, a hybrid algorithm proposed by Karahan et al. (2013), which consists of the combination of HS and the quasi-Newton, yielded rather stable results in their flood routing problem. However, this type of hybridization has not yet been operated for the calibration of CRR models except for a few exemplary applications (e.g., Qin et al. 2018). This approach has shed light on us to give a new perspective to the calibration step of CRR modeling. Instead of choosing the best one among many algorithms, the idea advocated in the study is to propose a hybrid algorithm, which is both fast convergent, and has an adaptable content at the point of relieving the uncertainties mentioned. On the other side, it is necessary not to ignore the issue discussed by Kavetski & Clark (2011) to avoid an inadequate selection of control parameters of a selected algorithm. In this respect, it was considered that coupling a computationally intensive algorithm like HS with a gradient-local search algorithm would not be a practical procedure. Therefore, it was decided that providing local search capability to an elementary algorithm such as PSO, which has a relatively minimal number of control constants, would be more consistent. In this respect, hybrid PSO (HPSO), which was based on the combined use of PSO and the Levenberg–Marquardt (LM) algorithm, was formulated and compared with several algorithms in terms of both robustness and convergence. To the best of our knowledge, the usage of HPSO has not yet been identified in the context of the CRR model calibration in the literature. Thus, an attempt was made to fill a gap in previous studies. The remainder of the presented study is arranged as follows. The details about methods employed are given under ‘Methods’. The study area and data are briefly explained under ‘Study area and data’, while the results derived from the adopted algorithms and the conclusions are presented under ‘Results’ and ‘Conclusions’, respectively.

In this study, two lumped CRR models, which are known as dynamic water balance model (dynwbm) and abcd model, were considered. Models have been widely used to simulate runoff under stationary and changed climatic conditions (e.g., Tekleab et al. 2011; Okkan & Kirdemir 2018; Shahid et al. 2018; Okkan & Kiymaz 2020). These models, which require only monthly total precipitation (P) and potential evapotranspiration (EPOT) as inputs, include a series of conceptual soil moisture and groundwater storage functions. Both models utilized in the study contain similar storage elements as seen in Figure 1. In the models, the sum of the two components obtained as direct runoff (Q_direct) and baseflow (Q_base) provides the modeled runoff (Q_modeled).

Among these two models, the dynwbm developed by Zhang et al. (2008) is an expanded analogous of the hypothesis introduced by Budyko (1958). Although there are four parameters in the original model, an additional parameter has been added to its groundwater storage function by Okkan & Kirdemir (2018) so as to improve runoff prediction performances. In the model, P is made up of two components which are Q_direct and catchment rainfall retention X, respectively. In this partition, the parameter α₁ plays an active role, and a larger α₁ value results in more rainfall retention and less direct runoff. The model has also a parameter termed as maximum soil moisture capacity (S_max) representing the soil and vegetation characteristics of the basin. Besides, the parameter α₂ controls evapotranspiration opportunity y, which is assumed to be composed of the sum of soil moisture content S and actual evapotranspiration E_act. During the month i, the available water content W_i can be expressed by the sum of the soil moisture remaining from the previous month (S_i−1) and X, as well as by the sum of S_i, E_act and the amount of recharge (Rec) draining into groundwater storage. After the Rec is taken from the budget calculation, the balance equation is made for the groundwater storage G, which is postulated to represent linear reservoir behavior, using the d and e parameters together, and the baseflow is then predicted (Okkan & Kirdemir 2018).

In the abcd model developed by Thomas (1981), the available water content in the basin during the relevant month is considered to be equal to the sum of the S_i−1 and P_i, with a coarser calculation compared with that of dynwbm. The evapotranspiration opportunity y is calculated depending on the parameters a and b, which represent the propensity of runoff occurrence and recharge before the soil being fully saturated and the soil saturation level, respectively (b in the abcd can be said to be the counterpart of S_max in dynwbm). While a larger a value results in less direct runoff and recharge amount, a larger b value brings about the contrary of this event. After a certain part of the water amount is reserved for y, the remaining portion is allotted to Q_direct and Rec components depending on the linear coefficient c. In the study, as in the previous model, the modified groundwater storage function was regarded in the abcd model, and it was transformed into a five-parameter abcde (see Figure 1).

For the steps following the joint implementation of Equations (1) and (2), the operating mechanism of the algorithms are briefly described under the following sub-headings. Basic operators and control parameters of algorithms are also summarized in Table 1.

After the implementation of all operators in the relevant iteration, new individuals are added to the old population, and then they are reordered so as to eliminate the chromosomes having higher SSE values (competitive exclusion). For both GA and other used metaheuristics, the individual having the best fitness in the current population is stored as a global solution. The above-mentioned operations of GA are continued until the maximum generation number is reached (Tayfur 2017). For the details about the real-coded GAs, the report that was edited by Peltokangas & Sorsa (2008) can be viewed.

In Equation (5), pb_i represents the best position that the ith particle has ever reached, while gb represents the global best solution among pb matrix with the N_pop × N_par dimension (where t = 0 to iter_max; i = 1 to N_pop; j = 1 to N_par). rand is the same variable that is previously defined in Equation (1). ω is the inertia weight controlling particle's velocity. Also, the acceleration coefficients c₁ and c₂ are both usually set to 2.0. pb(t = 0) is considered equal to the randomly generated solution matrix that is formulated in Equation (1). If the fitness value of the solution obtained after the position update for any particle i is improved compared with that of the previous iteration, the corresponding row of the pb is replaced with the new position. A similar approach is repeated iteratively for the gb as well. The operations described above are preceded until the iterations are completed (Tayfur 2017). In the operation of Equation (5), four inertia weight strategies, which are mentioned in Table 1, are examined referring to Rathore & Sharma (2017). For example, the variation of inertia weights over 500 iterations according to the random weight (RW), the linear decreasing weight (LDW), the nonlinear decreasing weight (NLDW) and the chaotic random weight (CRW) is shown in Figure 2(a). In the study, the inertia weight functions were requested to produce their values in a common range, and ω_max and ω_min were taken as 1.0 and 0.001, respectively.

Xu et al. (2012) have stated that F-values between 0.5 and 1.0 provide a similar contribution on the simulations. Besides, the effectiveness of the mutation is yielded by a non-uniform strategy applied during the crossover stage. In the crossover, following the random number (rand) derivation for each defined j gene in the present chromosomes, it is checked whether the rand_j ≤ C_r is ensured. In that case, the state of using the j gene of the vector obtained by Equation (6) occurs. Otherwise, the j gene of the corresponding chromosome is unchanged. Here, C_r is a crossover constant that gives plausible results with its values between 0.5 and 1.0 (Xu et al. 2012). Finally, a greedy criterion is implemented by means of Equation (2) to determine the chromosome that will be transmitted to the new generation. If the fitness function calculated from the candidate chromosome, which is consecutively subjected to mutation and crossover, improves compared with that of the old one, it is stored for the next generation, and the old solution is removed, or else the location of the old vector is preserved. The operations stated above are continued until the maximum generation number is reached (Tayfur 2017).

The seeds produced by Equation (7) are then randomly distributed to the field again. Depending on the number of Seed_i, a uniform random solution is generated in the neighborhood of the parental solution with the line number i. In this study, this process resembles the mutation operation defined in Equation (3), but is applied to the entire parameter vector. The time-related reduction of the distance of the newly produced seeds to their parents is iteratively adjusted by the NLDW stated in Table 1 (see Figure 2(b)). After seed production and spreading, it is checked whether the total population size exceeds P_max, the maximal population size of colony (P_max = N_pop in other algorithms). If such a situation occurs, only the weeds having better fitness values survive, depending on the competitive exclusion procedure. The operations summarized above are managed with regard to the maximum number of iterations. In the study, sensitivity analysis made by Asgari et al. (2016) is referred when assigning the control variables of IWA (see Table 1).

Once p_i > rand, the onlooker bee moves to this source and produces neighboring solutions, as performed in the phase of the employed bees. Similar to other phases, it is once again decided which solutions will be stored and the condition of the failure counter is. If the failure counter in a source exceeds a limit value set by the user, the employed bee using that source is assigned as a scout bee and search for another location in space by means of Equation (1). If the nectar quality after scout bee phase is improved compared with the old one, this solution derived is stored and the failure counter is zeroized again. All the above phases are iterated until the maximum generation number is achieved (Karaboga & Basturk 2008; Chen et al. 2015).

In this section, the integration process of the LM algorithm with PSO, which only requires c₁ and c₂ coefficients and inertial weight, is described. Essentially, the proposed HPSO is not very complicated. The hybrid approach is based on the introduction of the gb vector obtained by the standard operation of the PSO as an input to the LM in each iteration step. A similar procedure was previously operated by Zhang et al. (2007) and Nawi et al. (2014) in the training of artificial neural networks. Different hybrid versions of PSO were also proposed for job shop scheduling problem (Sha & Hsu 2006) and multi-modal functions (Kao & Zahara 2008). Both HPSO and other algorithms are coded in the MATLAB environment.

In Equation (10), λ is multiplied by a certain decay rate β when fitness decreases in a certain iteration step, and it is divided by β when fitness increases in another step. As suggested in the literature, λ_o, β and the maximum limit value of λ (λ_max) were set to 0.01, 0.1 and 10¹⁰, respectively. Since the formulations of RW, LDW, NLDW and CRW used as inertia functions in PSO are also performed for the hybrid approach, the hybrid ones are mentioned as HPSO1, HPSO2, HPSO3 and HPSO4, respectively. The flowchart of the proposed hybrid strategy is given in Figure 3.

In the study, Gediz Basin, which represents a major part of agricultural activities of the Aegean region in Turkey, was selected as a study area. CRR models were applied on three sub-basins that are differing by locations and hydro-meteorological properties. These sub-basins, which are represented by Acisu (E05A23), Hacihidir (D05A28) and Hacihaliller (D05A38) flow gauging stations, are located on the main stream of Gediz, Gordes creek and Nif creek, respectively (Figure 4). Station DO5A28 is also located on the branch that feeds the Gordes Dam reservoir. While the drainage area of the station is 808.2 km², the drainage area of the dam is 1,045.4 km². When switching to the reservoir inflow volumes, the flow volumes of the station DO5A28 are multiplied by 1,045.4/808.2. When topographically examined, it is seen that both Hacihidir and Acisu sub-basins fall within the mountainous area with higher altitudes than that of Hacihaliller sub-basin.

Since the 30-year window regarding the current climate normals is 1981–2010, which is also called a reference climate period, this period has been also considered in this study. Among the related sub-basins, the Hacihaliller sub-basin received the highest annual precipitation with the value of 733 mm, while the annual precipitation values of 573 and 505 mm were observed in Hacihidir and Acisu sub-basins, respectively. The highest mean annual temperature value was observed in Hacihaliller sub-basin as 16.8°C, while it was observed as 13.2°C in the rest of the study area. According to the data obtained from flow gauging stations, the highest flow depth was observed in the Hacihidir sub-basin (∼129 mm), such that the lowest EPOT is predicted in the related sub-basin with the value of 964 mm. Even if the Hacihaliller sub-basin stands for the wettest part among the sub-basins, its flow potential decreases on account of relatively higher temperature and EPOT values.

The all-flow measures obtained from The General Directorate of State Hydraulic Works (DSI) for the 1981–2010 water year period were then converted to millimeters to be used in modeling. Monthly time-scale Thiessen-weighted precipitation series (P) obtained from the data of Turkish State Meteorological Service (MGM) was prepared as a first model input. While monthly temperature and relative humidity data were compiled from the meteorological stations operated by MGM, ERA-Interim reanalysis data sets having 0.75° × 0.75° resolution were used for other variables (wind speed and solar radiation) needed for the Penman–Monteith equation to obtain areal mean EPOT estimations used as a second model input. These reanalysis data are served by the European Center for Medium-Range Weather Forecasts for several categories (http://apps.ecmwf.int/datasets/). When examined in terms of the drought index based on the P/EPOT ratio on the annual time scale, Acisu has semi-arid climate characteristics, while arid–semi-humid climate is prominent for the remaining sub-basins. For the land use/land cover classification of the region, readers can be directed to Droogers & Kite (1999). Further details about the study region and ERA-Interim grids that almost uniformly encompass the basin can be found in Okkan & Kirdemir (2018) and Okkan & Kiymaz (2020).

In the parameter optimization studies carried out with metaheuristics, the population size (N_pop) and the maximum iteration number (iter_max) to be used in the cycle of algorithms assuredly influence the calibration. Since the convergent and stable features of the proposed HPSO technique were advocated within the scope of the study, it was checked over whether the computational intensity diminishes with this approach by using fewer iterations and population sizes. For this purpose, in all algorithm experiments, N_pop and iter_max were fixed to 20 and 500, respectively. Additionally, a wide range of parameters have been defined for CRR models to make algorithms further force in the process of finding the global solution. The range of parameters S_max (in dynwbm) and b (in abcde) is 10–1000 mm, while the range 0.01–0.99 is assigned for parameters α₁, α₂, d (in dynwbm) and a, c, d (in abcde). The range of 0.5–0.99 is chosen for the common parameter e in both models. These parameter ranges are also within the physically possible limits, and there is no need to define a comprehensive penalty function within the fitness function. Not only the initial solutions generated by Equation (1) but also the iterative solutions within the algorithm loops are kept between these closed ranges. In the study, the split-sample procedure, where the observed runoff series were divided into two equal parts for calibration and validation, was implemented (Refsgaard & Knudsen 1996; Zhang et al. 2008). According to this, the data covering the water period of 1981–1995 compiled for sub-basins were evaluated during the parameter calibration stage of the CRR models, while data covering the water period of 1996–2010 were used in the validation.

In some daily hydrological model applications, initial soil moisture can also be defined as a parameter to be calibrated. One such application for the daily GR4 J model (modèle du Génie Rural à 4 paramètres Journalier) was performed by Mostafaie et al. (2018). As the models used in the study are monthly time-scaled ones, they are not overly sensitive to the initial soil moisture and the initial groundwater store values. In the first month of the water year period (at the end of the dry period), due to the low rainfall, the initial storage values in the conceptual reservoirs of CRR models (S₀ and G₀) have been chosen close to almost zero (5 mm was fitted for these initial values).

As mentioned in the previous sections, the algorithms were run 30 times according to their control variables specified in Table 1, and the CRR models' parameter ranges specified, and the results were then stored individually. For each CRR model, a total of 13 different algorithms, which are four different inertia weight variants of PSO and HPSO, and GA, DEA, IWA, ABC and LM, were operated. As the implementation of them on two CRR models for three sub-basins produced 78 graphical results, only some convergence graphs produced with dynwbm exerted to the Acisu sub-basin have been given as an example in Figure 5. From Figure 5, it can be clearly seen that some algorithms have relatively different convergence and stabilization properties. In this example, when focusing on the last 50 iterations of 30-run average values, it can be monitored that especially GA and PSO variants do not provide stable results and are possibly sensitive to initial solutions. On the other hand, DEA, IWA and ABC produced more stable results than other metaheuristics, but they displayed late convergence. The same example also proves that the HPSO shows both stable and rapid convergent performances. Instead of graphically interpreting the completed analyses for two CRR models applied to three sub-basins, descriptive statistics of the fitness values realised at the last iteration (the statistics concerning 30 fitness values generated under each algorithm) are extracted in Table 2. In order to quantify the number of iterations required and the convergence performance of the algorithms against different initial solutions, the formulas h and CR expressed in the section ‘Performance measures for algorithms’ were implemented. Following being stored of CR and h values obtained from the related cycle for each run, the mean statistics of these indices calculated from 30 runs are also stated in Table 2.

With the approach formulated as above, the performances of the algorithms in the corresponding model and sub-basin variations were rated from 0 to 1 and are stated in Figure 6. After the rank statistics of the p_geomean values were taken for each sub-basin and the CRR model, the general grading of utilized algorithms is introduced in Figure 7, in which p_geomean is inversely proportional to rank statistics. With the help of the mean statistics of the ranks compiled from the sub-basins, it is tested in Figure 7(a) whether the algorithms react distinctively depending on the CRR models. According to Figure 7(a), excluding a couple of PSO variants, there was no significant difference in the responses of the algorithms to the models. Here, the high determination coefficient detected in inter-model scattering has supported this remark (R² = 0.901). The detection made at this stage has also provided more comfortable comparison of the algorithms within themselves. Then, the mean statistics of the rank numbers arranged under all cases are indicated in Figure 7(b) for each algorithm.

Indeed, Figure 7 summarizes the indications presented in Table 2 and Figure 6. From Figure 7(b), it can be seen that the CRW inertia weight is more compatible with the PSO. Even though this inference is parallel to that advocated by Rathore & Sharma (2017), all PSO types have settled at the last ranks along with GA. The remaining metaheuristics (DEA, IWA and ABC) have produced relatively more reliable results than PSO and GA. At the same time, the most momentous diagnosis is that all HPSO variants appear in the first four and are conspicuously steady in all existing conditions. Also, HPSO4, operating with CRW, has notably come to the forefront compared with other ones since it has shown similar responses in both CRR models and has been treated as stable and well convergent. While LM and PSO algorithms individually give almost the worst results in all applications, it is an essential detection that the best solutions are achieved under each variation in case of their combining. This supports why we couple derivative algorithms and metaheuristics for an enhanced calibration.

Surely, it would be inevitable that LM combined with other metaheuristics would have the same impression. In the context of being exemplary, to what extent the HDEA approach made by the hybridization of DEA with LM has a resemblance to HPSO4 has been questioned. As the mean fitness values obtained from both hybrid algorithms after 30 runs are almost equal, the mean convergence rates and the standard deviations of the fitness values at the last iteration of completed runs are given in Figure 8. It can be seen from Figure 8 that there are no marked variations between the convergence performances of HPSO4 and HDEA on model calibrations. The standard deviations of the fitness values are in the range of 10⁻¹⁰–10⁻¹² in both hybrids as well. In those cases, HPSO4 with fewer control parameters has been proved to be a more reasonable method again. Furthermore, it will be preferred in various calibration practices since it has carried out its convergence with low population size and a moderate iteration number.

For sub-sampling formulations related to the terms in Equation (14), readers can go through the studies performed by Qi et al. (2016) and Yip et al. (2011). Similar to the procedure given by Qi et al. (2016), sub-combination groups belonging to the c₁ and c₂ coefficients ranging from 1.0 to 3.0 with an increment of 0.5 were formed. Together with the combinations generated between parameter value pairs, PSO4/HPSO4 algorithms were run 5 × 5 × 30 = 750 times in total. Figure 9(a) outlines the results of SSE collected from the runs by means of contour plotting for the Hacihidir sub-basin example. It can be clearly seen from Figure 9(a) that the coefficient combinations actuated in HPSO give more stable outputs and these are much closer to uniform, excluding the coefficients greater than 2.5. According to ANOVA analyses performed to analytically interpret Figure 9(a) and to decompose the sources of uncertainty, the contribution of internal variability over total variance was more dominant in both algorithms (see Figure 9(c)). This is, of course, a result of the stochastic structure in the algorithms. From Figure 9(b), it is noteworthy that the internal variability in HPSO is 90% less than that of PSO. (The vertical axis in Figure 9(b) is the logarithmic base.) Moreover, individual influences of parameters and their interactions do not have significant impacts on HPSO performance, while the same conditions add much more uncertainty to the PSO. It can be observed that individual parameter uncertainties (Varc1 and Varc2) and c₁–c₂ interaction uncertainties in HPSO are 98.3%, 98.9% and 91.4% less than those of PSO, respectively. Similar inferences were obtained for other sub-basins but not shared due to word limit. As a consequence, the fact that HPSO is not extremely sensitive to control parameters makes the algorithm more reliable again.

Although the main goal of this study is to scrutinize the performance of the algorithms and to recommend the coexistence of gradient approaches and swarm intelligence optimizations, for drawing a hydrological perspective, it would be useful to interpret the validation outputs of the CRR models, whose calibrated parameters are denoted in Table 3(a). Acisu and Hacihidir basins have a long-term average annual precipitation of about 550 mm, while this value is nearly 200 mm higher in the Hacihaliller basin. This was reflected in the estimated S_max and b parameters, and significant correlations were obtained between these soil moisture storage parameters and annual precipitation regimes. Additionally, since both models have similar evapotranspiration opportunity mechanisms, strong negative correlations were detected between the predicted α₂ and [1 − c] values and long-term annual EPOT values. Nash–Sutcliffe coefficient (NS), RSR, which standardizes root mean square error (RMSE) using the standard deviation of observed runoff, and the percentage of bias (Pbias), were also assessed in the performance check of abcde and dynwbm (see Table 3(b)). The details about the model ratings (very good, good, sufficient and insufficient) can be accessed from Moriasi et al. (2007). In their criterion, the model is rated very good if the NS is greater than 0.75 and the RSR is less than 0.50. In this regard, these two indices in both the calibration and validation period of CRR models refer to very good modeling. However, in general, dynwbm provides more appropriate results in terms of the indices, whereas abcde comes into prominence during the validation period. On the other hand, the Pbias shows that the dynwbm falls into the good model category in the validation period, while inter-model variability pertaining to this index is found to be more pronounced. (According to the same criteria, if the Pbias is less than ±10%, the model is rated very good, while the model is rated as good at the values of Pbias in the range of ±10–15%.) In this sense, there may be potential reasons why the model outputs display dissimilarities. For example, although the EPOT is based on the reference method Penman–Monteith, it has been determined by Okkan & Kiymaz (2020) that dynwbm in the same basin shows different responses to several empirical EPOT methods, and it is emphasized that the most appropriate EPOT equations should be decided locally at this stage. On the other hand, predictions obtained with the reference to a single CRR model may contain a set of uncertainties. Therefore, it has been emphasized in Okkan & Inan (2015) and Okkan & Kirdemir (2018) that multi-model ensemble approaches that evaluate multiple model outputs together are more consistent in the stage of making hydro-meteorological projections.

In this study, a hybrid scheme was offered to overcome some calibration handicaps, which are met in metaheuristics, such as slow convergence and failure to achieve global minimum solution in each individual run. In the application, algorithm experiments were carried out on two different CRR models by assessing three sub-basins having relatively different drought classes and runoff regimes, and it was questioned whether a common inference was derived from all practices made. First of all, the performances of exerted popular algorithms on the calibration process were stored. As a result of various criteria, the all variants of PSO and GA have not produced reasonable results, while DEA, IWA and ABC are more qualified among their companions. Although the DEA has been at the forefront of the implemented species, the main common matter in all these simulations is still characterized as late convergence. However, LM, which is often subjected to local minima drawback, has become surprisingly noteworthy when it is coupled with PSO algorithm. While the PSO having diverse inertia-weighted variants is the last in the overall ranking, the hybrid HPSO supplemented by LM is apparently much more prominent. Thus, the significance of the response obtained from the combination of PSO and LM algorithms has supported the original aspect of the study. Similar to the fact that the oxygen, which triggers the combustion, is transformed into the water with an extinguisher character after it is combined with hydrogen; the metaheuristics have become fairly robust structured when subjected to derivative algorithms. A similar deduction was also obtained with the HDEA and even approximate results were produced with the HPSO4 using CRW. Moreover, it is found more attractive to utilize the HPSO4 algorithm that is minimalist one in point of control parameters compared with the other hybrid approaches (e.g., Karahan et al. 2013; Liao et al. 2017) existed in the hydrological modeling literature. This outlook was supported by the implementation of ANOVA.

In summary, the advantages and possible limitation of the proposed hybrid algorithm are listed below:

• It provides fast convergence with less population and an iteration number.

• The variability in fitness values between runs is very minimal (its stability feature).

• It constitutes confidence in various modeling usages, as the uncertainty of the algorithm's control parameters is minimal.

• The hybridization style argued by HPSO can be easily accomplished with another metaheuristics like HDEA.

• However, due to the fact that use of the metaheuristic process and LM together in the same iteration step increase the number of function calls, and also the Jacobian matrix calculation is made, depending on the computer processor, the solution time may be relatively long.

On the other hand, as both examined models have almost similar structures, adaptation of this hybrid algorithm to other hydrology problems will be useful. The parameter estimation of Muskingum flood routing models, the calibration of rainfall–runoff models with intensive parameters (e.g., Zhang et al. 2009; Nourani et al. 2011; Nourani & Zanardo 2014) and geomorphological rainfall–runoff models (e.g., Nourani et al. 2009), reservoir operation optimization and training weights of artificial neural network-based streamflow prediction models are some of these applications.

Besides, since only the classic SSE objective function was used in the study, an alternative calibration scheme may also be needed to enclose the optimization of multiple objectives that quantify different aspects of the runoff series. For example, when the objective function was taken as [1 − NS], although algorithm grading did not change, there were minor alterations in the performance indexes of the CRR models (these results were not shared due to space limitations). In this sense, in the future works, a prospective automatic scheme with HPSO is planned to be formulated that regards the CRR model calibration issue in a multi-objective framework consisting of additional numerical performance statistics such as NS, and the average RMSE of peak and low runoff events (see Madsen 2000).

Furthermore, the uncommon use of some modern variants of evolutionary algorithms such as modified differential evolution with pb crossover, successful history-based DEA with linear population size reduction and genetic learning PSO in CRR modeling is also remarkable (e.g. Piotrowski et al. 2019). A wide comparison study including such methods will be one of our topics in the future as well.

This study was funded by the scientific research projects Department of Balikesir University under Grant No. BAP2017/145. Also, the authors are grateful to two anonymous reviewers for providing valuable comments and suggestions, which we considered to improve the quality of this paper.