Assessing the temporally dynamic parameters in hydrological models: dynamic operations and evolutionary processes

Abstract. The temporal dynamics of parameters can compensate for structural defects of hydrological models and improve the accuracy and robustness of the streamflow forecast. Given the parameters usually estimated by global optimization algorithms, a critical issue, however, which received little attention in the literature, is that the possible failure in finding the global optimum might lead to unreasonable parameter values. This may cause the poor response of the dynamic parameters to time-varying catchment characteristics (such as seasonal variations of land cover). In this regard, we propose a framework for identifying the difficulty of finding the global optimum for dynamic hydrological model parameters by investigating their evolutionary processes. Specifically, the probability distributions of the violin plots and the divergence measure of the polylines in the parallel coordinates are applied and developed to configure the evolutionary processes in the individual parameter spaces and multi-parameter space, respectively. Also, a complete solution for the dynamic operation of parameters is proposed. Furthermore, clustering operations, calibration scheme and correlation between parameters are further discussed. The results showed that the performance of the hydrological model with dynamic parameters achieves a significant improvement. However, the response of individual parameters (even high-sensitive parameter) to dynamic catchment characteristics is generally poor. The main reasons can be primarily attributed to the complexly linear and nonlinear correlation between parameters and poor ability in finding the global optimum. In this regard, the dynamic parameter set instead of individual dynamic parameters is suggested to extract dynamic catchment characteristics. Importantly, we found that the properties of hydrological-model parameters, including identifiability, sensitivity, correlation and the ability to find global optimum, interact with the response of parameters to the dynamic catchment characteristics. The ability to find global optimum has a significant influence on the hydrological model performance with dynamic parameters. Hence, the ability to find global optimum is suggested as one of the essential properties of the hydrological model parameters. The study provides a valuable benchmark for temporally dynamic parameters in hydrological models.



Introduction
This supporting information includes two sections that support the analysis. The introduction of the techniques is used to support the Methodology section in the main manuscript. The Results section is used to supplement the Results section in the main manuscript.

SCE-UA algorithm
The shuffled complex evolution approach (SCE-UA), as an effective global optimization method, is a commonly used algorithm, because it is open source and was the first algorithm aimed specifically at calibrating hydrological models (Khakbaz and Kazeminezhad, 2012;Eckhardt and Arnold, 2001;Duan et al., 1994;Sorooshian et al., 1993). The technical details about the SCE-UA can be shown in the flowchart (see Figure S1) (Duan et al., 1994). In the SCE-UA, the upper limit of the objective function evaluation is set to 10,000 times.
The stopping criteria are applied to prevent premature termination, while also avoiding unnecessary computations. There are three stopping criteria in SCE-UA. (1) The maximum number of function evaluations has been reached; this means that the optimization algorithm may not converge at the end of the run. (2) The population has prematurely converged to a pre-specified small geometric range, which indicates that the global convergence has failed to local optima. The method has failed to a zone that it thinks is not amongst the best. It cannot know that it is in a global optimum. (3) The improvement of the best point in the last loops is less than the specified threshold. It indicates that global convergence has been achieved. More details on the SCE-UA are provided in the supporting information (section 1.1). Certain stopping criteria of the global optimization algorithm lead to possible failure in finding the global optimum, which may lead to abnormal or unreasonable optimal parameters. Namely, the optimized parameters cannot effectively characterize the physical systems of the hydrological models (or catchment response behavior). Figure S1. The flowchart of the SCE-UA algorithm (Duan et al., 1992(Duan et al., , 1994.

Fitness landscape
A potent metaphor in global optimization is the fitness landscape (Aldrich, 1997). In evolutionary biology, fitness landscapes have been developed to visualize the relationship between the genotypes or phenotypes in a given population and their corresponding reproduction probability (Wright, 1932;Kauffman, 1993;Dawkins, 1997;Gavrilets, 2004). The idea of such visualizations goes Failure: The population has converged to a prespecified small geometric range.
Success: The best point has improved in last loops by less than the threshold.
back to Wright (1932), who used level contours diagrams in order to outline the effects of selection, mutation, and crossover on the capabilities of populations to escape local optimal configurations. The evolutionary algorithm research community has widely adopted the fitness landscapes as a relation between individuals and their objective values (Mitchell, 1998).

Violin plot
A violin plot is a combination of a Box Plot and a Density Plot showing more details of data distribution. As shown in Figure S2, the thick black bar in the center represents the interquartile range. The white dot represents the median. The thin black line is extended from the thick black bar and represents the 95% confidence intervals. On each side of the thin black line is a kernel density estimation to show the distribution shape of the data. Wider sections of the violin plot represent a higher probability that members of the population will take on the given value; the skinnier sections represent a lower probability (Hintze and Nelson, 1998). The violin plots can exactly show the kernel density distribution, avoiding the overlapping traditional density plot occur to become difficult to identify. Moreover, unlike bar graphs with means and error bars, violin plots contain all data points, which makes them an excellent tool to visualize samples of small sizes. Violin plots are perfectly appropriate even if the data do not conform to normal distribution. They work well to visualize both quantitative and qualitative data.

Parallel coordinates
The theory of parallel coordinates has been developed rigorously, and its point-line duality has been successively generalized to higher dimensions (Heer et al., 2010;Shenghui and Mueller, 2015). In order to show a set of points in an n-dimensional space, a backdrop is drawn consisting of n parallel lines, typically vertical and equally spaced. That is, a point in n-dimensional space is represented as a polyline with vertices on the parallel axes. Each axis corresponds to a variable, and each data item, having values for all variables, is represented as a series of line segments intersecting the axes at the corresponding values (Zhou and Weiskopf, 2018). An excellent complication of current research on parallel coordinates can be found in the state-of-the-art report by Heinrich and Weiskopf (2013).

Maximal information coefficient (MIC)
The MIC was proposed by Reshe et al. (2011) and did not rely on the distributional assumptions of the datasets. This measurement approach captures extensive mutual information of the variables for functional and non-functional relationships. For functional relationships, the MIC algorithm provides a score that is roughly equivalent to the coefficient of determination ( 2 ) of the datasets.
However, its statistical power reduced in detecting the associations in setting with a low sample size (Heller et al., 2012).

Scale
Density plot width = probability 95% confidence interval Interquartile range

Median
Higher probability Lower probability 5 1.6 HYMOD model structure For illustration purposes, the HYMOD model (Moore, 1985;Wagener et al., 2001;Vrugt et al., 2002;Yadav et al., 2007;de Vos et al., 2010;Pathiraja et al., 2018), with a simple and commonly used lumped model structure, is utilized. The HYMOD model consists of a simple rainfall excess model based on the probability-distributed moisture store, which characterizes the catchment storage as a Pareto distribution of buckets of varying depth as the soil moisture accounting component. It routes through three parallel tanks for quick flow and a tank for slow flow and required five adjustable parameters: , , , and . and are state variables characterizing the upper soil moisture content; is actual evapotranspiration which is calculated by linear correlations between the soil moisture state and the potential evapotranspiration; is effective precipitation; is excess precipitation to routing module generated from overflow of soil moisture accounting component; See Moore (1985) for a detailed description of the soil moisture accounting model; 1 , 2 , 3 and are the state variables of the individual tanks of the routing module; and are the flow values generated from the quick-and slow-flow tanks, respectively.   Figure S6. Evolutionary processes of dynamic parameters in multi-parameter space in the Xunhe basin.  Figure S7. Evolutionary processes of dynamic parameters with magnified details on the axes in multi-parameter space in the Mumahe basin.  Figure S8. Evolutionary processes of dynamic parameters with magnified details on the axes in multi-parameter space in the Xunhe basin.