Parameter sensitivity analysis for a physically based distributed hydrological model based on Morris' screening method

Achieving accurate parameter calibration is a relatively difficult problem in development and application of physically based distributed hydrological models. In the current research, Morris' screening method (one factor at a time [OAT]) and an error evaluation calibration were jointly adopted for parameter sensitivity analysis (PSA) for a distributed hydrological model developed for a single mountainous terrain. The results indicated that (1) the order of sensitivity of the parameters for rainfall−runoff calculation is thickness of the first soil layer, coefficient of permeability, infiltration velocity, porosity, and coefficient of roughness; (2) the order of the parameters for snowmelt calculation is bulk transfer coefficient, temperature correction factor, and albedo of snow cover layer surface; and (3) the value change of each parameter of the first layer in the soil vertical profile mainly impacts the numerical value of discharge of the surface flow, while the value change of each parameter of the second layer impacts the numerical value of discharge of the surface flow less, but results in the peak misplacement to some extent. The results are expected to provide methodology references for PSA and calibration for the physically based distributed hydrological model which developed for mountainous or otherwise single terrain.


| INTRODUCTION
Physically based distributed hydrological models are widely used for hydrological simulations in various environments. Achieving the accurate parameter calibration is a relatively difficult problem in the processes of development and application of a physically based distributed hydrological model (Sun et al., 2017). Parameter sensitivity analysis (PSA) is carried out to reveal model response caused by parameter changes, which is not only one of the important parts of model uncertainty analysis, but also a key link in model development and evaluation. PSA is the process of determining the rate of change in model output with respect to changes in model inputs. It is a necessary process to identify key parameters and the parameter precision required for calibration (Ma et al., 2000). PSA promotes understanding of model features and improves the structure of model stability (Wang, Xia, Liu, Ou, & Zhang, 2007). PSA of a hydrological model is the key to model uncertainty quantification; however, how to effectively validate the model and identify the dominant parameters for distributed hydrological models is an impediment to achieving parameter optimization (Song, Zhan, Xia, & Kong, 2012). The complexity of models and the large number of parameters make optimization of distributed hydrological models difficult. PSA speeds the calibration and validation process by identifying which parameters produce the largest changes in model output in response to perturbations (Demaria, Bart, & Thorsten, 2007;Lan, Thomas, & Alan, 2011;Scollo, Tarantola, Bonadonna, Coltelli, & Saltelli, 2008). Some available methods for conventional PSA have been proposed. Chen, Liang, and Chen (2011) used a perturbation analysis method to evaluate the variability of the parameter sensitivities of SWAT (Soil and Water Assessment Tool) as well as its influence on the modeling results of hydrological processes with different landscape cover types in Hailuogou Valley, China. Fang, Yang, Chen, Xu, and Philippe (2015) conducted PSA by using a sensitivity analysis approach called the SDP (State-Dependent Parameter) method for an application of SWAT to the Kaidu River Basin in China's Tianshan Mountains. Song et al. (2012) assessed the sensitivity of parameters of a distributed hydrological model (DTVGM, distributed time − variant gain model) by combining water balance (WB) coefficient, Nash−Sutcliffe (NS) coefficient, and correlation coefficient (RC). Zhan, Song, Xia, and Tong (2013) proposed an approach that integrates Morris' screening method (one factor at a time [OAT]) with a qualitative analysis method based on a statistical emulator to analyze the parameter sensitivity of a DTVGM. Jiang, Liu, Li, Liu, and Wang (2015) performed PSA and calibration for a distributed HIMS (Hydro-Information Modeling System) model applied to China's Luanhe River Basin by using OAT. From the above examples, two aspects of PSA of a distributed hydrological model can be identified. One is that OAT has been widely used and has been proven to be an useful approach (Gao, Sorooshian, & Gupta, 1996;Sun & Bosilovich, 1996); the other is that PSA have mostly been carried out for some typical distributed hydrological models, such as SWAT, DTVGM, and HIMS, etc. In order to facilitate the application of such models in various scale basins under different terrains, the number of parameters must be reduced while possessing good universality. Although the model parameters have physical meaning, the physical properties of some parameters were reduced to a certain extent, and, moreover, the model parameters also need to be re-calibrated when these models are applied to different basins.
The parameters of a physically based distributed hydrological model developed for a specified basin located in a single type of terrain maintain the original physical properties (Fu, Lu, & Huang, 2014;Huang, Hinokidani, Yasuda, & Kajikawa, 2008). To achieve efficient PSA for a physical distributed hydrological model developed for a specified basin in a single terrain type is significant for parameter calibration and model validation of such a model, and is significant for promoting the development and application of physically based distributed hydrological models. In the current state of research, OAT combined with an error evaluation calibration have been adopted and jointly used to analyze the sensitivity of parameters of an independently developed, physically based distributed hydrological model. Quantitative analysis was carried out for the parameter change impact on the numerical result of rainfall−runoff and snowmelt on the basis to which order of sensitivity each main parameter was assigned. The results are expected to provide the methodology reference for further research on parameter calibration and model validation of physically based distributed hydrological models for basins located in a single type of terrain.
2 | METHODOLOGY 2.1 | Overview of development of a distributed hydrological model An approach to distributed hydrological model development combining a GIS (Geographic Information System) with the kinematic wave model has been used widely in rainfall−runoff modeling for basins of different scales (Chua, Wong, & Sriramula, 2008;Yomoto & Islam, 1992). The kinematic wave model is based on the resistance law of river channels; it is a physically based model, as it uses physical parameters to characterize a basin (Tanaka, 2003;Tanaka, Fujita, & Kudo, 1999). The energy budget approach has been mostly used in snowmelt calculation (Herrero, Polo, Moñino, & Losada, 2009;Ohara & Kavvas, 2006). In the current state of research, a GISbased distributed hydrological model, developed by kinematic wave equations for rainfall−runoff calculation combined with the energy budget method for snowmelt calculation based on the physical conditions of the underlying surface and meteorology of Japan's Bukuro River Basin, was used to carry out a PSA.

| Study area
The Bukuro River Basin (134 12 0 -134 26 0 E, 35 25 0 -35 29 0 N) is located in the mountainous area of southwestern Honshu in Japan and it is a primary tributary of the Sendai River. The basin was majorly covered by natural shrub. Mean annual rainfall for this area is approximately 2,000 mm, snow cover and snowmelt generally occurs in the period from November to April. As a typical feature of a mountainous basin, the topography within the basin is very undulating, and the maximum height difference between the mountaintop and valley is over 900 m. The gradient of natural slope generally exceeds 0.20, while the mean riverbed gradient is greater than 0.005. In addition, there are more than 20 tributaries with a relatively concentrated distribution ( Figure 1) that generally cause flood events due to the short duration of high−intensity rainfall.

| Governing equations
Equations of kinematic wave model (Equations (1)-(2)) were used for the surface-flow calculation as follows: where t is the time factor (s), x is the spatial factor in flow direction (m), h is the depth of water (m), q is the discharge of unit width (m 2 /s), r is the effective rainfall (m/s), f 1 is the mean infiltration rate of the first soil layer (m/s), v is the flow velocity (m/s), R is the hydraulic radius (m), n is the coefficient of roughness (sÁm −1/3 ), and I is the gradient.
Equations (3)-(5) were used to calculate the infiltration flow in the first soil layer ( Figure 2): where λ is the effective porosity, h is the infiltration flow depth (m), q is the discharge of unit width of infiltration flow (m 2 /s), f 2 is the mean infiltration rate of the second soil layer (m/s), ET is the evapotranspiration (approximately replaced by a loss coefficient in the process of calculation) (m/s), v is the infiltration flow velocity (m/s), k is the coefficient of permeability (m/s), H is the hydraulic head of the infiltration flow (m), and the other factors correspond to those mentioned above.
Snowmelt calculation was achieved by calculating the heat budget process in the snow cover layer based on the energy budget approach using Equations (6)-(8) (Huang, Wang, Hinokidani, & Yamamoto, 2012): where Q g is the snowmelt heat obtained from the surface and bottom of the snow cover layer (W/m 2 ), R n is incident radiation (W/m 2 ), ε is the emissivity of the snow cover layer, σ is the Stefan-Boltzmann constant [5.67 × 10 −8 W/(m 2 ÁK 4 )], T s is the surface temperature of the snow cover layer (K), H is the sensible heat flux (W/m 2 ), lE is the latent heat flux (W/m 2 ), Q b is the heat exchange between the bottom of the snow cover layer and the soil (W/m 2 ), Q r is the heat contained in rain (W/m 2 ), S is the horizontal solar radiation (W/m 2 ), L is the downward atmospheric radiation (W/m 2 ), ref is the albedo of the snow cover layer surface, and T is the air temperature (K). Estimation of the sensible and latent heat flux is dependent on meteorological data, as seen in Equations (9) and (10): where c p is the specific heat of air at a constant pressure [J/(kgÁK)], q is the air density (kg/m 3 ), C H and C E represent the bulk coefficients corresponding to sensible and latent heat, respectively, U represents the wind speed (m/s), l is the latent heat of vaporization (J/kg), q sat (T) represents the saturation specific humidity corresponding to air temperature, Δ represents slope of the saturation vapor-pressure at air temperature (kPa C −1 ), and rh is the relative humidity; other factors were defined above.

| Modeling of soil vertical profile
Water movement is constrained by the actual physical conditions of the underlying topography in this model of the vertical soil profile. Modeling of the soil vertical profile of the Bukuro River Basin was developed based on the field survey results. The model was composed of the slope sector and the river channel sector, which involves the modeling of the thickness of the soil layers and the distribution of groundwater (phreatic water) ( Figure 2). Both the slope and river channel were separated into two soil layers in the vertical profile from the ground surface to bedrock.

| Method of PSA
OAT identifies the non-influential model factors through a limited number of model runs (Morris, 1991), and, as a screening method based on elementary effects, identifies a subset of inputs that have the greatest influence on outputs (Demaria et al., 2007;Scollo et al., 2008). By using this method for PSA, the value of a selected parameter was randomly changed; meanwhile, the values of other parameters were fixed. Parameter sensitivity was evaluated by analyzing the influence of parameter change on varying degrees of model output.
Order of importance for each parameter can be assigned according to the results of PSA (Francos, 2003;Saltelli et al., 2008). An error criterion (Equation (11)) recommended by the Japanese River Society (1997) was adopted for PSA combined with OAT. This error criterion was used to evaluate the error of the simulation results between the observed flow and the calculated results of runoff. Its allowable error was less than 0.03: where E r is the error, n is the number of calculations, Q o (i) is the observed discharge at moment i (m 3 s −1 ), Q c (i) is the calculated discharge at moment i (m 3 s −1 ), and Q op is the maximum discharge of the observed flow during the calculated time (m 3 s −1 ).

| Parameters for sensitivity analysis
Physical parameters such as the infiltration velocity (f ), coefficient of permeability (k), coefficient of roughness (n), layer thickness (H s ), and first-layer porosity (λ) of the soil vertical profile model ( Figure 2) have been adopted for sensitivity analysis of rainfall−runoff simulation ( Table 1). The main parameters of the snowmelt calculation adopted for sensitivity analysis of the snowmelt calculation were the bulk transfer coefficient (C H , C E ), the albedo of the snow cover layer surface (ref ), and the temperature correction factor (α T ) ( Table 2).

| Calculation periods
The accuracy of peak simulation is an important index for evaluating the function of the distributed hydrological model. Therefore, it is very necessary to evaluate the influence of parameter variation on the accuracy of peak simulation. A period from July 4 to July 24, 1997 was selected for rainfall−runoff calculation for PSA since two peaks occurred on July 12 and 17, respectively. Another period from February 1 to May 31, 1997 was selected for snowmelt calculation as this period covered the main snowmelt period in 1997 of the study area. Different numerical results were obtained due to multiple parameter-value changes, and it is difficult to show all the numerical results in the form of diagrams. Diagrams for each separate parameter-value change that significantly impacted the numerical results are given in the following discussions of each of the PSA conducted.

| RESULTS AND DISCUSSION
3.1 | PSA for main parameters of rainfall−runoff calculation

| Coefficient of permeability
The soil vertical profile model ( Figure 2) was assumed as homogenous as the topography and underlying surface experienced similar conditions within the range of the Bukuro River Basin. A same set of parameters value was thereby used for each sub-catchment, by which not only the amount of calculation can be significantly reduced, but the OAT-based results can be contrasted more clearly. The set of parameter values is shown in Table 1. Different values of the coefficient of permeability were chosen for the first and second layers in the soil vertical model (Figure 2) for the three cases designated Case-1, Case-2, and Case-3; the simulation results are shown in Figure 3. Figure 3 illustrates that the numerical results of Case-1 and Case-3 appropriately fit the observed flow discharge. The calculated peak flow of Case-2 is significantly larger than that of Case-1 and Case-3. The main reason is that coefficient of permeability of the first layer of Case-2 is smaller than that of Case-1 and Case-3, which caused groundwater runoff to the downstream in the first layer to be lower than that of Case-1 and Case-3. Under the same conditions of the infiltration rate from ground surface to the first layer and from the first layer to the second layer, the depth of groundwater in the first layer of Case-2 decreased more slowly than that of Cases-1 and Case-3, which led to the infiltration from the ground surface to the first layer to decrease in Case-2, and resulted in the peak flow of Case-2 being significantly larger than that of Case-1 and Case-3. Compared with Case-1 and Case-3, flow discharge in Case-2 was mostly less than that of Case-1 and Case-3 during the period except for the concentrated rainfall. The value of the coefficient of permeability of the second layer in Case-3 was larger than that of Case-1 and Case-2. Good fitting exists between runoff curves of Case-1 and Case-3 both in temporal process and in waveforms. No obvious difference exists in peak and non-peak periods between Case-1 and Case-3, while the peak of Case-3 was slightly lower than that of Case-1 on July 12, and a slight misplacement of peak flow appeared between Case-1 and Case-3 on July 17.  Figure 4 represents the decreased flow-discharge process of the three cases after the peak flow occurred on July 17. An obvious difference exists between the runoff curves of Case-2, Case-1, and Case-3. A visible difference also exists between the runoff curves of Case-1 and Case-3, mainly manifesting that flow-discharge decreased inconsistently in terms of the temporal process. Compared with Case-1, the flow-discharge decrease rate of Case-3 was faster. According to the above analysis that the following can be concluded: The value change of permeability of the first layer significantly impacted the flow scale, especially the peak flow, while the value change of the permeability of the second layer has almost no influence on flow scale. Rather, it mainly impacted the base flow, and impacted runoff occurrence in the temporal process to a certain degree, and caused a slight misplacement of runoff curves.
As the simulated results of the peaks, the peak flow of Case-2 were 1.75 times and 2.1 times those of Case-1 and Case-3, respectively (July 12), and 2.9 times those of Case-1 and Case-3 (July 17) (Figure 3).

| Thickness of layer
As the soil vertical profile (Figure 2) was assumed to be a homogenous model, the thicknesses of the first and second layers of the slope and river channel were assumed to be 0.35 and 2.0 m (Case-1), 1.0 and 2.0 m (Case-2), and 0.35 and 5.0 m (Case-3), respectively (Table 1). The simulated results are shown in Figure 5. Figure 5 depicts an obvious difference in the shape of the runoff curves between Case-2 and Case-1, the peak flow of Case-2 was significant lower than that of Case-1 and Case-3. With the same vertical infiltration conditions, the thickness of the first layer of Case-2 was larger than that of Case-1 and Case-3, which caused its water-holding capacity of the first layer to increase, and the saturation of the first layer needs more rainwater infiltration. Compared with Case-1 and Case-3, most of the rainfall infiltrated into the first layer in Case-2 in the process of model calculation, which resulted in a significant decrease of peak flow in Case-2. The simulated result of Case-3 was similar to that of Case-1 despite the value of the thickness of the second layer of Case-3 being greater than that of Case-1. The runoff curves of Case-1 and Case-3 represent good fitting in the temporal process, except for a slight misplacement appearing in the peck flow on July 17.
The peak flows of Case−2 were 5.7 and 2.7 times that of Case−1 during the intensive rainfall periods on July 12 and July 17, respectively. Compared with Case−3, the result of Case−1 more consistently fits the observed flow. From the above analysis, it could be concluded that the change of the thickness of the first layer in the soil vertical profile mainly affected the model result for flow scale, and the change of thickness of the second layer almost had no impact on flow scale, while causing a slight misplacement of peak flow.

| Infiltration velocity
Two cases were selected for a sensitivity analysis for the change of infiltration velocity of the first layer ( Figure 2). The values of infiltration velocity from the ground surface to the first layer were adopted as 6 × 10 −5 m/s for Case-1 and 6 × 10 −3 m/s for Case-2. The simulated results are shown in Figure 6. Figure 6 shows that the simulated results of the peak flows of Case-1 on July 12 and July 17 were much closer to the observed flow, while the simulated results of the peak flows of Case-2 were distinctly lower than the observed one. Equation (11) was used to calculate the errors of the simulated results of the intensive rainfall during periods in which peak flow occurred. The errors of Case-1 and Case-2 were 0.015 and 0.05 (10:00-13:00 on July 12), and 0.004 and 0.032 (8:00-11:00 on July 17), respectively. The errors of Case-1 were within the allowable range of criteria, while the errors of Case-2 exceeded the allowable range. Both the simulated results of Case-1 and Case-2 reproduced the observed flow inconsistently in the temporal process as well as no obvious misplacement of the peak flow. In the case of changing the infiltration velocity of the first layer from 6.0 × 10 −5 to 6.0 × 10 −3 m/s, the results indicated that with increasing infiltration velocity the flow discharge decreased, particularly for peak flow. However, there is no inverse relationship or linear relationship between the increase of infiltration velocity and decrease of peak flow. The result verifies that the change of infiltration velocity of the first layer mainly affected the simulated result of flow scale. As the simulated results, the peak flows of Case-1 were approximately 1.25 times of those of Case−2 on July 12 and July 17, 1997, respectively.

| Coefficient of roughness
The coefficient of roughness is a parameter that integratedly reflects river channel roughness affecting flow, and its value is generally obtained via experiment. It is one of the physical parameters of Manning 0 s averaged velocity formula used to calculate flow velocity. Values of the coefficient of roughness of slope and river channel were adopted as 0.15 and 0.06 for Case-1, 0.15 and 0.30 for Case-2, and 0.3 and 0.06 for Case-3; the simulated results of rainfall−runoff are shown in Figure 7. Figure 7 illustrates that the runoff curves of the three cases have good fitting in term of temporal process. A difference of peak flow exists among the three cases, despite the fact that difference is not obvious. The peaks of Case-1 were 1.2 times that of Case-3 on July 12 and July 17, 1997, respectively.

| Porosity
The porosity of the first layer exhibits a spatial difference within the range of the Bukuro River Basin, and it is difficult to achieve an accurate calibration of the spatial distribution of porosity corresponding to each sub-catchment on the river channel network of the study area ( Figure 1). As the soil vertical profile model (Figure 2) was considered homogenous, the porosity of the first layer was assumed to be 0.35 for Case-1 and 0.5 for Case-2. Figure 8 depicts the runoff curves of Case-1 and Case-2 have a relatively high similarity in shape and good consistency in the temporal process. However, the peak of Case-2 was lower than that of Case-1. The value of the porosity of the first layer of Case-2 was greater than that of Case-1, and with the other parameters being the same, the first layer of Case-2 could reserve more water than Case−1, which resulted in the peak of Case-2 being lower than that of Case-1. The peaks of Case-1 were 1.4 times and 1.3 times that of Case-2 on July 12 and July 17, 1997, respectively.

| Bulk transfer coefficients
The Bukuro River Basin is located in a snow zone, snow cover accompanied by snowmelt naturally occurs in Winter and Spring. The bulk transfer coefficient (C H , C E , C H ≈ C E ) is a parameter used to calculate the sensible heat flux (H) (Equation (9)) and latent heat flux (lE) (Equation (10)) (Niemeyer, Link, Seyfried, & Flerchinger, 2016). Different values for bulk transfer coefficients for the three cases (Case-1, 0.004; Case-2, 0.008; Case-3, 0.002) were adopted. Calculation of snowmelt for the three cases was carried out for the main snowmelt period from February 1 to May 31, 1997, and the simulated results are shown in Figure 9. Figure 9 illustrates that most of the calculated results of Case-2 were higher than those of Case-1 and Case-3 in February and March, whereas they were obviously lower than Case-1 and Case-3 since late April. As the value of the bulk transfer coefficient of Case-2 was the maximum among the three cases, compared with Case-1 and Case-3, the snowmelt of Case-2 earlier occurred in the process of model calculation. With the same model parameters, except for the bulk transfer coefficient, the quantity of snowmelt of Case-2 was the largest, and its flow discharge, mainly caused by snowmelt, was greater than that of Case-1 and Case-3 in February and March. Meanwhile, Case-2 needed a relatively shorter snowmelt time than Case-1 and Case-3, which resulted in its earlier end of snowmelt. The snowmelt time of Case-1 and Case-3 ended later than Case-2, as the flow discharge was greater than that of Case-2 since late April. The value of bulk transfer coefficient of Case-3 was the minimum among the three cases and its snowmelt time ended later than Case-1 and Case-2, as the flow discharge of Case-3 was generally higher than that of Case-1 and Case-2 since late April. Figure 10 depicts the change of the depth of snowcover layer of the three cases during the simulation period. No obvious snowmelt occurred before late February in all three cases, and with increased temperature since late February the snowmelt was faster, with the Case-2 snowmelt being faster than that of Case-1 and that of Case-1 being faster than that of Case-3. The simulated result of Case-1 reproduces the observed flow better than Case-2 and Case-3, which verifies the fact that the value of the bulk transfer coefficient of Case-1 was more reasonable.
Equation (11) was used to calculate the errors of the simulated results of the snowmelt; the errors are 0.004 (Case-1), 0.01 (Case-2), and 0.012 (Case-3). The errors of Case-2 and Case-3 are more than two times that of Case-1.

| Albedo of snow cover layer surface
The albedo of the snow cover layer surface is an essential parameter for calculating the heat budget in the snow  (7)). It is directly related to the altitude and azimuth of the Sun, quality of snow, water content in the snow cover layer, composition of dirt over the snow surface, etc. The albedo of the snow cover layer surface generally changes within a range of 0.4-0.9 and its value significantly decreases during snowmelt (Aoki, Fukabori, & Uchiyama, 1999;Kondo, Numata, & Yamazaki, 1988). The change of albedo of the snow cover layer surface impacts the incident radiation in the snow cover layer (Equation (7)) by impacting horizontal solar radiation. To evaluate the change of the albedo of the snow cover layer surface's impact on the model calculation of snowmelt, three cases of different values of the albedo of the snow cover layer surface were adopted 0.75 for Case-1, 0.50 for Case-2, and 0.30 for Case-3. The simulated results are shown in Figure 11. Figure 11 shows that the flow discharge of Case-1 was mostly lower than that of Case-3 in February and March, while the flow discharge of Case-2 mostly has a value between that of Case-1 and Case-3; from April to early May, the flow discharge of Case-1 and Case-2 was higher than that of Case-3. With an increasing value of the albedo of the snow cover layer surface, the reflected horizontal solar radiation increased, which resulted in a decrease of the incident radiation, and thus the snowmelt heat decreased. Thereupon, the flow discharge of Case-1 was lower than that of Case-2 and Case-3 in February and March. From the beginning of April, although the incident radiation of Case-1 was lower than that of Case-2 and Case-3, due to increased heat exchange between the snow cover layer and soil and to the fact that the snow cover layer absorbed the heat of rainfall (precipitation mainly in form of rainfall since late March), the snowmelt heat in the snow cover layer generally increased. Therefore, the runoff curves of the three cases exhibited different forms since the beginning of April. Figure 12 depicts the change of depth of snow cover layer of Case-1, Case-2, and Case-3. Since late February, although the decreased process of the depth of the snow cover layer was similar for the three cases, the decrease rate was different. Case-3 used the minimum value of albedo and its decreased rate of the depth of the snow cover layer was the maximum; in addition, snowmelt ended earlier in Case-3 than in Case-1 and Case-2. Case-1 adopted the maximum value of the albedo, and its decreased rate of the depth of the snow cover layer was the minimum; its snowmelt ended the latest among the three cases.
The errors were 0.004, 0.005, and 0.007 for Case-1, Case-2, and Case-3, respectively. The simulated result of Case-1 appropriately reproduced the observed runoff.

| Temperature correction factor
Air temperature is a dominant factor for multiple factors (Equations (8)-(10)) related to the snowmelt calculation. Air temperature naturally exhibits spatio-temporal variation due to the change of elevation. Elevation of the upstream in the Bukuro River Basin exceeds 1,000 m and a large difference in elevation generally exists between the mountaintop and valley floor, which leads to obvious temperature differences. Equation (12) was adopted to correct the air temperature when performing snowmelt calculations (Modeling of the snowmelt runoff for dam inflow predictions, 1990).
where α T is the air temperature correction factor ( C/m), h 0 is the elevation of the observation position (m), T 0 is the air temperature of the observation point ( C), T i is the corrected air temperature of each sub-catchment ( C)，i is the number of sub-catchment, and h i is the mean elevation of a sub-catchment (m). The values of the air temperature correction factor of the two cases were assumed to be −0.45 and − 0.65 C/ (100 m) for Case-1 and Case-2, respectively. The simulated results of the runoff process are shown in Figure 13. Figure 13 illustrates that the flow discharge of Case-2 was smaller than that of Case-1 before late February, while it was larger than that of Case-1 from March to the time of the snowmelt ended. From March to May, the flow discharge of Case-2 was mostly larger than the observed result. Because Case-2 used a relatively smaller value for the air temperature correction factor, its temperature was lower than Case-1 in Winter (before late February) and the precipitation in the form of snow was greater than Case-1, which caused Case-2's depth of snow cover layer to be higher than that of Case-1 ( Figure 14). When snowmelt occurred in Winter, the quantity of snowmelt and the flow-scale of Case-2 were lower than those of Case-1. However, with increasing air temperature since the beginning of March, the heat for snowmelt gradually increased in addition to the depth of the snow cover layer of Case-2 being higher than that of Case-1, which integratedly resulted in the flow discharge of Case-2 being larger than that of Case-1. The errors of the simulated results of Case-1 and Case-2 were 0.004 and 0.01, respectively.

| Numerical simulation
The parameters of the rainfall−runoff calculation (Table 1, Case-1), and those of the snowmelt calculation (Table 2, Case-1) were used to drive the model for numerical calculation of rainfall−runoff combined with snowmelt for the Bukuro River Basin. The period of the calculation was January 1, 1997 to December 31, 1999 (3 years). The calculated results ( Figure 15) illustrate that the calculated result appropriately reproduced the observed flow and that the two runoff curves fit well both in terms of flow scale and in the temporal process. The error was 0.005, which is significantly lower than the allowance value (0.03). This result indicates that highprecision simulation of rainfall−runoff combined with snowmelt can be achieved by using the calibrated parameters.

| Classification of parameter sensitivity
In the process of PSA using the OAT method, the initial value and the change range of each parameter was not snowmelt of the Bukuro River Basin randomly adopted. The value of each parameter was determined via field survey combined with referencing the relevant results of similar research. As the coefficient of permeability and infiltration velocity were determined by simple field infiltration experiments, the thickness of the first layer was determined by multi-point survey and the main parameters for snowmelt calculation were originally chosen by referencing related research (Aoki et al., 1999;Kondo et al., 1988;Niemeyer et al., 2016). The change range of each parameter was determined through error analysis comparing the simulated results and the observed flow with the different parameter values through numerous model calculations. The value change of each parameter was controlled within a reasonable range in the process of model calculation for PSA. For a few parameters, simulated results with non-obvious differences were obtained when an individual parameter value changed within a small range. To achieve simulated results with obvious differences, so as to evaluate the classification for parameter sensitivity, the change range of an individual parameter was adjusted, and the value exceeded its reasonable range, for example, the case of the thickness of the second layer was assumed to be 5.0 m (Case 3 in Figure 5), and the roughness of coefficient of the river channel (Case 2 in Figure 7) was assumed to be 0.30, etc. Parameter sensitivity for each parameter was evaluated by directly comparing the error of the peak flow between the simulated result and the observed flow when an obvious difference existed between the observed peak flow and the simulated result under a value change of a single parameter within the reasonable range. Meanwhile, error calibration (Equation (11)) was used to evaluate the error in the case of non-obvious differences between the simulated result and the observed value (particularly for peak flow). Based on the analysis results of parameter sensitivity as above mentioned, the classifications of parameter sensitivity have been separated and determined, with the results shown in Table 3.

| Limitation on applicability
The classification of each parameter's sensitivity (Table 3) were limited by the physical conditions of the underlying surface of the specified basin. In other words, the change range for each parameter used in the current study was merely applicable to the Bukuro River Basin or the basins whose topography and underlying surface conditions are similar to those of the Bukuro River Basin. In case of physical conditions of topography and the underlying Note: Values for a, b were 1.0-9.0; values for thickness of layer, permeability, infiltration velocity, and porosity were for the first layer in the soil vertical profile model.
surface are different from the Bukuro River Basin, such as a basin located in plain-or gully-type terrain. The results of the classification of PSA (Table 3) are inapplicable. However, the method for PSA used in this research, that is, OAT combined with error calibration (Equation (11)) can be referenced or adopted.

| Research deficiencies
For a physically distributed hydrological model, distributed physically based parameters must be assigned to each sub-catchment according to spatially separated results corresponding to the river channel networks. In the current research, the soil vertical profile model of the Bukuro River Basin was assumed as homogenous and a same set of values was adopted as the model parameters for each sub-catchment. When performing PSA for a large-scale basin which soil vertical profile model cannot be developed as a homogenous model, more efficient parameter calibration and sensitivity analysis methods must be explored and applied.

| CONCLUSIONS
In the current research, Morris' screen method (OAT) combined with error calibration was used to perform PSA for an independently developed distributed hydrological model for a small-scale mountainous basin, namely Japan's Bukuro River Basin. The main results are summarized as follows: 1. Classification of the parameter sensitivity of the main parameters influencing rainfall-runoff calculation was graded as thickness of the first layer of the soil vertical profile, coefficient of permeability, infiltration velocity, porosity, and coefficient of roughness. 2. Classification of the parameter sensitivity of the main parameters influencing snowmelt calculation was graded as coefficient of bulk transfer, temperature correction factor, and albedo of snow cover layer surface. 3. The value change of parameters of the first layer on the soil vertical profile model mainly impacts the model result for flow-discharge scale; the value change of parameters of the second layer on the soil vertical profile model has almost no impact on flowdischarge scale while causing a misplacement of peak flow.