Study on Applicability of Distributed Hydrological Model under Di ﬀ erent Terrain Conditions

: This research aimed to study the applicability and limitations of a distributed hydrological model under discontinuous steep topography and hydrogeological conditions. Based on GIS spatial analysis, typical cases of steep and gentle terrains were selected to construct the distributed hydrological model framework of the research areas (Qinhuangdao and Zhuanghe City, China). The observed runo ﬀ was used to test the applicability of the model in di ﬀ erent terrain watersheds and to analyze the versatility of the model structure and the relevant parameters of the core modules. The results show that: in the process of using a distributed hydrological model to build models for di ﬀ erent regions, problems such as a discontinuous dislocation of the empty area and poor connectivity of the water system will appear in the process of sub-basin division of a steep terrain. By determining the optimal threshold, selecting the best node, discontinuous dislocation, void fusion and other methods, we put forward the corresponding solutions to the problems in the division process and constructed the research area’s distributed hydrological model. The rainfall–runo ﬀ process in the study area was simulated accordingly, and the SUFI2 algorithm was used to calibrate the relevant parameters in the model. The relative error (Re), correlation coe ﬃ cient ( R 2) and Nash–Sutcli ﬀ e e ﬃ ciency ( NSE ), which meet the runo ﬀ accuracy in the study area, were obtained. The model veriﬁcation results show that the NSE of steep terrain is 0.90, and R 2 is 0.98; the NSE of gentle terrain is 0.91, and R 2 is 0.984: the simulation values ﬁt the measured values well, which makes the calibrated model suitable for both steep and gentle terrains. The results can provide a reference for the construction of a distributed hydrological model in watersheds with di ﬀ erent terrain.


Introduction
The continuous improvement of computer skills and the development of science and technology has concurrently led to the rapid development of digital elevation model (DEM), geographic information system (GIS) and remote sensing (RS) technology [1]. Data acquisition methods are becoming increasingly diverse with a continuous improvement in accuracy and coverage expansion [2]. Hamdi et al. [3] used Light Detection and Ranging (LiDAR) data and Google Earth Engine (GEE) to demonstrate a method on how to create accurate high-resolution riparian buffer maps. Sanzana et al. [4] used GIS to construct a Geo-PUMMA toolbox, which can be used to represent the terrain in distributed hydrological modeling in both urban and peri-urban scales. A distributed hydrological model is an ideal model with a physical basis [5] and has been rapidly developed. The distributed hydrologic model has been greatly improved in considering factors such as precipitation and the spatial variability of the underlying surface, describing the spatial and temporal variation of the hydrological cycle process in the basin. [6,7]. As powerful tools, the distributed hydrological models have gradually gained attention in modern hydrological simulation studies, through accurately and scientifically revealing hydrological variation regularity. Distributed hydrological models that are widely used currently include System Hydrological European (SHE) [8], Soil and Water Assessment Tool (SWAT) [9], Variable Infiltration Capacity (VIC) [10], and the Distributed Hydrological Soil Vegetation Model (DHSVM) [11]. Distributed hydrological model studies also pay more attention to model scale by thoroughly combining the functions of GIS spatial analysis and visibility [12,13]. At present, the distributed hydrological model mainly achieves the spatial dispersion of research objects through grids and sub-basin [14,15]. For applicability studies of distributed hydrological models, Fontaine et al. [16] improved the algorithm of snow melting process and extended the application range of the SWAT model from agricultural area to non-agricultural high-altitude mountainous area. Cambien et al. [17] tested the suitability of the Soil and Water Assessment Tool (SWAT) to simulate the dynamics of two different pesticides in the data-scarce Guayas River basin. Senent-Aparicio et al. [18] combined the Soil and Water Assessment Tool (SWAT) model and the chloride mass balance (CMB) method to improve the modeling of streamflow in high-permeability bedrock basins receiving inter-basin groundwater flow (IGF). Liu et al. [19] optimized the application of snow melting modules in the SWAT Model for the alpine regions of Northern China. Stone et al. [20] applied the modified SWAT model to Missouri catchment to simulate the effect of elevated CO 2 concentration on plant growth. Wang et al. [21] combined the advantages of lumped model-reservoir model and distributed hydrologic model-(SWAT) model to efficiently establish a simulation of the water cycle in a karst area. Li et al. [22] used the SWAT model to study the eco-hydrological characteristics of a plateau salt marsh wetland to help managers protect and restore water balance by using water salt dynamic model. Huang et al. [23] constructed a distributed hydrological model of different watersheds based on GIS and used the basic equation of motion wave theory combined with the topography (hydrogeological conditions) of the watersheds. They tested the practicability of the model in each basin by the numerical simulation of observed discharge. However, the distributed hydrological model had some limitations. The plain river network is flat, which led to the bottleneck problem of water system extraction and sub-basin division; a large scale area is challenging to accurately grasp the spatial scope of the study area and accurately determine the drainage outlet. To solve this problem, Turcotte et al. [24] proposed to use the digital river lake network (drln) as input to modify the flow direction of the model and network calculation. Liu et al. [25] proposed a watershed subdivision method using stem-branch topological codification. The subdivision method was proposed based on the automatic recognition of the outlets of the basin and the river network's fusion with variable catchment area thresholds.
Under the condition of steep terrain, there has been few studies on the applicability of a distributed hydrological model. Under a steep terrain, the application of a distributed hydrologic model have had some limitations: (1) the complex topography and geomorphology leads to great differences in the spatial distribution of water systems; (2) the existence of microtopography area with sharp elevation change results in the difference of catchment area between the basins without confluence connection; and (3) when the spatial scale of the underlying surface varies greatly, the steep terrain will generate discontinuous dislocation voids, watershed deformation and spatial dislocation. This makes the generalization of the distributed model for watershed analysis under disruption and steep terrain conditions challenging. In response to this practical technical problem, this study selected two research areas with large terrain differences as the research object and constructed a distributed hydrological model suitable for different regions. Finally, the rainfall runoff process is simulated, and the parameters are calibrated, and the applicability of the simulation results is evaluated.

Study Area
The typical research area for a steep terrain selected in this study is Qinhuangdao City, Hebei Province. The area has good primitive features, large terrain variation, obvious seasonal rainfall, large exposed rainfall-runoff, and small annual perennial flow. This area is a typical study area with no water and waterlogging during the rainy season. The area chosen for a gentle terrain in this study is Zhuanghe City in the south of Dalian. This area is a coastal city with relatively small topographic relief, high in the north and low in the south. The surface runoff flows slowly into the sea, which is a typical frequent urban waterlogging area. The underlying natural conditions of the two places are shown in Figure 1 and Table 1.

Study Area
The typical research area for a steep terrain selected in this study is Qinhuangdao City, Hebei Province. The area has good primitive features, large terrain variation, obvious seasonal rainfall, large exposed rainfall-runoff, and small annual perennial flow. This area is a typical study area with no water and waterlogging during the rainy season. The area chosen for a gentle terrain in this study is Zhuanghe City in the south of Dalian. This area is a coastal city with relatively small topographic relief, high in the north and low in the south. The surface runoff flows slowly into the sea, which is a typical frequent urban waterlogging area. The underlying natural conditions of the two places are shown in Figure 1 and Table 1.  A DEM is an important data source for watershed topography and feature recognition. Adequate surface information can be derived by using a digital elevation model, which includes the slope between cells in the watershed grid unit, the relationship between the aspect ratio and slope direction. The DEM data used in this study were obtained from the global digital elevation models (GDEM) 30 m resolution digital elevation data from the Geospatial Data Cloud (http://www.gscloud.cn/). We selected N39-E123, N39-E123, N40-E122, N40-E123 as the images of Zhuanghe City and N39-E118, N39-E119, N40-E118, N40-E119 as the images of Qinhuangdao city. After splicing and cutting the DEM data, the 3D function module of the ArcGIS software was used to extract the slope from the Slope function, and then the slope was reclassified (Figure 2), and the raster graphics attribute the table generated after the reclassification was exported. Finally, the  A DEM is an important data source for watershed topography and feature recognition. Adequate surface information can be derived by using a digital elevation model, which includes the slope between cells in the watershed grid unit, the relationship between the aspect ratio and slope direction. The DEM data used in this study were obtained from the global digital elevation models (GDEM) 30 m resolution digital elevation data from the Geospatial Data Cloud (http://www.gscloud.cn/). We selected N39-E123, N39-E123, N40-E122, N40-E123 as the images of Zhuanghe City and N39-E118, N39-E119, N40-E118, N40-E119 as the images of Qinhuangdao city. After splicing and cutting the DEM data, the 3D function module of the ArcGIS software was used to extract the slope from the Slope function, and then the slope was reclassified (Figure 2), and the raster graphics attribute the table generated after Sustainability 2020, 12, 9684 4 of 18 the reclassification was exported. Finally, the attribute data were statistically analyzed to calculate the slope composition of the study area in Table 2.
Sustainability 2020, 12, x FOR PEER REVIEW 4 of 18 attribute data were statistically analyzed to calculate the slope composition of the study area in Table  2.  The slope range of the study area was divided into six categories of ≤5°, 5°-10°, 10°-15°, 15°-30°, 30°-50°, 50°-70°, namely flat, gentle, oblique, steep, urgent, dangerous, respectively [26]. It can be found that the average slope of the study area 1 is 12.01° and the maximum slope reaches 68.14°, while the average slope of the study area 2 is 6.55° and the maximum slope reaches 54.43°. From the analysis of the large slope range in the study area, it can be found that Qinhuangdao City is mainly concentrated in the slopes between 15° and 30°, with widely distributed mountainous areas with large slopes and obvious elevation changes. It is a typical steep terrain, where erosion and block movement are more severe, and soil erosion is more serious. Zhuanghe City is mainly concentrated in slopes between 0° and 10°. The terrain gradually decreases from north to south. The terrain change is relatively gentle, the water body movement is stable, and the soil erosion is weak, which is a typical gentle terrain.

Land Use
For the simulation, the study used a land use distribution map of 2017 with a scale of 1:100,000. The data were downloaded from the website of the Resource and Environment Data Cloud Platform (http://www.resdc.cn/). The input data of the SWAT model include the land use distribution map and land use index table. The SWAT model converts the vector into a raster layer when loading the land use type and uses ArcGIS system tools to reclassify and project the land use type ( Figure 3 [26]. It can be found that the average slope of the study area 1 is 12.01 • and the maximum slope reaches 68.14 • , while the average slope of the study area 2 is 6.55 • and the maximum slope reaches 54.43 • . From the analysis of the large slope range in the study area, it can be found that Qinhuangdao City is mainly concentrated in the slopes between 15 • and 30 • , with widely distributed mountainous areas with large slopes and obvious elevation changes. It is a typical steep terrain, where erosion and block movement are more severe, and soil erosion is more serious. Zhuanghe City is mainly concentrated in slopes between 0 • and 10 • . The terrain gradually decreases from north to south. The terrain change is relatively gentle, the water body movement is stable, and the soil erosion is weak, which is a typical gentle terrain.

Land Use
For the simulation, the study used a land use distribution map of 2017 with a scale of 1:100,000. The data were downloaded from the website of the Resource and Environment Data Cloud Platform (http://www.resdc.cn/). The input data of the SWAT model include the land use distribution map and land use index table. The SWAT model converts the vector into a raster layer when loading the land use type and uses ArcGIS system tools to reclassify and project the land use type ( Figure 3). Then, it establishes a land use type index table by connecting the value of the land use raster map with the existing classification in the SWAT land cover/plant database. Finally, the construction of the land use database was completed, and it was divided into six basic land use types, including cultivated land, forest land, grassland, water area, construction land and unused land as shown in Table 3. with the existing classification in the SWAT land cover/plant database. Finally, the construction of the land use database was completed, and it was divided into six basic land use types, including cultivated land, forest land, grassland, water area, construction land and unused land as shown in Table 3.  The analysis of the land use types in the two research areas shows that the current land use distribution in the two research areas is similar. The proportion of various land use types differs by less than 10% and has good comparability. The most significant difference between the two is the different terrain: one is a steep terrain, and the other is a gentle terrain.

Soil
Soil texture and structural characteristics are the key factors affecting the hydrological cycle, which directly affect the redistribution of rainfall in the watershed. The soil type data collected in this study are the Harmonized World Soil Database version 1.1 (HWSD V1.1) (http://www.fao.org/soilsportal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/en/), which is in a grid format with a resolution of 1000 m. HWSD uses the USDA standard, which is consistent with the soil database standard included in the SWAT model, therefore it can be used directly to establish the SWAT model database without the need to convert particle size standards. With the help of SPAW [27], we calculated the available water capacity of the soil layer (SOL_AWC), saturated hydraulic conductivity (SOL_K), and the moist bulk density (SOL_BD). Soils with similar runoff generation capacity under the same external conditions such as the precipitation and surface water were classified into the same hydrological group (HYDGRP). This was done with reference to the Natural Resources Conservation Service (NRCS), and according to the permeability properties of soils. The hydrological group was then divided into four classes: A, B, C and D. Finally, a soil index table was established to call the attribute data with raster data, so the construction of the soil database was completed (see Tables SA1 and SA2 in Supplementary Materials). The results of the soil  The analysis of the land use types in the two research areas shows that the current land use distribution in the two research areas is similar. The proportion of various land use types differs by less than 10% and has good comparability. The most significant difference between the two is the different terrain: one is a steep terrain, and the other is a gentle terrain.

Soil
Soil texture and structural characteristics are the key factors affecting the hydrological cycle, which directly affect the redistribution of rainfall in the watershed.
The soil type data collected in this study are the Harmonized World Soil Database version 1.1 (HWSD V1.1) (http://www.fao.org/soils--portal/soil--survey/soil--maps--and--databases/harmonized--world--soil--database--v12/en/), which is in a grid format with a resolution of 1000 m. HWSD uses the USDA standard, which is consistent with the soil database standard included in the SWAT model, therefore it can be used directly to establish the SWAT model database without the need to convert particle size standards. With the help of SPAW [27], we calculated the available water capacity of the soil layer (SOL_AWC), saturated hydraulic conductivity (SOL_K), and the moist bulk density (SOL_BD). Soils with similar runoff generation capacity under the same external conditions such as the precipitation and surface water were classified into the same hydrological group (HYDGRP). This was done with reference to the Natural Resources Conservation Service (NRCS), and according to the permeability properties of soils. The hydrological group was then divided into four classes: A, B, C and D. Finally, a soil index table was established to call the attribute data with raster data, so the construction of the soil database was completed (see Tables S1 and S2 in Supplementary Materials). The results of the soil reclassification in the study area constructed by the above methods are shown in Figure 4. The proportion of different types of soil and HYDGRP are shown in Table 4.
Sustainability 2020, 12, x FOR PEER REVIEW 6 of 18 reclassification in the study area constructed by the above methods are shown in Figure 4. The proportion of different types of soil and HYDGRP are shown in Table 4.  The main soil types in Qinhuangdao are cambisols, fluvisols, luvisols and regosols. On the other hand, Zhuanghe is dominated by luvisols and phaeozems soil types. Both places have a large proportion of luvisols soils. (Table 4). Luvisols are a kind of soil with visible clay movement under moist soil moisture, which is fully leached by lime. In comparison with the HYDGRP of the two study areas, Qinhuangdao has a large portion of soils in class A and B soils. Zhuanghe comprises primarily class C and B soils. Class A soils have high permeability, mainly composed of sand and gravel, with poor water holding capacity. Most of the incoming water flows deep into the soil with low runoff generation capacity. Class B is a kind of soil with medium permeability, with medium drainage and   The main soil types in Qinhuangdao are cambisols, fluvisols, luvisols and regosols. On the other hand, Zhuanghe is dominated by luvisols and phaeozems soil types. Both places have a large proportion of luvisols soils. (Table 4). Luvisols are a kind of soil with visible clay movement under moist soil moisture, which is fully leached by lime. In comparison with the HYDGRP of the two study areas, Qinhuangdao has a large portion of soils in class A and B soils. Zhuanghe comprises primarily class C and B soils. Class A soils have high permeability, mainly composed of sand and gravel, with poor water holding capacity. Most of the incoming water flows deep into the soil with low runoff generation capacity. Class B is a kind of soil with medium permeability, with medium drainage and water diversion capacity and overall runoff generation capacity. Class C is the soil with low permeability, and there is a layer that hinders the water infiltration, the infiltration rate and hydraulic conductivity are low, and the runoff generation capacity is strong. It can be concluded that Zhuanghe soil has a stronger runoff producing capacity than Qinhuangdao, and the underlying surface has a better interception effect on precipitation.

Meteorology
The construction of this weather database is based on the CMADS (The China Meteorological Assimilation Driving Datasets for the SWAT model) dataset (http://www.cmads.org/) [28]. The dataset has been formatted and revised according to the SWAT model input data format and can be directly imported into the SWAT model for use [29]. The spatial location search function was used to screen meteorological stations around the study area. The stations were recorded based on the search results. We found the five essential elements of the meteorological data by the name of the recorded site, namely relative humidity, precipitation, solar radiation, temperature and wind, and read all the files in the SWAT model to complete the construction of the meteorological database. We water diversion capacity and overall runoff generation capacity. Class C is the soil with low permeability, and there is a layer that hinders the water infiltration, the infiltration rate and hydraulic conductivity are low, and the runoff generation capacity is strong. It can be concluded that Zhuanghe soil has a stronger runoff producing capacity than Qinhuangdao, and the underlying surface has a better interception effect on precipitation.

Meteorology
The construction of this weather database is based on the CMADS (The China Meteorological Assimilation Driving Datasets for the SWAT model) dataset (http://www.cmads.org/) [28]. The dataset has been formatted and revised according to the SWAT model input data format and can be directly imported into the SWAT model for use [29]. The spatial location search function was used to screen meteorological stations around the study area. The stations were recorded based on the search results. We found the five essential elements of the meteorological data by the name of the recorded site, namely relative humidity, precipitation, solar radiation, temperature and wind, and read all the files in the SWAT model to complete the construction of the meteorological database. We The comparative analysis shows that the average annual rainfall in Qinhuangdao City is 600.7 mm, and the rainfall is concentrated in July-August. The average annual rainfall in Zhuanghe City is 684.5 mm, and it is concentrated in June-September. It can be seen that the rainfall in the two The comparative analysis shows that the average annual rainfall in Qinhuangdao City is 600.7 mm, and the rainfall is concentrated in July-August. The average annual rainfall in Zhuanghe City is 684.5 mm, and it is concentrated in June-September. It can be seen that the rainfall in the two research areas has a similar trend. There are inter-annual fluctuations in rainfall in both cases, which are unevenly distributed in time, and the inter-annual distribution is generally balanced. The uneven distribution of water resources in time has led to drought in some seasons in some years, and its effective water collection and storage are of considerable significance to the balanced distribution of water resources in the region.

Modeling Technique Methodology
The hydrological analysis functions of ArcGIS 10.2 were used to analyze and process the DEM data. The optimal threshold was determined, and then the sub-basin was divided by SWAT. After dividing the whole watershed into multiple catchment units, the model calls the land use data, soil data, and meteorology data in the basic database through the index table. It then analyzes the results by performing simulation operations in sub-modules. The sensitivity analysis tool of the SWAT model was used to analyze the sensitivity of parameters. Then, the optimal value of parameters was determined by repeated operations using the SUFI2 algorithm [30] in SWAT-CUP software. Finally, the Re, R 2 and Nash Sutcliffe efficiency (NSE) were obtained, which were consistent with the runoff accuracy of the study area. The process is shown in Figure 6.

Hydrological Analysis
Accurate extraction of river networks is a key step in the division of the sub-basin. The main work includes: 1. Depression-free DEM generation: due to the resolution of the DEM data and the real terrain, there are always some depressions in the image data. To calculate the flow direction correctly, it is necessary to fill the original DEM first to get a non-depression DEM [31]. Extracting water flow direction: this paper used D8 algorithm [32] in GIS to calculate the water flow direction. The algorithm assumed that there are only 8 possible flow directions in a single grid, and the 8 directions are coded with different codes: 1-due east, 2-southeast, 4-due south, 8southwest, 16-due west, 32-northwest, 64-due north, 128-northeast ( Figure 7). The flow direction of a grid cell is defined as the cell with the steepest slope among the 8 grid points adjacent to it. The slope calculation is shown in Equation (1):

Hydrological Analysis
Accurate extraction of river networks is a key step in the division of the sub-basin. The main work includes:

1.
Depression-free DEM generation: due to the resolution of the DEM data and the real terrain, there are always some depressions in the image data. To calculate the flow direction correctly, it is necessary to fill the original DEM first to get a non-depression DEM [31]. Extracting water flow direction: this paper used D8 algorithm [32] in GIS to calculate the water flow direction. The algorithm assumed that there are only 8 possible flow directions in a single grid, and the 8 directions are coded with different codes: 1-due east, 2-southeast, 4-due south, 8-southwest, 16-due west, 32-northwest, 64-due north, 128-northeast ( Figure 7). The flow direction of a grid cell is defined as the cell with the steepest slope among the 8 grid points adjacent to it. The slope calculation is shown in Equation (1): where h i is the height of the grid unit, h j is the height of the adjacent grid unit, and D is the distance between the center points of the two grids. If the two grid cells are adjacent in the horizontal or vertical direction, D is the cell length; if it is in the diagonal direction, D is 2× cell length. According to this principle, the flow direction of all grid points in DEM was calculated one by one to obtain the flow direction raster map.
there are always some depressions in the image data. To calculate the flow direction correctly, it is necessary to fill the original DEM first to get a non-depression DEM [31]. Extracting water flow direction: this paper used D8 algorithm [32] in GIS to calculate the water flow direction. The algorithm assumed that there are only 8 possible flow directions in a single grid, and the 8 directions are coded with different codes: 1-due east, 2-southeast, 4-due south, 8southwest, 16-due west, 32-northwest, 64-due north, 128-northeast (Figure 7).

2.
Cumulative consolidation: assuming that DEM has one unit of water at each grid unit represented by a regular grid, the cumulative amount of confluence at each grid unit can be obtained by calculating the water flow through each grid unit according to the determined flow direction data of the grid unit. Grid units with high confluence accumulation values generally correspond to the rivers, and grid units with cumulative accumulation values of 0 generally correspond to watersheds.

3.
River network extraction: at present, the surface runoff model is mainly used to extract river networks. When the cumulative amount of confluence reaches a certain value, surface runoff will be formed. All the cumulative amount of confluence grid above the critical value will form a potential flow path. A river network is the network composed of these paths. The complex rivers with lines are amended to single rivers, the ring-shaped river networks are amended to dendritic drainage, and the rivers with more complex bends are amended to facilitate subsequent node selection and sub-basin division.

4.
Node generation: a node is a point where water flows out of an area. This point is usually the lowest point along the boundary of the basin, that is, the point where the accumulated maximum of the internal flow in the catchment area. The river nodes simulated in the model not only reflects the consistency of geographical data, but also provides a guarantee for the subsequent hydrological and geomorphic simulation process analysis. In this paper, the optimal node was selected by combining the distribution of nodes under different thresholds, and eliminating the useless node information, and then the reservoir point information was inputted into the river network node.

Model Calibration and Verification
Model calibration and validation are critical in SWAT model simulation as they determine the accuracy and reliability of the simulation results. Parameter calibration is to find the parameter that can make the simulated value and the observed value the most consistent. In this paper, the SUFI2 algorithm in the SWAT-CUP software was selected to determine the optimal value of parameters through repeated operations, and the best value of parameters was brought back to the model for simulation verification through the manual parameter adjustment function of the SWAT model [33]. The applicability of the distributed hydrological model is mainly determined by the following three indicators: relative error (Re), Nash Sutcliffe efficiency (NSE) [34] and correlation coefficient (R 2 ). These three indicators can better reflect the correlation between the simulated value and the measured value from different aspects. The three criteria are computed as follows in Equations (2)-(4): Sustainability 2020, 12, 9684 10 of 18 where W P is the simulated value; W o is the measured value; Q obs,i is the measured runoff (m 3 /s); Q sim,i is the simulated runoff (m 3 /s); Q oa is the measured runoff average (m 3 /s); Q sa is the average value of simulated runoff (m 3 /s); n is the number of measured data. When Re > 0, this means that the simulated value is greater than the measured value; when Re = 0, it means that the simulated value is consistent with the measured value; when Re < 0, it means that the simulated value is less than the measured value. R 2 is used to evaluate the degree of agreement between the simulated value curve and the measured value curve. The value range is 0 < R 2 < 1. The larger the R 2 value, the better the matching effect. When NSE = 1, this means that the simulation value is completely credible; when 0 < NSE < 1, the greater the value, the higher the credibility of the simulation value; when NSE < 0, it means that the credibility of using the simulated runoff value is less than that of directly using the average value of the measured runoff value. The closer the values of NSE and R 2 to 1, the closer the simulated values are to the measured values.

Threshold Determination
In the process of sub-basin generation, it is crucial to determine the threshold value of the minimum catchment area because the threshold value does not only affect the generation of the runoff path but also indirectly affects the accuracy of the water volume of the drainage basin. Therefore, it is necessary to determine a threshold method of a small catchment area that can accurately describe the study area's scope. At present, the scholars at home and abroad have explored a variety of methods to calculate the threshold value, including size-frequency method [35], river network density method [36], moderate index method [37], etc. Among them, the river network density method proves to be the most straightforward and widely used method, which is why it was selected to determine the minimum catchment area in the study area. A series of simulated river networks was generated by selecting different catchment area thresholds, as shown in Table 5. In this paper, 10 thresholds with a catchment area ranging from 1000 to 10,000 ha were set to generate a grid river network through the use of hydrology-processing tools in ArcGIS 10.2. As shown in Table 5, the extraction of different thresholds had a high impact on the river length and network density in the different study areas. An increase in the threshold values results in a sparser extracted river network, while a reduction in threshold value results in a denser extracted river network with numerous pseudo channels. Due to the manual setting of the river network threshold, the proximity between the 10 river network thresholds set above and the real river network water system is unknown. Therefore, this paper selected river network density as the quantitative analysis index for river network extraction in the study area and attempted to analyze the relationship between the river network threshold and river network density. Figure 8 demonstrates the power function relationship between the river network threshold and the river network density under a different regional environment, obtained using the formulas in Equations (5) and (6): where y 1 is Qinhuangdao drainage density, km/km 2 ; y 2 is Zhaunghe drainage density, km/km 2 ; x is the threshold value.
where 1 y is Qinhuangdao drainage density, km/km 2 ; 2 y is Zhaunghe drainage density, km/km 2 ; x is the threshold value.  Figure 9 demonstrates the power function relationship between the threshold and second derivative of drainage density under the different regional environment obtained using the formulas in Equations (7) and (8) Based on the analysis illustrated in Figure 10, it was revealed that the threshold value was 6000 ha, also, the simulated water system was more consistent with the actual water system. By comparing the two research areas, it was observed that the main factor affecting the characteristics of the river  Figure 9 demonstrates the power function relationship between the threshold and second derivative of drainage density under the different regional environment obtained using the formulas in Equations (7) and (8): where y 1 is the Qinhuangdao second derivative of drainage density, km/km 2 ; y 1 is the Zhaunghe second derivative of drainage density, km/km 2 ; and x is the threshold value.
where 1 y is Qinhuangdao drainage density, km/km 2 ; 2 y is Zhaunghe drainage density, km/km 2 ; x is the threshold value.
(a) (b) Figure 8. Relationship between the threshold and drainage density in Qinhuangdao (a) and Zhuanghe (b). Figure 9 demonstrates the power function relationship between the threshold and second derivative of drainage density under the different regional environment obtained using the formulas in Equations (7) and (8) Based on the analysis illustrated in Figure 10, it was revealed that the threshold value was 6000 ha, also, the simulated water system was more consistent with the actual water system. By comparing the two research areas, it was observed that the main factor affecting the characteristics of the river water system was the topography. Therefore, the distribution of water systems in the two study areas Based on the analysis illustrated in Figure 10, it was revealed that the threshold value was 6000 ha, also, the simulated water system was more consistent with the actual water system. By comparing the two research areas, it was observed that the main factor affecting the characteristics of the river water system was the topography. Therefore, the distribution of water systems in the two study areas obviously has different characteristics which are affected by terrain. Qinhuangdao city's catchment flow direction was multifarious; the water system was well developed, and the spatial distribution was not balanced. Due to the gentle terrain change, the number of nodes selected was correspondingly high. Zhuanghe City, on the other hand had a radial water system with regular confluence from northwest to southeast, and its spatial distribution was relatively balanced.
Sustainability 2020, 12, x FOR PEER REVIEW 12 of 18 was not balanced. Due to the gentle terrain change, the number of nodes selected was correspondingly high. Zhuanghe City, on the other hand had a radial water system with regular confluence from northwest to southeast, and its spatial distribution was relatively balanced.

Sub-Basin Division
A reasonable threshold value was used in the SWAT model to delineate the sub-watersheds. In the process of sub-basin division by the distributed hydrological model, a discontinuous dislocation gap appeared in the zoning process of Qinhuangdao (steep terrain). The discontinuous dislocation gap in this paper refers to the area not included in the sub-basin's catchment area in the generation process of sub-basin. The simulation results revealed that the blank area of Qinhuangdao mainly appeared near sub-basins 39, 40, 41, mainly due to the significant topographic changes in the region. The boundary line of the sub-basin surrounding the blank area was dragged to enable the sub-basin of the blank area merge with the surrounding area. Because the blank area proportion is small when the river network sub-basin is generated, and there will be a minor error caused to the sub-basin division.
Finally, the Qinhuangdao basin and Zhuanghe Basin were divided into 42 and 36 sub-basin, respectively ( Figure 11). Compared with the sub-basin division results of the two study areas, it was revealed that the division results of study area 2 were regular and continuous, and the degree of coincidence was high. Nonetheless, the results of the sub-basin division revealed discontinuous catchment areas in study areas 1, 4, 36, 40 and 41. The basin boundaries were distorted making it unable to completely match the geographical boundary of the actual features, hence, leading to the occurrence of the automatic extraction error of river network. The main reasons behind this are that the terrain of Qinhuangdao is uneven, the terrain changes abruptly and unevenly, there are abrupt changes, and showing step fault. Due to the complexity of underlying surface conditions and the uneven distribution of hydrometeorological elements, the runoff of Qinhuangdao basin shows a nonlinear state.

Sub-Basin Division
A reasonable threshold value was used in the SWAT model to delineate the sub-watersheds. In the process of sub-basin division by the distributed hydrological model, a discontinuous dislocation gap appeared in the zoning process of Qinhuangdao (steep terrain). The discontinuous dislocation gap in this paper refers to the area not included in the sub-basin's catchment area in the generation process of sub-basin. The simulation results revealed that the blank area of Qinhuangdao mainly appeared near sub-basins 39, 40, 41, mainly due to the significant topographic changes in the region. The boundary line of the sub-basin surrounding the blank area was dragged to enable the sub-basin of the blank area merge with the surrounding area. Because the blank area proportion is small when the river network sub-basin is generated, and there will be a minor error caused to the sub-basin division.
Finally, the Qinhuangdao basin and Zhuanghe Basin were divided into 42 and 36 sub-basin, respectively ( Figure 11). Compared with the sub-basin division results of the two study areas, it was revealed that the division results of study area 2 were regular and continuous, and the degree of coincidence was high. Nonetheless, the results of the sub-basin division revealed discontinuous catchment areas in study areas 1, 4, 36, 40 and 41. The basin boundaries were distorted making it unable to completely match the geographical boundary of the actual features, hence, leading to the occurrence of the automatic extraction error of river network. The main reasons behind this are that the terrain of Qinhuangdao is uneven, the terrain changes abruptly and unevenly, there are abrupt changes, and showing step fault. Due to the complexity of underlying surface conditions and the uneven distribution of hydrometeorological elements, the runoff of Qinhuangdao basin shows a non-linear state. Sustainability 2020, 12, x FOR PEER REVIEW 13 of 18 Figure 11. Sub-basin division map in the study area.

Simulation Calculus
After dividing the whole basin into several catchment units, the model utilized the daily measured precipitation data, soil composition data, land use type data and climate type data in a basic database through the index table, using this information to carry out simulation calculation and result analysis in modules.

Parameter Sensitivity Analysis and Model Calibration
In this paper, the sensitivity analysis tool built in the model was used for the parameter sensitivity analysis. This tool mainly adopts the global advantages of Latin hypercube sampling, which makes the parameter analysis more comprehensive. Secondly, the one factor at a time strategy method was used to ensure the accuracy of parameter analysis [38]. According to the data collected, 2000-2001 was set as the warm-up period, 2002-2009 was set as the rate period, and the verification period was 2010-2014. In this study, a total of five iterations and more than 2000 simulations were used to conduct the parameter calibration and model verification. After determining the coefficient, the parameter range and final calibration value of runoff simulation were obtained. The results of the parameter adjustment are shown in Table 6.

Simulation Calculus
After dividing the whole basin into several catchment units, the model utilized the daily measured precipitation data, soil composition data, land use type data and climate type data in a basic database through the index table, using this information to carry out simulation calculation and result analysis in modules.

Parameter Sensitivity Analysis and Model Calibration
In this paper, the sensitivity analysis tool built in the model was used for the parameter sensitivity analysis. This tool mainly adopts the global advantages of Latin hypercube sampling, which makes the parameter analysis more comprehensive. Secondly, the one factor at a time strategy method was used to ensure the accuracy of parameter analysis [38]. According to the data collected, 2000-2001 was set as the warm-up period, 2002-2009 was set as the rate period, and the verification period was 2010-2014. In this study, a total of five iterations and more than 2000 simulations were used to conduct the parameter calibration and model verification. After determining the coefficient, the parameter range and final calibration value of runoff simulation were obtained. The results of the parameter adjustment are shown in Table 6.
According to the parameter sensitivity analysis, CN2 and GW_DELAY, GW_REVAP were the main influencing factors in both study areas. The difference is that the sensitivity of SOL_AWC is stronger among the factors affecting the hydrological status of Zhuanghe city. This means that soil moisture has a greater impact on the runoff yield and concentration in Zhuanghe basin. However, the ALPHA_BF (base flow regression coefficient) was more sensitive to the hydrological distribution of Qinhuangdao City, which means that groundwater has an important impact on runoff yield and concentration in Qinhuangdao city.

Analysis and Evaluation of Simulation Results
According to the evaluation method of simulation results, the above conclusions were evaluated. When R 2 > 0.6, NSE > 0.5 and |Re| < 25% can be considered as the simulation results of the model. As shown in Figure 12, the simulation results of the regular and validation period of the model in the two study areas were generally consistent with the measured runoff, and all indicators within the standard range. When comparing the simulated data with the measured values, |Re| are all controlled within 5%, which indicates that the data credibility is high. It can be seen from Figures 13 and 14 that after the model calibration was completed, the simulation results of the monthly value in the validation period of the model were in good agreement with the actual runoff. R 2 ≥ 0.98, NSE > 0.9, and the simulation accuracy reaching the class A standard, indicating that the model has good applicability. The evaluation table of simulation results was shown in Table 7. In comparison with various indicators, it can be found that the simulation results of Zhuanghe City are better than that of Qinhuangdao city. Through the analysis of the runoff simulation value and measured value, we can see that there was a large uncertainty between the simulated value and the measured value, mainly concentrated from June to September. This was largely due to the precipitation, which mostly affects the runoff change, and the precipitation primarily occurs from June to September every year, reaching the peak in August, where the impact difference of topography on catchment can be amplified. Therefore, in general, the model is relatively accurate for the runoff simulation value in the dry season, but not fast enough for the peak value change in the wet season. The simulation effect of the distributed hydrological model is generally better on a gradually changing terrain than on a steep terrain. After adjusting the model parameters, the simulation results in both research areas can reach a better range during the validation period.

Conclusions
This study takes Qinhuangdao and Zhuanghe city as research units to compare and analyze the simulation process and the results of a continuous distributed hydrological model under the underlying surface conditions of a typical steep and gentle terrain. This study concludes that the simulation results of Zhuanghe City (gentle terrain) are better than that of Qinhuangdao City (steep terrain). With this study focusing on the problems of steep terrain, the process of constructing a

Conclusions
This study takes Qinhuangdao and Zhuanghe city as research units to compare and analyze the simulation process and the results of a continuous distributed hydrological model under the underlying surface conditions of a typical steep and gentle terrain. This study concludes that the simulation results of Zhuanghe City (gentle terrain) are better than that of Qinhuangdao City (steep terrain). With this study focusing on the problems of steep terrain, the process of constructing a hydrological database consisting of the application of slope, land use, soil, climate and other data were compared with that of the gentle terrain. It was revealed that these factors influence the runoff of river systems. In the process of the sub-basin delineation of steep terrain, there were some problems like the discontinuous dislocation of empty area and poor connectivity of the water system. These were mainly due to the complexity of underlying surface conditions and the uneven distribution of hydro-meteorological elements. To obtain an improved sub-basin division of steep terrain, the catchment area threshold was checked and calculated repeatedly, a reasonable threshold was set for confluence accumulation, and a sub-basin fusion processing on discontinuously dislocated areas was performed. Finally, in order to maintain the accuracy of the simulation results and reduce the uncertainty, the simulation results were evaluated and the parameter sensitivity analysis was carried out. The verification results showed that R 2 ≥ 0.98, NSE > 0.9, the simulation accuracy was high, and the simulation and measured value both had good line fittings. The results can provide reference for the construction of a distributed hydrological model for different terrain basins.