Sensitivity analysis of hydrological model parameters based on improved Morris method with the double-Latin hypercube sampling

WU

initial soil water content in the upper layer

WL

initial soil water content in the lower layer

WD

deep initial soil water content

WM

total initial soil water content

K

evaporation conversion factor

C

evaporation coefficient in the deep layer

B

square times of the storage capacity curve

IMP

impervious area as a percentage of total watershed area

SM

free water storage capacity

EX

free water storage capacity curve index

KG

groundwater outflow coefficient

KSS

loamy mid-stream outflow coefficient

KKG

subsurface runoff abatement factor

KKSS

abatement coefficient of mid-soil flow

ABC

Artificial Bee Colony Optimization Algorithm

SFLA

Shuffled Flog Leaping Algorithm

DC

Nash and Sutcliffe proposed deterministic coefficients

R²

The correlation coefficient squared

P_v

The peak error

D_v

The relative error

DC

Nash and Sutcliffe proposed deterministic coefficients

Hydrological models tend to use the simplified mathematical formulas or physical equations to describe conceptualized hydrological simulations and predict processes. In this way, it has more parameters and the complicated characteristics of uncertainty, high dimensionality and high nonlinearity influenced by many factors such as climate, meteorology, astronomy, subsurface, humanities, and so on, which leads to a very complex parameter rate determination process. Even small differences in parameters may lead to large differences in the simulation results of the model. As a matter of fact, only a small number of parameters play a key role in the model simulation.

In order to save computational resources, it is necessary to do sensitivity analysis of the hydrological model. The influence of each influencing factor on the model is quantitatively studied to identify the main factors that determine the variation of the simulation results of the hydrological model. In this way, the main effort can be concentrated on improving the observation accuracy of these major influencing factors, and even if those minor factors are less considered or ignored, a numerical model with simulation results close to the actual results can be obtained. Therefore, sensitivity analysis of the parameters is a crucial step in the hydrological simulation process.

In recent years, considerable efforts have been made to carry out sensitivity analysis of the parameters. The sensitivity analysis method is used to study the variability of the system at small uptake of the system influence factors, and the contribution of each influencing factor to the system can be quantified. The Morris method is one of the mainstream methods for sensitivity analysis of hydrological model parameters at present (Vanrolleghem et al. 2015; Jaros et al. 2019; Nielsen & Thorndahl 2019; Sreedevi & Eldho 2019; Huang et al. 2020; Pang et al. 2021; Wu et al. 2022), which is usually random sampling or using Latin sampling. Huang et al. (2018), Liu et al. (2019) and Sun et al. (2020) used the Morris method with a random sample to analyze the sensitivity of the hydrological model parameters. Meng (2012), Manache & Melching (2008) and Demirel et al. (2018) used the Morris method with Latin sampling to analyze the sensitivity of the hydrological model parameters.

Many mathematicians had improved the Morris method and applied it to hydrological model parameters. Nabi et al. (2021) analyzed the potential of Morris to achieve results tantamount to much cheaper computational volume, by using a new approach of increasing the number of replications for the Morris algorithm. Chen et al. (2012) conducted the parameter sensitivity analysis of the HEC-HMS model by modified Morris method using the independent variable in fixed steps variation. Wang et al. (2015) studied the sensitivity analysis of distributed hydrological models based on the vibrational method. Liu et al. (2015) quantitatively analyzed the local sensitivity of hydrological parameters of the ICM model by modified Morris screening method. Wang et al. (2016) analyzed the local sensitivity of SWMM model parameters by the modified Morris method. Chang et al. (2016) conducted the sensitivity analysis of InfoWorks ICM model parameters by the modified Morris method. Liu et al. (2021) analyzed the sensitivities of eight parameters in the study area by using the locally modified Morris screening method.

The modified Morris screening method mentioned in the above study is based on a percentage variation of the parameters by a fixed value, with the final sensitivity discriminant taking multiple averages of the Morris coefficients. In this paper, the improved Morris method also takes multiple averages of the Morris coefficients, but a different method is used in the parameter taking method.

Thus, there are some shortcomings in the current research on the application of the Morris method in hydrological models. For example, the sensitivity analysis of hydrological model parameters using the Morris method only considers the interaction between parameters due to the change of parameter values after fixing the parameter order, and does not consider the interaction between parameters due to the change of parameter order, which cannot reflect the parameter sensitivity completely and accurately, and then will have a greater impact on the model calibration.

In order to solve the shortcomings of the current Morris method, the optimal set of parameters of the Xin'anjiang model is first found by an artificial intelligence optimization algorithm (Li et al. 2013), and the model parameters are sampled within the specified range of value changes by using the Latin sampling method, and then the sequential arrangement of parameters is sampled by using the Latin sampling method to construct the input conditions of different combinations of the model and explore the sensitivity of the analyzed hydrological model parameters. This is the first time that the Morris method was improved using the double-Latin sampling method and applied to the sensitivity analysis of the Xin'anjiang model parameters.

Take Zhang River in Dayu County as a study area. Dayu County is located in the southwest end of Jiangxi Province, southwest of Ganzhou City, upstream of the Zhangjiang River, with a geographical location of 114°22′ ∼ 114°44′ east diameter and 25°15′ ∼ 25°37′ north latitude (Figure 2). Dayu County is in the upper reaches of the Zhang River, which joins with the Fujiang River to form an open basin area. And this entire basin is a key flood control area vulnerable to flooding. The Zhang River runs through the county from west to east, with the main tributaries Dangping River in the north and Wulishan River in the south. The two tributaries of Dangping River and Wulishan River converge into the Zhang River at Zhongqiao Bridge and Mudanting Park, respectively. The Zhang River basin covers an area of about 700 km², which is a subtropical Southeast Asian monsoon humid climate zone with mild climate, abundant sunshine, abundant rainfall, long frost-free period, and four distinct seasons. The average multi-year temperature is 18.4°C, the extreme maximum temperature is 39.8°C, and the extreme minimum temperature is −7.1°C.

There are 14 parameters in the three-source Xin'anjiang model, including five evaporation parameters, the evaporation conversion factor (K), the initial soil water content in the upper layer (WU), the initial soil water content in the lower layer (WL), deep initial soil water content (WD), and evaporation coefficient in the deep layer (C). Three flow-production parameters: total initial soil water content (WM), square times of the storage capacity curve (B), and impervious area as a percentage of total watershed area (IMP). Four sub-source parameters: free water storage capacity (SM), free water storage capacity curve index (EX), groundwater outflow coefficient (KG), and loamy mid-stream outflow coefficient (KSS). Two sink-in parameters: subsurface runoff abatement factor (KKG) and abatement coefficient of mid-soil flow (KKSS).

Due to the high degree of nonlinearity, it is difficult to obtain the global optimal parameters for the three-source Xin'anjiang model. Table 1 shows the range of values for the main parameters of the three-source Xin'anjiang model. The range of values of the parameters is based on a large number of relevant practical studies and manual experience. In particular, the value of total initial soil water content (WM) ranges from the sum of initial soil water content in the upper layer (WU), initial soil water content in the lower layer (WL), and deep initial soil water content (WD).

In the formula, O_s is the average value of simulated runoff and O_a is the average value of measured runoff.

Under the condition of other parameters keeping constant, local sensitivity analysis examines the extent to which model results are affected by a small change in the neighborhood around the initial value of a parameter.

In the formula, I represents the sensitivity discriminant factor. Y_i and Y_i+1 represent the output results of i and i+ 1, respectively, Y₀ represents the initial value of the parameter, and represent the change percentage of the parameter value relative to that of the Y₀ value when the model of i-th and i+ 1 run, respectively. N represents the number of model operations. The higher value of parameter sensitivity I, the more sensitive the model results influenced by the parameter.

The sensitivity index is a dimensionless index that reflects the sensitivity extent of model simulation results caused by parameter changes in order to verify the sensitivity of the hydrological model parameters. The paper classifies the sensitivity values into five levels (Lenhart et al. 2002) as shown in Table 2.

Although the local sensitivity analysis method is quite simple since it analyzes only one parameter at a time, such a method does not consider the situation that parameters interact with and influence each other. Meanwhile, when analyzing a single parameter, its sensitivity would be influenced by any value change of other parameters to some degree. Therefore, the results from the local sensitivity analysis of parameters often have limitations and do not accurately reflect the sensitivity of parameters.

The principle of global sensitivity analysis is firstly to vary the parameter value based on the given range of model parameters, then to analyze the comprehensive effect on the model when the parameters take different values and interact with each other.

The global sensitivity analysis method is an improvement on the local sensitivity analysis method so that the accuracy of the parameter sensitivity analysis results is higher. Compared with the local sensitivity analysis method, global sensitivity analysis has the following characteristics (Xu et al. 2004; Shu et al. 2007): (1) considering the influence of sensitivity value if parameter values change during the analysis process. (2) The model results are more objective and comprehensive because of being gained from the simultaneous change of all parameters being affected by each other. The global sensitivity analysis can be divided into qualitative sensitivity analysis and quantitative sensitivity analysis. Qualitative global sensitivity analysis is to analyze the effect of model results with any changes in model parameters, in order to obtain the sequence of sensitivity values for each parameter of the model with a small amount of computation. Quantitative global sensitivity analysis quantifies the impact frequency of changes on the model output when each parameter changes.

In the matrix, each column represents a parameter, element 0 represents the initial value of the parameter, and element 1 represents the parameter value after the change. Among them there is only one different parameter per two adjacent rows. Therefore, the corresponding model results could be obtained by substituting two adjacent rows into the model separately when calculating. The sensitivity value of that different parameter between two rows could be calculated by comparing the difference between the results of these two rows. Meanwhile, the sensitivity of m parameters is obtained by taking m sets of adjacent elements. The global sensitivity resulting from the combination of two or more parameters could be obtained by taking any two rows parameter of the matrix, then input them into the model and finally calculate the result difference between the corresponding models. The formula of the sensitivity discriminant factor is the same as that of the local sensitivity analysis.

The principle of the Morris method is to calculate the sensitivity value of each parameter in the form of matrix B by varying each parameter to a certain range and arrange them in a certain combination. However, since this method does not specify the parameters form and combination order, the number of parameters to be analyzed is ten at most, and the calculation workload would be enormous if considering all circumstances. Therefore, it is assumed that the parameter variation ranges from −20 to 20%, the fixed step is 5%, and a 9 parameter variation value including the initial value. Under this condition, it proposed the double-Latin hyper-geometric sampling method, namely to make samples from the form of parameter variation and the order of combination separately.

The Latin hyper-geometric sampling method is more effective than the common random sampling one (Iman 2008). It is to divide the vertical axis of the probability distribution function of each parameter into several equally spaced intervals without no overlap with each other, and each interval corresponds to an equally probable interval of the parameters of the horizontal axis. In short, it is to fix the order of all parameters, and take m samples when there are m change values. The practical method is to arrange the values of n parameters into n×m matrix and randomly change the order of the elements in each column, thus to obtain m inputs for all parameters.

The double-Latin hyper-geometric sampling uses the same method to change the order of parameters while taking different values of parameter variations. That is, it disrupts the order of each column element and each row element randomly and nonrepeatedly at the same time, then n×m samples can be obtained.

WU, WL, and WD are only relevant to the watershed, as a result, they are considered as nonsensitive parameters (Shu et al. 2008). The paper explores the local and global sensitivity of the first 10 parameters in Table 1.

In order to verify the superiority of the algorithm, three typical rainfall runoff processes are selected as test examples in this paper, the first two being for parametric rate analysis and the last one being for parametric verification analysis. The three rainfall runoff processes are all double flood peaks, with maximum runoff volumes of 471, 507, and 765 m³/s, respectively; the simulated time period is 49, 22, and 33 h, with a 1-h time interval; initial free water depths on average are 5, 10, and 5 mm, in order; initial loamy mid-stream flows are 5, 10, and 5 m³/s, respectively; initial subsurface runoff flow is 10, 10 , 10 m³/s; dimensionless unit hydrograph UH is taken for three time periods, 0.3, 0.6, 0.1, respectively.

In this paper, three artificial intelligence optimization methods are taken to optimize the rate determination of three water sources Xin'anjiang model. The above-mentioned three artificial intelligence optimization methods are Artificial Bee Colony Optimization Algorithm (ABC), Shuffled Flog Leaping Algorithm (SFLA), and SCE-UA, respectively. The model optimization rate determination results are shown in Table 3 (Li et al. 2013).

Nash and Sutcliffe proposed deterministic coefficients DC as an evaluation criterion for the effectiveness of hydrological model simulations, taking values in the range (−∞, 1] (the closer to one, the better the simulation effect). The correlation coefficient squared R² is used to describe the correlation between the simulated and measured values, taking values in the range [0,1], the closer to 1, the better the simulation effect. The peak error P_v is an index used to describe the peak error of simulated and measured runoff, and when P_v ≤ 10%, the simulation is quite good. The relative error D_v is an index used to describe the error of total simulated and measured runoff, and when D_v ≤ 10%, the simulation is very good.

The simulation results show that DC, R², P_v, and D_v obtained by the three optimization algorithms are all satisfactory. As can be seen from Table 3, the DC and R² after the first rainfall runoff rate determination are above 0.9 and 0.92, and P_v and D_v are within 9 and 11%, respectively, with very good simulation effect; the DC and R² after the second rainfall runoff rate determination reach above 0.81 and 0.83, and P_v and D_v are within 5 and 7%, respectively, with better simulation effect. The validation results of the third rainfall runoff had DC and R² above 0.83 and 0.85, respectively, while P_v and D_v were within 18 and 9%, respectively, which were better simulated.

In comparison, it is found that the ABC algorithm has the strongest optimization-seeking ability, fastest convergence speed, best overall performance, and highest stability among the above-mentioned three optimization algorithms. Therefore, a set of parameters with the best rate validation results during the rainfall runoff process were selected for three times as the benchmark for studying sensitivity analysis of parameters, with the results in detail shown in Table 4.

The rainfall runoff example can be found in Table 5 as follows, (1) the sensitivity values of parameters C, IMP, SM, EX, KG, and KKG are all less than 0.01, which are insensitive parameters. The sensitivity values of parameters K, B, and KSS are in the range of [0.01, 0.1], which are fewer sensitive parameters. The sensitivity value of parameter KKSS is greater than 1, which is the high-sensitivity parameter. (2) Among these parameters, it is shown that the soil mid-stream receding coefficient KKSS has a relatively strong influence on runoff and other parameters have little influence on runoff. In order of sensitivity, KKSS, B, KSS, K, KG, SM, IMP, EX, KG, C are in descending order. (3) Since the local sensitivity analysis does not analyze the interaction between parameters, the calculation results still have limitations and cannot accurately represent the sensitivity of parameters, so the analysis results only provide a reference basis for the global sensitivity analysis.

Based on local sensitivity analysis, the effect of a single parameter change on the total runoff simulation is linearly correlated. Therefore, to conduct global sensitivity analysis using the Morris method, it would be better to assume that there is only 10% increase in each parameter. In this way, we can not only get credible analysis results but also simplify the calculation process.

By substituting each combination of matrix B parameter into the three water Xin'anjiang model, we could find the differences between the simulated runoff results, and obtain the sensitivity of each parameter and the interaction between the parameters. We should put the order in the optimized parameters, without considering any change in the order of the parameter values (namely the single-Latin sampling method), the variability of the model runoff results for 55 parameter combinations was obtained by Morris method, as shown in Table 6. The magnitude (%) of runoff changes for various combinations of methods would be shown in the bar graphs, as shown in Figure 4.

The condition of model input was determined by the double-Latin hyper-geometric sampling method, and then substituted in the order of the B matrix by Morris method. In this way the results of the parameters can be obtained as shown in Table 7 and Figure 5.

For this example of rainfall runoff, from Table 7, we can see that: (1) if the sensitivity values of parameters C, SM, IMP, EX, KG, and KKG are all less than 0.01, they belong to insensitive parameters; if the sensitivity values of parameters K and B are in the range of [0.01, 0.1], they belong to less sensitive parameters; if the sensitivity values of parameter KSS is in the range of [0.1, 0.2], it belongs to low-sensitivity parameters; if the sensitivity value of parameter KKSS is greater than 1, it belongs to the high-sensitivity parameter. (2) among these parameters, it reflects that parameter KSS has an influence on runoff to some degree, and parameter KKSS has a relatively strong influence on runoff. Other parameters have less influence on runoff; in order of sensitivity values from high to low, they are KKSS, B, KSS, K, KG, SM, IMP, EX, KG, and C. (3) The local sensitivity values of the Xin'anjiang model parameters (K, B, IMP, KKG and KKSS) are 0.0422, 0.0747, 0.0004, 0.0062, and 2.3348, respectively, when only individual parameter variations are considered. After considering the interaction between the parameters, the global sensitivities of the Xin'anjiang model parameters (K, B, IMP, KKG and KKSS) are 0.0328, 0.0659, 0.0001, 0.0052, and 1.9634, respectively. It can be seen that the sensitivity values of parameters are reduced by the interaction between parameters. (4) The local sensitivity values of the Xin'anjiang model parameters (SM, EX, KG and KSS) are 0.0009, 0.0004, 0.0004, and 0.0727, respectively. The global sensitivities of the Xin'anjiang model parameters (SM, EX, KG and KSS) are 0.0061, 0.0024, 0.0007 and 0.1017, respectively. It can be seen that the sensitivity values of parameters are increased by the interaction between parameters. (5) The results of the improved Morris global sensitivity analysis more fully reflect the positive and negative correlation effects of the interaction between parameters. Those results can be used as the final reference basis for model parameter sensitivity analysis.

In the paper, we propose the Morris global sensitivity analysis method based on double-Latin sampling, and conduct sensitivity analysis by local sampling, single-Latin and double-Latin sampling. The main conclusions drawn from this paper can be summarized below:

• (1)

  The parameters C, SM, IMP, EX, KG, and KKG are insensitive parameters. The parameters K and B are fewer sensitive parameters. The parameter KSS is the low-sensitivity parameter, and KKSS is the high-sensitivity parameter.

• (2)

  Compared with the local sensitivity method and the modified Morris method with double-Latin sampling, the results show that the sensitivity values of Xin'anjiang model parameters (K, B, IMP, KKG, and KKSS) are reduced by the interactions among parameters, and the sensitivity values of parameters (SM, EX, KG, and KSS) are increased by the interactions between parameters. The results of the improved Morris method with double-Latin sampling more fully reflect the positive and negative correlation effects of the interaction between parameters.

• (3)

  Compared with the local sensitivity analysis method and Morris method with single Latin sampling, the improved Morris method with double-Latin sampling fully considers the influence of the convenience of parameter taking values and the sequential arrangement of parameter combinations on the interaction between parameters, thus its calculation results are superior and more accurate. None of the other global sensitivity analysis methods consider the independence between parameters.

• (4)

  In this paper, when we conduct global sensitivity analysis by Morris method, the variation in the values of model parameters has been simplified, but such simplification is feasible to find the interaction effects between the parameters.

• (5)

  In this study, the number of selected flood process fields is small, and the results after the optimization of the model parameters have some limitations, and more typical field flood processes will be selected for further study in the next step. However, the Morris method based on double-Latin sampling is proposed in this paper, which shows a better superiority in theory and results, and can be extended.

This research was funded by the Key R&D Program of Jiangxi (Grant Nos. 20212BBG71008, 20203BBGL73224, and 20203BBGL73221), Water Resources Department Science and Technology Project of Jiangxi (Grant No. 202123ZDKT16), National Natural Science Foundation of China (Grant Nos. 52179013, U2240217), and the Belt and Road Special Foundation of the State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering (Grant No. 2020491511).

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.