Investigating the Dynamic Influence of Hydrological Model Parameters on Runoff Simulation Using Sequential Uncertainty Fitting-2-Based Multilevel-Factorial-Analysis Method

Abstract

Hydrological model parameters are generally considered to be simplified representations that characterize hydrologic processes. Therefore, their influence on runoff simulations varies with climate and catchment conditions. To investigate the influence, a three-step framework is proposed, i.e., a Latin hypercube sampling (LHS-OAT) method multivariate regression model is used to conduct parametric sensitivity analysis; then, the multilevel-factorial-analysis method is used to quantitatively evaluate the individual and interactive effects of parameters on the hydrologic model output. Finally, analysis of the reasons for dynamic parameter changes is performed. Results suggest that the difference in parameter sensitivity for different periods is significant. The soil bulk density (SOL_BD) is significant at all times, and the parameter Soil Convention Service (SCS) runoff curve number (CN2) is the strongest during the flood period, and the other parameters are weaker in different periods. The interaction effects of CN2 and SOL_BD, as well as effective hydraulic channel conditions (CH_K2) and SOL_BD, are obvious, indicating that soil bulk density can impact the amount of loss generated by surface runoff and river recharge to groundwater. These findings help produce the best parameter inputs and improve the applicability of the model.

1. Introduction

Hydrological models play a crucial role in simulating the hydrological process of river basins. These hydrological models are generally composed of several parameters, whose values cannot be directly determined by field observations but which can be calibrated through input/output records, which inevitably contain the basin error response [1]. Among many hydrological models, the Soil and Water Assessment Tool (SWAT) model has been widely used in many countries for its ability to completely reflect the influence of spatiotemporal heterogeneity, such as topography, soil, and land use, on the water cycle of the river basin [2,3,4]. However, the uncertainty of the SWAT model parameter is difficult to evaluate, because the model parameters are numerous and difficult to obtain. It also brings difficulties to decision-makers when the hydrological process is accurately described, as well as the regional relationship between the model parameters and the watershed characteristics. Thus, more effort is required to quantify the uncertainty in the hydrological simulation.

To date, a variety of optimization algorithms have been developed for calibration and uncertainty analysis, and good results have been achieved [5,6,7,8]. Kouchi et al. [9] use three different optimization algorithms (sequential uncertainty fitting 2 (SUFI-2), particle swarm optimization (PSO), and generalized likelihood uncertainty estimation (GLUE)), as well as eight evaluation indexes in a SWAT model for emphasizing that the combination of target functions of each optimization algorithm may lead to different optimal parameter sets, and have the same performance at the same time. Trudel et al. [10] use seven different objective functions to study hydrological model calibration and the structural uncertainty in low-flow simulations under climate change conditions. Muleta et al. [11] use comparative analysis of parameter sensitivity in the high- and low-flow period, and found that there was a significant difference in the parameter sensitivity between the high- and the low-flow periods, and the same parameters were different in different periods. However, the algorithm limitations, such as with GLUE and SUFI-2, lie in the quantitative analysis of the impact of their parameters on system performance. In fact, model parameters describing different hydrological processes have different individual effects on the model output. Factorial analysis can help study the individual and interaction effects of the parameters [12]. The factorial-analysis-of-variance method is used to diagnose the curve relationship between the parameters and the response [13,14,15]. Nevertheless, no previous study has been conducted to investigate the dynamic influence of hydrological model parameters on runoff simulation using the SUFI-2-based multilevel-factorial-analysis method.

As a case study, this study regarded the source region of the Yellow River, which is known as the “water tower” of the Yellow River basin and contributes 35% of total annual runoff from about 16.2% of the basin area [16]. More importantly, the Yellow River plays a key role in the water supply for 107 million people and for about 13% of the agricultural production of the country’s total cultivated area.

The objective of this study was to develop a SUFI-2-based multilevel-factorial-analysis method to address the dynamic influence of hydrological model parameters on runoff simulation. The main steps of the study included (i) the China Meteorological Assimilation Driving Datasets (CMADS), which were used to drive the SWAT model; (ii) a multivariate regression model to conduct parametric sensitivity analysis, which was based on the results of the Latin hypercube sampling (LHS) method; (iii) the confidence interval of each sensitive parameter, found using the SUFI-2 algorithm; and (iv) using the SUFI-2-based multilevel-factorial-analysis method, we quantitatively evaluated the individual and interactive effects of parameters on the hydrologic model output. The results of the study are helpful for improving the simulation and prediction ability of the hydrologic model for water resources.

2. Study Area

The source region of the Yellow River (Figure 1) is located in the northeast Qinghai–Tibet Plateau between longitudes 95°50′ E and 103°30′ E and latitudes 30°30′ E and 35°0′ E, covering 12.19 × 104 km² and occupying 16.2% of the entire Yellow River basin (75.24 × 104 km²) [17,18]. The average temperature is about 5 °C, and the temperature here varies greatly between day and night. The average annual precipitation varies between 320 and 750 mm. Precipitation in June to September accounts for 80% of the total year. Alpine vegetation and alpine meadows are the major vegetation types, accounting for the total area of 70% in 2010. The major soil type in the watershed is loam, and most of the soil has poor water retention and low fertility.

3. Methodology

This paper uses the CMADS data set to drive the SWAT model. The constructed SWAT model employs parameters to depict the characteristics of the hydrological process. To identify the parameter sensitivity to the runoff simulation of the SWAT model, we applied the LHS-based multivariate regression model. Then, the SUFI-2 algorithm was used to explore the confidence interval of these identified sensitive parameters. We employed the multilevel-factorial-analysis method to quantitatively evaluate the individual and interactive effects of the parameters on the hydrological model output, i.e., the simulated discharge. Finally, the physical mechanism of the parameters involved in the hydrological process was traced.

3.1. Construction of the Soil and Water Assessment Tool Model

The inputs required for the SWAT model include a DEM (digital elevation model), land use, and soil data sets. The DEM is the Shuttle Radar Topographic Mission (SRTM) (90) DEM, which comes from the geospatial data cloud (http://www.gscloud.cn). The soil data were obtained from the China Soil Data Set (v1.0), based on the World Soil Database (HSDW). The land use data (LCC2010) of the study region were derived from the dry area scientific data center in the cold region. The China Meteorological Assimilation Driving Datasets for the SWAT model Version 1.1 (CMADS V1.1, http://www.cmads.org), which was developed by Dr. Xianyong Meng from the China Agricultural University (CAU), has received worldwide attention [19]. This data set is widely used by countries throughout the world [20,21,22,23,24,25,26]. Meng et al. [27] chose the Manas River Basin (MRB) in China as a research area, in order to verify the adaptability of the China Meteorological Assimilation Driving Datasets for the Soil and Water Assessment Tool model (CMADS); the results showed that the SWAT model could reproduce the runoff process of two stations (Kenswat and Hongshanzui) in the research area well using data from CMADS. Zhang et al. [28] used CMADS to drive the SWAT model for a runoff simulation in the Hunhe River Basin, and the results showed that both Nash-Sutcliffe efficiency coefficient (NSE) values and R² were found to be greater than 0.74 in calibration, and are greater than 0.58 in validation. Related results indicate that CMADS performs particularly well in runoff simulation.

3.2. Parameter Sensitivity Analysis

There are many parameters in the SWAT model. Different parameters have different sensitivities to the model simulation [29,30,31]. The model can eliminate the parameters that have little influence on the model results from sensitivity analysis, and then reduce the influence of the uncertain transmissions. This study used NSEs as an objective function, defined as

$$\mathrm{NSE}=1-\frac{{\displaystyle \sum _{i=1}^n{(Q_{obs,i}-Q_{sim,i})}^2}}{{\displaystyle \sum _{i=1}^n{(Q_{obs,i}-{\overline{Q}}_{obs})}^2}}$$

(1)

where Q_sim,i is the ith simulated discharge, Q_obs,i is the ith observed discharge, ${\overline{Q}}_{obs}$ is the mean of the observed data, and n is the simulation period.

According to the physical meaning of each parameter, we calibrated the SWAT model using the same initial parameter ranges, according to the calibration protocol presented by Abbaspour [32]. The sensitivity analysis of the parameter values is generated by LHS sampling, and the value of the target function by the multiple regression model. The calculation expression of the parameter sensitivity is written as

$$g=\alpha +{\displaystyle \sum _{i=1}^m{\beta }_i{b}_i}$$

(2)

where g is the objective function value, α is the regression constant, β is the coefficient of parameters, b_i is the parameter value, and the m is the number of parameters. The t-test method is used to determine the sensitivity of each parameter.

The following method is used to deduce the confidence interval of each parameter. First, the sensitivity matrix calculation formula for the objective function is

$$J_{ij}=\mathsf{\Delta }{g}_i/\mathsf{\Delta }{b}_{j}i=1,\cdots ,C_2^n,j=1,\cdots ,m$$

(3)

where $C_2^n$ indicates the number of rows in the sensitivity matrix, j represents the number of parameters, i means the group number, Δb_j is the parameter of the j rate, and Δg_i represents the parameter sensitivity. The Hessian matrix calculation formula for the objective function is

$$H=J^{T}J$$

(4)

where H indicates the haessen matrix, J represents the matrix of the number of parameter columns.

According to Kramer’s theorem, the covariance matrix C for estimating the lower limit of the parameter is calculated as

$$C=s_g^2{\left(J^{T}J\right)}^{-1}$$

(5)

where $s_g^2$ is the deviation of the result of the n simulation of the target function.

The standard variance of parameter b_i and its 95% confidence interval (CI) are calculated by the diagonal elements in C, as follows:

$$s_j=\sqrt{C_{jj}}$$

(6)

$$b_{j,lower}=b_j^{\ast }-t_{\nu ,0.025}\cdot S_j$$

(7)

$$b_{j,upper}=b_j^{\ast }+t_{\nu ,0.025}\cdot S_j$$

(8)

where $b_j^{\ast }$ is the optimal solution of the parameter b, v is the degree of freedom (n − m), b_j,lower is the lower limit of the confidence interval, and b_j,upper is the upper limit of confidence interval.

3.3. Parameter Uncertainty Evaluation Index

In order to judge the influence of parameter uncertainty on runoff simulation under different levels, this paper uses the two indexes of variation rate (VR) and the relative length of the confidence interval (RL) as the model uncertainty evaluation index. These can be expressed as

$$VR=\frac{Abs\left|Q_{sim}-Q_{obs}\right|}{Q_{obs}}$$

(9)

$$RL=\frac{Q_{upper}-Q_{lower}}{Q_{obs}}$$

(10)

where VR is the variation rate, RL is the relative length of confidence interval, Q_upper is the upper limit of the runoff simulation under the 95% confidence interval, and Q_lower is the lower limit of the runoff simulation under the 95% confidence interval.

3.4. Multilevel Factorial Analysis

The factorial analysis is a multivariable reasoning method. It performs excellently in testing the effects of individual variables and their interactions on the dependent variable [33,34,35]. Factors A and B have m and n levels, respectively. Thus, a full-factor factorial design contains all possible factor combinations. The factor model for this factorial experiment can be expressed as

$$Y_{ijk}=\mu +{\tau }_i+{\beta }_j+{\left(\tau \beta \right)}_{ij}+{\epsilon }_{ijk}\{\begin{array}{c}i=1,2\dots a\\ j=1,2\dots b\\ k=1,2\dots n\end{array}$$

(11)

where μ is the total average effect, ε_ijk is a random error effect, τ_i is the effect of factor A at the ith level, β_j is the effect of factor B at the jth level, and ${\left(\tau \beta \right)}_{ij}$ is the interaction effect when A is at the ith level and B is at the jth level. There is a total of abn experiments, where n is the number of repeated experiments. In order to test the influence of the parameter main effect and the interaction effect on the runoff simulation, the F-statistic can be used as follows:

$$F_A=\frac{M{S}_A}{M{S}_E}=\frac{S{S}_A/a-1}{S{S}_E/ab(n-1)}$$

(12)

$$F_B=\frac{M{S}_B}{M{S}_E}=\frac{S{S}_B/b-1}{S{S}_E/ab(n-1)}$$

(13)

$$F_{AB}=\frac{M{S}_{AB}}{M{S}_E}=\frac{S{S}_{AB}/(a-1)(b-1)}{S{S}_E/ab(n-1)}$$

(14)

where MS_A, MS_B, MS_AB, and MS_E are the mean squares for factors A, B, their interaction with each other, and the error component, respectively. The SS_A, SS_B, SS_AB, and SS_E are the sum of squares for factors A and B, their interaction, and the error component, respectively. Each mean square deviation is the squared sum of the corresponding effects, divided by its degree of freedom. SS_T is the sum of the total effect square. This can be calculated by

$$S{S}_A=\frac{1}{bn}{\displaystyle \sum _{i=1}^a{y}_{i..}^2-\frac{y_{...}^2}{abn}}$$

(15)

$$S{S}_B=\frac{1}{an}{\displaystyle \sum _{j=1}^b{y}_{.j.}^2-\frac{y_{...}^2}{abn}}$$

(16)

$$S{S}_{AB}=\frac{1}{n}{\displaystyle \sum _{i=1}^a{\displaystyle \sum _{j=1}^b{y}_{ij.}^2}-\frac{y_{...}^2}{abn}}-S{S}_A-S{S}_B$$

(17)

$$S{S}_T={\displaystyle \sum _{i=1}^a{\displaystyle \sum _{j=1}^b{\displaystyle \sum _{k=1}^n{y}_{ijk}^2}}}-\frac{y_{...}^2}{abn}$$

(18)

$$S{S}_E=S{S}_T-S{S}_{AB}-S{S}_A-S{S}_B$$

(19)

where y_i.., y_.i., and y_ij. represent the ith level of the factor A, the jth level of the factor B, and the ijth interaction between factors A and B, respectively.

In particular, there are k factors in the 3^k factorial design. Each factor has three levels, so there is a total of 3k factor level combinations and 3k degrees of freedom. There is a total of k main effects, each of which is two degrees of freedom. There is an interaction effect of k factors, and the degree of freedom is 2k − 1. If n repeated tests are performed for each factor level combination, the total degree of freedom is n3k – 1, and the degree of freedom of the error is 3k(n − 1). The sum of squares for the main effects and interaction effects is usually obtained by the factorial analysis method.

4. Results and Discussion

4.1. Parameter Sensitivity Analysis, Calibration, and Verification of Model

Parameter sensitivity analysis is an indispensable part of the evaluation model, and it is helpful in developing a deep understanding of the model characteristics. Hence, this paper takes NSE as the objective function, and uses LHS-based multiple regression models to analyze the parameter sensitivity. The t-test method is used to determine the parameter sensitivity—the higher the absolute value of t, the stronger the sensitivity of the parameters. The related CI of the parameters was determined by the SUFI-2 algorithm. Parameters’ descriptions and their CIs are listed in

Table 1. Final value range of sensitivity parameters.

Parameter	Description	CI		Calibrated Value	t Value
		Min	Max
GWQMN	Threshold depth of water in the shallow required for return flow to occur (mm)	700.61	780.61	742.93	1.17
HRU_SLP	Average slope steepness (m/m)	0.24	0.26	0.25	1.5
ALPHA_BF	Baseflow regression constant (days)	0.00	0.20	0.02	−1.52
SFTMP	Snow temperature (°C)	3.43	3.63	3.45	−2
SOL_AWC	Effective water capacity of soil layer (mmH2O/mm soil)	0.15	0.17	0.16	2.19
SOL_K	Soil hydraulic conductivity (mm·hr−1)	−0.32	−0.29	−0.32	−4.04
CH_K2	Effective hydraulic conductivity of channel (mm/h)	110.97	125.97	121.19	4.59
ESCO	Soil evaporation compensation coefficient (mm/h)	0.73	0.75	0.73	−7.36
CN2	Initial SCS runoff curve number to moisture conditions II	0.26	0.28	0.26	−15.28
SOL_BD	Soil bulk density (g/cm3)	0.45	0.47	0.46	−58.26

SWAT parameter ranges (and allowable percentage changes): 0 ≤ GWQMN ≤ 5000; 0 ≤ HRU_SLP ≤ 0.6; 0 ≤ ALPHA_BF ≤ 1; −5 ≤ SFTMP ≤ 5; 0 ≤ SOL_AWC ≤ 1 (−50% ≤ SOL_AWC ≤ 50%); 0 ≤ SOL_K ≤ 2000 (−80% ≤ SOL_K ≤ 80%); 0.01 ≤ CH_K2 ≤ 150; 0 ≤ ESCO ≤ 1; 35 ≤ CN2 ≤ 98 (−50% ≤ CN2 ≤ 50%); 0.9 ≤ SOL_BD ≤ 2.5 (−50% ≤ SOL_BD ≤ 50%); reference [36].

. For the most part, specific values are given for the parameters. However, some parameters, such as the runoff curve number, can take on a variety of values. During calibrations, these parameters may be changed by increasing or reducing them by a certain percentage, until the calibration objective function is met. Thus, values of these parameters have been reported as a percentage change from a specified value. The Latin Hypercube sampling (LH-OAT) method was performed for 10 parameters, in order to screen out the most sensitive parameters of the model, and the results are illustrated in Figure 2.

As we can see from Figure 2 and Table 1, the most sensitive parameters for runoff simulation are followed by SOL_BD, CN2, ESCO, CH_K2, SOL_K, SOL_AWC, SFTMP, ALPHA_BF, HRU_SLP, and GWQMN. In order to design the experimental scheme, the four most sensitive parameters are selected: The soil bulk density (SOL_BD) was the most sensitive parameter, followed by the Soil Convention Service (SCS) runoff curve number for moisture conditions II (CN2), the soil evaporation compensation coefficient (ESCO), and effective hydraulic channel conductivity (CH_K2). A detailed interpretation of these parameters will be provided in the following section.

On the basis of CMADS (2008–2015), this paper selected years 2008–2009, 2010–2013, and 2014–2015 as warm-up, calibration, and validation periods, respectively. The NSE, R², and the absolute value of relative error (|Re|) indicators were used to evaluate the model for calibration and validation periods. The results show that the NSE values were 0.73 and 0.81, respectively, R² was 0.82 and 0.87, respectively, and |Re| was less than 10% for both periods. The results indicate a good performance of SWAT in describing the runoff simulation, based on the CMADS data in the source region of the Yellow River. Che [37] investigated the source area of the Yellow River, and used the SWAT model to simulate the daily runoff. The results of the study showed that both NSE values and R² were less than 0.74 in calibration (validation), which means that CMADS data is superior to other data in watershed runoff simulation.

Figure 3 shows the simulated daily runoff from 2010–2015 (i.e., the calibration and validation periods) for the Tangnaihai hydrological station. Specifically, this figure indicates that the 95% CI width of the daily runoff simulation varies with the flow amount. To further examine the above results, three different frequencies are set, in order to identify the influence of parameter uncertainty on flow of different levels (i.e., more than 75%, between 25% and 75%, and less than 75%).

Figure 4 shows that when the flow at the high level, the variation rate and the relative length of confidence interval are small; on the contrary, as the flow reaches the lower level, the variation rate and the relative length of confidence interval are larger. Given this, we can preliminarily infer that the physical mechanism of runoff generation dynamically changes during different periods.

4.2. Multilevel Factorial Analysis and Dynamic Changes in Parameter Sensitivity

As mentioned, the physical mechanism of runoff generation varies at different periods. In the SWAT model, parameters reflect the characteristics of the hydrological processes. Thus, based on the above identified sensitive parameters, we used the analysis of variance method (ANOVA) to quantify the response of individual and interactive parameters in the runoff simulation.

Four parameters, CN2, ESCO, CH_K2, and SOL_BD, were selected as factorial experimental factors (denoted by A~D, respectively). Subsequently, we designed the 3⁴ factorial design scheme shown in

Table 2. The 3⁴ factorial design scheme.

Parameter	Level
	Low	Medium	High
CN2	0.257	0.267	0.277
ESCO	0.726	0.736	0.746
CH_K2	110.97	120.96	130.97
SOL_BD	0.431	0.452	0.471

.

Figure 5 shows the parameter contributions for the hydrological responses. The results show that parameter D is significant at different times, especially in the non-flood period (November to March), when its contribution reached 0.98. The contribution of parameter A in the pre-flood period (April to May) increased, but the effect was relatively small, and D was still dominant. The contributions of A, B, and C during the flood period (June to September) gradually increased, whereas D gradually weakened. The linear individual effects of A, B, and D, as well as the AC and CD interaction effects, are thus significant for modelling runoff simulation. During the post-flood period (October), the contributions of A, B, and C weaken, while parameter D increases.

4.3. The Individual and Interactive Effects of Parameters on the Hydrologic Model Output in Different Periods

From the previous section, we find that the parameter sensitivity in different years shows similar characteristics in the non-flood, pre-flood, flood, and post-flood periods. To quantitatively evaluate the dynamic effect of the parameters on the runoff simulation, we selected 2012 as a typical year, and three consecutive days of non-flood, pre-flood, flood, and post-flood periods in a typical year. The objective is to explore the dynamic parameter individual and interaction effects on the runoff simulation. For simplicity, this article only presents the first-day results of the four periods in a typical year, and the other results are attached to Supplementary Materials.

4.3.1. The Statistically Significant Individual and Interaction Effects on Runoff Simulation in Non-Flood Period

Table 3. Results of ANOVA for the runoff simulation in the non-flood period.

Model Term	Sum of Squares	F Value	p-Value	Significance
A	28.47	16,998.37	0.00	**
B	0.76	456.16	0.00	*
C	0.00	0.00	1.00
D	12,890.57	7,695,860.81	0.00	***
AB	0.12	34.59	0.00	*
AC	0.00	0.00	1.00
AD	2.73	814.23	0.00	*
BC	0.00	0.00	1.00
BD	0.12	34.83	0.00	*
CD	0.00	0.00	1.00
Error	0.04
Total	12,922.80

presents the results of variance analysis for the simulation of daily runoff in the non-flood period. In this paper, a p-value of less than 0.05 indicates that parameters A, B, and D have significant effects on the model output in the non-flood period. In particular, the influence of D is more significant. Moreover, the interaction effects of AB, AD, and BD have statistical significance.

For the non-flood period, the main effects of these parameters on the runoff simulation are shown in Figure 6. From this, we see that changes in the levels of A, B, and C have a weak effect on the runoff simulation, while changes in the level of parameter D have significant negative effects on the runoff simulation.

In essence, during the non-flood period, D is an important parameter affecting the model simulation, which is closely related to the high elevation and cold temperatures of the source region of the Yellow River, characterized by low rainfall and low temperature. Moreover, lower flow values relate to lower compensation for glaciers and snowmelt in soil water, resulting in decreased streamflow [38,39].

Figure 7 presents the parameter interaction effects in the non-flood period of the runoff simulation. Results show that the AD, BD, and CD interaction effects are significant in this period, while the others have less influence on the runoff simulation. The main reason for this is that the soil bulk density (i.e., parameter D) is the parameter that has the greatest impact on the runoff simulation.

Note that for the AD interaction plot in Figure 7 (bottom left), the red, blue, and black lines represent parameter D at low, medium, and high levels, respectively. This plot discloses that the changes differ across the three levels of parameter D, depending on the level of parameter A. When D is at the high level, A has obvious negative effects. However, when D is at the middle level, the negative effect of A is weakened. This result further indicates that A has a significant impact on the runoff simulation.

4.3.2. The Statistically Significant Individual and Interaction Effects on the Runoff Simulation in the Pre-Flood Period

The parameters and their interactions have statistically significant effects on the runoff simulation (

Table 4. Results of ANOVA for the pre-flood period runoff simulation.

Model Term	Sum of Squares	F Value	p-Value	Significance
A	145.81	591,251.63	0.00	***
B	2.36	9586.61	0.00	**
C	0.00	5.42	0.01	*
D	420.67	1,705,827.37	0.00	***
AB	0.00	7.40	0.00	*
AC	0.00	0.94	0.45
AD	0.03	66.12	0.00	*
BC	0.00	1.39	0.25
BD	0.00	1.96	0.12
CD	0.00	0.46	0.76
Error	0.01
Total	568.88

). The linear individual effects of A, B, and D on the runoff simulation are significant in the pre-flood period, while other parameter effects (individual and interaction) have little impact on the response.

Figure 8 presents the individual effect plot for the four parameters in the pre-flood period. Results show that A has a positive effect on the runoff simulation, while D has a negative effect. Nevertheless, parameters B and C have a slight effect on the flow.

The flow value reaches 112.12 m³/s when A is at the high level. In contrast, the simulated flow value is only 96.20 m³/s. Furthermore, compared with other levels of runoff simulation, the three levels of D lead to smaller runoff simulation values. These results are attributed to a higher compensation for snow melt, and thus increased flow [40]. However, the slope of D is irregular, implying an obvious nonlinear effect on the runoff simulation, owing to the complex topography and soil characteristics of the basin.

Figure 9 depicts the parameter interaction effects in the pre-flood period. If the slope of one curve differs from another, there are interactive effects between the parameters. The parallel lines in the interaction plots of parameters A, B, C, and D indicate tiny parameter interactions. Specifically, the interaction of A and B, as well as B and C, are not apparent in Figure 9.

In addition, the results indicate that parameter A has the largest nonlinear effect on the runoff simulation, while the other nonlinear effects are weak. The identical results of these plots can be used to confirm the findings in the above numerical value model.

4.3.3. The Statistically Significant Individual and Interaction Effects on the Runoff Simulation in the Flood Period

The results of ANOVA for the runoff simulation during the flood period are shown in

Table 5. Results of ANOVA for the flood period runoff simulation.

Model Term	Sum of Squares	F Value	p-Value	Significance
A	31,490.77	3587.56	0.00	***
B	3800.32	432.95	0.00	**
C	745.95	84.98	0.00	*
D	15,926.84	1814.45	0.00	***
AB	1.53	0.09	0.99
AC	287.46	16.37	0.00	*
AD	177.23	10.10	0.00	*
BC	37.23	2.12	0.09
BD	26.79	1.53	0.21
CD	263.16	14.99	0.00	*
Error	210.67
Total	52,967.95

. Results indicate that parameters A, B, C, and D have significant linear main effects on the runoff simulation. Among them, the significance of A is the strongest, followed by D, B, and C. Meanwhile, parameters AC, AD, and CD interaction effects are of important significance.

The influence of the main parameter effects on the runoff simulation during the flood period is shown in Figure 10. It shows the influence of the parameters on the runoff simulation and compares the relative size of the influence. Parameter A has the largest positive individual effect when A is at low, middle, and high levels, when the simulated flow value is 2375.24, 2398.45, and 2420.22 m³/s, respectively. Parameter A is a comprehensive reaction of the underlying surface characteristics, which directly determines the size of the flow. However, parameter D has the greatest nonlinear effect.

Therefore, the smaller the soil bulk density, the weaker the ability to resist precipitation and thus increased flow. However, with the increase of D, the ability to resist precipitation is enhanced, which aids infiltration; however, the flow is then increased, due to a sharp decrease in soil porosity [41,42].

Figure 11 plots three level interaction effects of the flood period for the four parameters. Results disclose that AC, AD, and CD have noticeable interaction effects, while the others have almost no interaction effect.

Take the plot on the bottom left, with the interaction of parameters A and D as an example. The plot depicts the changes in D at its low, middle, and high levels, depending on the level of A, as well as the interaction between the low level and the middle level curve of D. This reveals the interaction of A and D has a significant influence on the runoff simulation. Therefore, the parameter interaction must be emphasized at calibration.

4.3.4. The Statistically Significant Individual and Interaction Effects on the Runoff Simulation in the Post-Flood Period

Table 6. Results of ANOVA for the post-flood period runoff simulation.

Model Term	Sum of Squares	F Value	p-Value	Significance
A	1267.75	258,334.54	0.00	***
B	662.22	134,942.54	0.00	**
C	0.03	5.48	0.01	*
D	66,212.33	13,492,324.58	0.00	***
AB	0.05	4.88	0.00	*
AC	0.00	0.35	0.84
AD	27.18	2769.41	0.00	*
BC	0.00	0.50	0.73
BD	4.98	507.67	0.00	*
CD	0.00	0.50	0.73
Error	0.12
Total	68,174.67

shows the ANOVA results for the runoff simulation in the post-flood period. Results reveal that the selected parameters, A, B, and D, have a significant (p < 0.05) effect on the runoff simulation; among them, the significance of D is the strongest. Similarly, the interaction effects of A and C have an important impact on the runoff simulation (Table 6).

The individual parameter effects for the runoff simulation in the post-flood period are shown in Figure 12. The result shows that the level change of D has the greatest negative effect on the runoff simulation. In essence, the above results are attributable to the continuous precipitation decrease, and the river runoff mainly depends on the recharge of the interflow, e.g., the greater the A value, the less recharge in the soil. Moreover, a decrease in soil evaporation is caused by lower temperature. The nonlinear variation of the D curve is mainly due to the complex physical characteristics of soil. When D is at the low level, the soil viscosity is larger, and the interflow has less recharge for the runoff. As D reaches the high level, less precipitation and lower porosity combine to reduce runoff in the post-flood period [43].

The interaction effects of parameters in the post-flood period are shown in Figure 13. The results disclose that the interaction effects of D and other parameters are negative for the runoff simulation, particularly the interaction effects of A and C.

4.3.5. Contributions of Parameter Individual and Interaction Effects for the Runoff Simulation in Four Periods

Figure 14 accurately quantifies the contribution of individual and interaction parameter effects to runoff in the four periods. The results show that the parameter contribution to the runoff simulation is significant for the different periods.

In detail, the parameter D have the greatest impact on the runoff simulation in the non-flood and post-flood period, contributing 0.99 and 0.97, respectively. The contribution of D to the runoff simulation is reduced while others increase in the pre-flood period—especially A, which contributes 0.26 (i.e., 26%). However, the contribution of A to the runoff simulation is the most significant in the flood period (0.60), implying that the infiltration-excess runoff production caused by concentrated precipitation has a great influence on this parameter, while soil porosity and soil moisture have great influence on D. The influence of temperature on B is also significant.

5. Conclusions

In this study, using CMADS to drive the SWAT model, we developed a SUFI-2-based multilevel-factorial-analysis method in order to disclose the effect of parameters in the SWAT model on runoff simulations. Here, the SUFI-2 was used to explore the CI of identified sensitive parameters.

Subsequently, we applied the multilevel-factorial-analysis method to explore the individual and interactive effects of parameters on the runoff simulation in different periods (i.e., non-flood, pre-flood, flood, and post-flood). The developed method was exemplarily applied to the source region of the Yellow River, due to its key role in water resource supply. Important conclusions drawn from this study are as follows:

(1)

The influence of parameters CN2, ESCO, CH_K2, and SOL_BD (i.e., A, B, C, and D, respectively) on the runoff simulation is significant in different periods (Figure 2, Table 3, Table 4, Table 5 and Table 6). In general, the linear individual effects of factors A, B, and D, as well as the AD interaction effects, are thus significant, while the others have little influence on the response.

(2)

The contributions of different parameters to the runoff simulation are different in different periods (Figure 14). The effect of soil bulk density (D) on the runoff simulation is significant in four periods, contributing 0.99, 0.73, 0.30, and 0.97, respectively. The effect of the initial SCS runoff curve number (A) on the runoff simulation is significant in the non-flood and flood periods, contributing 0.26 and 0.60, respectively.

(3)

The interaction effects of parameters on runoff simulation are significant in the flood period. Take parameters A and D as an example: The changes differ across the three levels of parameter D, depending on the level of parameter A. The slope curve is distinctly different between parameters A and D. This reveals the interaction of A and D has a significant influence on the runoff simulation. Therefore, the parameter interaction must be emphasized in flood periods.

(4)

In essence, the soil bulk density moisture content and infiltration-excess runoff production are important water inputs for the hydrological system in the source region of the Yellow River. It is further explained that soil bulk density will affect the loss of surface runoff and river recharge groundwater.