Simulating the Hydrological Processes of a Meso-Scale Watershed on the Loess Plateau, China

: The Soil and Water Assessment Tool (SWAT) model is widely used to simulate watershed streamﬂow by integrating complex interactions between climate, geography, soil, vegetation, land use / land cover and other human activities. Although there have been many studies involving sensitivity analysis, uncertainty ﬁtting, and performance evaluation of SWAT model all over the world, identifying dominant parameters and conﬁrming actual hydrological processes still remain essential for studying the e ﬀ ect of climate and land use change on the hydrological regime in some water-limited regions. We used hydro-climate and spatial geographical data of a watershed with an area of 3919 km 2 , located on the Loess Plateau of China, to explore the suitable criterion to select parameters for running the model, and to elucidate the dominant ones that govern the hydrological processes for achieving the sound streamﬂow simulation. Our sensitivity analysis results showed that parameters not passing the sensitive check ( p -value < 0.05) could play a signiﬁcant role in hydrological simulation rather than only the parameters with p -value lower than 0.05, indicating that the common protocol is not appropriate for selecting parameters by sensitivity screening only. Superior performance of the rarely used parameter SOL_BD was likely caused by a combination of lateral and vertical movement of water in the loess soils due to the run-on inﬁltration process that occurred for meso-scale watershed monthly streamﬂow modeling, contrasting with traditionally held inﬁltration excessive overland ﬂow dominated runo ﬀ generation mechanisms that prevail on the Loess Plateau. Overall, the hydrological processes of meso-scale watershed in the region could be well simulated by the model though underestimates of monthly streamﬂow could occur. Simulated water balance results indicated that the evapotranspiration in the region was the main component leaving the watershed, accounting for 88.9% of annual precipitation. Surface runo ﬀ contributed to 63.2% of the streamﬂow, followed by lateral ﬂow (36.6%) and groundwater (0.2%). Our research highlights the importance for selecting more appropriate parameters for distributed hydrological models, which could help modelers to better comprehend the meso-scale watershed runo ﬀ generation mechanism of the Loess Plateau and provide policy makers robust tool for developing sustainable watershed management planning in water-limited regions.


Introduction
The physical process-based distributed parameters hydrological models that incorporate watershed heterogeneity and spatial distribution of climate, terrain, soil, vegetation, and land uses are widely used for watershed management planning, water resources allocation, flood protection, climate and land use change assessment, and environmental and pollution control evaluation [1][2][3][4][5][6]. However, applying hydrological models in the real world is often challenged by calibration, validation, and predictability issues [7][8][9][10]. During the calibrated process, parameter sensitivity and uncertainty analysis often become a hurdle, even though hydrological models are normally calibrated to find optimal parameters set with the optimum objective functions [11][12][13]. Unless there is a careful sensitivity and uncertainty analysis, overestimation or underestimation of hydrologic regimes can cause over-design of mitigation measures or insufficient preparation for potential situation [13,14].
The Soil and Water Assessment Tool (SWAT) model is one of the most widely used physical processes based hydrological models worldwide due to its capability to model the full spectrum of hydrological processes and its user-friendly interface [15][16][17][18]. Moreover, capability of SWAT-CUP program that comprises a semi-automated approach (SUFI-2) can help users conduct sensitivity and uncertainty analysis in a more rigorous way [10,19]. Most of the previous studies pointed out that less than ten parameters were sufficient to represent simple rainfall-runoff processes through models [20][21][22]. However, the result of acceptable simulation of streamflow in a watershed does not signify a correct performance of runoff generation processes [23]. Because of the large number of parameters in the comprehensive watershed model, parameterization and calibration of the model become complicated [10]. All parameters governing the hydrologic processes can be used for calibrated SWAT model, and many parameters have the impact on multiple processes [10]. However, correct parameterization cannot only make model calibration be faster, be more accurate, and have lower prediction uncertainty [10], but also describe the complete hydrological cycle [6,24].
Some of the frequently used calibrated parameters that have been considered as the most sensitive parameters in controlling hydrological processes are SCS runoff curve number (CN2), soil evaporation compensation factor (ESCO), available water capacity (SOL_AWC), groundwater delay (GW_DELAY), snowfall temperature (SFTMP), baseflow alpha factor (ALPHA_BF), and surface runoff lag time (SURLAG) [6,14,15,17,23,[25][26][27][28]. Most of the previous studies suggested that CN2 is the most sensitive parameter for simulating streamflow on the Loess Plateau as the infiltration excessive runoff generation mechanisms prevail for small watersheds in the region [3,20,28,29]. Nevertheless, recent studies have noted that rarely used parameters, such as moist bulk density (SOL_BD), may also be more sensitive in streamflow simulation, not only for the small watershed in humid region [30], but also for large-scale watersheds on the Loess Plateau. For example, Gong et al. (2017) identified that SOL_BD was the most sensitive parameter rather than CN2 for runoff simulation in Yanhe watershed (7725 km 2 ) [31]. Such a finding was attributed to the higher permeability of loess soil in the Loess Plateau by ecological restoration, the lower rainfall intensity since the 1990s, and the deep-seated tunnel systems that may reach the channels [32][33][34]. However, it still remains unclear if and in what way such rarely used parameters should be considered for meso-scale watersheds in the region.
Although there have been many studies involving sensitivity analysis, uncertainty fitting and performance evaluation of SWAT model all over the world [2,11,14,16,18,35], identifying dominant parameters, reducing the magnitude of uncertainties for streamflow simulation and confirming actual hydrological processes still remain essential for some distinct regions. Uncertainties in hydrological modeling come from input data, model structure, and parameters [2,21]. However, there are many key parameters in a watershed that describe the bio-geophysical characteristics and hydrological processes of the watershed, among which some parameters are hard to measure directly and can be gained only by empirical estimation and literature references, resulting in the increased uncertainties of the modeling system [11,17,21,36]. Several studies showed that using more informative data, modified model structures, or more parameter sets could reduce uncertainties [2,37,38]. In addition, the SUFI-2 method is capable of assessing the uncertainties from all sources for the SWAT model [21]. Moreover, the parameters cannot be identified easily due to the phenomenon of equifinality [13]. Therefore, a scientific understanding of model uncertainties effect on the model performance and capacity of model to capture hydrologic processes can help users apply the model in a more rigorous way.
Most notably, the implementation of a large-scale land-conversion program, e.g., converting from agriculturally used slope land into forest or/and grassland in the Loess Plateau region of China has raised wide concerns over the water yield decrease, water carrying capacity of forest coverage, and balance between afforestation and food production [39,40]. However, changes in hydrological processes and watershed balance components associated with land use and land cover change are integrated elements for assessing the ecosystem services and healthy status at watershed scale [20,29,41,42]. Clearly, it is the physical-processes-based distributed-parameter hydrological models that can achieve such objectives. However, the models require a considerable number of parameters for their performance evaluation, including parameterization and calibration. Here, we used long-term hydro-climate data combined with land use/land cover change data interpreted from remote-sensing images, soil distribution, and landform data of Xinshui watershed located on the Loess Plateau, one of the main tributaries in the middle reach of the Yellow River, China, to investigate the applicability of SWAT model, hydrological processes, and water balance dynamics. Our specific aims are to (1) explore the suitable criterion to select parameters for running the model by challenging general protocol of sensitivity screening; (2) elucidate the dominant parameters that govern the hydrological processes for achieving the sound streamflow simulation at meso-scale watersheds that are traditionally assumed infiltration excessive overland runoff generation mechanisms prevail; and (3) quantify the long-term water balances for developing sustainable watershed management plan in the region.

Study Area
The Xinshui watershed lies in the southeastern part of the Loess Plateau, Shanxi Province, Northern China. The river originated from the Lvliang Mountains and is the secondary contributor to the Yellow River. The drainage area is approximately 3919 km 2 (36 • 36 -36 • 57 N, 110 • 30 -111 • 27 E) ( Figure 1). Temperate continental monsoon climate dominates the region. The annual average temperature ranges from 7.9 to 10.7°C, and rainfall is from 289 to 778 mm, during the past 60 years. The watershed terrain is dominated by mountain ranges, hills, and valleys. The watershed elevation varies from 699 to 2010 m ( Figure 1). Fifty-five percent of the watershed is covered by Calcaric Cambisols, while other soil types include Calcic Luvisols, Calcaric Fluvisols, Eutric Leptosols, Rendzic Leptosols, Eutric Regosols, Calcaric Regosols, Haplic Luvisols, and Eutric Cambisols. The texture of these soils is sandy loam soil composed by the sand and silt, with minor clay. Land-use types and percentage of each land use within the watershed are 7.48% forestland, 26.20% shrubland, 37.81% grassland, 28.31% cropland, 0.18% residential areas and built-in land, and 0.02% water body. Grassland is the most dominant land-use type. The major agricultural crops are corn, soybean, wheat, potato, and sorghum. Major forest types are Quercus liaotungensis, Betula platyphylla Suk, Populus davidiana, Pinus tabulaeformis, Platycladus orientalis, and Robinia pseudoacacia Linn.

SWAT Structure
The SWAT model is a typical semi-distributed, process-based hydrologic model [43]. It has a strong physical foundation for runoff generation and can simulate water, sediments, and nutrients from individual sub-watersheds to its outlet in a basin [43]. During the construction of the SWAT project, the watershed separates to sub-basins by watershed delineation according to DEM data, and some Hydrological Response Units (HRUs) constitute one sub-basin. The HRU is unique and consists of one slope type, one soil type, and one land-use type [43]. The water-balance equation in SWAT is as follows: where SWt and SW0 means the final and the initial soil water content, respectively; t means the time; Pday means precipitation; Rsurf means surface runoff; ETa means the actual evapotranspiration; wseep means the water entering the vadose zone; and Rgw means the return flow.
The hydrologic cycle is climate driven [10]. The primary elements of hydrologic processes of the SWAT model are evapotranspiration, infiltration, surface runoff, return flow, lateral flow, tile drainage, water stored in the soil profile, and transmission losses [10]. The plant-growth model within the SWAT model simulates all land-cover types and is capable of distinguishing annual and perennial vegetation, in order to evaluate biomass production and transpiration [10]. Curve number method is used for predicting surface runoff [10]. Channel routing is estimated by the Muskingum method [15].
The water-balance equation for a closed drainage basin is as follows [44]: where P = precipitation (mm); ET = evapotranspiration (mm); R = water yield (mm); ∆S = Sj − Si = soil water storage change (mm); and ∆S= 0 at longer time intervals.
In this study, 1986-1999 climatic inputs, land-use data and soil data were used to rerun the calibrated model and combined with the water-balance equation, to identify the proportions of hydrological elements in the watershed.

SWAT Structure
The SWAT model is a typical semi-distributed, process-based hydrologic model [43]. It has a strong physical foundation for runoff generation and can simulate water, sediments, and nutrients from individual sub-watersheds to its outlet in a basin [43]. During the construction of the SWAT project, the watershed separates to sub-basins by watershed delineation according to DEM data, and some Hydrological Response Units (HRUs) constitute one sub-basin. The HRU is unique and consists of one slope type, one soil type, and one land-use type [43]. The water-balance equation in SWAT is as follows: where SW t and SW 0 means the final and the initial soil water content, respectively; t means the time; P day means precipitation; R surf means surface runoff; ET a means the actual evapotranspiration; w seep means the water entering the vadose zone; and R gw means the return flow. The hydrologic cycle is climate driven [10]. The primary elements of hydrologic processes of the SWAT model are evapotranspiration, infiltration, surface runoff, return flow, lateral flow, tile drainage, water stored in the soil profile, and transmission losses [10]. The plant-growth model within the SWAT model simulates all land-cover types and is capable of distinguishing annual and perennial vegetation, in order to evaluate biomass production and transpiration [10]. Curve number method is used for predicting surface runoff [10]. Channel routing is estimated by the Muskingum method [15].
The water-balance equation for a closed drainage basin is as follows [44]: where P = precipitation (mm); ET = evapotranspiration (mm); R = water yield (mm); ∆S = S j − S i = soil water storage change (mm); and ∆S= 0 at longer time intervals. In this study, 1986-1999 climatic inputs, land-use data and soil data were used to rerun the calibrated model and combined with the water-balance equation, to identify the proportions of hydrological elements in the watershed.

Data Preparation and Model Construction
The basic data required to build the SWAT model are weather data, digital elevation model (DEM) data, soil data, and land-use data [43]. The soil database and weather generator database need to be site-specific, and these data are introduced in detail below.
The DEM (Figure 2A) data with a resolution of 30 m × 30 m were downloaded from Geospatial Data Cloud site of Chinese Academy of Sciences, which was used for watershed delineation, stream creation, slope gradient, and channel slope extraction ( Figure 2B).
Land-use data in 1995 were obtained from Landsat5 TM satellite remote sensing image product and processed by ENVI software (The Environment for Visualizing Images, version 5.1, Exelis Visual Information Solutions, Herndon, VA USA), at a 30 m resolution. The accuracy of the classification was then verified, yielding 1995 Kappa coefficient of 0.771 [20]. Six categories were obtained by classifying land use: urban and built-up areas, water, forest, shrub, grassland, and cropland. These six categories were transformed into URML, WATR, FRST, RNGB, PAST, AGRL, respectively, in the SWAT project ( Figure 2C). Each land use had different surface vegetation types which influence runoff, evapotranspiration, and groundwater in the watershed. Before the large-scale Land Conversion Program (e.g., Grain for Green) was implemented since 1999, there had been no very significant land-use and land-cover changes in the region from 1986 to 1999, due to the progressive implementation of Three North Shelterbelt Programs started in 1978 [45]; therefore, we used this land-use data for our simulation, both for calibration and validation periods.
The soil database (1 km resolution) was obtained from the Harmonized World Soil Database (HWSD) and was cut out by the grid map of the study area [2]. Since SWAT used FAO soil grading, we did not convert the soil particle size content. We used SPAW software (Washington State University, WA, USA) to calculate available water content, saturated hydraulic conductivity, and matric bulk density for different layers of each soil type [43]. Other values in the soil database were calculated by the equations from SWAT documentation. Therefore, we reclassified our study soil types into 11 subtypes ( Figure 2D).
Climatic inputs for running SWAT included daily wind speed, solar radiation, relative humidity, maximum and minimum temperature, and precipitation. Because the ArcSWAT tool (USDA-ARS Grass-land, Soil and Water Research Laboratory and the Texas AgriLife Blackland Research Center in Temple, Texas, USA) allows users to load weather station locations into the current project and assign weather data to the sub-watersheds, these data were then transformed into appropriate text format of SWAT project. We used three meteorological stations (from Xi County, Ji County, and Houma City) and six rainfall stations (from Daning, Majiayao, Sanduo, Xiali, Guxian, and Jiaokou) in the project ( Figure 1). Sub-watershed climate data were furnished by the station nearest to the centroid of the sub-basin that achieves the point weather data to watershed weather data. The Penman-Monteith method was applied to calculate potential evapotranspiration.
Our SWAT model generated 31 sub-basins and 845 HRUs, by setting a threshold of 2%, 5%, and 10% for land use/land cover, slope, and soil types, respectively. The aim of setting thresholds is to ignore unneeded large number of HRUs [43]. Water 2020, 12, x 6 of 18

Sensitivity Analysis
The first step of model calibration and validation is sensitivity analysis. Its primary objective is to identify parameters that have a great influence on model outputs (e.g., streamflow). We selected the global sensitivity analysis method to determine sensitive parameters, using the SUFI-2 method of SWAT-CUP [19]. Global sensitivity analysis distinguishes the sensitive rank of whole considered parameters related to streamflow by Latin hypercube regression analysis [19].
The t-stat and p-value, two statistical measurements, could assess the sensitive rank of each parameter. The t-stat represents a range of sensitivity, while the p-value identifies the significance of sensitivity. Higher absolute value of t-stat and lower value of p-value (<0.05) indicate that a parameter is sensitive [17].
Twenty-four input parameters for running SWAT model were selected for sensitivity analysis according to the literature review from watersheds with similar characteristic to Xinshui watershed of the Loess Plateau (Table 1) [10,20,28,30,35].

Sensitivity Analysis
The first step of model calibration and validation is sensitivity analysis. Its primary objective is to identify parameters that have a great influence on model outputs (e.g., streamflow). We selected the global sensitivity analysis method to determine sensitive parameters, using the SUFI-2 method of SWAT-CUP [19]. Global sensitivity analysis distinguishes the sensitive rank of whole considered parameters related to streamflow by Latin hypercube regression analysis [19].
The t-stat and p-value, two statistical measurements, could assess the sensitive rank of each parameter. The t-stat represents a range of sensitivity, while the p-value identifies the significance of sensitivity. Higher absolute value of t-stat and lower value of p-value (<0.05) indicate that a parameter is sensitive [17].
Twenty-four input parameters for running SWAT model were selected for sensitivity analysis according to the literature review from watersheds with similar characteristic to Xinshui watershed of the Loess Plateau (Table 1) [10,20,28,30,35].

Uncertainty Analysis
The values of P-factor and R-factor in SUFI-2 algorithm of SWAT-CUP, were applied to assess the SWAT model simulation uncertainties [19]. P-factor is the percentage of measured data enveloped by model simulation results, 95PPU. R-factor is 95PPU envelope thickness [19]. All the uncertainty factors (including driving variables, parameters, observed data, conceptual model, etc.) in the whole process of simulation can be accounted by parameter uncertainty [15]. Therefore, for streamflow, if the P-factor is higher than 70% and the R-factor is lower than 1, that means the simulation results are acceptable [19].

Calibration and Validation
Our SWAT monthly streamflow was calibrated in comparison with the monthly observed streamflow at Daning Hydrological Station of Xinshui watershed. The model was implemented for 17 years, from 1983 to 1999, by applying the first 3 years for warm-up period, allowing the hydrologic cycle to be fully represented in the model. Calibration data were from 1986 to 1992, whereas validation data were from 1993 to 1999. Land-use data in 1995 were applied to ensure high accuracy of land use during the whole calibrated process [46].

Model Performance Evaluation
We also used the NS (Nash-Sutcliffe coefficient) [47], R 2 (coefficient of determination), PBIAS (percent bias), and RSR (standardize the RMSE using observation standard deviation) to determine the model performance [19].
where Q is streamflow; m is observed data; s is simulated data; the bar is represented for average; and j is the j th measured or simulated data. The performance evaluation criteria (NS, R 2 ≥ 0.5, PBIAS < ±25% and RSR < 0.7) means satisfactory model calibration [48].

Sensitivity Analysis
Twenty-four hydrological parameters were tested for sensitivity analysis via the SUFI-2 method for streamflow simulation in the Xinshui River Watershed. The absolute value of t-stat ranged from 0.05 to 14.26, and the p-value ranged from 0 to 0.96 (Table 2). According to the parameter sensitivity evaluation criterion, 11 parameters passing the sensitivity check (p-value < 0.05) were SOL_BD, RCHRG_DP, GWQMN, ESCO, CN2, SLSUBBSN, OV_N, GW_DELAY, ALPHA_BNK, CH_N2, and SOL_K. Hydrologically, the SOL_BD parameter controls the horizontal and vertical movement of soil water; the RCHRG_DP and GWQMN parameters regulate the occurrence of groundwater; the ESCO parameter controls the contribution of soil water to evaporation; the CN2, SLSUBBSN, and OV_N parameters control the formation of surface runoff; the GW_DELAY and ALPHA_BNK parameters regulate the retention time of groundwater or baseflow. The CH_N2 parameter plays an important role in regard to channel flow; the SOL_K parameter controls the movement of soil water [19].

Model Calibration and Validation
Model calibration was first conducted by using the 11 parameters passing the sensitivity check (Table 2). However, model performance criteria R 2 , NS, PBIAS, and RSR were 0.7, 0.66, 32.3% and 0.58, respectively, indicating the simulation was unsatisfactory and erroneous. The SWAT model could satisfactorily simulate the monthly streamflow for the watershed, until the first nineteen parameters were all used for calibration (Table 3). Model performance indictors NS, RSR, and PBIAS were 0.83, 0.41, and 20% for the calibration period and 0.88, 0.35, and 17.6% for the validation period, respectively (Table 4), achieving the satisfactory level. The linear regression coefficients in calibration and validation less than 1 indicated that monthly streamflow was likely underestimated in the Xinshui River Watershed (Figure 3). V indicates that the parameter value will be replaced by the given value; R indicates that the parameter value is multiplied by (1 plus the given value).

Uncertainty Analysis
The SUFI-2 method could capture the monthly streamflow regime during the whole calibration period. First, most of the observed monthly streamflow data were within the 95PPU bracket (77% in

Uncertainty Analysis
The SUFI-2 method could capture the monthly streamflow regime during the whole calibration period. First, most of the observed monthly streamflow data were within the 95PPU bracket (77% in calibration period, 87% in validation period). Second, the R-factor equaled 0.85 and 0.9, respectively, for calibration and validation periods. The simulated annual averaged streamflow during the calibration period was 2.14 m 3 /s, which was 20% lower than the observed data (2.68 m 3 /s). Similarly, the observed and simulated annual averaged streamflow during the validation period were 2.69 and 2.21 m 3 /s, respectively, with −17.6% simulation errors. Both results showed that some streamflow were underestimated, especially in the winter and during the flood periods in summer in 1990, 1991, and 1992 ( Figure 4). However, they remained within the boundaries of 95PPU of simulated streamflow, indicating that the uncertainty of the SWAT model was within the allowable range. calibration period was 2.14 m 3 /s, which was 20% lower than the observed data (2.68 m 3 /s). Similarly, the observed and simulated annual averaged streamflow during the validation period were 2.69 and 2.21 m 3 /s, respectively, with −17.6% simulation errors. Both results showed that some streamflow were underestimated, especially in the winter and during the flood periods in summer in 1990, 1991, and 1992 ( Figure 4). However, they remained within the boundaries of 95PPU of simulated streamflow, indicating that the uncertainty of the SWAT model was within the allowable range.

Water Balances
For a watershed, precipitation, evapotranspiration, surface runoff, lateral flow, and groundwater are the most important water balance components. The long-term mean annual water balance components of the Xinshui River Watershed during 1986 to 1999 are shown in Figure 5. The distribution proportion of precipitation for water balance components was obtained by SWAT-Check tool (USDA-ARS Grass-land, Soil and Water Research Laboratory and the Texas AgriLife Blackland Research Center in Temple, Texas, USA). The incoming precipitation flux was balanced by evaporation of soils and plants (88.9%), followed by water yield (3.8%) and other components stored in soils and/or lost from channel by transmission (7.3%) ( Table 5).
indicated that the amount of precipitation was lost mainly by evapotranspiration and surface runoff. Due to the dominated infiltration excessive overland runoff generation mechanism in the region, surface runoff contributed 63.2% to the total streamflow, followed by lateral flow (36.6%) and groundwater (0.2%). The baseflow, which is constituted by lateral and groundwater flow together, contributed to 36.8% of the total streamflow ( Figure 6). The seasonal trend in potential evapotranspiration (PET) showed an increase trend from January to June, and then a decrease trend from July to December ( Figure 6). The trend of seasonal precipitation (P), evapotranspiration (ET), water yield (WYLD), and surface runoff (SURQ) showed similar patterns during the simulation period.    Similarly, the seasonal water balance from 1986 to 1999 is shown in Figure 6. The results also indicated that the amount of precipitation was lost mainly by evapotranspiration and surface runoff. Due to the dominated infiltration excessive overland runoff generation mechanism in the region, surface runoff contributed 63.2% to the total streamflow, followed by lateral flow (36.6%) and groundwater (0.2%). The baseflow, which is constituted by lateral and groundwater flow together, contributed to 36.8% of the total streamflow ( Figure 6). The seasonal trend in potential evapotranspiration (PET) showed an increase trend from January to June, and then a decrease trend from July to December ( Figure 6). The trend of seasonal precipitation (P), evapotranspiration (ET), water yield (WYLD), and surface runoff (SURQ) showed similar patterns during the simulation period. Relative error (%) equals to (P -SURQ -PERC -LATQ -ET − ∆S) / P x 100%, representing the water balance errors of the watershed.

Discussion
Sensitivity analysis usually provides a parameter selection scheme for running a physical-processes-based distributed-parameter hydrological model. All the sensitive parameters with higher absolute value of t-state and lower p-value than 0.05 are used for calibrating the SWAT model, while other parameters are treated as insensitive and not used [10,17]. However, the model performance could achieve the satisfactory level until 19 parameters calibrated (Table 3 and Table 4). This result is consistent with previous studies indicating that the good results in prediction would be achieved based on a large number of parameter sets rather than one "best" parameter set [37,49,50]. Clearly, our study suggested that the insensitive parameters could also play a role in hydrological simulation according to other model-performance evaluation criteria (R 2 , NS, PBIAS, P-factor, and R-factor).

Discussion
Sensitivity analysis usually provides a parameter selection scheme for running a physical-processes-based distributed-parameter hydrological model. All the sensitive parameters with higher absolute value of t-state and lower p-value than 0.05 are used for calibrating the SWAT model, while other parameters are treated as insensitive and not used [10,17]. However, the model performance could achieve the satisfactory level until 19 parameters calibrated (Tables 3 and 4). This result is consistent with previous studies indicating that the good results in prediction would be achieved based on a large number of parameter sets rather than one "best" parameter set [37,49,50]. Clearly, our study suggested that the insensitive parameters could also play a role in hydrological simulation according to other model-performance evaluation criteria (R 2 , NS, PBIAS, P-factor, and R-factor).
In our study, the sensitivity rank of CN2 was lower than SOL_BD, RCHRG_DP, GWQMN and ESCO, suggesting that more complex hydrological processes may occur in the watershed. As the most sensitive parameter among all 19 parameters for calibrating SWAT model in the Loess Plateau, SOL_BD could significantly improve the performance of the model for simulating monthly streamflow. The performance evaluation criteria R 2 , NS, and RSR, without considering SOL_BD, were 0.75, 0.72, and 0.53 for calibration and 0.82, 0.81, and 0.43 for validation, accordingly. When this parameter was involved, R 2 , NS, and RSR were improved to 0.85, 0.83, and 0.41 for calibration and 0.89, 0.88, and 0.34 for validation, respectively (Table 4). Our sensitivity analysis results showed that SOL_BD was the most sensitive parameter, a finding that differs from the previous conclusion that CN2 is the most sensitive parameter for streamflow simulation supported by the long hold infiltration excessive overland flow runoff generation mechanism in the Loess region [3,20,28,51]. Superior performance of SOL_BD was likely caused by a combination of lateral and vertical movement of water in the loess soils, as it was also observed in the Yanhe River Basin of the Loess Plateau [31]. Hortonian overland flow (run-on infiltration process) is the main infiltration mechanism at the hillslope scale, resulting from high rainfall intensity but low duration [8]. As the Hortonian flow moves to the river channel, the movement of water via the Dunne flow (saturation excess) determines the hydrological cycle at the watershed scale [33]. Furthermore, the main soil type in the Xinshui River Watershed is Calcaric Cambisols, so low-intensity rainfall would infiltrate the soil [33]. Parameter SOL_BD is closely related with soil infiltration in that the initial infiltration capacity of soil decreases against the increase in bulk density, and the attenuation rate of infiltration capacity increase with the increase of bulk density [51,52]. These results indicated that the lateral and vertical redistribution movement of water in the soil layers into the watershed outlet had a significant effect on hydrological processes, such as later flow, groundwater, and evapotranspiration. Our study suggested that, in addition to the frequently used parameters, the rarely used parameter SOL_BD should be used for model calibration, to improve the model performance and help modelers understand the main hydrological process in the meso-scale watershed of the Loess Plateau and other similar regions in the world.
Understanding uncertainty can improve on the deeper comprehending of both hydrologic cycle and model predictions [37]. The larger the P-factor value, the greater the contribution of parameter uncertainty to the uncertainty of the simulation [2]. The lower P-factor means that input or structure may dominate model simulation uncertainty [2,38]. Our uncertainty analysis results during the validation period (P-factor = 87%, R-factor = 0.9) were better than those during the calibration period (P-factor = 77%, R-factor = 0.85), implying that parameter uncertainty contributed the most. These results suggested that using appropriate parameters makes it possible to represent model simulation uncertainty and reduce input and/or structure uncertainty. In addition, there were some underestimations of peak flows in summer, particularly during the calibration period. Monthly streamflow in this temperate monsoon climate region is usually contributed to by several thunderstorms characterized with quick response. However, simulated monthly streamflow could be tamped down by the time step of 30 days, leading to the underestimation for summer peak flows [8,33,53]. Moreover, a simplified approach in the SWAT model is employed to reckon soil temperature and confirm the frozen and thawed state of soils, and both also increase a small part of simulation uncertainty [2,40]. Uncertainties of the SWAT model were within acceptable limits, meaning that all source uncertainties were captured by parameter uncertainty in our study [15,17,21,23,54]. However, this does not neglect the significance of the spatiotemporal variability of precipitation [38,53] and model structures [2,55,56] during the whole calibration processes.
The hydrological processes of the Xinshui River Watershed were satisfactorily simulated by SWAT model in our study. The average simulated ET from 1986 to 1999 was 447.3 mm (Figure 5), and the average measurement ET (which equals to averaged P minus averaged streamflow) was 457.8 mm. Compared with measured ET, the relative error of simulated ET was only −2.3%. Our simulated results showed that the 88.9% of P returned to atmosphere through ET, and a similar result was found in Yanhe River Basin [57]. Only 3.8% of P was in form of streamflow through the watershed outlet, and the rest of P was stored in soils and/or lost from channel by transmission (7.3%). The Xinshui River Watershed, located in the Loess Plateau, under temperate continental monsoon climate, which is known to have over 80% of P, appears from May to October, and the result is consistent with [8]. SURQ and WYLD had a similar trend with P in seasonal scale because the dominated hydrologic processes were overland flow and water yield generated from thunder storms [8]. However, LATQ showed a smoother upward trend due to lateral redistribution movement of water [33]. The mean annual P (503.1 mm) was equal to the sum of ET (447.3 mm), R (18.98 mm) and (38.2 mm) was 504.17 mm, the absolute error of water balance was −1.38 mm, and the relative error was 0.3%. The reason for simulated water balance error may be the poor simulation of snowmelt process of the SWAT model [56]. Although the relative error of watershed water balance from 1986 to 1999 ranged from −2.26% to 2.17% (lower than 5%) (Table 5), the hydrologic processes simulated by SWAT model can used for further analysis. In addition, the simulated total annual water yield for 1986-1999 was 18.98 mm, which is constituted by surface runoff (63.2%) and baseflow (36.8%). The result is in accordance with Kang et al.'s (2019) finding that the average annual baseflow index in Xinshui River Watershed is 0.368 [58]. Hence, knowing the watershed hydrologic processes is helpful for mangers to make decisions on the water resources' utilization.

Conclusions
In this study, the values of performance evaluation criteria (R 2 , NS, PBIAS, RSR, P-factor, and R-factor) during the calibration period were 0.85, 0.83, 0.41, 20.0, 77% and 0.85, while during validation period, those values were 0.89, 0.88, 0.34, 17.6, 87% and 0.9, showing that streamflow and actual hydrological processes of the Xinshui River Watershed (~4000 km 2 ) could be modeled by using SWAT with the SUFI-2 method. However, unlike most published works, we found that more parameters used could significantly improve the model performance and sensitivity analysis may not be applicable universally for determining parameter sets. We also concluded that the inclusion of SOL_BD, the most sensitive parameter controlling the horizontal and vertical movement of soil water, could improve the model performance greatly, challenging the long-held ideas that the CN curve number is the most sensitive parameter relevant to infiltration excessive overland flow generation mechanisms in the region. It was anticipated that model uncertainty could be captured by parameters uncertainty. In addition, the watershed water-balance results suggested that most of the precipitation was evaporated into the atmosphere; therefore, only a small percentage of precipitation formed streamflow. Surface runoff is the primary component of the streamflow within the study area; the significance of baseflow in hydrological processes cannot be neglected. The calibrated SWAT model in Xinshui River Watershed can be used to explore the effect of climate and land-use-change scenarios on hydrologic cycle.