## Abstract

Hydrological models are simplified imitations of natural and man-made water systems, and because of this simplification, always deal with inherent uncertainty. To develop more rigorous modeling procedures and to provide more reliable results, it is inevitable to consider and estimate this uncertainty. Although there are different approaches in the literature to assess the parametric uncertainty of hydrological models, their structures and results have rarely been compared systematically. In this research, two different approaches to analyze parametric uncertainty, namely *direct* and *inverse* methods are compared and contrasted. While the direct method employs a sampling simulation procedure to generate posterior distributions of parameters, the inverse method utilizes an optimization-based approach to optimize parameter sets of an interval-based hydrological model. Two different hydrological models and case studies are employed, and the models are set by two distinct mathematical operations of interval mathematics. Findings of this research show that while the choice of the interval mathematic method can affect the final results, generally, the inverse method cannot be counted on as a reliable tool to analyze the parametric uncertainty of hydrological models, and the direct method provides more accurate results.

## INTRODUCTION

Hydrological modeling, like any other physical and empirical process-based representation of a natural phenomenon, is a simplification of real-world systems; therefore, it always deals with inevitable uncertainty (Solomatine *et al.* 2009). Since an important task of hydrological models is to predict future characteristics of water cycle systems, it is crucial to consider and evaluate their uncertainty. Different sources of uncertainty are as below (Beven & Binley 1992; Lall & Sharma 1996; Wagener & Gupta 2005; Kavetski *et al.* 2006; McMillan *et al.* 2011; Hattermann *et al.* 2018):

Uncertainty of measurements (input and output data),

Uncertainty of model structure, and

Uncertainty of model parameters.

However, generally speaking, there are three main sources for modeling uncertainty: errors related to data to calibrate model (input and output data), flaws and imprecision in model structure, and uncertainty in model parameters (Refsgaard & Storm 1990). The successful application of a hydrological model chiefly depends on how well its parameters are selected and calibrated (Jiang *et al.* 2010). Therefore, among all different sources of uncertainty, the uncertainty of model parameters is significant and has received more attention; hence, there are a variety of methods and procedures aiming at its evaluation and interpretation (Kuczera & Parent 1998; Vrugt *et al.* 2003; Moradkhani *et al.* 2005; Benke *et al.* 2008; Nasseri *et al.* 2014a; Ahmadi *et al.* 2019). Nasseri *et al.* (2013) classified this wide range of uncertainty assessment methods into probabilistic, possibilistic, and hybrid classes. However, from another standpoint, these methods can be classified into two well-known and widely used categories, namely direct and inverse methods.

Direct approaches usually employ probabilistic sampling simulation procedures, receiving the model efficiency of likelihood feedbacks (e.g. Monte Carlo and Markov Chain Monte Carlo simulations) to generate accepted samples and to determine their statistical characteristics (Kuczera & Parent 1998; Bates & Campbell 2001; Marshall *et al.* 2004; Vrugt *et al.* 2008; Abbaszadeh *et al.* 2018). By analyzing posterior distributions of parameters and model output (e.g. runoff), their uncertainty can be evaluated and eventually quantified (Jin *et al.* 2010). In this class, the focus is on parameter sampling methods and their efficiencies employing a feedforward flow of parameter generation and its feedback in the form of likelihood values. Among different direct methods to estimate the parametric uncertainty of hydrological models, the generalized likelihood uncertainty estimation (GLUE), developed by Beven & Binley (1992), is extensively used in many different studies (Freer *et al.* 1996; Beven & Freer 2001; Blasone *et al.* 2008; Yu *et al.* 2015; Cu Thi *et al.* 2018).

Contrary to direct approaches, inverse methods do not provide posterior distributions of model parameters. Basically, they optimize the hydrological model's parameters to yield output as similar to the observed data as possible. In other words, inverse approaches try to minimize the model error, which is the difference between the observed and computed outputs. In the literature, there are various optimization algorithms used to calibrate hydrological models, e.g. particle swarm algorithm (Jiang *et al.* 2010), multi-objective particle swarm (Gill *et al.* 2006), and genetic algorithm (GA) (Seibert 2000; Nasseri *et al.* 2008). Generally, optimization algorithms are used to calibrate hydrological models based on the available observed data. This means that the outcome of the calibration is a single set of parameters which provides the best similarity between the observed and computed data. However, in the case of uncertainty analysis, we are interested in the parameters' intervals (i.e. lower and upper bounds of parameters) to generate the uncertain bound of model response (i.e. output) instead of a deterministic answer. To deal with this issue, usually, interval uncertainty assessment based on fuzzy mathematics or standard interval mathematics (SIM) is employed (Alvisi *et al.* 2006; Alvisi & Franchini 2011; Nasseri *et al.* 2013, 2014a; Ahmadi *et al.* 2019). However, for this purpose, there is another alternative in the literature named grey mathematics with similar mathematical formulation (Alvisi & Franchini 2012; Alvisi *et al.* 2012).

According to the equifinality thesis, proposed by Beven (2006), when assessing the uncertainty associated with prediction, rather than one correct answer, there are many acceptable representations of the system's reality that cannot be easily rejected and should be taken into account. In other words, there are various ways of modeling to employ, not a single correct representation. Having this point in mind, various studies are trying to pool together the output of multiple models to generate an ensemble of the results (e.g. hydrographs) to represent uncertainty (Carpenter & Georgakakos 2004; Georgakakos *et al.* 2004; Marshall *et al.* 2007). Like providing an ensemble of the results, attempts to compare different methods to assess the uncertainty of hydrological processes can cast light on the modeling procedure and, therefore, have merit. Despite the considerable variation in the literature of the parametric uncertainty assessment, different methods have rarely been compared and contrasted; therefore, there is a critical research gap.

Trying to fill this gap, in the current article, two different interval mathematic formulations which come from fuzzy mathematics, namely SIM and modified interval mathematics (MIM) are applied to develop interval-based hydrological models. Then, by the means of GA, the best bounds of the models' parameters are calibrated as the output of the inverse models. In addition to the inverse method, GLUE is used to analyze the parametric uncertainty of the same models and to evaluate the inverse method's results. So, two specific objectives have been followed in the current article: to analyze the parametric uncertainty of two conceptual hydrologic models, each assigned to an appropriate case study, using direct and inverse methods; and to compare and contrast the results of these two methods and to examine the dissimilarities in their results and discuss the reasons for them.

In the following parts of the paper, first, the implemented direct and inverse methods are presented, followed by an introduction to the case studies and hydrological models. Then, the comparative results are presented and subsequently discussed and concluded.

## HYDROLOGICAL MODELS AND CASE STUDIES

Figure 1 shows the flowchart depicting the methodology presented in the current research. The steps of the methodology are discussed in the following sections. The first step is the selection of the hydrological models and case studies that are presented in the following subsections (Step 1: Preprocessing).

### Hydrological models

Two monthly water balance models are used in this research. The first model, which is referred to as the Guo model, was originally developed by Guo *et al.* (2002); it is a water balance model appropriate to model wet climate areas. The model consists of a surface water storage, without groundwater and subsurface flow. Precipitation and evapotranspiration are its input variables, and its output is total streamflow. The model has three parameters to calibrate, namely evapotranspiration coefficient (*C*), initial soil water content (*S*1), and soil moisture scale factor (SC). The model structure is depicted in Figure 2(a).

The second hydrological model, which is referred to as the enhanced Guo model with snowpack, is another monthly water balance model, first proposed by Guo *et al.* (2005) (without considering snowpack) and developed further by Nasseri *et al.* (2014b, 2019). Like the first model, this hydrological model employs water continuity and mass balance equations, but it also includes portions of snowmelt, soil water content, and groundwater on simulated streamflow. The model has ten parameters to calibrate, namely precipitation scale factor (SF), initial snowpack (SP(1)), minimum and maximum temperature criteria (*T*_{s} and *T*_{m}, respectively), surface runoff coefficient (*K*_{s}), snowmelt coefficients (*K*_{sn} and *K*_{sn1}), base flow coefficient (*K*_{g}), soil moisture scale factor (*S*_{max}), and evapotranspiration parameter (*P*_{1}). The model structure is illustrated in Figure 2(b).

### Case studies

The first case study of the paper is the Adour–Garrone basin (at Pont de Maussac) in southwest France, which is assigned to the first hydrological model because, in this basin, the snowmelt contribution to the total runoff is very low. The data of this case study are from January 1954 to December 1994 and have been used in some earlier studies (Perrin *et al.* 2001; Nasseri *et al.* 2013; Ahmadi *et al.* 2019). The location map of the Adour–Garrone basin is illustrated in Figure 3.

The second case study is the Roodzard basin, which is a middle-sized sub-basin of the Zohre–Jarrahi watershed, in the Khuzestan province, southwest of Iran. The climatic and hydrological characteristics of this basin differ from the first case study very widely. Since this basin is located in a mountainous area and snowmelt has a considerable contribution to its total runoff, it is assigned to the second hydrological model. The data relating to this basin are from October 1974 to April 2003. The data show substantial seasonality, with significant fluctuations in long-term average monthly discharges, from 2.31 m^{3}/s in August to 13.83 m^{3}/s in February (Nasseri *et al.* 2013, 2014b, 2019; Ahmadi *et al.* 2019). The location map of the Roodzard basin is depicted in Figure 4. Table 1 reports some hydrological information on the case studies.

Basin . | Area (km^{2})
. | Rainfall (mm/km^{2}). | Evaporation (mm/km^{2}). | Runoff (mm/km^{2}). | |||
---|---|---|---|---|---|---|---|

Minimum . | Maximum . | Minimum . | Maximum . | Minimum . | Maximum . | ||

Adour–Garrone | 82.3 | 2.8 | 439.5 | 12.9 | 131.2 | 0.3 | 257.0 |

Roodzard | 900 | 0 | 310.0 | 46.9 | 538.6 | 0 | 215.0 |

Basin . | Area (km^{2})
. | Rainfall (mm/km^{2}). | Evaporation (mm/km^{2}). | Runoff (mm/km^{2}). | |||
---|---|---|---|---|---|---|---|

Minimum . | Maximum . | Minimum . | Maximum . | Minimum . | Maximum . | ||

Adour–Garrone | 82.3 | 2.8 | 439.5 | 12.9 | 131.2 | 0.3 | 257.0 |

Roodzard | 900 | 0 | 310.0 | 46.9 | 538.6 | 0 | 215.0 |

## MODELING PROCEDURE

As depicted in Figure 1, the first step for implementing both direct and inverse methods is setting the range of model parameters (Step 2: Uncertainty Assessment). It is important to choose the same range for each parameter in different methods, which is indispensable for an impartial comparison of the approaches' results. The range of the parameters of each hydrological model is chosen based on the physical and mathematical characteristics and also the parameters' nature. The models' parameters and their ranges are presented in Table 2. In the following sections, direct and inverse uncertainty assessment methods are described.

Simple Guo model . | Enhanced Guo model . | ||
---|---|---|---|

Parameter . | Range . | Parameter . | Range . |

Evapotranspiration coefficient (C) | (0.5,2) | Precipitation scale factor (SF) | (0.5,2) |

Soil moisture scale factor (SC) | (200,6000) | Initial snowpack (SP(1)) | (0,100) |

Initial soil water content (S1) | (100,250) | Lower temperature criterion (T_{s}) | (−5,10) |

Upper temperature criterion (T_{m}) | (0,40) | ||

Runoff coefficient (K_{s}) | (0,2.5) | ||

Snowmelt coefficient (K_{sn}) | (0.5,9) | ||

Base flow lag coefficient (K_{g}) | (0,0.4) | ||

Snowmelt coefficient (K_{sn1}) | (0,10) | ||

Soil moisture scale factor (S_{max}) | (0,1500) | ||

Evapotranspiration parameter (P_{1}) | (0,2.5) |

Simple Guo model . | Enhanced Guo model . | ||
---|---|---|---|

Parameter . | Range . | Parameter . | Range . |

Evapotranspiration coefficient (C) | (0.5,2) | Precipitation scale factor (SF) | (0.5,2) |

Soil moisture scale factor (SC) | (200,6000) | Initial snowpack (SP(1)) | (0,100) |

Initial soil water content (S1) | (100,250) | Lower temperature criterion (T_{s}) | (−5,10) |

Upper temperature criterion (T_{m}) | (0,40) | ||

Runoff coefficient (K_{s}) | (0,2.5) | ||

Snowmelt coefficient (K_{sn}) | (0.5,9) | ||

Base flow lag coefficient (K_{g}) | (0,0.4) | ||

Snowmelt coefficient (K_{sn1}) | (0,10) | ||

Soil moisture scale factor (S_{max}) | (0,1500) | ||

Evapotranspiration parameter (P_{1}) | (0,2.5) |

### Direct uncertainty assessment method

*N*is the number of simulation's time steps. In the direct method, three different acceptance rates () are considered; and for each, the GLUE algorithm is run to generate 50,000 accepted samples. Basically, an accepted sample consists of a set of parameters that are used in the hydrological model under simulation and generates the time series of streamflow that has an NS value more than the acceptance rate. All 50,000 accepted parameter samples along with their computed streamflow time series will be processed to achieve their posterior distributions.

### Inverse uncertainty assessment method

The results of the inverse method to analyze parametric uncertainty are explained in an interval form for each parameter. In other words, it provides the optimized upper and lower bounds of the parameters instead of their distributions. In this research, a real-coded GA is used to optimize the hydrological models compatible with the interval mathematics. For all optimizations, the ranges of parameters are the same and equal to the direct method. Moreover, GA parameters are the same for all optimizations, namely mutation probability, crossover probability, population size, and maximum number of generations which are 0.09, 0.8, 100, and 2,500, respectively. By setting a penalty for the fitness function, the GA optimization is conditioned to generate results with specific characteristics. It should be noted that contrary to the direct method that generates a single deterministic output, the inverse method generates an uncertain interval. Therefore, deterministic criteria like NS are not suitable to evaluate the inverse method's performance.

*et al.*2009):where and are the upper and lower limits of streamflow, which are computed by the optimized hydrological model in the time step by applying the optimized set of parameters to the hydrological models, and is the observed streamflow in that time step. Again,

*N*is the number of observations.

*et al.*2013, 2014a). In the formulation of NUE, POC along with average relative interval length (ARIL) is used. The formulations of ARIL and NUE are as follows:where

*N*, , , and are the same as the prior equation, and

*w*is the scale factor of POC versus ARIL and it is considered to be equal to 1 in this study. As can be understood from Equations (2)–(4), setting the fitness function of GA equal to NUE implies that the optimized set of parameters is the set resulting in the narrowest possible streamflow band that contains as many observed instances as possible. Having the constraint of the optimization, this set of parameters should also have a POC amount more than or equal to a pre-specified value.

*et al.*2010; Nasseri

*et al.*2014a). For example, addition operation () in SIM (as in SFM) is as follows:where and represent the lower values and and represent the upper values of the interval numbers and (as the fuzzy interval number), respectively.

*et al.*(2013, 2014a). For example, the formulation of addition operation () in MIM is as follows:

## RESULTS OF THE PARAMETRIC UNCERTAINTY ASSESSMENTS

As mentioned earlier, in both direct and inverse methods, all parameters are selected from a pre-specified range, which is the same for both methods. The parameters and their allocated ranges are presented in Table 2.

In Figures 5 and 6, the histograms of 50,000 accepted samples of the parameters generated by GLUE as the direct method for three different acceptance criteria (NS = 60%, 65%, and 70%) are presented. In each figure, the lower and upper limits of the parameter sets achieved by the inverse method for three different optimization constraints are highlighted as well. The number above each panel is the NS values for the direct method, GLUE, which is equal to the POC values (represented as percent, e.g. ) for the inverse method.

By a glance at Figures 5 and 6, it can be seen that there is not a meaningful correlation between the limits of the histograms generated by the GLUE method and the limits of the parameter sets optimized by GA, for neither SIM nor MIM. This is an important result, implying that optimized sets of parameters cannot be generally counted on as the limits of the parameters' distributions of the probabilistic uncertainty assessment method.

Moreover, we know that the shape of the distribution indicates the degree of uncertainty. In other words, peaked, convex, and sharp distributions are associated with well-detectable parameters (i.e. lower uncertainty), while flat/semi-uniform posterior distributions indicate more vagueness about the results (i.e. more uncertainty). Therefore, it can be inferred from Figures 5 and 6 that the inverse method is not able to identify not only parameters' distributions but also parameters' limits. On the other hand, the direct method can accurately provide the meaningful distributions of the model's parameters, causing a better understanding of the model's parametric uncertainty. It can be seen in the figures that in some cases (e.g. parameter *S*1 in the first hydrological model, Figure 5(c)), the inverse method does not even calculate two limits for the parameter set, but it results in just one number. In other words, the lower and upper limits are the same. This happens for both standard and modified operations (Figures 5 and 6).

To investigate the uncertain parameter bounds resulted from direct and inverse methods and to compare them more accurately, Table 3 shows the POC. This amount represents the fraction of samples generated by the direct method that falls into the bound calculated by the inverse method (Figures 5 and 6). As can be seen in Table 3, there is not any specific and concrete correlation between the bounds calculated by the inverse method and the histograms or parameters' distributions resulting from the direct method.

Simple Guo model . | Enhanced Guo model . | ||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Parameter . | MIM . | SIM . | Parameter . | MIM . | SIM . | ||||||||

60% . | 65% . | 70% . | 60% . | 65% . | 70% . | 60% . | 65% . | 70% . | 60% . | 65% . | 70% . | ||

C | 45 | 60 | 39 | 100 | 100 | 100 | SF | 100 | 0 | 0 | 100 | 100 | 0 |

SC | 17 | 48 | 63 | 100 | 100 | 100 | SP(1) | 17 | 0 | 60 | 0 | 40 | 0 |

S1 | 40 | 0 | 45 | 0 | 0 | 0 | T_{s} | 100 | 0 | 61 | 100 | 100 | 100 |

T_{m} | 0 | 95 | 100 | 10 | 12 | 9 | |||||||

K_{s} | 44 | 16 | 2 | 0 | 0 | 0 | |||||||

K_{sn} | 47 | 0 | 74 | 25 | 11 | 15 | |||||||

K_{g} | 37 | 8 | 57 | 23 | 0 | 0 | |||||||

K_{sn1} | 0 | 0 | 0 | 77 | 80 | 74 | |||||||

S_{max} | 0 | 0 | 0 | 71 | 0 | 88 | |||||||

P_{1} | 16 | 61 | 91 | 0 | 0 | 98 |

Simple Guo model . | Enhanced Guo model . | ||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Parameter . | MIM . | SIM . | Parameter . | MIM . | SIM . | ||||||||

60% . | 65% . | 70% . | 60% . | 65% . | 70% . | 60% . | 65% . | 70% . | 60% . | 65% . | 70% . | ||

C | 45 | 60 | 39 | 100 | 100 | 100 | SF | 100 | 0 | 0 | 100 | 100 | 0 |

SC | 17 | 48 | 63 | 100 | 100 | 100 | SP(1) | 17 | 0 | 60 | 0 | 40 | 0 |

S1 | 40 | 0 | 45 | 0 | 0 | 0 | T_{s} | 100 | 0 | 61 | 100 | 100 | 100 |

T_{m} | 0 | 95 | 100 | 10 | 12 | 9 | |||||||

K_{s} | 44 | 16 | 2 | 0 | 0 | 0 | |||||||

K_{sn} | 47 | 0 | 74 | 25 | 11 | 15 | |||||||

K_{g} | 37 | 8 | 57 | 23 | 0 | 0 | |||||||

K_{sn1} | 0 | 0 | 0 | 77 | 80 | 74 | |||||||

S_{max} | 0 | 0 | 0 | 71 | 0 | 88 | |||||||

P_{1} | 16 | 61 | 91 | 0 | 0 | 98 |

In Table 4, the results of implementing parameter sets optimized by the inverse method are presented. Three different metrics, namely ARIL, POC, and NUE, are employed to examine the fitness and accuracy of uncertainty bounds generated by the inverse method.

. | . | Simple Guo model . | Enhanced Guo model . | ||
---|---|---|---|---|---|

. | . | SIM . | MIM . | SIM . | MIM . |

60% | ARIL | 0.0029 | 0.8056 | 1.7337 | 0.7447 |

POC | 0.0020 | 0.6044 | 0.7278 | 0.7361 | |

NUE | 0.6729 | 0.7502 | 0.4198 | 0.9884 | |

65% | ARIL | 0.0029 | 0.8085 | 1.6426 | 0.5422 |

POC | 0.0020 | 0.6501 | 0.7000 | 0.6639 | |

NUE | 0.6729 | 0.8041 | 0.4261 | 1.2244 | |

70% | ARIL | 0.0029 | 0.9594 | 1.7050 | 0.6755 |

POC | 0.0020 | 0.7236 | 0.7167 | 0.7472 | |

NUE | 0.6729 | 0.7542 | 0.4203 | 1.1061 |

. | . | Simple Guo model . | Enhanced Guo model . | ||
---|---|---|---|---|---|

. | . | SIM . | MIM . | SIM . | MIM . |

60% | ARIL | 0.0029 | 0.8056 | 1.7337 | 0.7447 |

POC | 0.0020 | 0.6044 | 0.7278 | 0.7361 | |

NUE | 0.6729 | 0.7502 | 0.4198 | 0.9884 | |

65% | ARIL | 0.0029 | 0.8085 | 1.6426 | 0.5422 |

POC | 0.0020 | 0.6501 | 0.7000 | 0.6639 | |

NUE | 0.6729 | 0.8041 | 0.4261 | 1.2244 | |

70% | ARIL | 0.0029 | 0.9594 | 1.7050 | 0.6755 |

POC | 0.0020 | 0.7236 | 0.7167 | 0.7472 | |

NUE | 0.6729 | 0.7542 | 0.4203 | 1.1061 |

As mentioned in the methodology section, the inverse method implements GA to optimize the set of parameters with a specific threshold of POC (0.6, 0.65, or 0.7) as its optimization constraint; while its fitness function is NUE. Therefore, the algorithm tries to find the best set of parameters that maximizes the amount of NUE, while its POC metric is greater than or equal to the threshold. Table 4 shows that for both hydrological models reformulated with SIM, the GA is not able to find a suitable set of parameters and to satisfy the models' POC constraint. For the simple Guo, as the amounts of ARIL show, GA results in a very narrow bound, which cannot satisfy the constraint. The low value of ARIL and narrow bound can mathematically provide a relatively high amount of NUE (Equation (4)), but the result is not reliable at all because the lower and upper limits of streamflow uncertainty bound are equal or their difference is very close to zero. In the second hydrological model reformulated with SIM, GA is able to satisfy the constraints, but the amount of ARIL is very high. Actually, the lower limit of streamflow bound is equal to zero (Figure 7(a)). Therefore, it can be inferred that in the case of SIM, the inverse method cannot provide an accurate and reliable uncertainty bound.

Looking at the results of both hydrological models reformulated by MIM, it can be understood that they can be counted on as reliable uncertainty bounds. In all cases, optimization constraint is satisfied and the amount of NUE is maximized with the best set of parameters (Figure 7(b)).

## DISCUSSION AND CONCLUSION

In the current paper, to fill a research gap of comparative studies in the realm of uncertainty assessment methods, two frameworks of parametric uncertainty simulation, namely direct and inverse methods, have been described and evaluated.

To conduct this comparison, two hydrological models and two case studies have been employed. In addition, the GLUE uncertainty assessment method has been selected as the direct framework because of its feedforward information flow and its selection of the behavioral parameter sets considering the pre-specified acceptance rates.

To assess the parametric uncertainty using the inverse framework, an optimization procedure reformulated by interval mathematics has been used. To achieve the optimized lower and upper bounds of the parameters, a real-coded GA is used with different POC values as the fitness functions.

GLUE, like most of the sampling simulation methods, is based on the key hypothesis of equifinality, and in these models, the generation of behavioral parameters and their posterior distributions are based on this concept. However, as the inverse method results in a single optimized set of parameters, its main difference with the direct method is the absence of equifinality assumption.

Nasseri *et al.* (2013, 2014a) showed that the uncertainty analyzed by fuzzified direct frameworks are comparable with GLUE and other deterministic methods. Though the mathematical operators (SFM and MFM) are similar to what is used in the current research, the parameters' bounds have been extracted from their fuzzy membership function and uncertainty functions have been inferred from their posterior distributions based on the GLUE uncertainty assessment. So, the point of compatibilities in their output uncertainties was the similarity in their parametric uncertainties. This is the exact point of the contrast of the parametric uncertainty between the presented direct and inverse methods in the current analytical/critical research.

The current inverse framework resulting optimized pair of the lower/upper limits of parameter sets transferred the uncertainty assessment method from a probabilistic/possibilistic topic to a single-objective optimization problem. So, based on the results and statistical metrics, it can be concluded that the direct method can generally provide more reliable insight in different hydrological models and with different datasets than the current inverse uncertainty assessment framework.

As one possible future research idea, it seems that using multi-objective optimization to create a set of acceptable behavioral parameter sets using the indirect framework may be a solution to improve its statistical performance. Moreover, testing the findings of this research with other hydrological models and under higher temporal resolutions (e.g. daily) can further verify their accuracy and generality.

## ACKNOWLEDGEMENTS

The authors acknowledge the anonymous reviewers for their invaluable, constructive, and helpful comments.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/hydro.2020.190.