Review and comparison of performance indices for automatic model induction

Increasing demand for integrated environmental modelling and diversity of hydrological modelling approaches motivates the development of unified modelling frameworks (Clark et al. 2015b) that integrate, compare and evaluate multiple model representations. Unified hydrological modelling frameworks can be implemented by defining a general set of conservation equations and their numerical implementation, model structure units/reservoirs, flux parameterizations, model complexity that governs granularity of process representations, spatial discretizations and hydrologic connectivities. In this context, flexible multi-model frameworks such as Rainfall-Runoff Modelling Toolbox (RRMT) (Wagener et al. 1999), Framework for Understanding Structural Errors (FUSE) (Clark et al. 2008) and multi-component frameworks such as FARM (Euser et al. 2013), Superflex (Fenicia et al. 2011, 2014; Kavetski & Fenicia 2011; van Esse et al. 2013) and SUMMA (Clark et al. 2015a) have been developed. In earlier applications, selection of the optimal model configuration was carried out based on hypothesizing model representations using prior knowledge of interactions among observed data and field observations, followed by evaluation of all possible combinations of the selected hypotheses using Bayesian framework, which is computationally intensive. On the other side, with advancements in computing power, use of machine learning tools is gaining popularity and data are becoming drivers for automatic optimal model selection, formulating modelling strategies. One recent study has presented the application of data mining algorithm Automatic Model Configuration Algorithm (AMCA) (Vitolo 2015) to identify the most suitable model configuration using minimum data and FUSE as sample model inventory. The success of an ideal data-driven multi-component modelling framework is not only to perform automatic selection from certain hypotheses, but also to search an entire model space extensively to evolve the most suitable combinations of conceptual model components and generate novel modelling philosophies to fill knowledge gaps. Evolutionary computation (EC) techniques (Goldberg 1994; Babovic et al. 1995, 2017; Babovic & Abbott 1997; Babovic & Keijzer 2000) mimic the mechanisms of natural evolution to provide global or near global solutions to real world optimization problems. Literature shows the successful application of EC techniques in modelling hydrological processes of complex and chaotic nature in application areas, namely, rainfall–runoff prediction (Khu et al. 2001; Babovic & Keijzer 2002), settling velocity estimation (Babovic et al. 2001), water quality (Jeong et al. 2003), sediment transport (Kizhisseri et al. 2006), model error correction (Zechman & Ranjithan 2007), reservoir operation (Rani & Moreira 2010), water distribution systems (Savic & Walters 1997; Bi et al. 2016), etc. Earlier studies involve the use of EC techniques for evolving empirical, grey box models for short-term prediction rather than harnessing the idea of automated knowledge evolution and improving our current understanding of the system. Recent studies have attempted to incorporate hydrological concepts into the framework of an evolutionary data-driven algorithm, say, genetic programming (GP) (Koza 1992) to automatically estimate the input (rainfall) history for the prediction of runoff (Havlíček et al. 2013) and to automatically evolve Sugawara tank (Sugawara 1979) model configurations (Chadalawada et al. 2017) for the catchments of interest.

In this study, GP is used in conjunction with Superflex titled Evolutionary Superflex framework (http://scholarbank.nus.edu.sg/handle/10635/138678) for automatic evolution of optimal model configuration for the given catchment. As mentioned earlier, Superflex (Fenicia et al. 2011; Kavetski & Fenicia 2011) is a multi-component hydrological modelling framework that consists of building blocks that are customized into spatially lumped flow network configurations. The building blocks are reservoirs, flux and lag functions, connection elements. Twelve submodels (M1 to M12) (Figure 1) from Superflex framework (Fenicia et al. 2014) of varying complexities along with basic algebraic functions form the function set of GP. The selected submodels consist of five different types of reservoir units (interception (IR), unsaturated (UR), riparian (RR), fast reservoir and slow reservoir (FR and SR)) with serial, linear and parallel connections as illustrated in Figure 1. The constitutive functions (e.g., Monod, Modified logistic, Linear, Power, Reflected Hyperbolic, Tessier) governing storage–discharge relationships and evaporation losses, shape of lag functions (e.g., Triangle lag scheme), connection elements and ranges of associated parameters, which include, correction factors to inputs, splitters, coefficients of storages and outflows, smoothing factors, lag time are defined based on literature (Fenicia et al. 2011). The selected submodels represent a subset of possible conceptual model constructs of Superflex framework. Submodels M1 and M2 contain a single reservoir unit varying in terms of functions governing storage–discharge relation. M2 to M6 contain reservoir units with serial connections differing in terms of number and type of reservoir units, lag components, flux functions and associated parameters. M8 has two reservoir units in parallel with a splitter (D) dividing the incoming flux precipitation (P) between them. M9 to M12 represent the complex structures with either one (UR) or two reservoirs (IR and UR) preceding the parallelly connected FR and SR. The unsaturated zone reservoir is initialized to a percentage of its maximum capacity. All other reservoirs are assumed to be empty initially. The inputs to GP framework include observed forcings Precipitation (P) and Evaporation (E), Discharge (Q) and random constants. Other GP settings are listed in the Appendix (Table A1, available with the online version of this paper). For more details on GP framework consult (http://scholarbank.nus.edu.sg/handle/10635/138678). In Figure 1, rectangular boxes represent storage units and variables within parentheses represent parameters of the constitutive functions.

The model identification and its performance evaluation should be based on ability to reproduce hydrological behaviour and capture properties of catchments rather than merely matching observed and simulated responses. There are a variety of performance indices that measure the ability of a hydrological model to reproduce observed streamflow. It is important for the simulated time series to preserve critical aspects of observed streamflow hydrograph (Vis et al. 2015) so that the governing model structure can also be used as a catchment classifier for predictions in ungauged basins. One of the recent studies (Shafii & Tolson 2015) highlights the point that the results are more sensitive to objective function formulation as opposed to the choice of optimization algorithm. Therefore, in this paper, a range of entropy, statistical and hydrological signatures listed in Table 1 are used as optimization objectives of Evolutionary Superflex approach.

In this review, a single objective optimization scheme is adopted to clearly assess each of the 14 objective functions highlighted in Table 1. The objective functions defined consist of one criterion or combined criteria representing different aspects of streamflow. This includes three one criterion (NS0, MD0 and RD0) and the rest as combined criteria objective functions.

The main characteristics of the metrics used in this study can be summarized as follows:

1. Nash–Sutcliffe efficiency (NSE) and mean absolute relative error (MARE) are sensitive to peak and low flows, respectively.

2. Modified NSE and logarithmic NSE emphasize low flows.

3. Mean absolute error (MAE) measures the agreement between the average simulated and observed catchment runoff volume and root mean squared error (RMSE) measures the overall agreement of hydrograph shape.

4. Volumetric efficiency (VE) measures the accuracy in prediction of overall runoff volume over the entire simulation period.

5. Correlation coefficient (r) measures the success in replicating overall timing and magnitude of discharge but disregards differences in absolute values.

6. Kling–Gupta efficiency (KGE) combines correlation, bias and variability into one objective metric and is proven to have overcome the problems associated with NSE.

7. Persistence index (PI) aims at comparing the performance of a model against a simple model using the observed value of the previous day as the prediction for the current day as opposed to a constant mean. This is based on the assumption that the process is a Wiener process, i.e., variance increases linearly with time.

8. Index of agreement is based on a two-part squared distance measurement in which the absolute difference between model and simulated value and observed mean is added to the absolute difference between observed record and observed mean. Modified index agreement (MD) involves replacing the squared (power = 2) differences that are sensitive to extreme values by different values of power. On the other hand, relative index of agreement (RD) quantifies the differences between observed and predicted values as relative deviations which reduce the influence of absolute differences during high flows and enhance the influence of absolute low flow differences.

9. Signature measures of variability of runoff generation mechanisms (runoff coefficients and recession curves) are believed to capture significant catchment information that is not conveyed regularly using least squares or conventional variance-based measures of fit. Flow duration curve (FDC) can be considered a hydrologically informative signature of catchment behaviour. Entropy-based descriptors are more sensitive to subtle differences in streamflow (Pechlivanidis et al. 2014). Conditioned entropy difference (CED) is a diagnostic FDC signature representing static information of flow frequency distribution. CED is temporally insensitive and is said to result in better performance when combined with measures that characterize temporal dynamics such as KGE.

The main objective of this work is to evaluate the ability of different objective functions in bringing out the best performance of the Evolutionary Superflex approach thereby resulting in simulated response that adequately preserves observed flow characteristics. The outline of the remainder of this paper is as follows. The next section provides a description of the synthetic dataset and Evolutionary Superflex simulations. The results and analysis of the simulations follow and finally the key conclusions are summarized.

Data used in this study are collected at headwater catchments Attert Basin, Luxembourg (Fenicia et al. 2014), over a period of five years from 1 September 2004 to 31 August 2009, covering a wide range of runoff events. Data consist of daily time steps of precipitation (P) and potential evaporation (E) in mm/h. Two synthetic streamflow time series are generated using the observed forcings (P and E) of two catchments, namely, Huewelerbach and Weierbach. Huewelerbach and Weierbach catchments have contrasting geology dominated by sandstone-marly and schists, respectively. Huewelerbach is a comparatively larger catchment (2.7 km²) with landuse consisting of forests and agricultural land, while Weierbach (0.42 km²) is a fully forested catchment. Runoff behaviour varies between stable, low reactivity of Huewelerbach catchment with runoff coefficient (R_c) of 0.3 and strong seasonality of Weierbach catchment with R_c of 0.98 in wet and 0.095 in dry periods, respectively.

Sum of Standardized Indices (SIS) is a combined measure of signature indices derived from FDC (Yilmaz et al. 2008) that distinguishes the merits of the evolved model alternatives. FDC signature indices considered in the formulation of SIS include high flows corresponding to extreme rainfall events (FHV (FDC high-segment volume): probability of flow exceedance <2%), intermediate flows between high and medium runoff (FMV (FDC intermediate-segment volume): probability of flow exceedance between 2% and 20%), medium flows covering the mid slope of FDC corresponding to vertical moisture distribution (FMS (FDC mid-segment slope): flow exceedance probability between 20% and 70%) and low flows corresponding to base flows (FLV (FDC low-segment volume) flow exceedance probability >70%). These indices identify the portion of hydrograph where bias is the highest. SIS compares the deviation of simulated FDCs from observed FDC weighting all parts of FDC equally. Negative SIS values indicate above average performance and the model configuration with the lowest SIS value is the one that performs the best for the given catchment.

The equations of resultant model configurations induced by Evolutionary Superflex simulations using two synthetic case studies are presented in Tables 2 and 3, respectively. It is evident that the respective underlying model structures are induced correctly for all objective functions except for CED|KG10 under Q_MIV.

M3 and M4 submodels consist of two serially connected reservoirs (UR and FR) in which precipitation enters UR, whereas its excess storage is routed through FR (Figure 1). The outflow from UR is a power function of storage in M4 rather than a threshold function in M3 and CED|KG10 fails to capture this difference.

Other resultant model configurations show differences in terms of parameter values. Further investigation is required to identify the objective functions and other settings that enable Evolutionary Superflex to efficiently search model parameter space and find the optimum.

In terms of identifying the underlying model structure, the objective functions that are able to find the right structures (termed true positives) in at least five out of ten independent runs are noted. Comparing Tables 2 and 3, the combined objective function Borsanyi has performed consistently well in evolving the correct model structure as the one with minimum fitness in 50% of the total runs considering both synthetic case studies of varying complexities.

The most commonly occurring false positives (wrong model structures established as the ones with minimum fitness in respective runs) in the case of simulations using Q_MIV and Q_MXII are M9, M3 and M9, M10, respectively. The objective functions resulting in less than five true positives and more than two different false positives in the case of Q_MIV and Q_MXII simulations are Dawson, Mai0, NS0, Vis_C2 and CED|KG10, RD0, respectively. The top ten model configurations (in terms of training fitness) resulting from 50 generations of ten independent runs (500 generations) using each of the 14 objectives are analysed. It is found that the objective functions which identified the right model structure all ten times are Dawson, Mai0, Madsen, Price, RD0, Vis_C1, Vis_C2 (7 out of 14) and Borsanyi, Dawson, KG10, Mai0, MD0, Madsen, NS0, Price, RD0, Vis_C1, Vis_C2, Vis_C3 (12 out of 14) in the case of Q_MIV and Q_MXII simulations, respectively. The objectives other than those mentioned above fail to adequately distinguish between hypotheses (submodels) of the model structure space in the respective case studies. It is evident that greater efficiency in searching model structure space is observed in the case of M12 as compared to M4. This might be due to the fact that M12 is the most unique submodel of GP function set consisting of four reservoirs with both serial and parallel connections, whereas M4 is reflected as not significantly different from M3 (differs from M4 only in terms of UR outflow function) and M9 (differs from M4 in terms of UR outflow function and excess UR storage being routed through two parallel reservoirs FR and SR).

The variability of parameter values for the best model configurations simulating Q_MIV and Q_MXII are analysed using boxplots shown in Figures 2 and 3. For ease of comparison, all associated parameter values are normalized between 0 and 1. Dots indicate the spread of parameter values and the lines connect the parameters of multiple feasible descriptions (model configurations with corresponding best training fitnesses and negative SIS values) of the target. From Figure 2, it can be observed that the least and the most uncertain parameters of evolved M4 model configurations are Smax_UR (maximum capacity of UR) and β_E_UR (smoothing factor for function governing evaporation loss from UR), respectively. In resultant M4 configurations whose structure consists of serially connected UR and FR, greater variation is observed in the parameters of flux functions of UR, namely, SiniFr_UR (Intial UR storage = SiniFr_UR*Smax_UR), β_E_UR and β_Qq_UR (outflow function parameter).

Figure 3 shows that the least and the most uncertain parameters of the resultant M12 configurations whose structure consists of serial network of IR and UR before parallelly connected FR and SR are K_Qq_FR and D_S, respectively. K_Qq_FR denotes hydraulic conductivity of outflow from FR and D_S is the splitter that diverts excess storage of UR to FR and SR, which implicitly governs the outflows of FR and SR. Greater variation is visually evident in parameters of IR, UR and SR of M12 configurations.

In GP runs simulating Q_MIV and Q_XII, the lowest variation is observed in the values of evaporation multiplication factor (Ce) and FR parameters (α_Qq_FR, K_Qq_FR).

Evolutionary Superflex approach shows greater uncertainty in filtering out combinations of parameters from model parameter space given that inherent model structure is accurately found at most instances. Beven (2006) suggests that this uncertainty should be viewed as a problem of decidability between different probable descriptions of the working of hydrological systems. Sources of uncertainty include model structure (process representations), numerical solution approximations and parameter values. Many different parameter sets within a chosen model structure may be acceptable in reproducing the observed behaviour of that system and this concept is called ‘Equifinality’.

Literature suggests that uncertainty estimation techniques allowing for equifinality can be broadly classified as formal (e.g., Bayesian) limited by requirements on prior knowledge and informal (e.g., generalized likelihood uncertainty estimation, GLUE) limited by subjectivity but simple (Vrugt et al. 2009). Here, GLUE, a useful paradigm that avoids overconditioning and estimates prediction quantiles is chosen for uncertainty analysis. All Evolutionary Superflex configurations generated using 14 objective functions that have rightly induced the inherent model structures M4 and M12 for Q_MIV and Q_MXII, respectively, are listed and the corresponding parameter sets are chosen for uncertainty analysis using GLUE. From the chosen parameter sets belonging to behavioural model space, elimination of those that do not perform adequately is carried out based on SIS criterion (likelihood measure) and the threshold is set as 0. More weight is assigned to parameter sets that result in negative SIS values and the prediction boundaries are constructed as shown in Figures 4 and 5. The percentage of observed values falling outside the prediction limits are observed as 1.49% and 3.58% for Q_MIV and Q_MXII, respectively. The width of the prediction boundary for Q_MIV (0.67) is smaller than for Q_MXII (0.99), which can be attributed to higher complexity of inherent model (M12) of Q_MXII. It can be seen in Figures 4 and 5 that the observed values closely follow the 50% quantile (median) of the simulated values. It has also been verified that uncertainty bands derived using parameter sets given by Evolutionary Superflex are narrower than using parameter sets randomly drawn from uniform distributions (with user specified parameter ranges) for the same likelihood measure (SIS).

In order to select the most appropriate, the best model configurations (Tables 2 and 3) of both the synthetic case studies are evaluated with respect to biases in FDC signatures. The FDCs (Figures 6 and 7) reflect the differences in the behaviour among the best configurations of Q_MIV and Q_MXII with negative SIS values.

Figure 6 illustrates that M4 model configuration evolved based on Vis_C1 with the least SIS value (−3.1) matches most of the Q_MIV characteristics showing bias in FHV and negligible bias (less than ⁠) with respect to other FDC signatures (FMV, FMS and FLV) (see subplot of Figure 6 entitled ‘Vis_C1’). KG20 (SIS = −2.69) and Madsen (SIS = −2.16) also resulted in good overall approximation of the target with slightly greater deviations in FMS and FMS, FLV respectively (see subplots of Figure 6 entitled ‘KG20’ and ‘Madsen’), when compared to the best objective Vis_C1. Although the wrong model structure M3 (differs from M4 in terms of UR outflow function) identified as the best by CED|KG10 has a negative SIS value, it results in the poorest performance with respect to low flows (see subplot of Figure 6 entitled ‘CED_new’).

Figure 7 shows that M12 model configurations evolved based on Madsen (SIS = −3.05), Vis_C3 (SIS = −2.83) and Vis_C2 (−2.49) perform well in comparison to others. The FDC signatures of Madsen, Vis_C2 and Vis_C3 show that high and medium flows are better approximated by M12 model configuration evolved based on Vis_C3 in comparison to Madsen and Vis_C2. Intermediate flows are well captured by M12 model configuration evolved based on Madsen. Considering low flows, the lowest value of FLV bias is observed in the case of Vis_C2.

Comparing Figures 6 and 7, it is evident that the optimal model reported by Evolutionary Superflex approach in the case of simulations using Q_MIV is superior (close match to the target both in terms of structure and parameters) to those using Q_MXII. This can be attributed to higher complexity of M12 submodel both in terms of structure and number of parameters to be optimized. It can be concluded that combined objective functions, say, Vis_C1, KG20 and Madsen for Q_MIV simulations and Madsen, Vis_C3 and Vis_C2 for Q_MXII simulations resulted in Evolutionary Superflex performing better compared to others. Madsen fitness metric is observed as the common best objective irrespective of the complexity of model search space.

Tables 4 and 5 show the performance evaluation of the best model configurations of Evolutionary Superflex simulations using Q_MIV (Table 2) and Q_MXII (Table 3), respectively, based on hydrological efficiency criteria KGE1, KGE2, NSE and VE, correlation (r) and numerical accuracy metrics MAE, MD, RMSE and PI. No considerable variation with respect to SUSE, KGE1, KGE2, NSE, MD and r is observed (except in the case of RD0 for Q_MXII simulations) making them unfit to decide between different system representations. Very high values of Pearson correlation coefficient (r) (see Tables 4 and 5) associated with all the best configurations indicate the absence of phase errors and good agreement of observed and simulated means and variances.

In Table 4, the low RMSE values of 0.055 and 0.066 correspond to the best M4 model configurations induced based on Vis_C1 and Madsen, respectively, which is in agreement with the decision on the most suitable configurations made using SIS criterion. On the other hand, the highest RMSE and the lowest VE (see Table 4) observed in the case of model configuration induced based on CED|KG10 for Q_MIV is not to be considered, although it has a negative SIS value. KG20 is associated with better MAE (1.6 × 10⁻³) than Madsen (0.013) while RMSE of Madsen (0.066) is superior to KG20 (0.167). Both SIS criterion and performance evaluation based on a wide range of metrics indicate that M4 configurations evolved based on KG20 (QGP_KG20), Madsen (QGP_Madsen) and Vis_C1 (QGP_Vis_C1) can be considered the most suitable for Weierbach catchment synthetic dataset used in this study (see Figure 8).

Table 5 shows that the best M12 model configurations evolved based on CED|KG10 (SIS = 1.25) and RD0 (SIS = 6.18) with low PI values of 0.63 and −0.16, high RMSE values of 0.2 and 0.34, low VE values of 0.76 and 0.66, respectively, can be ignored. The best configurations evolved using the objective functions Borsanyi, Dawson, MD0, Madsen, NS0, Vis_C2 and Vis_C3 have better fitness in comparison to the rest (Table 5). Among the above mentioned, all configurations have negative SIS values except for those evolved using Dawson (SIS = 1.03) and MD0 (SIS = 1.45). As the complexity of the system increases, the model selection decisions based on FDC signature indices and other performance evaluation metrics do not match each other completely. In such cases, the choice based on SIS is considered superior as it equally weights different flow aspects. Thus, M12 model configurations evolved using Madsen (QGP_Madsen), Vis_C2 (QGP_Vis_C2) and Vis_C3 (QGP_Vis_C3) are considered the most suitable for Huewelerbach catchment synthetic dataset used in this study (see Figure 9).

This review provides an assessment of abilities of different objective functions that enable Evolutionary Superflex modelling framework to identify the most suitable model configurations (structure + parameters set) that reproduce several target flow regime signatures simultaneously. As expected, model configurations evolved based on one criterion objective function do not ensure good agreement with respect to others. NS0 has performed better than the other one criterion objective functions used in this study. Combined objective functions have delivered good performance across a range of flow characteristics. The proposed approach performs well with respect to refinement of model structure space using objective functions defined in this study. Greater variation is observed in filtering out feasible parameter combinations from model parameter space. The uncertainty bands of simulated model parameters are quantified using GLUE approach and it is found that observed values closely follow the median of simulated values. The uncertainty in the prediction of parameter values increases as the complexity of the structure increases. Among the five types of reservoirs (IR, RR, UR, FR and SR) that constitute the structures of submodels employed in this study, the parameters of FR show the least variation. The best configurations of Evolutionary Superflex approach evolved using each of the 14 objective functions can be considered as multiple descriptions of the system under study. Deciding on the most suitable configurations among them is done using SIS criterion and other selected performance evaluation metrics. Each performance evaluation metric measures one aspect of the difference between observed and simulated responses. Some of the metrics, say, SUSE, KGE, NSE, MD and r assign similar values to different hydrographs and are not suitable for selection of the most suitable configurations. Difference is observed between selection of suitable configurations based on SIS criterion and other performance evaluation metrics, as complexity of the underlying system increases. The model configuration selected based on SIS criterion overrules in case of disparity as SIS equally weights (different portions of the hydrograph). The combined objective functions, say, Madsen, KG20, Vis_C1 and Madsen, Vis_C2, Vis_C3 are found to be the best for induction of inherent models of Weierbach and Huewelerbach catchments synthetic datasets, respectively. Madsen metric that equally weights agreement of overall volume, overall shape of system response, high and low flows is found as the common best objective for both case studies, irrespective of the complexity of search space. It remains a challenge to pinpoint an objective function that brings out the best performance of Evolutionary Superflex approach in effectively searching the model space to evolve optimal conceptual model constructs and parameters for given datasets. This study aims to run the evolutionary modelling framework using different, strategically aggregated objective functions and to select the most suitable model based on the collective analysis of results obtained.

This research will be extended to capture hydrological processes of larger and complex catchments in order to reinforce the current findings. An attempt can be made to constrain equifinal parameters to ranges associated with observed catchment behaviour using identifiability analysis which can result in narrow uncertainty bands. Also, the efficiency of the proposed modelling approach will be tested in a multiobjective optimization context.

The authors thank Dr Laurent Pfister, Luxembourg Institute of Science and Technology (LIST) for granting the permission to use Luxembourg catchments dataset for this research study. Additionally, we are grateful to Dr Fabrizio Fenicia for his assistance with Superflex framework.