Sensitivity and Interaction Analysis Based on Sobol' Method and Its Application in a Distributed Flood Forecasting Model

Sensitivity analysis is a fundamental approach to identify the most significant and sensitive parameters, helping us to understand complex hydrological models, particularly for time-consuming distributed flood forecasting models based on complicated theory with numerous parameters. Based on Sobol' method, this study compared the sensitivity and interactions of distributed flood forecasting model parameters with and without accounting for correlation. Four objective functions: (1) Nash–Sutcliffe efficiency (ENS); (2) water balance coefficient (WB); (3) peak discharge efficiency (EP); and (4) time to peak efficiency (ETP) were implemented to the Liuxihe model with hourly rainfall-runoff data collected in the Nanhua Creek catchment, Pearl River, China. Contrastive results for the sensitivity and interaction analysis were also illustrated among small, medium, and large flood magnitudes. Results demonstrated that the choice of objective functions had no effect on the sensitivity classification, while it had great influence on the sensitivity ranking for both uncorrelated and correlated cases. The Liuxihe model behaved and responded uniquely to various flood conditions. The results also indicated that the pairwise parameters interactions revealed a non-ignorable contribution to the model output variance. Parameters with high first or total order sensitivity indices presented a corresponding high second order sensitivity indices and 2925 correlation coefficients with other parameters. Without considering parameter correlations, the variance contributions of highly sensitive parameters might be underestimated and those of normally sensitive parameters might be overestimated. This research laid a basic foundation to improve the understanding of complex model behavior.


Introduction
A distributed hydrological model is an essential tool and has been widely used in hydrology and water resources management due to its good ability in simulating temporal and spatial variation of hydrological variables [1,2].In practice, the performance of a hydrological model is highly correlated to the setting of model parameters.Because of the existing phenomenon called "equifinality" [3,4], i.e., different parameter sets resulting in the same model simulations, there is large uncertainty in the application of hydrological models.Various selection of parameters will induce a large variety of simulation results; on the other hand, considering that the distributed hydrological models always contain complex structures with large number of parameters, the optimization choice of parameters is considered to be a difficult and time-consuming task.
However, the sensitivity analysis (SA), which involves testing the relative impact of model parameters on the runoff prediction, has been proved to be an efficient way to identify the parameters in hydrological models and has been applied in many studies [5,6].Some studies have stated that SA can improve parameter calibration efficiency, avoid over-parameterization [7], reduce model uncertainty [8][9][10], aid in understanding model structure [11,12], and generate more accurate model outputs [6,13].
In general, SA methods can be broadly grouped into two types: local methods and global methods [6,7,14,15].The local methods only make one parameter change at one time, while the other parameters are kept as constant [7,16].The advantage of these methods is that only a small computation is needed during their application, while the drawback is that they cannot determine the interaction between parameters, thus may result in a large bias of SA.On the contrary, the global methods overcome the drawbacks of local methods by changing all of the model parameters at the same time [6,7,14,15].This approach can consider the linear or nonlinear interaction between model parameters, and it is suggested to be more appropriate to complex and physically-based distributed hydrological models than local methods.Many global methods have been proposed and used, such as Regional Sensitivity Analysis (RSA) [17], Parameter Estimation Software (PEST) [18], Analysis of Variance (ANOVA) [14,19], Fourier amplitude sensitivity test (FAST) [15,20], Sobol' method [6,21], etc.Among these methods, Sobol' is proved to be one of the most robust and effective tools to describe individual and cooperative sensitivities [6,14,15,[20][21][22][23], and it is widely applied not only in hydrological models [14,21,[23][24][25], but also in environmental models [26], chemical models [27], biological models [28], ecological models [29], etc.
Tang et al. [14] and Gan et al. [15] compared Sobol' method with other methods to analyze the sensitivities of the SAC-SMA model.Tang et al. [14] conducted a contradistinctive study between Sobol' method and three other methods for two watersheds with two evaluation metrics including non-transformed root mean square error (RMSE) and transformed RMSE (TRMSE).Gan et al. [15] also compared Sobol' method with eight other methods for first, second, and total order sensitivities based on the mean absolute error (MAE) of the simulated and observed flow measurement, in the South Branch Potomac River basin in the U.S. Reusser et al. [20] examined the dominant parameter and parameter interaction of TOPMODEL and a complex WaSiM-ETH model by the comparison of the Sobol' method with FAST and extend FAST.Zhan et al. [6] proposed a highly efficient integrated sensitivity analysis tool on the basis of Sobol' method for a DTVGM model in the Huaihe River Basin, China.The findings revealed the parameters' sensitivities varied with regard to three different evaluation criteria, namely, water balance coefficient, Nash-Sutcliffe coefficient, and correlation coefficient.More recently, Zhang et al. [23] improved the understanding of the SWAT model through Sobol' sensitivity analysis method at the Yichun River Basin in China, by exploring model performance based on four metrics-RMSE, Nash-Sutcliffe, runoff coefficient error, and slope of the flow duration curve error under dry, moderate, and wet climate conditions.
Previous studies have validated the robustness and ease of implementation of Sobol' method to generate global sensitivity analysis (including first, second, and total order) in numerous hydrological models.Little attention has been paid to the independent assumption of parameters in Sobol' method.Most researchers [6,14,15,23] did not consider the correlation underlying the parameters, even though great second order sensitivity indices were displayed.For this case, a generalized Sobol' method [30] that considers correlation is more suitable for correlated parameters.Besides, as Wagener et al. [24] pointed out the significance of multiple objective functions during the model performance assessment, researchers realized the importance of reducing the uncertainty of sensitivity analysis by multi-metrics, various characterized catchments, and climate conditions.However, most distributed models were carried out with daily or monthly data [6,15,23], seldom with hourly data.In addition, even though multiple objective functions were strengthened, little attention was paid to the peak discharge timing and magnitude.
In this study, the Morris screening method and the extended [31] and generalized [30] Sobol' method for uncorrelated and correlated parameters are implemented to the distributed flood forecasting model named the Liuxihe model by four objective indicators that can explicitly characterize the flood processes: (1) Nash-Sutcliffe efficiency (ENS); (2) water balance coefficient (WB); (3) peak discharge efficiency (EP); and (4) time to peak efficiency (ETP).Different types of flood magnitude (small, medium, and large) are considered to compare the sensitivity and interaction of parameters in the Liuxihe model.The Nanhua Creek catchment, which is located in the Pearl River basin, China, is selected as the study area, and the hourly rainfall-runoff data is used.

Methodology
The sensitivity and interaction analysis based on Sobol' method can be carried out by the following steps: (1) screening ineffective parameters out by Morris method [5]; (2) identifying the parameters correlations through Generalized Likelihood Uncertainty Estimation (GLUE) [3,4]; (3) sampling parameter sets with and without considering the identified parameters correlations; (4) evaluation of the sensitivity indices of model correlated and uncorrelated parameters.

Morris Screening
Morris [5] proposed a "one-at-a-time" (OAT) design to screen the most influential parameters easily and effectively with a small amount of computation, which is now named the Morris screening method or Elementary Effect (EE) method.Although defined as a local SA method, the Morris screening method is broadly taken as a preliminary step to separate out the negligible parameters and decrease the huge computation demand for dozens of model parameters [6,[32][33][34].The detail of the Morris screening method is given below.
Let a model with n parameters , the elementary effect of the ith parameter can be estimated by merely changing the ith parameter with a step i Δ while the rest of the parameters are kept unchanged.This process can be commonly defined as: where i Δ is a value in the set {1 ( ) are a random and the designed sample set, respectively, both of them within the predefined n-dimensional hypercube space.Two sensitivity measures are proposed to investigate the parameters [5]: μ and σ , which evaluate the overall effect of each parameter on the model output, and the parameters' higher order effects caused by nonlinear or interactive relationship between parameters, respectively [35].As suggested by Morris, R randomly sampled different trajectories are preferred to assess the elementary effects.Statistically speaking, μ and σ are the mean and the standard deviation of the R × (n + 1) elementary effects.When the model is non-monotonic, Campolongo et al. [35] improved the estimation of μ by * μ with the absolute value of ( ) A large * i μ underlines the fact that the ith parameter is a relatively influential parameter and a large i σ underlines a strong nonlinear or interactive relationship between the ith parameter and the others [6].

Sobol' Method for Uncorrelated Parameters
Sobol' method is one of the widely used global sensitivity analysis methods based on variance decomposition [6,14,21,22].Assume that the model parameters are independent or uncorrelated, the total variance of model output is composed of the variance from individual parameters, and the variances from cooperative parameters.The proportion of variance resulting from individual and cooperative parameters to the total variance is applied to estimate the first order and interactive sensitivity indices of parameters.The variance decomposition can be expressed as: where V(y) denotes the total variance of model output; Vi denotes the first order contribution of the ith model parameter; Vij denotes the pairwise interactive effect of the ith and jth parameters; and n denotes the number of model parameters.Obviously, the first, second, and total order of sensitivity indices can be respectively expressed as: where V~i denotes the contribution of all except the ith parameter.The variances in Equation ( 5) can be estimated by approximate Monte Carlo numerical integrations [14].The computation requirement of the original Sobol' method to evaluate the first and total order sensitivity indices is m × (2n + 1), where m denotes size of sample sets, and the second order sensitivity indice cannot be calculated in such a case with tremendous model runs.Saltelli [31] proposed an extended method that generated more robust and satisfactory sensitivity outputs by double calculations of the first, second, and total order sensitivity indices with m × (2n + 2) model runs.The estimation for Vi, Vij, V~i, and V(y) is described as follows: where m denotes the sample size; Θ presents all the sampled parameters selected from sample 1 except the ith one from sample 2. The combinations of sampled parameters input from sample 1 and sample 2 produce double times calculation of model outputs (see Equations ( 11), (14), and ( 18)).Equations ( 9)-( 11) can be used to estimate the first order sensitivity indices according to Equation (6).Similarly, Equations ( 12)-( 14) and Equations ( 16)-( 18) can be used to estimate the second and total order sensitivity indices according to Equations ( 6) and (7), respectively.Before the sensitivity analysis is implemented, the sampling scheme must be carefully selected since it has an important effect on the convergence efficiency of Monte Carlo approximation.As suggested by Zhan et al. [6], Tang et al. [14], and Saltelli [31], the Latin Hypercube sampling (LHS) [7,21,36,37], based on Monte Carlo, which generates an efficient estimation of model output through a stratified sampling method, is adopted in this study.

Sobol' Method for Correlated Parameters
The original Sobol' method can obtain good efficiency in the model of independent or uncorrelated parameters, which might be unsuitable for complex hydrological models with high nonlinearity and interaction.Kucherenko et al. [30] proposed a new method based on Sobol' method to estimate the global sensitivity indices with dependent parameters, and made it appropriate for complex hydrological models.

GLUE Method for Correlation Identification
Before the implementation of the Sobol' method for correlated parameters, correlation in parameters must be investigated in advance.In this paper, the Generalized Likelihood Uncertainty Estimation (GLUE) method [3,4] is employed to detect the correlated parameters [38].The detailed procedures are listed as follows: (1) randomly sample a huge number of parameter sets from the n-dimensional hypercube space, the same as Sobol' method.Because no supposition about the a priori distribution of these parameters is essential, a uniform distribution is a common choice; (2) calculate corresponding model outputs and likelihood measures with these sample sets; (3) define a likelihood measure threshold to separate behavioral and non-behavioral parameters sets, termed as b Θ and nb Θ , respectively; (4) estimate the mean i μ , standard derivation i σ , and pairwise correlation coefficient ij ρ of each parameter vector from b Θ .The pairwise correlation coefficients then constitute the correlation matrix c M .
The first three steps are abstracted from the GLUE method, which is applied to generate the behavioral parameters.The subsequent uncertainty analysis with these behavioral parameter sets is beyond the scope of this paper.Details about the GLUE method can be seen in Beven and Binley [3,4].The fourth step shows the in-depth analysis of these behavioral parameters sets and relevant works referred to are Hornberger and Spear [17] and Ratto et al. [38].

Generation of Correlated Sample Sets
Based on the identified parameters correlations by the GLUE method, correlated parameter sets can be generated through a general procedure introduced by Kucherenko et al. [30].The procedures are implemented as follows: (1) randomly generate a parameters set from standard normal distribution, termed Φ by LHS.If the sample size is m, and parameters number is n, then the dimension of Φ is m n × .Since Φ is produced aimlessly, the parameters vector in Φ is definitely independent and uncorrelated; (2) decompose the correlation matrix c M by an appropriate decomposition method.Cholesky is preferred when the correlation matrix is positive definite, whereas this is not easy to promise.Therefore, orthogonal decomposition method, which is more universal, is employed: where Q is the orthogonal matrix and R is the upper triangular matrix.Both dimensions are n n × ; (3) transform the uncorrelated standard normal sample sets Φ to correlated standard normal sample sets Ψ : where the dimension of Ψ is still m n × ; (4) derive the correlated parameters sets ( ) where ci θ and i ψ are the ith parameter vector from c Θ and Ψ , respectively; and i μ and i σ are the mean and standard derivation of the ith parameter vector, respectively.

Generalization of Sobol' Method
Concluded from Kucherenko's method, the Sobol' method to estimate the first and total order sensitivity indices for pairwise correlated parameters is slightly different from that for uncorrelated parameters, as Zhao et al. summarized in [39].
For uncorrelated parameters, the first order sensitivity indice can be obtained from Equations ( 6) and ( 9)- (11).Equation (11) estimates the variances of model outputs from the ith parameter with the explicit meaning that the ith parameter in sample 1 is replaced by the corresponding ith parameter in sample 2 or vice versa.Similarly, the first order sensitivity indice for correlated parameters can be derived from replacing several vectors which must contain the ith parameter in sample 1 from those in sample 2. If the pairwise correlation is considered, then the dimension of replaced vectors is equal to 2, with the following equation: where 1

~i,k
Θ presents all the sampled parameters are selected from sample 1 except the ith and kth ones from sample 2. The others keep the same meaning as introduced above.Attention has to be paid to the fact that for correlated parameters, samples 1 and 2 are generated from procedures introduced in Section 2.3.2.Analogously, Equation ( 18) for uncorrelated parameters expressed the connotation that the estimated variance from all except the ith parameter is the product of model output from sample 1 or 2 and the ith parameters in the samples replaced by those from sample 2 or 1.Therefore, the equation to estimate the variance for all except the ith parameter for correlated parameters can be presented as:

Liuxihe Model
In this study, the Liuxihe model, which is a physically based distributed watershed model proposed by Chen et al. [40] that has been applied successfully to numerous catchments in China [41][42][43][44], is selected as the used distributed hydrological model.This model comprises five sub-models: basin digitization, evapotranspiration, runoff generation, runoff routing, and parameter derivation.In this model, the distributed model parameters are acquired directly from the DEM and channel cross-sectional properties from available free public data (also see the details in Chen et al. [40]).

Basin Digitization
The Liuxihe model divides the catchment horizontally into three topographic cell types, namely hillslope, channel, and reservoir.The catchment is also divided vertically into three layers, namely canopy, soil, and sub-surface.Runoff generation and flow routing is carried out independently for each grid cell.

Evapotranspiration
In the Liuxihe model, evapotranspiration is determined by the triplex pattern: where E is the actual evaporation (mm); v is the evaporation coefficient relevant to vegetation type; FC is the field water capacity (%); Wp is the soil water content at wilting point (%); θ is the temporary soil water content (%); and Ep is the potential evaporation.Cells that are entirely open water (reservoirs and channels) have an evapotranspiration coefficient of 1, while other cells have a value less than 1.

Runoff Generation
The model calculates runoff generation using the saturation-excess mechanism.Rainfall excess, or net rainfall, is derived from the difference between actual rainfall and evapotranspiration.Positive rainfall excess results in infiltration to the soil.Rainfall excess becomes surface runoff only after the soil becomes saturated.
Interflow in the soil layer is calculated according to water mass conservation: where h is the soil thickness (mm); L is the cell width (m); Qlat is the interflow per unit width (m 2 /s); r is the rainfall excess vertically and adjacent interflow input per unit width (m 3 /s); vlat is the interflow velocity (m/s); and Qper is the percolation per unit width (m 2 /s).
According to Darcy's law, the interflow velocity vlat is: where α is the hillslope cell slope in degree (°); 0 S is the hillslope cell slope in radians; K is hydraulic conductivity at water content (mm/h); and θ is derived according to Campbell [45]: where SC is the saturation water content (%); Ks is the saturated hydraulic conductivity (mm/h); and b is the soil porosity characteristic coefficient.

Runoff Routing
The Liuxihe model assumes three types of flow routing according to the three water sources: surface flow, interflow, and groundwater.Also, there are three different flow routing methods corresponding to the three different grid cell types of hillslope, channel, and reservoir.Surface flow occurs in the hillslope cells; there are total five different routing categories.
For hillslope routing, the kinematical wave approximation (the 1-D Saint-Venant equations) is adopted with the assumption that the hydraulic cell slope denotes the hillslope: where Q is the surface flow (m 3 /s); h is the rainfall excess (mm); and q is the lateral flow per unit width (m 2 /s).In addition, according to the Manning equation, Q can be derived by: where For river routing, the 1-D diffusive wave method is used.The model assumes a trapezoid cross-section shape of area A (m 2 ), water depth h (m), river bottom width w (m), and side slope β.The continuity equation is: According to the Manning equation and continuity equation, river routing can be calculated by: where and χ is the wetted perimeter of the cross-section (m).
For reservoir routing, inflow to the reservoir is calculated as the sum of all inflow cells to the reservoir based on Chen's idea that discharge on such cells can be regarded as inflow without lapsing.
For interflow routing, the model determines the interflow by Equation ( 25) from adjacent soil cells since the runoff generation process already takes the interflow routing into consideration.
For groundwater routing, the linear reservoir method is calculated as: (1 ) where Qg,t and Qg,t + ∆t are groundwater flow at time t and t + ∆t, respectively (m 3 /s); and Kg is the groundwater recession coefficient.

Parameters Derivation
Based on the evapotranspiration, runoff generation, and routing of the Liuxihe model introduced above, there are 12 relevant parameters in this model (see Table 1).These 12 parameters can be categorized into five different types related to various aspects, including climate, land use, soil type, terrain, and others.A vital and challenging task before the sensitivity and correlation analysis is the determination of appropriate ranges and distribution for the parameters [46,47].Tong and Graziani [48] and Zhan et al. [6] suggested that the ranges of parameters should be decided in advance considering their significant effect on the sensitivity results.Fortunately, the Liuxihe model, as one of the physically based distributed hydrological models, developed a reasonable method to derive parameters from terrain data.The parameters can be derived from the following two procedures, i.e., initial parameter determination and parameter adjustment.Soil related parameters (Table 1) were pre-fixed by Soil Water Characteristics Hydraulic Properties Calculator [49].Other parameters were defined according to references.After the initial parameters values were pre-decided, the parameters were adjusted mildly in turn by multiplying a change factor f (f > 0).The best-performing parameters set was then selected judging from the model performance.It is noteworthy that the Liuxihe model is physically based, so the predetermined parameters can perform well to some degree.A slight adjustment is aimed at improving the performance.Details about Liuxihe model parameters derivation can be obtained from Chen et al. [40].
After the calibration through the above predesign and slight adjustment processes, the large sample sets for sensitivity analysis based on Sobol' method are generated by multiplying a change factor matrix with the calibrated parameters.The change factor matrix is sampled from the normal distribution and ensures the corresponding sampled parameters are within each ranges listed in Table 1.

Objective Functions
The assessment of parameters sensitivity depended greatly on the choice of objective functions.In this study, four widely used metrics were used to test model performance: (1) Nash-Sutcliffe efficiency (ENS) [50]; (2) water balance coefficient (WB); (3) peak discharge efficiency (EP); and (4) time to peak efficiency (ETP).These four metrics were calculated according to the following equations: where Q(t) was observed or measured flow; ( ) Q t was simulated or forecast flow; Q was the mean of observed flow; p T and p T were the simulated and observed time to peak, respectively; and p Q and p Q were simulated and observed peak discharge, respectively.The ENS is widely used to assess overall hydrograph goodness-of-fit [6,7,23,24,51,52].The WB shows the gap between the total simulated flood volume and total measured volume.This objective function is merely influenced by runoff generation process and is extensively applied to the water resources management [6,24,53].The EP and ETP are the other two essential elements in flood forecasting, and these two elements are also crucial factors in flood prevention, which is adopted in the "Forecasting norm for hydrology intelligence (GB/T 22482-2008) of China".

Study Area
The Nanhua Creek catchment (Figure 1) is a minor tributary of the upper Pearl River in China.It is located between 112°37′ E ~ 113°05′ E and 24°53′ N ~ 25°16′ N with a drainage area of 120 km 2 and a length of 82 km.The Nanhua Creek catchment belongs to the monsoon climatic zone, with cold and dry winters and hot and humid summers.The long-term mean annual precipitation in the catchment is about 1300-1500 mm, the majority of which occurs in summer.The steep slope of the catchment, combined with concentrated rainstorms, makes this area prone to flash runoff response.

Data
The Liuxihe model requires elevation (DEM), land use, and soil type data.In this study, DEM was obtained from NASA Shuttle Radar Topographic Mission (SRTM, http://srtm.csi.cgiar.org/)at a resolution of 90 m × 90 m.Land use data came from the USGS (http://landcover.usgs.gov/)and soil type data were from the FAO/UNESCO (http://www.isric.org/).Eighteen storm floods with hourly measured rainfall data of three rain gauge stations and discharge data of the catchment outlet hydrological station (Figure 1) were collected from 1980 to 1987 to calibrate the model parameters.Ten storm floods from 1990 to 1995 were applied to validate the model.
In order to investigate the response of Liuxihe model to various flood condition, three different types of floods with small, medium, and large magnitudes were selected from the available data set.These typical flood events with different magnitudes were selected with flood frequency of 75%, 50%, and 25%, respectively.The characteristics of these floods were listed in Table 2.

Results of Calibration and Validation
The average of the four objective functions of the 18 storms in the calibration period were 0.82, 1.08, 0.10, and 0.7 h respectively (see Table 3).According to the "Forecasting norm for hydrology intelligence (GB/T 22482-2008) of China", these satisfactory results met the second grade of the norm, which requires 0.7 < ENS < 0.9, 0.8 < WB < 1.2 and |EP| < 0.2.The calibrated results demonstrated that the Liuxihe model was suitable for flood forecasting in the Nanhua Creek catchment and the parameters were reasonable.Similar results were also obtained from the model validation period.The average of the four objective functions of the 10 storms in the verification period were 0.74, 1.06, 0.12, and 1.7 h, respectively.Even though the performance of validation was slightly worse than the calibration period, the former still met the second grade of the norm, which strengthened the confidence of the Liuxihe model and calibrated parameters.A comparison of observed and calculated flow for the three representative storms is illustrated in Figure 2. It can be seen from these hydrographs that the Liuxihe model generated good results in the rising limbs, flood peaks, and recession limbs.

Results of Morris Screening
The Morris screening method, or Elementary Effect (EE) method, was employed to qualitatively identify the unimportant parameters.The identified relatively important parameters were then investigated for the following global sensitivity analysis.Two parameters should be predefined before the application of Morris screening, i.e., level number p and trajectories number R. As suggested by Campolongo et al. [32,35], p is preferentially even and R is in general between 10 and 50.In this paper, level number p was determined to be 10 and R was 40.Subsequently, 520 parameters sets were generated randomly by uniform distribution within the rational ranges.
The Morris screening results can be presented in the form of ranking table [33,34], Tornado diagrams [33], and scatter plots [6].Scatter plots are confirmed to be the most efficient and clearest with modified means * μ as the x-axis and standard deviations σ as the y-axis, while the former two failed to illustrate the characteristics of standard deviations.The scatter plots of the 12 parameters in the Liuxihe model were compared in three different flood conditions as well as four different objective functions, as showed in Figure 3. On the basis of the fact that a larger * μ value reveals higher sensitivity, it is interesting to see that even under various conditions and objective functions of interest, the sensitive rankings could all be categorized into three rough classes: highly sensitive, sensitive, and insensitive; the classified parameters also gave consistent results.Among all the parameters, SC and FC were highly sensitive; b, Mn, h, n, and Ks were sensitive; and Ep, Kg, Wp, v, and S were insensitive.Explicit ranking among the 12 parameters differed within various situations.The top sensitive parameters were found among SC, FC and Mn from the most to the least frequent.The least sensitive parameters covered were S, v, and Wp, from the most to the least frequent.Our results were in accordance with those of Xu et al. [54], who stated that SC and FC were first-rate; h and b were second-rate; Wp, n, and Ks were third-rate; and S, Ep, Mn, Kg, and v were insensitive.

Precipitation
Observed Flow Simulated Flow Along with these sensitivity rankings reflected by * μ , the nonlinear and interactive effects between parameters reported by σ were also remarkable.All the subplots in Figure 3 indicated that SC and FC had a larger σ than the others, which inferred that the nonlinear relationship between parameter values and model outputs was more noticeable or that their interaction with other parameters was more evident.Similar results can also be obtained from b when taking WB as the objective function and Mn when ETP is the objective function of interest.It is plausible to suspect that parameters with higher modified mean values ( * μ ) had a higher possibility of presenting a correspondingly larger standard derivation value ( σ ).Zhan's work [6] supported suspicion with another distributed hydrological model (DTVGM).Therefore, attention must be paid to these high ranking parameters with large σ since the high nonlinearity and interaction with other parameters decreased the reliability of high sensitivity classification.This was also one of the dominating disadvantages of the Morris screening method: the nonlinearity and interaction between parameters cannot be differentiated.On the other hand, the μ were fortunately able to display a correspondingly lower standard derivation value σ (see Figure 3), which achieved the original goal of screening out insensitive or non-behavioral parameters with high reliability.
Therefore, the insensitive parameters of the Liuxihe model can doubtless be determined even under diverse flood conditions and objective functions.These were S, v, Wp, Ep, and Kg.These parameters were then kept constant in the subsequent GSA with negligible influences on model output.Parameters passing the validation screening, including SC, FC, Ks, b, h, n, and Mn, were chose for the further investigation of GSA based on Sobol' method.The original 12 parameters in the Liuxihe model were reduced to seven for GSA, reducing the computation by almost half.That the reason for the broad application of the Morris screening method as a first step to lower computation demand for GSA [6,[32][33][34].

Results of Correlation Identification
Since the Sobol' method is a variance-based decomposition method with the assumption that all parameters are independent, the validation of independence is definitely a prerequisite.However, this procedure was often ignored.In this paper, the GLUE method was employed to characterize the correlation relationship between the Liuxihe model parameters before the implementation of the Sobol' global sensitivity analysis method.The GLUE method must define the likelihood measure and its threshold for dividing behavioral and non-behavioral.As Beven and Binley concluded [4], the Nash-Sutcliffe efficiency (ENS) (see Equation ( 33)) was the most extensively used likelihood measure.For the sake of simplicity without loss of generality, this paper also set ENS as the likelihood measure.In addition, according to the "Forecasting norm for hydrology intelligence (GB/T 22482-2008) of China", the threshold was decided to be 0.7 in order to meet the requirement for the second grade.Based on the undemanding hypothesis about the prior distributions of the model parameters, a sample of 5000 sets was generated by the uniform distribution within the ranges of parameters.
The Pearson correlation analysis and significance level between the two parameters are listed in Table 4.The higher the correlation between the parameters, the darker gray in Table 4.The poorly correlated parameters are showed with no color.Highly correlating parameters were suggested to cut parameters' hypercube space.Positive values designated that the relevant parameters had the same effect on model performance and worked as a quotient or difference, whereas negative values designated opposite effects from two relevant parameters on model performance and worked as a product or sum [38].The same or opposite effects from parameters can be easily obtained from the principles of the model in Section 2.4.Interestingly, the SC and FC displayed the highest positive Pearson correlation, which was extremely close to 1 at the significance level of 0.01.Other parameters, like n and h, Mn and h, h and b, Mn and SC, and Ks and FC with Pearson correlation coefficients of 0.32, −0.34, 0.30, 0.26, and −0.28, respectively, were also relatively correlated at the significance level of 0.01.These highly or relatively correlated parameters were all classified as the highly or relatively sensitive parameters in the Morris screening design in Section 4.2.Individual parameters demonstrating a larger standard deviation in the Morris screening scatter plots in Figure 3 exhibited a higher correlation relationship with other parameters (e.g., SC and FC).These two methods strengthened the reality that the correlated parameters in the Liuxihe model were attentive and the straightforward application of the Sobol' method was inappropriate.Different from the figurative illustration by the Morris screening method, the GLUE method highlighted the correlation analysis by directly tabling the correlation coefficients, which laid a fundamental basis for the subsequent sensitivity analysis accounting for correlated parameters.

Results of First and Total Order Sensitivity Analysis
In order to compare the difference of the sensitivity analysis with and without considering parameters' correlation, the generalized and extended Sobol' method [30,31] were implemented to estimate the first and total order sensitivity indices of the seven parameters (i.e., SC, FC, Ks, b, h, n, and Mn) in the Liuxihe model.In this study, as suggested by Gan et al. [15], 2000 parameter sets were generated from normal distribution by the Latin Hypercube Sampling (LHS) method [7,21,36,37] for the extended Sobol' method, ignoring correlations.As for the generalized Sobol' method accounting for parameter dependency, for the sake of consistency, 2000 parameter sets were generated from normal distribution by LHS as well.The difference with the latter was that the parameters were correlated with corresponding coefficients identified from Section 4.3 and were generated through the procedures in Section 2.3.2.Therefore, the comparison analysis required a total of 32,000 model runs, half for uncorrelated parameters and half for the correlated.
Independent and dependent first and total order sensitivity indices based on the four different objective functions under the three different flood magnitudes are shown comparatively in Figure 4.The red and blue bars presented the first order sensitivity indices for uncorrelated and correlated parameters, labeled with UCS1 and CS1, respectively.The green and yellow bars exhibited the total order sensitivity indices for uncorrelated and correlated parameters, labeled with UCST and CST, respectively.The difference between total order and first order indices can indicate the interactive indices between the parameters of interest and other parameters.The different height of the bars revealed the various levels of sensitivity.The higher the bar, the more sensitive the parameter.For the uncorrelated and correlated cases, the overall contribution of SC to the four objective functions under the three flood conditions varied from 0.23 to 0.5, which was the highest and indicative of the most sensitive.Then, the partial effect of FC was the second largest with the range of 0.23-0.48.The other parameters demonstrated sensitivity indices within the range of 0.01-0.1.As Tang et al. [14], Zhan et al. [6], and Zhang et al. [21] mentioned, parameters with sensitivity indices larger than 0.1 can be regarded as highly sensitive; sensitivity indices between 0.01 and 0.1 can be classified into sensitive.Those parameters with little contribution to model outputs' variance (e.g., sensitivity indices less than 0.01) can be considered as insensitive.Therefore, both SC and FC were indicated to be highly sensitive with the largest first and total order indices for uncorrelated or correlated cases, whereas Ks, n, Mn, h, and b were regarded as sensitive parameters.The Sobol' results did not classify any parameters into insensitive, which revealed that the parameters fixed as constant by the Morris screen in Section 4.2 were reasonable.The Sobol' method strengthened the Morris screening results and quantified the sensitivity indices of each parameter with detailed values.In addition, our results were in accordance with Liao et al. [44], who investigated the Liuxihe model parameters based on the E-FAST algorithm using ENS as objective function by 1560 model runs.
Despite these similar outcomes obtained by uncorrelated and correlated parameters, it is useful to distinguish the differences between the two cases as well.For the uncorrelated case, the general order of average first order sensitivity indices of the seven parameters studied, from the highest to lowest, was 0.300, 0.292, 0.047, 0.039, 0.036, 0.035, and 0.029 for FC, SC, Mn, h, n, Ks, and b, respectively.For the correlated case, the order of first order sensitivity indices, from the highest to lowest, was 0.380, 0.377, 0.054, 0.045, 0.034, 0.032, 0.029 for FC, SC, Mn, h, n, b, and Ks, respectively.A large increase (larger than 0.08) of the sensitivity indice was observed from the uncorrelated case to the correlated, referred to FC and SC.A slight increment (about 0.006) was indicated from the uncorrelated to the correlated when it came to Mn and h.Finally, n, b, and Ks present a small decrease (about 0.003) for the correlated case.Similar results can also be obtained for the total sensitivity indices comparing the uncorrelated and correlated cases.These results were in agreement with the results given by Pan et al. [55], who revealed that the sensitivity of van Genuchten n to the percolation flux of interest improved when the van Genuchten n was correlated with van Genuchten α and the sensitivity indice of permeability declined when no interaction was considered.Furthermore, the improved sensitivity of van Genuchten n increased as the correlation coefficient went higher, which meant that parameters' correlations played a significant effect on the sensitivity analysis and should not be overlooked.In this study, a nearly linear regression between SC and FC was exhibited in Section 4.3.Correspondingly, the variation of correlated analysis compared with the uncorrelated was very large (even larger than the sensitivity indices of the five sensitive parameters) in the case of highly correlated parameters.For the poorly correlated parameters (Mn and h), the increments were indistinct and when it came to non-significantly correlated parameters (n, b, and Ks), the sensitivity indices eventually decreased.This may be caused by the sample sets, considering the correlation, diminishes the parameters' space [38].In other words, the correlated samples are more concentrated.Therefore, the fractional contributions of correlated parameters to model output variance are highlighted compared with that of uncorrelated parameters.In addition, highly correlated parameters cast a more profound influence on the highlighting than the poorly correlated.From this point of view, as indicated by the results from Morris screening and the GLUE method, it is very likely that the highly sensitive parameters with large modified means reflected high correlation coefficients.Hence, the sensitivity analysis of these parameters may be underestimated if the interaction is not considered.On the contrary, the sensitivity analysis of sensitive parameters with non-significant correlation may be overestimated when other highly correlated parameters were not considered.Another interesting result should also be noted: the differences between the total and first order sensitivity indices, indicative of interaction effect, shrank compared with the uncorrelated case.Since only pairwise correlations were taken into account, the shrunken yet retained gaps between the total and first order sensitivity indices can be considered to be caused by the multivariable interactions and should be investigated in the future.
The sensitivity results also varied among different flood magnitudes.In both uncorrelated and correlated cases, the sensitivity indices decreased as the magnitudes increased.For example, the total order sensitivity indice of SC decreased from 0.468 for small floods to 0.353 for medium ones and 0.317 for large ones, based on the ENS measure for the uncorrelated.In addition, the gap between highly sensitive and sensitive parameters decreased as magnitudes increased.This can be explained by the fact that most small floods occurred after a long drought, during which the soil moisture was kept at a very low level [56] and a hydrograph was vulnerable to the disturbance and attenuation of the underlying surface.Therefore, a slight fluctuation of model parameters would play a noteworthy role in determining the hydrograph of a watershed outlet.In this way, the sensitivity indices and their variance for small floods were more significant.On the contrary, large floods mainly happened after long precipitation, during which the soil was moist enough and the disturbance of the underlying surface cast a less important effect on outlet hydrographs.Therefore, the sensitivity indices and their variation for large floods were less significant.
Additionally, attention should also be paid to the variations of sensitivity indices related to different objective functions.Taking ENS as the objective function, the most sensitive parameters, SC and FC, were significant and five sensitive parameters were at about the same level.While taking WB as the objective function, the indices of n and Mn declined compared with other sensitive parameters.However, this situation was reversed when taking EP and ETP as the objection functions.These differences came from the diverse emphasized hydrograph process, magnitude, and timing from the four objective functions.ENS focused on the variation between simulated hydrograph and the average measured one.WB concentrated only on the gap between the simulated and measured volume.Moreover, EP and ETP merely paid attention to the efficiency of the simulated peak discharge and the time to peak.Therefore, flow concentration parameters such as n and Mn that determine the flow velocity would have a more momentous impact on the peak discharge presented by EP and ETP, while little impact on the flood volume presented by WB.On the contrary, runoff generation parameters such as Ks, h, and b would play a more important role in influencing WB than EP and ETP compared with n and Mn.Both runoff generation and flow concentration parameters affect the shape of a hydrograph; their competing consequences result in the same sensitivity indices level for ENS.These results implied that the first and total order sensitivity indices of model parameters were greatly affected by the selection of objective functions as well as the choice of flood magnitudes.Previous studies [6,7,14,23,24,51,53] that took ENS or RMSE as the objective function would not classify the runoff generation and flow concentration parameters and obtain deep insight into the model behavior.

Results of Second Order and Interactive Sensitivity Analysis
In this study, further insight of the second order indices between the seven parameters without consideration of correlation relationship was also illustrated.Table 5 showed that the sum of the first order sensitivity indices of these seven parameters was less than 1.Zhan et al. [6] explained that the sum of first order indices less than 1 was caused by the interactive effects among parameters, which could not be considered in the single parameter sensitivity analysis.In addition, this outcome indicated that the interactive impacts among parameters cannot be neglected.This indication was reinforced by the gaps between the first and total order sensitivity indices (see Figure 4 and Table 5).On the other hand, the sum of total order sensitivity indices was slightly larger than 1.This phenomenon was explained by Tang et al. [14]: the truncation and Monte Carlo integration implemented in the Sobol' method can result in slight numerical errors.However, Tang et al. [14] acknowledged that this unimportant error would not affect the sensitivity ranking of parameters.Baroni and Tarantola [57] also investigated the robustness of sensitivity indices under the change of sample size.The total order sensitivity indices proved to be more stable than the first order ones.Therefore, the total order sensitivity indices of the most and least sensitive parameters (i.e., FC and Ks) were iteratively explored by increasing the sample size by 200 at a time from 200 to 2000.Results in Figure 5 showed that both total order sensitivity indices of Ks (left) and FC (right) fluctuated within a wide range when the sample size was less than 1200.Furthermore, the variation of highly sensitive parameters was more significant than that of sensitive parameters.The total sensitivity indices of Ks with EP and ETP as objective functions were seen to converge at a very low computation cost.Even though the total sensitivity indices changed more significantly when taking ENS and WB as objective functions, the differences between the maximum and minimum were 0.006 and 0.015 for Ks, respectively.However, from the perspective of FC, only when taking EP as objective function did the total sensitivity indice vary slightly, and the differences between the maximum and minimum values were 0.047, 0.018, and 0.035 for ENS, WB, and ETP, respectively, which were much larger than those of Ks.These results were in consistent with the research of Baroni and Tarantola [57].However, as long as the sample size increased to a large value, e.g., 1600, the total sensitivity analysis began to converge and generate steady results for both highly and generally sensitive parameters considering various objective functions.These results underlined that the sample size of 2000 set in the Sobol' analysis in Section 4.4 was reasonable.
As a matter of fact, it is obvious that the sum of the first and second order sensitivity indices accounted for nearly all the total order (see Table 5), which indicated that the interactive effects among parameters can be mainly explained by pairwise interactions.In other words, the second order sensitivity indices of the Liuxihe Model can denote most parameter interactions.More precisely, as given in Table 5, the average second order sensitivity indices were 0.32, 0.25, 0.30, and 0.27, respectively, based on ENS, WB, EP, and ETP measures, with corresponding first order indices of 0.75, 0.77, 0.73, and 0.76, respectively.In order to better analyze the two-parameter interactions, the pairwise variation-induced second order sensitivity indices (see Figure 6) provided a more exhaustive understanding of the second order indices listed in Table 5.These figures displayed the lower triangle matrix of pairwise parameters interactions where dark gray shadow denoted higher interactions and light shadow denoted lower interactions.Remarkably, the pairwise interactions became more significant as the parameters grew more sensitive.For instance, the interactive sensitivity indice between SC and FC was as high as 0.09 based on the WB measure.Furthermore, the pairwise interactions between SC and the remaining sensitive parameters, and between FC and the remaining sensitive parameters, were also relatively higher for all objective functions under the three flood magnitudes.However, the two-way interactions between sensitive parameters were relatively lower (e.g., the second order sensitivity indice between Mn and n was less than 0.01 even though both parameters described the process for flow routing).
Our results were in good agreement with the results provided by Tang et al. [14], Zhan et al. [6], and Zhang et al. [21], even though the hydrological models were different.Tang et al. [14] pointed out that the lower zone storage parameter and percolation factor parameter in the model combined by SAC-SMA and SNOW-17 presented the highest number of interactions, while the two parameters both also showed the highest influence on model output in the first order sensitivity analysis.Zhan et al. [6] concluded that the second order indice between soil moisture parameter and thickness of upper soil layer in the model named DTVGM was 0.121 high for the objective function WB and their corresponding first order indices were 0.288 and 0.139, respectively, ranking inside the top three.Based on the Xinanjiang model, Zhang et al. [21] concluded that the pairwise interaction between areal mean free water storage capacity parameter (Sm) and recession constant for channel routing (Cs) had a profound effect on model performance, while Sm and Cs were the two most sensitive parameters in the total order sensitivity evaluation in the meantime.The application of Sobol' method into various models' sensitivity analysis indicated that the parameters that showed high sensitivity indices in the first and total order analysis were prone to presenting a correspondingly significant pairwise interaction with other parameters in the second order analysis.This was an expected result, as the changes in the sensitive parameter, coupled or not coupled with other parameters' changes, would definitely contribute to the great variance of model  Additionally, pairwise interaction investigation among different flood magnitudes revealed that the second order parameter sensitivities were largely influenced by the selected flood magnitudes.Small floods presented the highest second order sensitivity while large floods presented the lowest sensitivity.The reason was similar to the differences caused in the first and total order sensitivity analysis.A dry and low soil moisture content condition before a small flood kept the watershed susceptible to the complex perturbations (e.g., changing pairwise parameters) of the underlying surface, while high soil moisture content before a large flood made the watershed more tolerable to the multipart disturbances.This result supported the above interaction analysis.Coupling parameters with higher sensitivity would generate higher interactive sensitivity indices (e.g., for small magnitude) than coupling lower sensitive parameters (e.g., for large magnitude).

Conclusions
In this study, the sensitivity and interaction of the distributed hydrological model parameters was investigated to help present and understand a model's behavior.The Morris screening method was employed to detect the insensitive parameters and the GLUE method was applied to identify the parameters correlations.Subsequently, the extended and generalized Sobol' methods were implemented to compare the sensitivity indices of the Liuxihe model parameters in the case of uncorrelated and correlated parameters based on four objective functions: (1) Nash-Sutcliffe efficiency (ENS); (2) water balance coefficient (WB); (3) peak discharge efficiency (EP); and (4) time to peak efficiency (ETP).Hourly rainfall-runoff data in the Nanhua Creek catchment, Pearl River, China, were used.Contrasting outcomes for the sensitivity and interaction analysis were illustrated with three different flood magnitudes, i.e., small, medium, and large.The main conclusions drawn are as follows: 1.The Morris screening method is helpful for computationally demanding distributed models with dozens of parameters.For the 12 parameters in the Liuxihe model, S, v, Wp, Ep, and Kg are detected as insensitive parameters no matter which objective functions or flood magnitudes are selected.The other seven parameters were retained for further global sensitivity analysis.2. The GLUE method provides a convenient approach to identify the correlation relationship between parameters.SC and FC are observed to be nearly linear regressed, with a correlation coefficient of about 0.96 at the significance level of 0.01.Other parameters demonstrated a certain correlated relationship as well.3.For both uncorrelated and correlated parameters, the classification of parameters as first or total order sensitivity indices was affected neither by the choice of objective functions nor by the selection of different flood magnitudes.FC and SC are highly sensitive parameters and the remaining five parameters are sensitive.Meanwhile, the ranking of parameters sensitivity indices was greatly determined by the various objective functions specified and particular flood magnitudes stated.Due to the distinctive antecedent watershed conditions inducing different tolerance to watershed disturbances, the sensitivity indices for small floods were relatively higher than those for large floods.With regard to specific objective functions, flow concentration parameters such as n and Mn had a greater effect on EP and ETP than WB, whereas runoff generation parameters such as Ks, h, and b display totally reversed consequences.Hydrological models behave and respond uniquely to various evaluation functions and application conditions.4. Highly correlated parameters present relatively larger sensitivity indices in the correlated case compared with the uncorrelated, whereas poorly correlated parameters displayed inverse results.Considering that highly sensitive parameters are prone to correlate with each other and commonly sensitive parameters are inclined to be distinguished from other parameters, the sensitivity indices of highly sensitive parameters might be underestimated and that of sensitive parameters might be overestimated without accounting for correlations.Further in-depth second order sensitivity analysis indicated that pairwise interactions, which cannot be neglected, accounted for most interactions.
In order to understand how model parameters contribute to model uncertainty, the preliminary research on Sobol' sensitivity analysis was carried out accordingly.Further investigation of the GLUE uncertainty analysis can be conducted instead of simple calibration/validation analysis since the GLUE method also provides an estimation of uncertainty in the performance of the hydrological models.

Figure 1 .
Figure 1.The location, elevation, river networks, and distribution of rain gauge and hydrologic stations of the Nanhua Creek catchment.

Figure 3 .
Figure 3. Scatter plots of the Morris screening results based on different objective functions measures under small, medium, and large flood magnitudes.
lower modified mean value *

Figure 4 .
Figure 4. Comparisons between the first and total order sensitivity indices of seven parameters in the Liuxihe model with and without accounting for correlations based on different objective functions measures under three different flood magnitudes (UCS1 and UCST designated the uncorrelated first and total order sensitivity indices; CS1 and CST designated the correlated first and total order sensitivity indices).

Figure 5 .
Figure 5. Variations of the total order sensitivity indices along with increasing sample size.(a) For the least sensitive parameter Ks; (b) for the most sensitive parameter FC.

Figure 6 .
Figure 6.Second order sensitivity indices between seven parameters based on different objective functions measures under three different flood magnitudes.Dark grey shadow denoted high parameter interactions.Light gray shadow denoted low parameter interactions.

Table 1 .
Summary of the Liuxihe model parameters.

Table 2 .
Characteristics of the selected various flood conditions.W0 denoted the soil moisture content before the flood event began.

Table 3 .
Goodness of fit statistics between Liuxihe predicted and observed hydrographs in the Nanhua Creak catchment.

Table 4 .
Correlation analysis between effective parameters.Dark grey shadow denoted highly correlated parameters.Light gray shadow denoted relatively correlated parameters.White cells denoted poorly correlated parameters.Values labeled "**" meant its significance level was 0.01.Values labeled "*" had a significance level of 0.05.

Table 5 .
Summations of the first, second, and total order sensitivity indices based on four different objective functions.