In this application, we demonstrate the application of the UNEEC method to a rainfall-runoff model of the Bagmati catchment in Nepal. The HBV model is also used as a rainfall-runoff model. Compared to the previous case studies, the Bagmati catchment has the following characteristics: the size of the catchment is bigger, the length of the data is larger, the resolution of the data is daily, and the quality of the data is comparatively poorer. In this study we also extend the methodology to estimate the probability distribution function of the model output.

Description of case study

Bagmati catchment lies in the central mountainous region of Nepal. The catchment encompasses nearly 3700 km² within Nepal and reaches the Ganges River in India. The catchment area draining to the gauging station at Pandheradobhan is about 2900 km² (see Figure 7.20) and it covers the Kathmandu valley, including the source of the Bagmati River at Shivapuri and the surrounding Mahabharat mountain ranges. The catchment covers eight districts of Nepal and is a perennial water body of Kathmandu. The length of the main channel is about 195 km within Nepal and 134 km above the gauging station.

 

Time series data for rainfall at three stations (Kathmandu, Hariharpurgadhi, and Daman) within the basin with a daily resolution for eight years (1988 to 1995) were collected. The mean areal rainfall is calculated using Thiessen polygons. Although this method is not recommended for mountainous regions, the mean rainfall is consistent with the long-term average annual rainfall computed with the isohyetal method (Chalise et al., 2003). The long-term mean annual rainfall of the catchment is about 1,500 mm with 90% of the rainfall occurring during the four months of the monsoon season (June to September). Daily flows are recorded from only one station at Pandheradobhan. Long-term mean annual discharge of the river at the station is 151 m³/s but the annual discharge varies from 96.8 m³/s in 1977 to 252.3 m³/s in 1987 (DHM, 1998). The daily potential evapotranspiration is computed using the modified Penman method recommended by the FAO (Allen et al., 1998).

 

Two thousand daily records from 1988/01/01 to 1993/06/22 are selected for calibration of the process model (in this study, the HBV hydrological model as described in section 5.2) and data from 1993/06/23 to 1995/12/31 are used for the validation (verification) of the process model. The first two months of the calibration data are considered as a warming-up period and are hence excluded from the study. This separation of the 8 years of data into calibration and validation data is done on the basis

Figure 7.20. Location map of Bagmati catchment considered in this study. Triangles denote the rainfall stations and circles denote the discharge gauging stations.

Table 7.9. Statistical properties of the runoff data sets

 

of previous studies (Shrestha and Solomatine, 2005; Shrestha et al., 2006). These data sets are used with all the uncertainty analysis methods considered in this study. The some statistical properties of runoff data sets are presented in Table 7.9.

Calibration of model parameters

The HBV model is calibrated automatically using the global optimisation routine – adaptive cluster covering algorithm, ACCO (Solomatine et al., 1999) and followed by manual adjustments of the model parameters as done in the previous case study (see section 5.3). The HBV model parameters ranges used in the calibration procedure and their calibrated values are given in Table 7.10.

 

The Nash-Sutcliffe efficiency (CoE) value of 0.83 was obtained for the calibration period; this value corresponds to a root mean squared error (RMSE) value of 92.31 m³/s. The model was subsequently validated by simulating the flows for the independent validation data set. The CoE was 0.87 for this period with a root mean squared error value of 127.6 m³/s. Note that the standard deviation of the observed discharge in the validation period is 54% higher than that in the calibration period, and this apparently affects the increased performance in the validation period with respect to the CoE. The simulated and observed hydrographs, along with rainfall and simulation error, are shown in Figure 7.21.

Table 7.10. Ranges and calibrated values of the HBV model parameters. The uniform ranges of parameters are used for calibration of the HBV model using the ACCO algorithm and for analysis of the parameter uncertainty of the HBV model by the GLUE method.

 

Figure 7.21. Simulated discharge for the Bagmati catchment for a validation period. The top and bottom plots show the precipitation and model errors respectively during the same period.

Analysis of model residuals

Figure 7.22. Normal probability plot of (a) model residuals in the calibration period and (b) model errors in the validation (or verification) period.

The analysis of the model residuals in the calibration period shows that the model residuals are highly correlated with the observed flows. Most of the high flows have relatively high residuals whereas the low flows have small residuals. The presence of heteroscedasticity in the residuals is observed as well. These results are consistent with those obtained in the Brue case study (section 7.3). Figure 7.22 shows the distribution of the residuals in the calibration and in the validation periods, and it can be seen that the residuals are not normally distributed. The Kolmogorov-Smirnov and Lilliefors (Lilliefors, 1967) tests support this hypothesis.

Clustering

The clustering is performed using Fuzzy C-means algorithm based on previous experience (Shrestha and Solomatine, 2006a, 2008). Selection of the input variables and the optimal number of clusters are discussed in this section.

Selection of input variables

Several approaches have been reported in the literature (Guyon and Elisseeff, 2003, Bowden et al., 2005) to select the model input variables; we follow the approach used in the previous case study (section 7.3). The input variables X_c used in clustering are constructed from the rainfall, the potential evapotranspiration, and the observed discharge. Several structures of the input data including the lagged variables are considered following the analysis of the correlation and the AMI between the rainfall, runoff and evapotranspiration with the model residuals. It appears that the inclusion of the potential evapotranspiration does not improve the results obtained for the cross-validation data set, and it can be said that its inclusion would introduce “confusion” into the model. For example, during the low flow season (i.e. the dry season of April and May) there is a very high potential evapotranspiration due to the hot weather in this period. The calibration of the hydrological model shows that the model captures the low flow reasonably well. However, this hydrological condition (low flow, negligible or zero rainfall, and very high potential evapotranspiration) is not identified as a low flow season in the clustering. So it was decided not to include the potential evapotranspiration as a separate variable, but rather to use the effective rainfall (see equation (5.10)).

 

After some trials with the different input combinations the following variables are selected for clustering: RE_t and Q_t (see Table 7.11). The principle of parsimony is followed by avoiding the use of a large number of inputs. In addition to this, the absolute values of the model residuals are used, which explicitly force the input data having similar values of model residuals to be in one group. Since the rainfall and discharge have different units with different orders of magnitude, and their values do not represent the same quantities, all input variables are normalized to the interval from zero to one. This scheme can prevent the model from being dominated by variables with large values, and is commonly used in data-driven modelling.

Table 7.11 Summary of input variables selected for the models/processes.

Note: P is the partition matrix that contains the membership values of data for each cluster and V is the matrix of the cluster centres.

Selection of number of clusters

Figure 7.23. Finding the optimal number of clusters using partition index SC (top figure), separation index S (middle figure), and Xie-Beni index XB (bottom figure).

Figure 7.24. Sensitivity of the statistics of uncertainty measures with the number of clusters for the calibration period.

As previously mentioned, an important issue in clustering is to identify the optimal number of clusters. Three validation indices are used to select the optimal number of clusters: partition index (SC), separation index (S), and Xie-Beni index (XB) (Xie and Beni, 1991; Bensaid et al., 1996). The results are shown in Figure 7.23. The SC and S indices hardly decrease when the number of clusters c equals 6. However, the XB index reaches this local minimum at c = 5. It must be mentioned again, that the use of a single validation index is not reliable, so three indices are used, and the optimum can only be determined through a comparison of all the results. The partitions with fewer clusters are preferred when the differences between the values of a validation index are minor. So the optimal number of clusters is chosen to be 5, which is also confirmed by estimating prediction bounds for the calibration data (see Figure 7.24).

 

Figure 7.24 depicts the sensitivity of uncertainty measures to the number of clusters c. It is observed that MPI decreases with the increase of c. However, in the case of PICP there is no obvious pattern. The MPI fluctuates around the value 97.5% after c=5. At c=5, MPI and deviation of PICP from the desired confidence level (i.e., 90%) is smaller as compared to those with c=6. This value is also consistent with previous case studies, which have shown that this value is reasonable to represent different situations related to the runoff generation process.

Analysis of clusters

Figure 7.25 shows fuzzy clustering of input examples. The input variables, effective rainfall RE_t (Figure 7.25a) and discharge Q_t (Figure 7.25b), are on the abscissa and the model residuals are on the ordinate. Note that each input data belongs to all 5 clusters with different degrees of fuzziness (see Figure 7.26). However, in the plot the cluster which has the maximum membership function is shown. It is observed that there is a well defined pattern of the model residuals with the input variables such as defined by RE_t and Q_t. One can see that high flows and high (effective) rainfall generally have higher values of model residuals which are well identified (Cluster C3) by the clustering process. On the other hand, the conditions characterized by low flows are also separated into one cluster (Cluster C1) which has very low values of the model residuals. These results are also consistent with the results of the previous case studies.

Figure 7.25. Fuzzy clustering of the input data in the calibration period (from 1988/01/01 to 1993/06/23) showing correlation of (a) effective rainfall and (b) discharge with model residuals. The labels C1 through C5 indicate the cluster ID.

Figure 7.26. Membership functions of fuzzy clustering to the data. The top most plot shows the observed hydrograph. Each data point belongs to all Clusters C1, C2, C3, C4, and C5 with different membership function values.

Figure 7.27 presents the separation of the hydrograph with respect to the different clusters. One may notice that the majority of the flows have the highest membership in Cluster C1, which can be interpreted as the base flow. The peaks and most of the high flows are attributed to Cluster C3. Some high flows, especially on the recession of the hydrograph, are attributed to Cluster C5, because these examples have less or no rainfall

Figure 7.27. Fuzzy clustering of the input data in the calibration period (from 1988/01/01 to 1993/06/23). The bottom plot shows the enlarged view of the monsoon event of 1990.

and cannot be grouped into C3, which has high values of both flow and rainfall. It can be said that the Fuzzy C-means clustering is able to identify the clusters corresponding to the various mechanisms of the runoff generation process, such as peak flow with high rainfall, high flow with no or less rainfall (recession of the hydrograph) and base flow, etc.

Selection of input variables for uncertainty model

In order to select the most important influencing variables for the uncertainty model U, an approach similar to the one used for clustering is followed. Correlation analysis and AMI analysis between input variables RE_t and Q_t (including lags) and the quantiles of the model error are conducted Figure 7.28a shows the correlation coefficient and AMI of RE_t and its lagged variables up to 7 days, i.e., RE_t-1, RE_t-2,…, RE_t-7 with the 5% and 95% quantiles. It is observed that the variables RE_t and RE_t-1 are strongly correlated with both quantiles; so these two variables are included in the input vector x_u.

 

The correlation and AMI analyses between Q_t and the quantiles of the model error are presented in Figure 7.28b. It is also observed that Q_t and Q_t-1 are strongly correlated with the quantiles. Although the lag 2 variable Q_t-2 also has a high correlation, only Q_t and Q_t-1 are included in the input vector. The reason is that the flow Q_t is highly auto-correlated and the inclusion of too many lagged variables of Q_t may lead to the redundancy of the model structure. Note that during the model application, Q_t is not available and we use its approximation made by model M. Although the simulated Q_t may bring additional uncertainty to the model U, our experiments have shown that this approach results in the more accurate model U (in terms of PICP and MPI).

Figure 7.28. Average mutual information (AMI) and correlation coefficient (Corr. Coef.) of 5% and 95% quantiles of the model errors with (a) effective rainfall and (b) discharge. The thin dark line shows the 5% quantile, the thick grey line shows the 95% quantile.

Selection and validation of uncertainty model

M5 model tree (MT) is used as an uncertainty prediction model U. There are certain advantages of using MT compared to other machine learning methods; it is simple, easy, and fast to train. The results are interpretable, understandable, and reproducible. Solomatine and Dulal (2003) have shown that MT can be used as an alternative to ANN in rainfall-runoff modelling. There is only one parameter in MT, the pruning factor (or, alternatively, the minimum number of data allowed in each linear model component), which controls the complexity of the model. The following shows the structure of the input data for the model U to predict 5% and 95% quantiles (see also Table 7.11):

where e⁵ and e⁹⁵ are the 5% and 95% quantiles of the model error respectively, and pf is the pruning factor. The following is the example of a generated model tree by U⁵for e⁵ with pf = 4:

 

 

There are five linear models (namely LM1, LM2, LM3, LM4, and LM5) generated for various intervals of Q_t, RE_t, and Q_t-1. Note that numbers inside the parenthesis are the numbers of the data vectors in the sub-sets. Similar structures of the linear models are obtained for e⁹⁵. Table 7.12 shows the mean and standard deviation statistics of the generated quantiles of the model error and performance of U⁵ and U⁹⁵ in the training, cross-validation and calibration data sets.

 

The mean and the standard deviation of the quantiles on the training and cross-validation data sets are consistent with the high variability in the original data. The performances of the uncertainty models U as measured by RMSE and CoE for 5% and 95% quantiles are quite reasonable for both the training and cross-validation data sets, in spite of the high variability of the quantiles. The RMSE and CoE values of the models U on the calibration data set are consistent with those for the training and cross-validation sets, and this ensures the predictability of the models U. Figure 7.29 depicts the scatter plot for the generated and predicted quantiles in the cross-validation data. It is observed that both models are quite good for approximating the relationship between the input space variables and the quantiles of the model error.

Table 7.12. Mean and standard deviation statistics of the 5% and 95% quantiles of the model error, and the performance of uncertainty prediction models U.

Note: Performances are RMSE and CoE. The training and cross-validation data constitutes 67% and 33%, respectively, of the calibration data (data used to calibrate process model). Values in parentheses correspond to statistics of 95% quantiles of the model error and the performance of uncertainty prediction models U⁹⁵. Std. dev., RMSE, and CoE are standard deviation of model errors, root mean squared error, the Nash-Sutcliffe between predicted and target values of the quantiles, respectively.

 

Figure 7.29. Scatter plot of the predicted and the target quantiles of the model errors for (a) 5% and (b) 95% quantiles in cross-validation data (part of the calibration data).

Analysis of model uncertainty

Figure 7.30. Cumulative distribution of the model residuals for each cluster in the calibration period. The last figure (second row and third column in the figure panel) shows the cumulative distribution of the model residuals for the calibration data.

Table 7.13. Cluster centres and computed quantiles of the model error for each cluster.

Note: Cluster C1 is identified by low values of flows and rainfall. Cluster C3 is identified by high flows and high rainfalls.

Figure 7.30 shows the empirical cumulative distribution of the model residuals for each cluster. In this plot sets S1, S2, S3, S4 and S5 belong to the sets of the input data which has the highest membership grade to the Cluster C1, C2, C3, C4 and C5, respectively. Obviously the distributions are quite different among them.

Figure 7.31 A comparison of 90% prediction bounds (darker shaded region) estimated with (a and b) the UNEEC method, (c and d) the GLUE method, (e and f) the meta-Gaussian method and (g and h) the quantile regression method in the validation period. The figures on the left shows the monsoon period of 1993 and the figures on the right show the monsoon of 1995. The dots show the observed discharge, the line shows the simulated discharge.

5% and 95% quantiles of the model residuals are computed for each of the clusters. The results are presented in Table 7.13. Investigating the clusters with their centres and quantiles reveals that the clusters with high values of rainfall and runoff have large values of the quantiles, while the clusters with low values of rainfall and runoff have small values of the quantiles. The asymmetrical values of these quantiles are obvious due to the fact that the HBV hydrological model calibration result is biased. It is observed that most of the time the hydrological model underestimates the river flows in the calibration period.

 

Figure 7.31 shows the observed discharge, the 90% hydrograph prediction uncertainty and comparison to the other three methods. The details of the comparison follow later. Figure 7.31a highlights the flood event that occurred during the monsoon of 1993. Interestingly enough the HBV model captures the highest peak flow very well and consequently the peak flow is bracketed by the predicted PIs. Figure 7.31b focuses on another monsoon event during 1995. This event was underestimated by the HBV model. One can see that the estimated uncertainty bound fails to enclose the highest peak discharge of the 1995 monsoon.

 

It is observed that 88.07% of the observed data points are enclosed within the computed PIs. 6.4% of the validation data points fall below the lower PI, whereas 5.53% data points fall above the upper PI. The average width of the uncertainty bounds i.e. MPI is 165 m³/s. This value is reasonable if compared with the order of magnitude of the model error in the validation data. The further analysis reveals that the distribution of the observed discharge below the lower PI is relatively consistent with the observed discharge. As far as the upper PI is concerned, less data are outside in the low flow (range of 0-250 m³/s). This means that the upper PIs are unnecessarily overestimated. However, the width of the upper PI in the intermediate flows (range of 250-750 m³/s) is considerably narrower.

1.1.6      Comparison of uncertainty results

In this section the results are compared with the widely-used Monte Carlo method GLUE, the meta-Gaussian, and the quantile regression (Koenker and Bassett, 1978) method. Note that the comparison with GLUE is performed purely for illustration since it analyses the parametric uncertainty, and other methods – the uncertainty based on the analysis of the optimal (calibrated) model residuals.

 

The experiment for the GLUE method (Beven and Binley, 1992) is set up similar to the Brue case study except, the number of behavioural parameter sets is increased to 25,000. The comparison results are reported in Figure 7.31c, d. One may notice the differences between the prediction bounds estimated by the UNEEC and GLUE method, but, again, these two techniques are different in nature, which was discussed in the previous case study (section 7.3).

 

It can be noticed that only 63.9% of the observed discharge values in the validation data fall inside the 90% prediction bounds estimated by the GLUE method. As expected, the width of the prediction bounds is smaller than those obtained with the other methods. The average value of the prediction bound width is 120.35 m³/s (see Figure 7.32). The further detailed analysis reveals that only 6.51% of the observed discharges are below lower PIs. The majority of the observed flows (29.61%) fall above the upper PIs.

 

The 90% prediction bounds estimated with the meta-Gaussian approach are shown in Figure 7.31e, f. Interestingly, it is found that 90% (more accurately, 90.02%) of the observed discharge values in validation data fall inside the estimated 90% prediction bounds (see also Figure 7.32). Further analysis of the results reveals that 3.7% of the

Figure 7.32. A Comparison of the statistics of uncertainty estimated with the UNEEC, meta-Gaussian, GLUE, and quantile regression methods in the validation period.

observed data are below the lower PI whereas 6.3% of data are above the upper PI. However, one can see that the bounds width is consistently larger and the average width of the prediction bounds is 225.79 m³/s, which is about 35% larger than that estimated by the UNEEC method.

 

Quantile regression method was used to compute the 5% and 95% error quantiles. The input variables x (see equation (2.26)) used in this method are the same as those used in the UNEEC method. We have also made several experiments with other input variables combination, but the results are not so good. The best results are reported in Figure 7.31g, h. It is observed that 86.12% of the observed discharge values in validation data fall inside the estimated 90% prediction bounds (see Figure 7.32). Further analysis of the results reveals that 4.01% of the observed data are below the lower PI, whereas 9.87% of data are above the upper PI. The average width of the prediction bounds is 216.86 m³/s.

Estimation of probability distribution of model errors

In the previous sections, only 90% prediction intervals, i.e., 5% and 95% quantiles are estimated. In this section the extension of the method to derive a more accurate estimate of the pdf is demonstrated. Several quantiles (such as 2.5, 5, 10:10:90, 95, 97.5%) are computed for the calibration data after clustering of the input space. Then regression models are trained for each quantile independently:

where p = 2.5, 5, 10:10:90, 95, 97.5 %. In this experiment, a total of 13 regression models are trained. Once the models U ^p are trained on the training data, they are used to predict the quantiles for the new input data (e.g., for the validation data set). Note that since MTs are used as regression models, it takes only a couple of seconds to train a single model, so the computational cost to estimate the full distribution is not a major concern. However, this could be an issue for computationally extensive algorithms, such as support vector machine or ANNs with long records of input data.

 

The cumulative probability distribution (cdf) for the peak discharge of the flood dated 21 July 1993 (the highest peak event of the 1993 monsoon shown in Figure 7.31) is shown in Figure 7.33a. The cdf computed from the GLUE, meta-Gaussian, and quantile regression methods are also presented for comparison. One can notice that the cdf computed from the UNEEC method is relatively steep. The GLUE and quantile regression methods produce a comparatively flat cdf. For this particular flood event this means that the uncertainties estimated by the GLUE and quantile regression methods are very high and the uncertainty estimated by the UNEEC method is lower. The meta-Gaussian method gives intermediate results. Additional analysis of the cdf for the flood event of 14 August 1995 supports the finding that the uncertainty estimated with the UNEEC method is consistently lower for the flood events (Figure 7.33b).

 

Figure 7.33. A Comparison of estimation of cumulative probability distribution for peak discharges of the monsoon period of (a) 1993 and (b) 1995.

Figure 7.34. A Comparison of the statistics of uncertainty measures. (a) Prediction interval coverage probability (PICP). In an ideal case, the plot between PICP and confidence level follows the thick grey line. (b) Mean prediction interval (MPI) for different values of the confidence level.

Further analysis is performed in order to compute the percentage of the observed discharge data falling within the estimated uncertainty bounds (i.e. PICP) and the average width of the uncertainty bounds (i.e. MPI) for various specified confidence levels (ranging from 20% to 95%) used to produce these uncertainty bounds. The results are presented in Figure 7.34. As far as the PICP is concerned, the ideal would be to follow the thick grey line (Figure 7.34a). Points below this line indicate that less data are bracketed by the uncertainty bounds. On the other hand, if more data are enclosed in the uncertainty bounds, the PICP line would be above the ideal line. It can be seen that the PICPs computed with the meta-Gaussian and QR methods are very close to the desired confidence levels. In the UNEEC method more data are enclosed at lower values of the confidence levels. The GLUE method produces consistently lower values of PICPs for all confidence levels. The results with the GLUE method are consistent with the results reported in the literature, namely that the percentages of the observed discharge falling within the uncertainty bounds estimated by the GLUE is normally much lower than the specified confidence levels (see, e. g., Montanari, 2005; Xiong and O’Connor, 2008). The low values of PICP of the GLUE uncertainty bounds can be attributed to many factors which are discussed in section 5.6.2.

 

As far as the MPI is concerned, the GLUE method generates relatively narrower uncertainty bounds (Figure 7.34b). The MPI values estimated by UNEEC and meta-Gaussian methods are comparable at the lower values of the confidence levels. However, uncertainty bounds estimated by the meta-Gaussian method increase faster as compared to those obtained with the UNEEC method after 60% confidence levels. The MPI estimated with the QR method is similar to that obtained with the meta-Gaussian method.

Conclusions

This case study presents the application of the UNEEC method for the prediction of uncertainty of the simulated river flows referring to the case study of the Bagmati catchment in Nepal. In this case study, we also observe that the model residuals are neither normally distributed nor homoscedastic. Compared to the case study of the Brue catchment, the uncertainty of the model prediction is higher; this is obvious due to the following reasons: (i) the size of catchment of this case study is larger and thus basin average rainfall, temperature data may not be representative; (ii) the quality of data is poorer; and (iii) the resolution of data is higher.

 

The results show that the estimated uncertainty bounds are consistent with the order of magnitude of the model errors in validation period. We also compare the uncertainty bounds with those estimated by GLUE and meta-Gaussian and quantile regression methods. The comparison results shows that the UNEEC method gives reasonable estimate of the uncertainty in the model prediction; and the percentage of observed data falling inside the uncertainty bounds estimated by the UNEEC method is very close to the desired degree of confidence level used to derive these bounds. It has been also demonstrated the extension of the UNEEC method to approximate the probability distribution function of the model output.

Multiobjective calibration and uncertainty

In sections 7.1-7.4, uncertainty analysis of the optimal model is carried out. The model is calibrated using the single objective function value. This section presents some results with multiobjective calibration and uncertainty.

 

Practical experience with the calibration of the hydrological model suggests that a single objective function value is often inadequate to measure properly the simulation of all the important characteristics of the system that are reflected in the observations. Gupta et al. (1998) pointed out that there may not exist an objective “statistically correct” choice for the objective function, and therefore no statistically correct “optimal choice for the model parameters. Furthermore, recent advances in computational power have led to more complex hydrological models, often predicting multiple hydrological fluxes simultaneously. These issues have led to an increasing interest in the multi -objective calibration of hydrological model parameters (Gupta et al., 1998; Yapo et al., 1998; Madsen, 2000; Khu and Madsen, 2005).

 

Multiobjective calibration problem can be stated as follows:

where m is the number of objective functions, f_i(q), i = 1, ..., m are the individual objective functions, and q  is the set of model parameters to be calibrated. Due to trade-offs between the different objectives, the solution to equation (7.15) will no longer, in general, be a single unique parameter set. Instead, there exist several solutions that constitute a so-called Pareto optimal set or non-dominated solutions. Any solution q^* belongs to the Pareto set when there is no feasible solution q that will improve some objective values without degrading performance in at least one other objective. Mathematically, the solution q^*  is Pareto-optimal (i) if and only if f_i(q^*) £ f_i(q)  for all k = 1, ..., m and (ii) f_j(q^*) < f_j(q) for some j = 1, ..., m. According to these two statements the Pareto-optimal solution q^*  has at least one smaller objective value compared to any other feasible solution q  in the decision space, while performing as well or worse than q  in all remaining objectives.

 

The multiobjective optimisation algorithm used in this study is the non-dominated sorting genetic algorithm (NSGA-II) developed by Deb et al. (2002). NSGA-II is capable of handling a large number of objective functions and provides an approximate representation of the Pareto set with a single optimisation run. The NSGA-II algorithm is outlined briefly as follows:

1.       Generate an initial population of size p randomly in the domain of the feasible range;

2.       Evaluate the population and sort the population based on non-domination using a bookkeeping procedure to reduce the order of the computation;

3.       Classify the population into several Pareto fronts based on the non-domination level. Individuals belonging to the first Pareto front are assigned with rank 1, individuals belonging to the second Pareto front (second Pareto front is the Pareto front after removing the individuals from the first front) are assigned with rank 2, and so on;

4.       Form a new population by generating an offspring population from the parents and combining them with the parents (following standard GA procedure);

5.       Compute the crowding distance for each individual;

6.       Select the individuals based on a non-domination rank. The crowding distance comparison is made if the individual belongs to the same rank;

7.       Repeat the steps 3 – 6 until the stopping criterion is satisfied. The stopping criterion may be a specified number of generations, maximum number of function evaluations or computation time.

Preference ordering of the Pareto-optimal points

It is observed that the number of the Pareto-optimal solution grows with the dimensionality of the objective functions used in the multiobjective optimisation. It was reported that the qualitatively of the solutions are indistinguishable from each other. However, it is often required to select a Pareto-optimal point from a set of many and it becomes highly intractable when the number of solutions becomes large. For example, the number of Pareto-optimal solutions in calibrating hydrological model parameters in four dimension objective functions reaches easily above 100. Since the decision maker generally prefers a small numbers of solutions, modelers have to face the task of selecting a set of suitable model parameters from the numerous Pareto-optimal sets.

 

The concept of the preference order of efficiency, which was first introduced by Das (1999), can be applied to select the most preferable solutions by setting up a hierarchical ordering among the various solutions obtained from the multiobjective optimizations. The preference ordering (PO) uses the following two theorems (Das, 1999), which are strong domination criteria compared to Pareto domination:

 

The first theorem, efficiency of order k (k-Pareto-optimal points), considers all the possible k-dimensional subspace of the original m-dimensional objective function space (1£ k £ m). A solution is defined as being efficient of order k, if this solution is not dominated by any other solution in any of the k-dimensional subspaces. The second theorem is the efficiency of order k with degree p (denoted as [k, p]). A solution is defined as being efficient of order k with degree p, if it is not dominated by any other solutions for exactly p out of the possible _mC_k  k-dimensional subspaces.

 

The condition of efficiency of order can thus be used to help reduce the number of solutions from the Pareto-optimal set by retaining only those that are regarded as the best compromises. When the number of points selected is still considerably large, a more stringent condition of efficiency with degrees is required to sort out better solution. Thus by reducing the efficiency of order and increasing the degree of order in a sequential manner, the criteria of Pareto-optimality is tightened, which enables the modeller to sieve through the large number of Pareto-optimal sets and determine the preferred solution (Khu and Madsen, 2005).

Results and discussions

The NSGA-II algorithm was applied to Bagmati catchment for multiobjective calibration of HBV model parameters. Four different model performance measures that emphasise different aspects of the hydrograph are considered as objective functions for calibration. These objective functions are as follows:

where Q_i is the observed discharge at time i, is the simulated discharge, N is the number of time steps in the calibration period, N_l is the number of time steps in low flow events, N_h is the number of time steps in high flow events. In this study, low flow events are defined as periods with flow below a threshold value of 50 m³/s and high flow events are defined as periods with flow above 500 m³/s.

 

Since GA is a stochastic process, 10 independent calibration runs are carried out to compare the quality of the solutions of each run. It is observed that each independent run is comparable to the other. The maximum number of the generation and the population size is set at 50 and 80, respectively. The probability of crossover and mutation are set to 0.9 and 0.1, respectively. We also use a single objective GA algorithm to calibrate each of the single objective functions, (equations (7.16) - (7.19)), i.e., optimising for volume error f₁(q), RMSE f₂(q), low flow RMSE f₃(q), high flow RMSE f₄(q).The following results are from one of the ten runs.

 

Table 7.14 presents the number of points on the Pareto front for all possible combinations of the four objective functions (equations (7.16)-(7.19)). Since there are four objective functions (m=4), the possible combinations are given by _mC_k , where K =

Table 7.14. All possible combinations of the four performance measures, denoted as 1-2-3-4 and the resulting number of Pareto-optimal solutions.

Runs	1-2-3-4	1-2-3	1-2-4	1-3-4	2-3-4	1-2	1-3	1-4	2-3	2-4	3-4
1	413	349	33	362	88	17	67	18	33	11	39
2	347	260	33	289	73	13	46	8	30	12	41
3	392	315	41	318	70	17	50	12	31	12	27
4	300	240	28	255	62	13	48	8	27	7	35
5	425	327	46	328	81	17	44	12	37	8	41
6	342	296	32	275	56	13	39	9	29	4	24
7	371	303	47	321	64	12	52	11	34	14	36
8	387	340	27	331	47	11	46	10	32	4	34
9	294	231	41	229	79	14	55	16	36	12	41
10	382	335	38	339	43	14	54	14	28	10	37
Avg.	365	300	37	305	66	14	50	12	32	9	36

 

{4, 3, 2}. In total there are 11 combinations – (i) one combination of all four objective functions (denoted as (1-2-3-4); (ii)  four combinations of any three objective functions (i.e., (1-2-3), (1-2-4), (1-3-4), (2-3-4)); (iii) six combinations of any two objective functions (i.e., (1-2), (1-3), (1-4), (2-3), (2-4), (3-4)). It can be seen from Table 7.14 that the numbers of solutions for each of the 10 runs are not significantly different from each other. Although this does not guarantee that the good solutions are present, it does indicate that the quality of each run is comparable to the others.

 

It is observed that the resulting number of Pareto-optimal points in 4-dimensional space is, on average, about 365. With such a large number of Pareto-optimal solutions, it is not feasible to examine each and every one of the solutions to identify and select the most appropriate calibration parameter set. One must therefore be able to sieve through these solutions to obtain a smaller number of representative solutions. If we examine the four possible combinations of any three objective functions (columns 3– 6), the resulting number of Pareto-optimal points ranges between 37 and 305. If we examine the six possible combinations of any two objective functions, the resulting number of Pareto-optimal points ranges between 9 and 50. However, it should be noted that these Pareto-optimal points for the two objective functions may not be Pareto-optimal in the other objective function spaces.

 

Preference ordering has been applied to select a small number of solutions from the higher dimension Pareto-optimal solutions. The results are presented in Table 7.15. As the efficiency of order (k) decreases, the number of preferred solutions decreases; and as the degree (p) of that order increases, the number of preferred solutions decreases. From Table 7.14 we see that the number of solutions reduces from an average of 361 to 2 when the degree of order 3 increases from 1 to 4. Thus the problem becomes highly manageable when the number of solutions becomes small while the qualities of these solutions in terms of Pareto-optimality are good (they are of a high degree and high order of efficiency, i.e., [3,4]).

 

Taking the example of run 1 in Table 7.15, it can be seen that if all four objective functions are considered, there are a total of 413 solutions. However, a closer examination of these solutions indicates that not all of them are Pareto-optimal if the

Table 7.15. Results of applying preference ordering to determine the number of preferred solutions in the 4th, 3rd, and 2nd order of efficiency (k) for different degrees (p) denoted by [k, p].

Runs	[4,1]	[3,1]	[3,2]	[3,3]	[3,4]	[2,1]	[2,2]	[2,3]	[2,4]	[2,5]	[2,6]
1	413	410	360	61	1	146	33	6	0	0	0
2	347	342	275	35	3	123	23	4	0	0	0
3	392	389	304	50	1	113	28	8	0	0	0
4	300	298	246	41	0	105	26	7	0	0	0
5	425	418	309	55	0	118	34	7	0	0	0
6	342	334	281	42	2	87	24	7	0	0	0
7	371	367	323	43	2	121	29	9	0	0	0
8	387	382	318	43	2	103	30	4	0	0	0
9	294	289	244	45	2	142	28	4	0	0	0
10	382	377	327	48	3	116	34	6	1	0	0
Avg.	365	361	299	46	2	117	29	6	0	0	0

 

objective function combinations are reduced and changed. When the combinations of any three out of the four objective functions are considered, (i) 410 of the original 413 Pareto-optimal solutions are present in at least one of the four possible combinations; (ii) 360 of the original 413 solutions are present in at least two of the four combinations; (iii) 61 of the 413 solutions are present in at least three of the four combinations; and (iv) one of the 413 solutions are present in all four combinations. Thus the preference ordering techniques has selected one preferred solution out of the original 413 Pareto-optimal solutions. If the number of solutions has not been reduced sufficiently, the order of efficiency should be reduced to 2 instead of 3 and the degree of efficiency be considered in turn. The remaining columns in Table 7.15 show such a process ([2,1], [2,2], [2,3], [2,4], [2,5], [2,6]). For run 1 no solutions are present in degrees higher than 3 (i.e., four or more of the six two-objective combinations). Thus, in this case the one [3,4] Pareto-optimal solutions are selected for further analysis.

 

Table 7.16 shows the comparison of the parameter values of the preferred [3,4] Pareto-optimal set and the optimal values of the parameter for each of the single objective functions. As expected, the values of each parameter vary depending on the selected objective functions used for calibration. Note that the values of the parameters of the preferred Pareto-optimal set are within the ranges of the parameter values of the single objective calibrations except first two parameters. The variation of the optimal model parameters sets of four single objective calibrations along the Pareto front is shown in Figure 7.35.

 

Table 7.17 presents the comparison of the objective function values of the preferred Pareto-optimal set and the objective function values for each of the single objective optimisation. The widely used Nash-Sutcliffe coefficient (CoE) is also calculated as an additional performance measure. As expected, the best calibration and validation results for each objective functions are obtained using single objective functions optimisation. For example, the lowest RMSE value of 7.00 m³/s on low flow events is obtained when the parameters are optimised for single objective function f₃(q). Similarly, the lowest RMSE value of 265.78 m³/s on high flow events is obtained via calibrating objective function f₄(q). The results of the Pareto-optimal set are a good compromise between all

Table 7.16. Comparison of the preferred model parameters with the optimal parameters corresponding to the four single objective function values.

 

Figure 7.35. Normalised range of parameter values along the Pareto front. Thick solid magenta, black, green, and red lines indicate parameter sets calibrated (single objective) for volume error, RMSE, low flow RMSE, and high flow RMSE, respectively.

four objective functions and are within the best and worst solutions from the single objective function optimisations.

 

Figure 7.36 shows the Pareto sets (best 80 solutions) of various two-objective function combinations. Note that all these solutions are from Pareto sets of optimisation results using all four objective functions. For comparison, the solutions obtained from single objective optimisations and the preferred Pareto-optimal set are also shown in the figure. Further analysis of correlation between the objective functions shows that the objective function f₂(q) (i.e. RMSE) is highly correlated (linear correlation coefficient of 0.99) with the objective function f₄(q) (i.e. high flow RMSE). There is little positive correlation (correlation coefficient of about 0.38) between objective function f₁(q) (volume error) and f₂(q). As expected, the objective functions f₁(q) and f₄(q) are also

Table 7.17. Comparison of the ranges of model performances corresponding to the Pareto font with the four single objective function values in calibration and verification period

Note: Bold type indicates the best value.

Figure 7.36. Various two-objective function combination plots of Pareto-optimal sets. Circle indicates Pareto solutions in 4-D objective functions space, plus marks in circles indicate Pareto solutions in 2-D space. Asterisk and cross indicates solutions corresponding to the single objective calibration and square indicates the preferred Pareto-optimal set.

correlated with an almost equal value of the correlation coefficient (0.39). However, as seen in the figure, other two-objective functions combinations such as f₁(q) and f₃(q) or f₂(q) and f₃(q) or f₃(q) and f₄(q) have negative correlation (correlation coefficient ranges from 0.53-0.60). The negative correlation means the two objective functions are conflicting.

 

The simulated hydrographs obtained by multi and single objective optimisations in both calibration and validation periods are shown in Figure 7.37. The dark shaded region is the range of the simulation obtained by the best 80 Pareto-optimal sets. It is observed that simulations obtained by single objective optimisations are within the range of hydrographs obtained by the Pareto-optimal sets. Furthermore, as expected, the simulation obtained by the parameters optimised for low flow RMSE reproduces low flow accurately, while the simulation obtained by the parameters optimised for high flow RMSE is close to the high flow events. The preferred Pareto-optimal set gives the compromise simulation.

Uncertainty assessment

Equifinality principle of Beven and Binley (1992) is restated here that no unique optimum parameter set exists in calibration of hydrological model, but rather there will always be several different models that mimic equally well an observed natural process. In a multiobjective context, there is a multitude of parameter combinations that are equally good in reproducing the observed responses; thus multiobjective calibration accords with the equifinality principle. These equally good parameters sets are the Pareto-optimal set.

 

An illustration of the interpretation of Pareto-optimal sets of parameter within equifinality principle is shown in Figure 7.35. It is observed that parameter sets cover a large range of parameter values, but produce virtually equal good simulations according to a specified objective function criterion. Consequently these parameter sets produce a range of simulated hydrographs as shown in Figure 7.37. These figures illustrate that there are considerable ranges of the model simulations which are equally good according to a specified objective function criterion. In a way, this range can be interpreted as uncertainty of the model prediction corresponding to the Pareto-optimal parameter sets, and it provides additional information necessary for risk based decision making. With respect to the UNEEC method, as mentioned in Chapter 6, we could apply the UNEEC method for each solution of the Pareto-optimal sets, thus producing the several uncertainty bounds corresponding to Pareto-optimal sets. More research is anticipated in the future in this direction.

Figure 7.37. Ranges of simulated hydrographs corresponding to the Pareto parameter sets. The ranges are also compared with the solutions obtained with four single objective functions and preferred solution. (a and b) Part of the calibration period and (c and d) part of the verification period.

Conclusions

This chapter presents the application of UNEEC method to estimate uncertainty in the model prediction for the three different catchments, which have highly varied hydro-climatic characteristics. The uncertainty bound of the model prediction estimated with the UNEEC method is compared with the observed data. The results show that the percentage of the observed discharge falling within the estimated uncertainty bound is very close to the specified confidence level used to produce these bounds. It has been demonstrated that uncertainty bounds are consistent with the order of the magnitude of the model errors.

 

A comparative assessment of the uncertainty bounds with those obtained by other existing uncertainty methods including GLUE, the meta-Gaussian, the quantile regression indicates that the uncertainty estimation with the UNEEC method is reasonably accurate. The UNEEC method is computationally efficient and can be applied to any kind of models. This application concludes that the machine learning techniques can be used to analyse, model and predict the residual uncertainty of the optimal model output. Note that residual uncertainty considers all sources of uncertainty in aggregated form through model errors.