The third-level catchment (S-catchment) represents the smallest element of a hydrological model, and the S-catchment division accuracy directly determines the accuracy of subsequent hydrologic analyses. To ensure that the catchment areas are consistent with land types, the division of S-catchments is based on topographic flow directions and fine-resolution land-use types. First, the D8 algorithm is applied to perform flow analysis based on an amended DEM to determine the main flow directions. The SL catchments are then roughly partitioned into direction-based (D-B) catchments based on these flow directions. An iterative S-catchment division algorithm is finally proposed to finely partition these D-B catchments based on their fine-resolution land-use types to produce S-catchments.

#### 3.3.1. Partitioning of D-B Catchments Based on Flow Directions

- (1)
DEM Amendment

The surface runoff processes occurring within an urban area are controlled not only by the natural topography but also by buildings, roads, and canals. In the method proposed, the original DEM is amended to account for the influences of urban underlying surfaces on runoff; this is performed by incorporating information about land objects that influence the overland flow routing into the DEM. The elevation values at DEM grid points containing land object information are then increased or decreased to alter their catchment capacities. It should be noted that to avoid the artificial depressions caused by DEM amendment, the building data are preprocessed to avoid the building with holes. Since the influences of roads and canals on runoff have already been accounted for during the partitioning of SL catchments, when amending the DEM, it is necessary only to account for the impacts of buildings in obstructing runoff. The procedure and equation for DEM amendment are shown in

Figure 7 and Equation (1).

where E

_{a} is the amended DEM elevation value, E

_{i} is the original DEM elevation value, and H is the corresponding building’s height.

- (2)
D-B Catchment Division Based on Flow Directions

Step 1: Flow analysis. The hydrologic analysis tool in ArcGIS 10.2 is used to fill sinks in the amended DEM. It should be noted that if there exist artificial depressions caused by DEM amendment, this paper will first judge whether it is a depression surrounded by buildings. If so, the continuous iteration will be applied to fill the depression until there is no depression caused by buildings.

The D8 algorithm is then used to calculate the DEM elevation differences in eight directions from the central grid (north, west, east, south, southeast, northeast, southwest, and northwest), and the descending direction with the greatest elevation difference is defined as the flow direction (purple line in

Figure 8). The basic D8 algorithm is probably the most popular method for automated drainage recognition and catchment area determinations [

14]. During the extraction of flow directions, the catchment threshold will affect the selection of trunk flows and tributaries. Based on multiple experiments, the optimal threshold was determined to be 1000 in this paper.

Noted that the threshold can affect the results of flow directions. The larger the threshold, the less tributary at the end of the stream is extracted, leading to the main flow direction is clearer, whereas the details of tributaries are missing. The threshold for the formation of flow in each region is different. Since this paper is mainly for plain urban areas, the difference in surface slope is small, the threshold should be as small as possible to describe the surface runoff in detail, but it is necessary to avoid excessive fragmentation. Considering the fine division of catchment is based on land use data, the threshold whose flow results are more consistent with the land use distribution is selected as the optimal threshold in this paper. Taking about 5.2 km

^{2} as the sample area, the thresholds of 500, 1000, 2000, 5000, and 8000 are used to generate the flow direction, and the number of main flow directions generated by each threshold is obtained. Then, the ratio of the number of flows to land use patches are calculated, as illustrated in

Figure 8. It can be found that as the threshold increases, this ratio gradually decreases. When the threshold is 1000, the ratio is the closest to 1, suggesting the flow results are more consistent with the land use distribution. Therefore, the threshold of 1000 in this paper is selected as the optimal threshold of the study area. It should be noted that the flows under different thresholds do not denote true geomorphological significance, so it is necessary to select the optimal threshold according to regional characteristics and data conditions.

Step 2: Extraction of catchment points. The flow accumulation of each grid is determined by calculating the number of grids upstream whose flows pass through this grid. The grid elements with high flow accumulation levels will be used to form catchment ridge lines (green lines in

Figure 9), and the intersections between catchment ridge lines are catchment points (blue points in

Figure 9).

Step 3: Preliminary division of catchment areas. Based on the catchment points obtained in Step 2, all the grid elements with the same flow direction are partitioned into the same independent catchment area. These catchment areas are the third-level D-B catchments.

#### 3.3.2. The Smallest Catchment Division Based on Land Type Data

To accurately determine urban catchments, an iterative fine division algorithm was proposed to partition D-B catchments based on high-precision land-use type data and catch basin data, thus obtaining the smallest catchments (S-catchments). It should be noted that due to that different land uses have different infiltration intensity, generally, the higher the precision of land uses data, the more accurate the catchment area will be. However, high precision land uses data are usually fragmented, which will lead to the difficulty of automatic determination of catchments in large-scale urban area. The proposed algorithm can address this issue. To guarantee the accuracy of the sub-catchments, the land use data used in this paper is derived from the Department of Natural Resources of Shandong Province, China. The data accuracy and quality conformed the standard specification DB37T 2761.2-2016 (Technical specification of geographic information public service platform, Part 2: Framework Data) issued by Shandong Bureau of Quality and Technical Supervision. The land use patch with area larger than 100 m^{2} can be obtained, which resolution is about 10m. The data used in this paper include lawn, forest, cropland, unused land, buildings, and plazas/roads.

- (1)
Calculation of Land-use patch Areas

The area of each land-use patch,

Area_{i} (

i-th patch), is calculated by superimposing the land-use type data onto the D-B catchment data (

Figure 10).

- (2)
Determination of a Single Catchment

The ratio between the area of each land-use patch and the total area of their D-B catchment,

AreaShare_{i}, is calculated as follows:

where

m is the total number of land-use patches and

Area is the total area of the D-B catchment in which the

i-th land-use patch is located.

The single catchment areas and mixed catchment areas are determined based on the AreaShare_{i} results. To prevent subcatchments from being overly partitioned and fragmented, which would affect subsequent hydrologic analyses, an iteration algorithm was proposed in this paper.

The threshold value δ for the land-use patch area is defined as follows: when

$AreaShar{e}_{i}\ge \delta $ and the number of catch basins in the patch is greater than 0, land-use patch

i is deemed a single subcatchment. For instance, as illustrated in

Figure 10, due to that the area share of patch A

$\text{}\ge \delta $, the land use patch A is determined as a single sub-catchment. While the area share of patch B

$\text{}\delta $, it is determined a mixed sub-catchment. Based on the characteristics of the study area and a large number of tests, the value of the threshold,

δ, was determined to be 0.12.

- (3)
Iterative algorithm for the determination of mixed catchment areas

When

AreaShare_{i} ≤

δ, land-use patch

i is topologically associated and merged with its adjacent land-use patches (

j) to form a mixed catchment area (Panel B in

Figure 10). The function for the topological association and merging of land-use patch

i with its adjacent patches may be expressed as:

where

n is the total number of patches adjacent to land-use patch

i,

AreaSNei_{ij} is the ratio between the area of patch

i and that of (adjacent) patch

j,

Direction_{ij} is the angle between the flow directions of patch

i and patch

j,

Rainwater_{ij} is the number of catch basins and catchment points between patch

i and its adjacent patch, and

SemNei_{ij} is the semantic similarity between patch

i and patch

j. The mathematical model for semantic similarity can be calculated based on Li et al. [

26] with the following expression:

where

L_{i} and

L_{j} are the land-use types of land-use patch

i and its adjacent patch

j, respectively; m is the number of land-use types that are semantically adjacent to

L_{i}, for the land use data used in this paper, m is 6; and Distance (

L_{i},

L_{j}) is the distance between the positions of

i and

j in the set of semantically adjacent land types. The semantic positions of the land types are based on the arrangement of permeating coefficient. Patches with similar permeating coefficients are considered more compatible and have priority to merge. The similarity is reflected by a value between 0 and 1, where 0 means not compatible at all and 1 means totally compatible. For instance, the semantic positions are plazas/roads, buildings, unused land, cropland, lawn and forest, their permeating coefficients are 0.011, 0.012, 0.06, 0.17, 0.20, and 0.50, respectively. According to Li et al. [

26], define the semantic distance between adjacent land uses as 1 unit, if L

_{i} and L

_{j} are unused land and plaza/roads, the Distance (

L_{i},

L_{j}) is 2 unit, so the

$SemNe{i}_{ij}$ equals 2/3 (1-2/6).

The merging function for land-use patch

i and its adjacent patch

j is

where

w_{1} and

w_{2} are the weights of each index, and the weights are given by the criteria importance through intercriteria correlation (CRITIC) method.

Next, all of the land use patches are iteratively processed, and a merging priority is used to determine the optimal adjacent patch j to merge with land-use patch i. This process continues until all S-catchments have areas greater than δ, or the number of catch basins in each individual land patch is greater than 0, or the flow directions between all adjacent patches are different. The end of this iterative process marks the completion of the S-catchment division process.

Finally, the flow directions of the S-catchments are used to add upstream/downstream relationships between the S-catchments, which will be used in subsequent hydrologic analyses (

Figure 11). This process helps to improve the catchment division accuracy, and the thresholds and iterative processes of this method also ensure that the catchment areas will not be so fragmented as to preclude their application in large-scale hydrological models.