Automated calibration of a two-dimensional overland flow model by estimating Manning's roughness coefficient using genetic algorithm

Hydrological modelling results vary greatly depending on the values assigned for the parameters chosen for model simulation. Due to complexities of physical processes in hydrology, the exact values of empirical parameters in such models are often uncertain. Manning's roughness coefficient, denoted by ‘n’, is one of the empirical parameters, whose value is inevitable for flow modelling in open channels. It represents the amount of frictional resistance applied by the channel to the flow. Besides channel flow, surface runoff estimation due to single rainfall events, a key factor in drainage system design, also needs the estimation of Manning's roughness coefficient. In general, tools that are commonly used for determination of runoff hydrographs from rainfall excess also require appropriate values of ‘n’. Conventionally, the roughness coefficient is selected from the available literature for various channel conditions based upon the number of controlling factors (Chow 1959; Yen 1991). An average value for Manning's roughness coefficient is assumed from the specified range.

Parameterization is done through optimal control theories that are capable of minimizing the discrepancies between the measurements and computations. Minimization of these discrepancies is mathematically expressed as an objective function. The parameter identification then turns out to be an inverse problem for searching the minima or maxima of an objective function, by considering the value of the parameter as being dependent on the governing equations. In the modern era, intelligent techniques such as fuzzy logic, artificial neural networks, neuro-fuzzy and genetic algorithms (GA) are being used effectively in many such optimization processes. The advantage of using such methods is that the whole processes of model simulation and calibration can be automated through integrated approaches.

In hydrology, models may have numerous local optima on their objective function surface. Hence, it is well and good to use a global optimization method instead of the local method (Sorooshian et al. 1993) for their automatic calibration. Global optimization methods can efficiently search optimum parameter solutions that can minimize or maximize objective functions from a population. Various population-based evolutionary global optimization algorithms applied in hydrology are GA (Wang 1991), shuffled complex evolution algorithm (Duan et al. 1992) and simulated annealing (Sumner et al. 1997). Jeon et al. (2014) recommended GA for quicker convergence of model parameters. Also, numerous studies have observed the utility of GA for calibration of rainfall–runoff events (Franchini 1996; Cooper et al. 1997; Franchini & Galeati 1997; Wang 1997; Ndiritu & Daniell 2001; Cheng et al. 2006; Liu et al. 2007; Lim et al. 2010; Kamali et al. 2013; Shinma & Reis 2014; Xu et al. 2016). Apart from calibration of hydrological models, the technique is widely applied in pipe flow optimization processes (Dandy et al. 1996; Chamani et al. 2013; Olszewski 2016), pollution source identification processes (Ritzel et al. 1994; Yang et al. 2008; Bu et al. 2013; Cantelli et al. 2015; Hu et al. 2015), groundwater monitoring network design (Babbar-Sebens & Minsker 2010) and optimal management of coastal ground water (Ketabchi & Ataie-Ashtiani 2015). The main advantage of a GA is that it takes the optimal control theories for parameter estimations and can provide realistic parameter sets. Zhang et al. (2015) performed optimization process using multi-objective calibration of the SHETRAN hydrological model. Shafii et al. (2015) also worked on hydrologic model calibration via model pre-emption. VanGriensven & Bauwens (2003) presented a probability scale approach for calibration of water quality catchment models. Xu et al. (2013) compared three global optimization algorithms for model calibration for a large watershed.

Literature review suggests the possibility of estimating optimal value for Manning's roughness parameter through a recursive parameter updating approach by using a GA. A reliable calibration process must include: (1) a criterion to be defined to evaluate model performance; (2) an appropriate algorithm to minimize or maximize the objective function; (3) model validation using calibrated parameter against new datasets. Thus, starting from these considerations, the article presents an automated calibration procedure for a two-dimensional overland flow routing model. The specific objective of this paper is estimating Manning's roughness coefficient using an automated single event model integrated calibration framework by applying techniques of the GA. The calibration problem is formulated as an optimization problem and the GA tool in MATLAB is applied. A Python subroutine links the model and GA to allow the data transfer between the two programs. The methodology adopted can be used as a viable alternative to manual approaches in any parameter estimations and can be explained as a best tactical approach in the calibration of single event hydrological models.

This section explains the details of the model used, datasets for the model, GA details, objective function formulation and two-stage model optimisation frameworks using model integrated GA. The two-stage optimization framework includes model integrated single objective genetic algorithm (SOGA) and model integrated multi-objective genetic algorithm (MOGA).

The hydrodynamics of overland flows are modelled as two-dimensional shallow flows using diffusive waveform of St. Venant's equations. The linearized form of the partial differential equation is solved numerically by applying weighted implicit finite volume method (FVM). FVM techniques are preferred because they enforce volume conservation across all control surfaces and the discretized governing equations retain their physical interpretation. A finite volume-alternating direction implicit scheme proposed by Lal (1998) is used in this model code. The model code written in C language has inputs of rainfall intensity, time of rainfall, slopes in x, y directions and time steps to simulate the runoff. This model is selected for the present study because the only parameter that requires calibration is Manning's roughness coefficient.

Datasets obtained from the published works of Maksimovic & Radojkovic (1986) and Xiong & Melching (2005) are used for the present study. Both datasets compiled by the above researchers are based on the Watershed Experimentation System (WES) of Chow & Yen (1974). WES is a laboratory urban catchment at the University of Illinois and is designed to study only the surface runoff on an impervious surface. Longitudinal and lateral slopes are adjustable in the WES. Figure 1 represents the simplified watershed laboratory of WES, where S_x and S_y are overland flow slopes along x and y directions, respectively. L represents length and W represents the width of the WES watershed. WES used a watershed size 12.2 m by 12.2 m.

From the obtained datasets, the storm events that reached equilibrium discharge are used for calibration and storms that did not achieve equilibrium discharge are used for validation in the present study. Thus, the experiments having code US2L005, US2L006, US2L009, US2L010, 613, 634, 806 and 816 are selected for parameter estimation. The remaining four events with code 821, 624, 615 and 822 are used for validation purposes. Table 1 shows the geometric characteristics of the watershed and rainfall parameters for selected experiments. The time offset bias corrections (McCuen et al. 2006) based on equilibrium discharge are applied on the experimental data for calibration.

GA, the search and optimization tool based on the ‘survival of the fittest’ scenarios, gained momentum due to the flexibility and adaptability in different environments and conditions (Holland 1992). Bäck & Schwefel (1993) provided the conceptual algorithms for the canonical GA. Enormous resources are available to understand the essential components of GA (Goldberg 1989; Reed et al. 2000; Haupt & Haupt 2004; Nicklow et al. 2009). The main difference between GA and most of the traditional optimization methods is that GA uses a population of points at one time in contrast to the single point approach by traditional optimization methods. GA solver in MATLAB is used as the tool for parameter estimation in the current study.

Hydrological and hydraulic models use different efficiency criteria as a fitness function for optimization which quantifies how much goodness of fit for each potential solution. Green & Stephenson (1986) gave criteria for comparison for single event models. Accordingly, three different performance measures selected for the present study are percent bias (PBIAS), root mean square error (RMSE) observations standard deviation ratio (RSR) and Nash–Sutcliffe efficiency (NSE). Krause & Boyle (2005) suggested that a combination of dimensionless model evaluation statistics and error index evaluation statistics are necessary for a scientific sound model calibration. Here, dimensionless model evaluation statistics is represented by NSE and error index evaluation statistics is represented by PBIAS and RSR. Accordingly, three fitness functions used are described below.

First, a single objective GA with minimization of PBIAS as objective function is evaluated. For this, a set of single event experimental data is chosen. The input parameters like rainfall intensity, time of rainfall, slopes in x, y directions and time steps are prepared separately for the event. The Manning's roughness coefficient is the parameter to be estimated through GA. It is separately fed through another input file, since it needs to be automatically updated every time as the GA gets initiated. The overland flow routing model code is then called via a subroutine in Python. The model runs through these inputs and simulated runoff is saved into another file. Data files containing observed runoff datasets are called for comparison.

Once fitness values are evaluated, GA will retain the parameter with best fitness value for the current generation. If the fitness value satisfies the convergence criteria, then the process is stopped and the best Manning's value is updated. Otherwise, the next generation is produced by genetic operators and the whole process is continued. Thus, the entire calibration framework gets automated. Later, model integrated SOGA is performed by minimizing RSR and maximizing NSE for the same event. Default GA parameter values in the MATLAB solver are chosen for the present study. The entire methodology is explained in Figure 2. The process is repeated for eight single events.

In a multi-objective context, model calibration can be performed by different criteria. In hydrological modelling, multi-site, multi-variable and multi-response modes are the different kinds of criteria to be considered while using multiple objectives for calibration (Madsen 2003). In this study, since a single parameter is the matter of concern and spatial variability of the experimental watershed is very minimal, only multi-response modes need to be analysed. Hence, multi-response modes, i.e., objective functions that measure various responses of the hydrological processes (PBIAS, RSR and NSE) are considered for the MOGA.

If bias is an objective function, one may find a set of parameters that provides a very good simulation of volume, but a poor simulation of the hydrograph shape or peak flow, and vice versa. RMSE-based observations and NSE tend to emphasize the high flows, and consequently, are oversensitive to extreme values and outliers (Legates & McCabe 1999). On the opposite, the mean absolute percent error tends to emphasize the low flows (Yu & Yang 2000).

are different objective functions. The optimization problem is constrained in the sense that is restricted to the feasible parameter space ⁠. However, since second and third objective functions are not conflicting, it is not necessary to consider all the objective functions simultaneously. Thus, two separate sets of model integrated MOGA are performed using these objective functions. Minimization of PBIAS, RSR (Equation (6a) and (6b)) make the first set and minimization of PBIAS and maximization of NSE (Equation (6a) and (6b)) constitute the second set. MATLAB toolbox (Mathworks 2014) is used for MOGA analysis.

The solution for conflicting multi-objectives will not be a unique set of parameters but will consist of the so-called Pareto set of solutions according to the trade-offs between the different objectives. Pareto-fronts are generated individually for each event. The process is repeated for another set of objective functions, and ⁠. This kind of separate evaluation is performed due to the non-conflicting nature of and ⁠. The flowchart for the entire process is presented in Figure 3.

All these processes are automated in the present study. The parameter obtained by integrating all the events is the balanced measure of that set.

Model integrated SOGA and MOGA frameworks have been autocalibrated for finding Manning's roughness coefficient value for the overland flow model. The simulated results are then compared with Xiong & Melching (2005) for the same rainfall events.

Series of tests are carried out to optimize the roughness coefficient through GA integrated model framework. To optimize the parameters, eight single rainfall events are used for calibration purposes. Three fitness functions are evaluated separately. The obtained Manning's roughness value is close to around 0.01 (when rounding off to two digits) from a specified lower limit of 0.001 to an upper limit of 0.2. Table 2 shows the results of such parameterization process obtained after model integrated SOGA for each eight events. Corresponding objective function values obtained are presented within parentheses. The results show that all the calibrated ‘n’ values are in the range of Chow's predictions. To have a visual comparison, hydrographs are plotted (Figure 4) for all the calibrated events using the single objective optimized average value. It shows a good overall agreement of the hydrograph shape.

The results showed a variation in parameter value when minimizing PBIAS or RSR. When RSR is minimized or NSE is maximized, the optimum parameter value is almost the same for each event. This indicates the non-conflicting nature of RSR and NSE. It is also noticed that there is a difficulty in obtaining the best global maximum or minimum with respect to all objectives. For example, with respect to the first objective function minimization, event 816 shows the best value. At the same time, it is not giving the best result when compared to second or third objective functions. Another event (US2L010) has the minimum value of RSR, whereas a higher value with respect to PBIAS minimization. Moreover, the parameter value also varies depending on the objective function used. In short, no single set of parameter can be obtained through model integrated SOGA. This shows the existence of a trade-off solution which needs a multi-objective optimization process.

For multi-objective GA, two objective functions are evaluated simultaneously using the first level of aggregation. Simulations are carried out as per the pooled statistics of Equation (7), to evaluate the objective functions simultaneously. A trade-off solution between the objectives is obtained and a Pareto-front is generated. Comparative results between PBIAS and RSR as well as between PBAIS and NSE are presented in Figure 5.

If for every member x in a set P there exists no solution y dominating any member of the set P, then the solutions belonging to set P constitute a locally Pareto-optimal set. For a given Pareto-optimal set, the corresponding objective function value in the objective space is called the Pareto-front (Deb 2001). From model integrated MOGA results, there exists no feasible vector of decision variables which would decrease PBIAS without causing a simultaneous increase in RSR. Similarly, there is no parameter which decreases PBIAS without causing a simultaneous increase in NSE. Thus, a single parameter value is not yielded, but rather a set of solutions called the Pareto-optimal set is yielded.

The results of MOGA indicate that the solutions on the Pareto-front are equally good with respect to the first objective, but may not be suitable for the second objective. Another notable feature is that the tails of the Pareto-front represent the results of single objective optimization. The estimated Pareto-front plots represent a trade-off between the objectives. A very good calibration of PBIAS provides a bad calibration of RSR and vice versa. In other words, when PBIAS is high, RSR will be low, and vice versa. Higher NSE solution points have higher PBIAS too.

By analyzing Pareto-front results of multi-objective optimization, it is observed that the Manning's roughness values are different for each single event experiment for a single set itself. Hence, to achieve a single optimum parameter for all the events, one more calibration run is performed by integrating all the events simultaneously. All the eight calibrated events are invoked at the same time through model integrated GA. Figure 6(a) represents the convergence of aggregated function values obtained through model integrated MOGA by considering eight events simultaneously for different sets of optimization. It is noted that after the 15th iteration, the function value converged. 0.0132 is the final Manning's roughness coefficient value obtained during both the optimization processes. The balanced optimum value of the parameter obtained by minimum distance to the origin is 0.0132. To validate this parameter value, one more run is performed for similarly trended events randomly chosen from the experiment datasets. For each Manning's roughness coefficient, corresponding PBIAS and NSE values are plotted in Figure 6(b). From the figure, it may be noted that the plots of PBIAS and NSE cross near 0.0132, the optimal Manning's roughness value stated earlier.

Error analysis is important for quantifying the significant differences between the events. Here, runoff depth, peak runoff and time to peak are calculated for all the events with multi-objective optimized parameter value.

Table 3 denotes the percentage change in the derived parameters from hydrographs to observe the improvement in using the automated calibrated parameter in simulations. Bold values indicate the lowest and highest error observed for the calibration process. When automated calibration integrated model is used, runoff depth is simulated for a minimum error percentage of 0.175%. It is observed for the event US2L005. The maximum error observed in runoff volume is 1.36% (event 634). The highest percentage bias observed in peak runoff and time to peak is 6.45% and 25.22%, respectively. Bold values indicate the lowest and highest error observed for the calibration process.

Average percentage error is the ratio between the sum of absolute error percentage in a criteria obtained for all events to the total number of events in the watershed. Xiong & Melching (2005) simulated the events 613, 634, 806 and 816 using Manning's roughness coefficient obtained by trial and error method. The average percentage errors in runoff depth, peak runoff and time to peak analysed from their study are 1.29%, 2.22% and 20.59%, whereas, the average percentage errors obtained after automated multi-objective calibration process for these parameters are 0.69%, 3.4%, and 13.91%, respectively. Thus, the overall error has reduced considerably when the parameter is obtained through the automated calibration process.

To validate the entire processes, it is necessary to check the model performances with datasets that are not used for calibration. Four separate storm events are simulated by applying the roughness value previously obtained. Figure 7 shows the hydrographs plotted for validation events. PBIAS values of rainfall events 821, 624, 615 and 822, respectively, are 3.91, 7.75, 2.42 and 8.72. The range of RSR values are 0.248, 1.3657, 0.226 and 1.406, respectively. The NSE values observed for these validation events are 0.938, 0.865, 0.977 and 0.874, respectively.

Implementation of automated calibration techniques is one of the major concerns of all hydrological models. The result produced by hydrodynamic models largely depends on the accuracy of flow resistance parameters. In general, these parameters are mainly determined by applying trial and error techniques and are assumed to be more or less the same for model simulations. However, with the advent of artificial intelligent techniques, automated parameter estimation techniques can be performed by integrating model simulation environments. GA is one of the techniques which increases the chances of finding better solutions from a population-based evolutionary algorithm. Hence, the study is performed to estimate Manning's roughness coefficient by integrating GA with the model.

A multi-step approach (model integrated SOGA and model integrated MOGA) is used in the automation of calibration processes by applying GA. A two-dimensional overland event-based rainfall–runoff model is integrated with GA. The model has Manning's roughness coefficient as the only parameter which needs to be calibrated. The model is calibrated on eight single event datasets of the model watershed. The calibration phase results are compared based on their performance with regard to different objective functions. Three objective functions (PBIAS, RSR and NSE) have been tried individually for single stage optimization processes. Visual agreement test is performed for assessing the graphical differences in the simulated hydrographs. The result shows that the calibrated parameter values are dependent on the type of the objective function used. Conflicts are observed between the different objectives and no single set of the parameter can optimize all objectives simultaneously. However, parameter obtained by minimizing RSR and maximizing NSE is found to be the same. Therefore, model integrated MOGA is performed for individual storm events by considering two objective functions simultaneously. The results of MOGA indicate that the solutions that lie on the Pareto-front are equally good with respect to the single objective considered, but may not be suitable for other objectives. Later, a balanced aggregated objective function is formed to combine all the events to obtain single optimal parameter value. Manning's roughness value is found to be 0.0132 for the final calibration integrating all the events in the MOGA. The quantitative error analysis of hydrograph parameters also shows significant improvement when optimized parameters are used. The error in runoff depth simulation is reduced to 0.175%.

Further, validation is performed on four single event datasets using the optimum parameter value. The validation results show the same level of accuracy as the calibration results, confirming that the balanced optimal parameter value determined is deemed to be good for all the events. The study thus highlights the solution of the calibration problem by recursive parameter updating in the model and checking for the best model predicted solutions in a hydrological environment. In real-world applications, especially in operational frameworks and for real-time predictions, this type of robust calibration becomes highly essential.

The authors wish to thank the financial support from INSPIRE component, Department of Science and Technology, Govt. of India for the present study.