A Method for Dynamical Sub-Watershed Delimitating by No-Fill Digital Elevation Model and Defined Precipitation: A Case Study of Wuhan, China

Abstract

Watershed delimitation is very important in flood control management. The traditional sub-watersheds delimitated by a filling digital elevation model (DEM) may change the real sink area, such that it may not be the best choice in studies sensitive to sub-watershed storage. This paper proposes a dynamical watershed delimitation method using a no-fill DEM and precipitation. It considers a closed sink area containing cells that fully flow into a large special cell, which can flow out when its water level is “higher than outlet”. We took Wuhan City as a study area and defined the precipitation in return periods of 1, 5, 20, or 100 years to derive the sub-watersheds. It is found that, in the four delimitations, the ratio of isolated basic units which could not flow outside were 27%, 9%, 5%, and 1%, respectively, as the precipitation increased. The results show that the provided method satisfies the assumption that the sink area might overflow with increased precipitation. The sub-watershed delimitated by the proposed method has higher correlation with the distribution of waterlogging points than those delimitated according to the D8 algorithm. These findings indicate that the proposed method can derive reasonable sub-watershed delimitation and that it may be helpful in the practice of urban flood control management.

1. Introduction

China has recently suffered serious flooding and waterlogging threats. Waterlogging is often caused by local rainfall that far exceeds an area’s drainage capacity, resulting in a certain degree of surface rainwater. Among the current 654 cities in China, 641 have been exposed to frequent flooding, and more than one hundred cities face flooding and waterlogging risks to different degrees [1]. Urban flooding has a significant impact on the health of the entire city and its residents [2]. Flood disasters may cause public transportation, water, and electricity to be blocked, which seriously disrupt the economic order. They may also cause direct damage to the property of residents, such as houses, cars, and warehouses storing goods being submerged. Furthermore, for those who need to travel through the waterlogged areas, there exist potential risks, such as being trapped in deep waters, falling into an open rainwater well, or getting an electric shock, which may pose a direct threat to their health and even to their lives.

Watershed delimitation is very important in the subject of water or sedimentary transport in disaster and environmental fields [3,4], and basin morphometric parameters largely control a catchment’s hydrology response [5]. In flood control management studies, different watershed division schemes can produce different simulation quantity results in the distribution of rainwater and sediment at the catchment scale [6,7,8]. In urban flooding carrying capacity estimation studies, the sub-watershed is the basic unit used to estimate the storage ability, the receiving upstream water ability, the potential of drainage to main outlets, and the emergence of pumping. These studies are sensitive to the original storage ability of the sub-watersheds.

Digital elevation models (DEMs) are the main data in watershed delimitation. The traditional watershed delimitation method is comprised of the single flow direction (SFD) and multiple flow direction (MFD) algorithms [9,10], both of which need to obtain continuous flow directions in the depression areas. As errors may result from the DEM production process, such as in the survey, smoothing, and transferring (the data type changes from float to integer) steps [11], most of these algorithms deal with construct depressions by filling in the DEM. Martz and Garbrecht (1992) elaborated the Jenson and Domingue (1988) [12] and Martz and Jong (1988) [13] of treating sinks in an effective manner by “filling” each depression in the DEM to the elevation of the lowest overflow point out of the sink, and the difference between the two methods was that the Jenson and Domingue modified flow directions from sink bottom to outlet to obtain continuous flow path by assuming that sink was primarily data error or artifacts [11]. Anold (2010) reviewed the historical progress of removing depression algorithms in [14], and the “priority-flood” theory has been widely adopted to determine the flow outlets of closed sink areas [15,16,17,18]. However, these “filling” DEMs will change the elevation of cells in depression areas (such as lakes, ponds, and reservoirs or wicket pits) [19].

In urbanized catchments, low-lying and flat areas may be the main terrain features. Whether the exact drainage topology can be extracted from a filled DEM has attracted attention in the literature. For example, Kenny et al. and Matthews et al. (2005) described that using verified hydrology stream data to burn in the DEM could indeed improve the watershed delimitation [20]; however, they also noted that the stream burning method still yielded error in predictable area within a wetland, particularly at complex stream junctions [20,21]. Arnold (2010) figured out that finding depressions which received flows from several upstream catchments but only passed to a single catchment may be significant in determining the routing path of a catchment [14]. As stereo remote sensing and Light Detection and Range (LiDAR) images have been widely used in the study of land type classification [22,23,24] and change detection [25,26], more and more refined digital terrain models can describe the real structure of the underlying surface. For example, the classical DTMs produced by the ALOS World 3D (AW3D) and the TanDEM-X methods can portray high resolution terrain features [27]. Therefore, some studies have begun to face the problem that depressions are a real type of terrain feature.

In [28], the suggestion that adopting an unfilled DEM to address the depressions can keep the original DEM data unchanged was reviewed. In 2004, Chou et al. adopted the preference ranking organization method for the enrichment evaluations (PROMETHEE) technique to deal with depressions based on an unfilled DEM dataset [29]. However, even the authors acknowledged that the evaluation matrix was too complex. In 2010, Arnold introduced a new “fill and overflow” algorithm to calculate flow accumulation using unfilled DEM, which transformed the depression problem into identifying and finding the outlets of depressions, and then calculated the network topology by tree-ordering the depressions to route the area [14]. In 2015, Jongmin and Yeong proposed the maximum depth tracing algorithm (MDTA), which used “saddle” to mark sub-depression flow outlets and “true outlet” to mark the final outlet of a large compound depression to delimitate watersheds based on unfilled DEM [30]. All of the former algorithms considered a depression to only pass down its area to one catchment. However, in some cases, the depression area might be able to pass to more than one of its neighboring catchments; for example, if there are several suitable outlets on the edge.

The essential process of distributing the flow direction in a sink area is inundation analysis, which can be grouped into “source flood” and “non-source flood” analysis [31]. “Source flood” analyzes the inundation with external water, which is determined by the inflowing water and the terrain features, while “non-source flood” analyzes the inundation according to a certain elevation, which is considered by flooded cells below a defined water level. Most watershed delimitation algorithms have utilized “non-source flood” analysis. For example, all fill-DEM algorithms based on the “priority-flood” concept derive the flow path in a depression by tracing the path from the minimum outlet to the center, considering all cells under the water level of the outlet to be overflowed, such that the flow path can reach the outlet [15,16,17,18]. The existing no-fill DEM algorithms, such as those in [14,29,30], consider each depression to have potential outlets in a neighborhood of the depression, which would directly consist as a compound depression and, so, belong to the “non-source flood” analysis methods.

The physical sense of “source flood” agrees with the idea of convergence. In a certain rainfall scenario, whether the neighborhood of a depression could consist of a new compound depression is determined by the rainwater. Only those connected neighborhood depressions can form a larger compound depression when the inside water level reaching the outlet achieves a certain level.

Therefore, in this paper, we propose a watershed delimitation method through the concept of multiple flow directions and “source flood” analysis, where the sink area can flow out when its water level is high enough to reach the outlet. The sub-watersheds are delimitated based on the unfilled DEM and the defined precipitation. The sub-watershed delimitated by this method can retain the real terrain characteristics, with low-lying and flat lands left unchanged. Therefore, the rainwater storage ability of the sub-watershed is kept as same as that in the original DEM dataset. This character indicates that the proposed method can be used in the analysis of rainwater and sedimentary carrying capacity [32], which is sensitive to the original sub-watershed storage ability.

2. Materials and Methods

2.1. Study Area and Materials

Wuhan is located in the east of Jianghan Plain. There are 13 districts in Wuhan, among which Jiang’an, Jianghan, Qiaokou, Hanyang, Wuchang, Hongshan, and Qingshan Districts are the seven downtown districts. The six new towns are the Dongxihu, Caidian, Jiangxia, Huangpi, Xinzhou, and Hannan Districts. We chose the seven downtown districts in Wuhan as the study area, as shown in Figure 1.

The average elevation in the study area is about 23.3 m, and the main range of elevation is between 18–874 m. The city has a humid subtropical monsoon climate with abundant rainfall; its annual precipitation is between 1150–1450 mm, with 40% of the precipitation concentrated from June to August every year [33]. Huang et al. calculated its average cumulative rainfall threshold over 12 h as 196 mm, according to the waterlogging and rainfall observation data from 2011–2016 [34]. Hong et al. analyzed the cumulative rainfall in 12 h in the 10-year and 20-year return periods to be 168.1 mm and 205.6 mm, respectively [35]. There are 166 lakes in the study area, 38 of which are located in the centers of the cities. The Yangtze River meets its largest branch (the Han River) in the center of this area, which divides Wuhan into three parts: Hankou, Wuchang, and Hanyang. In the range of Wuhan, the common rainy season water level of the Yangtze River is 23 m, and the historical water level record once reached 28–29 m.

Due to its low-lying land elevation and unique geological and hydrological features, Wuhan has the potential of water overflowing from inner lakes and drainage depression from the crossed river, which makes this a typical city threatened by the risk of flood and waterlogging disasters. For example, during 30 June 2016 to 6 July 2016, more than 200 waterlogging events occurred in the downtown area of Wuhan, which triggered concerns in society and the Chinese government [36]. Therefore, taking Wuhan as an example to study the sub-watershed delimitation is meaningful in the practice of flood control management.

In this study, the main material we used can be grouped into three types: the DEM dataset, the precipitation data, and the waterlogging points.

• DEM dataset. The DEM dataset is the main material used to delimitate watersheds during the process of calculating flow direction and identifying sink areas. According to the fluvial landscape contrast results of STRM30, ASTER, AW3D30, and TanDEM-X in paper [27], Boulton et al. (2018) found that AW3D30 generated the best derivative river network and was highly accurate, even in areas of steep topography and high relief, compared to the others. Therefore, we adopted the AW3D30 dataset for watershed delimitation in this research. The AW3D30 dataset was retrieved from the website of the Japan Aerospace Exploration Agency (JAXA) (available from http://www.eorc.jaxa.jp/ALOS/en/aw3d30/index.htm).
• Precipitation. Precipitation data is used to analyze the rainwater which could be collected over a certain rainfall event. It is a very important input element, which is used to derive the merged basic units and, finally, the watershed delimitation according to the provided method. In the regional hydrology study field, the precipitation in a certain area can be described by a statistical curve model. In China, the extreme value distribution curve, the negative exponential distribution curve, and the Pearson-III type distribution are the three empirical frequency distribution curve models which have been used for city design rainfall intensity frequency adjustments [37]. In this paper, we quoted four kinds of 12 h rainfall identities to be the input precipitation; that is, in the 1-year, 5-year, 20-year, and 100-year return periods, according to the results from [35], which were calculated by Guoping Hong et al. (2018) [35] based on the extreme value distribution curve. The values are listed in

  Table 1. The precipitation of short-duration rainfall in the urban area of Wuhan (mm) [35].

  Duration (h)	Return Periods in Wuhan Urban Area (mm)
  	1-year	5-year	20-year	100-year
  12	73.6	138.1	205.6	336.2

  .
• Waterlogging points. Waterlogging points are used to verify the corresponding relationships between the delimitated watershed and the real waterlogging points. In this paper, we collected 40 typical waterlogging points on the main road inside the third ring road of Wuhan during the flood and waterlogging events during 30 June 2016 to 6 July 2016. The original waterlogging road cross-information have been reported in the Flood Control and Rapid Report (6 July 2017) by the Wuhan Water Bureau (www.whwater.gov.cn/water/fxkht/7958.jhtml, visited: 15 November 2017). Based on the report, we acquired the coordinates of the waterlogging points using the system for obtaining points from POI (http://api.map.baidu.com/lbsapi/getpoint/index.html) provided by the Baidu Company.

2.2. Research Method

The provided watershed delimitation method using no-fill DEM and precipitation considered the sink area to be part of the real terrain and, so, could not be filled. All the closed low-lying or plateau areas were isolated from the outside unless the water level collected in certain rainfall conditions could reach its outlets; in which case, the flow direction of the sink area was analyzed using the concept of “source flood”. This means that, in certain conditions, the sink area could be connected with other cells, and the whole sink area was considered as a special “big cell” (this means all the cells in the sink area will be logically considered as a whole), which only flowed out through the potential flow direction of the original outlets of the sink area. According to the flow paths of the ordinary cells and the special cells in the whole catchment, the watershed scheme could be delimitated in the same way.

As well as we know, the sub-watershed shall be determined by the terrain characteristics and the distribution of the drainage system. In this paper, the precipitation is used to determine the basic unit, which is used to derive the sub-watershed. We assumed that the rainwater will be finally drained into the sink area by the drainage system, and the sink area can overflow through outlets while its water level reaches to its outlet to form a basic unit with the connected neighbor sink areas. Only the overflowing basic units can trace flow accumulation with connected cells to be identified as part of a larger sub-watershed. Those basic units which cannot overflow will be considered as isolated sub-watersheds. Therefore, the sub-watershed delimited by the proposed method is a logical component of the catchment that is divided in terms of the rainwater storage of flat and sink areas. Additionally, the sub-watershed delimitated by the proposed method is supposed to correspond with the rainwater distribution in certain precipitations. The provided method is comprised of a flow calculation and watershed delimitation, as shown in Figure 2.

• Flow calculation (the first step). The main process of this part consists of flow direction calculation, original sink areas identification, and flow accumulation analysis. The input to this part is the DEM dataset, and its outputs are the flow direction dataset, the original sink area, and the flow accumulation dataset. First, the flow directions are calculated based on the multiple flow concept, where all the lower cells are considered as potential outflow directions. Then, the original sink area can be identified by labeling the cells having no flow direction, which are grouped into consolidated units, in terms of spatially adjacent cells forming continuous regions. Finally, the basic units are consolidated into a special big cell, adopting the flow accumulation method (as for an ordinary cell) to obtain the flow accumulation dataset.
• Watershed delimitation (the second step). The watershed can be delimitated by the input precipitation and the datasets of merged basic units and flow direction. For each of the original sink areas, the limitation precipitation is determined by its area and storage, which is related to the outlets. When the input precipitation is higher than its limitation precipitation, it can flow out from its outlets, either merged with those of its neighborhood units or directly pouring through the directions of its outlets. The watershed can be delimitated in terms of the merged basic units and flow direction datasets. Additionally, the sub-watershed delimited by the proposed method is a logical component in terms of the rainwater storage of flat and sink areas, which are supposed to be used in analyzing the rainwater distribution in certain precipitations in the geological view.

2.2.1. Part 1: Flow Calculation

(1) Calculate the Flow Direction by No-Fill DEM

The traditional single-flow algorithm D8 uses the direction with the greatest elevation drop as the flow direction; in the case where there are several cells which have the same lowest elevation, it chooses one as the flow direction, which introduces some randomness into the flow direction choice. The D8 algorithm fills the sink areas to obtain a continuous flow path, which changes the real terrain features.

In this paper, based on no-fill DEM, we adopt the flow direction calculating method of adopting all of the cells having a lower elevation than the central one as the potential outlets. We adopted C# to develop the program, and used Microsoft Visual Studio 2010 as the development environment. Referring to the flow matrix in the D8 algorithm, we used a move window size of 3 cells × 3 cells to calculate the flow path around the eight adjacent positions and recorded their numeric values in the set {1, 2, 4, 8, 16, 32, 64, 128}, which respectively indicate the directions of east, south-east, south, south-west, west, north-west, north, and north-east. Furthermore, we added a numeric value of “0” to record that all the eight neighborhood cells could not flow out. We used a continuous value between 1–255 to record the potential multiple flow directions by accumulating the numeric value of the potential flow directions together, as shown in Equation (1):

$$Di{r}_{flow}+=\{\begin{array}{l}{2}^i,if\left(Z_0-Z_i\right)>0\\ 0,if\left(Z_0-Z_i\right)\le 0\end{array},\left(0\le i\le 8\right),$$

(1)

where $Di{r}_{flow}$ is the flow direction of the current cell; $i$ is the index of the eight adjacent cells, starting from the east, south-east, west, and so on, in a clockwise order; $Z_0$ is the elevation of the center cell; and $Z_i$ is the elevation of the ith adjacent cell.

(2) Identify the Original Basic Units

The original basic units are those cells which have a flow direction recorded as “0” and all of their adjacent cells can flow into them, as shown in Figure 3a.

The main pseudo-code for identifying a basic unit is shown in Algorithm 1, which identifies sink areas and their flows into cells. We adopt the flow direction dataset as the main input data and label basic units by a unique number starting from 256. Those cells which can fully flow into the sink area (which are labeled by the unique value “number”) are labeled as “number + 1”. All of the group consisting of cells labeled by “number” and “number+1” are marked as in the same basic unit to derive the basic unit dataset (named “mSink”). At the same time, we list the unique basic units stored in the list (named “B_list”).

Then, we calculate the main parameters of the basic units, including Edge, Outlets, Next, Area, and Storage. Their meanings are listed as follows:

✧

Edge: The cells on the border of the basic unit.

✧

Outlet: The cells having the lowest elevation on the edge. There may exist one or more outlets.

✧

Next: The neighboring basic units whose limitation precipitations are lower than the considered cell. If the outlet is into another neighborhood basic unit, the unique number of the basic units will be recorded in the parameter. This parameter is permitted to be null.

✧

Area: The amount of the cells belonging to the basic unit.

✧

Storage: The rainwater volume of the basic unit, in the condition that the water level is lower than the outlet.

(3) Tracing the Flow Accumulation among the Ordinary Cells and the Basic Unit

In the flow accumulation step, the flow accumulation is calculated among the ordinary cells and the basic unit, as shown in Figure 3b. We consider the basic unit as a special big cell, which collects rainwater inside of its area as a whole, flowing outside from its outlets evenly when its water level is higher than the elevation of the outlets. Therefore, both the ordinary cells and the special big cells will be distributed a default unit of water, as shown in Equation (2).

$$Wate{r}_{orig}=\{\begin{array}{l}1,if(orig=cell)\\ Uni{t}_{area},if(orig=B_{unit})\end{array},All\left(cell\right)\cup All\left(B_{unit}\right)=Catchment,$$

(2)

where $Wate{r}_{orig}$ represents the original water for the accumulation flow in each cell (it is recorded by $orig$), $cell$ represents a normal cell, and $B_{unit}$ represents a basic unit (special cell).

Considering the complex runoff process, including rainfall, infiltration, and runoff production, the water collected inside a basic unit can be calculated by hydrological numerical simulation models, such as those based on the traditional model—the comprehensive runoff coefficient model. The comprehensive runoff coefficient can be calculated by the runoff coefficient of the land type, according to weighting by the area of the underlying surface structure. Then, the precipitation limitation can be calculated based on the storage, area, and comprehensive runoff coefficient, as shown in Equation (3):

$$P_{\lim }=\frac{Uni{t}_{storage}}{{\psi }_{av}\ast Uni{t}_{area}},{\psi }_{av}=\frac{{\displaystyle \sum A_i{\psi }_i}}{{\displaystyle \sum A_i}},$$

(3)

where $P_{\lim }$ is the limitation precipitation (mm); $Uni{t}_{storage}$ and $Uni{t}_{area}$ are the storage and the area of the basic unit, respectively; and ${\psi }_{av}$ is the regional comprehensive runoff coefficient in the basic unit (i.e., ${\psi }_{av}$ represents the integrated runoff coefficient). In this paper, we set the default runoff coefficient as ${\psi }_{av}=1$, rather than adopting the vegetation-impervious surface-soil-water (VIS-W) underlying surface model to calculate the runoff coefficient.

Finally, all of these basic units can be considered as a big special cell. If the defined precipitation is higher than the limitation precipitation, the basic unit can converge with the outside through its outlets, and the big special cell can flow out from its outlet. To avoid circulation of water flow between the basic unit and its outlets, the flow direction of outlets which can flow into the basic unit are removed from its original direction. Those special cells that cannot flow out are considered as isolate areas that can only receive flow. The flow accumulation of all the converged special cells and the ordinary cells can be obtained based on the flow traced by unit water through each cell, which is distributed to each potential direction evenly. As every cell may flow into each of its eight neighboring directions, we set the default basic water of each cell to 8.4 in order to avoid the unit water being unable to be divided by the directions (the number of flow directions is between 1–8). The main pseudo-code of flow accumulation is demonstrated in Algorithm 2.

2.2.2. Part 2: Watershed Delimitation

In the concept of basic unit overflow in the condition of “water level reaching outlet”, the defined precipitation determines whether the basic unit merged with these basic units contains its outlet. The basic units are taken from the original basic unit dataset.

If the defined precipitation ($P$) is higher than the limitation precipitation ($P_{\lim }<P$) of the current basic unit, the neighborhood of basic units adjacent to its outlet will be merged with the current visited basic unit by updating its unique number to the visited number of the basic unit, following which, the parameters Edge, Outlets, Next, Area, and Storage will be updated. Until all of the basic units have been visited, the merged basic unit dataset will be used. Based on the basic units and flow direction datasets, the watershed can be delimitated by tracing the flow path through the whole catchment, as shown in Figure 3c.

The main pseudo-code of watershed delimitation is shown in Algorithm 3.

3. Results and Discussion

3.1. The Sink Area Identified as A Basic Unit Related to Precipitation

The original sink areas are identified by the flow direction dataset, while the final sink areas are marked as a basic unit according to the precipitation, as shown in Figure 4.

The cells with flow directions recorded as “0” are identified as sink areas, and the adjacent and continuous cells which can flow into the sink area are mapped as the same original basic units, as shown in Figure 4a. We find that the cells belonging to a large continuous area in a vector map (e.g., in the water layer, the cells in lakes, rivers, and reservoirs) cannot be definitely marked as in the same sink area. The reason for this might be that these elements in the vector map might represent one entity but their elevations are different, and so, they will derive different areas of flow range to derive different sink areas. It is obvious that the sink areas in low-lying areas seem to be larger than those sink areas in common land areas, as large, closed sink areas have more chances to occur in the low-lying areas.

The basic units merged by the defined precipitation are shown in Figure 4b–e. The precipitation in the four figures are 73.6 mm, 138.1 mm, 205.6 mm, and 336.8 mm, respectively. We find that the basic units merged by higher precipitation produced a larger area of basic units and that the number of basic units with large areas increased slowly from Figure 4b–e. For example, in Figure 4b, the precipitation was 73.6 mm. There were three basic units covered by the rectangle located in the Yangtze River, while in Figure 4c–e, only two basic units were left, as the two basic units at the top had been merged together. In the square area located at Lake Dong, the number of basic units showed the same performance; however, the number of basic units in Figure 4b was less than the number of basic units in Figure 4c–e, where the number of basic units were the same in Figure 4c–e. This phenomenon indicated the basic units merged in the scheme were stable within a certain phase of precipitation.

We analyzed the amount and area changes of the basic units in eight aspects, as shown in

Table 2. The statistical information of basic units in four types of return periods of Wuhan City.

Index	Statistical Items		Precipitation (mm)
			73.6 (1 year)	138.1 (5 years)	205.6 (20 years)	336.2 (100 years)
1	Amount of basic units		51,811	47,844	46,962	46,293
2	Amount of isolated units		13,794	4311	2175	325
3	Area of basic units (cells)	Average	30.55	33.08	33.70	34.19
4		Minimum	1	1	1	1
5		First quartile	7	6	6	6
6		Median quartile	11	11	11	11
7		Third quartile	21	23	24	24
8		Maximum	44,960	45,535	45,675	45,791

. The first element in row 1 (the amount of basic units) is the amount of all the basic units. The second element in row 2 is the amount of isolated basic units whose limitation precipitation were larger than the defined precipitation. The six elements in rows 3–8 are statistical information about the area of the basic units. The element in row 3 is the average area. The elements in rows 4–8 are the minimum, the first quartile (25%), the second quartile (50%), the third quartile (75%), and the maximum value of area.

We find that the average area of the four basic units was about 30 cells, where the change was very small as the precipitation increased, and most of the basic units were very small. Among the four kinds of precipitation, more than 75% of the basic units had an area of less than 25 cells, lower than the average value. For example, when the precipitation was 73.6 mm, the average area of a basic unit was 30.55 cells, while the minimum and maximum areas were 1 and 44,960 cells, respectively. The 25%, 50%, and 75% quartiles of these basic units were no more than 7, 11, and 21 cells, respectively. We also found that the total amount of the four basic units decreased as the precipitation increased. Among the basic units, the isolated basic units were the areas whose precipitation was higher than the defined precipitation in each scheme.

From Table 2, row 2, we can see that the amount of isolated units decreased as the precipitation increased. This meant that, although the number of basic units was large and there were a lot of small areas, most of the basic units could converge with their neighborhoods through their outlets.

We constructed a histogram based on the amount of isolated basic units, as shown in Figure 5. We found that the ratio of isolated basic units to total units in the four schemes were 27%, 9%, 5%, and 1%, respectively. This result was consistent with our expectations of the proposed algorithm: with the increment of precipitation, more and more of the isolated basic units will flow into their outlets, either to their neighborhood basic unit or to ordinary cells having certain flow directions.

We analyzed the unique area value of the basic units, which is shown in Figure 6. All the histograms in Figure 6a–d adopted the unique area as the abscissa axis, and all maps adopted the count of the unique area value to draw the histogram in blue, which were mapped according to the abscissa axis and the main vertical axis. All of them used the cumulative ratio of the unique area value to draw the cumulative line graph in yellow, which were mapped according to the main abscissa axis and the secondary vertical axis.

From Figure 6a–d, we found that, in the four schemes, 99% of the basic units contained less than 176 cells. Compared to the DEM original resolution of 30 m × 30 m, 176 cells represented about 0.16 km², which is close to a square with a side length of 400 m. This mapping result also indicated that, at the catchment scale, most of the sink areas were small no matter which input precipitation was used.

From Figure 6a–d, we can observe that the cumulative ratio of basic units increases rapidly at the beginning, before leveling-off for the remainder of the cumulative ratio curve. This indicated that most of the basic units were small in size. In

Table 3. The cumulative ratio of area in the four schemes of precipitation to divide basic units.

Scheme	Cumulative Ratio	Area (Cells)	Count	Ratio	Scheme	Cumulative Ratio	Area (Cells)	Count	Ratio
73.6	3%	1	1478	2.83%	138.1	3%	1	1469	3.07%
	6%	2	1583	3.06%		6%	2	1574	3.29%
	15%	4	2838	5.48%		10%	3	1786	3.73%
	24%	6	2489	4.80%		24%	6	2434	5.09%
	47%	9	7900	15.25%		48%	10	1201	2.51%
	71%	19	842	1.63%		72%	20	697	1.46%
	80%	26	412	0.80%		80%	28	371	0.78%
	90%	47	151	0.29%		90%	52	113	0.24%
	99%	221	4	0.01%		99%	241	3	0.01%
	100%	44,960	1	0.00%		100%	45,535	1	0.00%
205.6	3%	1	1469	3.13%	336.2	3%	1	1470	3.18%
	6%	2	1575	3.35%		7%	2	1573	3.40%
	16%	4	2811	5.99%		12%	4	2802	6.05%
	26%	6	2399	5.11%		24%	6	2391	5.17%
	49%	10	1165	2.48%		49%	10	1148	2.48%
	72%	21	626	1.33%		72%	21	616	1.33%
	80%	29	355	0.76%		80%	30	346	0.75%
	90%	54	102	0.22%		90%	55	131	0.28%
	99%	251	6	0.01%		99%	257	1	0.00%
	100%	45,675	1	0.00%		100%	45,791	1	0.00%

, we list the basic unit area distribution in the four schemes of precipitation.

From Table 3, we can observe that when the cumulative ratio maintains a certain value—as the precipitation increases from 73.6 mm to 336.2 mm—the size of the basic unit increases by a small amount. For example, when the cumulative ratio is 80%, the areas at 73.6 mm, 138.1 mm, 205.6 mm, and 336.2 mm of precipitation are 26, 28, 29, and 30 cells, respectively. In the four schemes, the areas of 90% of the basic units are less than 60, and of 99% of the basic units, are less than 260. The areas of the four schemes differ little. As a result, the precipitation has little effect on the size of the basic unit.

3.2. The Flow Accumulation Reflects the Flowing Relations Inside the Catchment

Flow accumulation data is the basic information used to trace flow paths. Since the flow accumulation of certain cells is calculated by its flow-in cells and flow-out directions, either a common cell or special “big cell” can achieve a flow-in condition in which its flow-in cells and flow-out direction are assured. For a basic unit, as it is considered as a special big cell, which can only flow-out through the outlets, this means that the flow-out direction can be assured. As all the cells can fully flow into basic units that have been identified as part of the basic unit, these outside cells play a part in the flow-in direction of a “big cell”. Therefore, the “big cell” successfully dictates the flow-in cells and flow-out directions. Finally, either the common cells or the special “big cells” can derive the flow accumulation in the same way, by summing the values of the flow-in basic water and its original accumulation water. The final flow accumulation dataset is shown in Figure 7.

As the number of the flow direction for the original cell was between 0–8, we allocated the basic water of each cell as 8.4, such that it could be fully divided in the situation of the highest count of flow directions. For example, as the amount of flow directions is {1, 2, 3, 4, 5, 6, 7, 8}, the average distributed water according to the flow direction will be {8.4, 4.2, 2.8, 2.1, 1.68, 1.4, 1.2, 1.05}, respectively. After tracing the flow path, the accumulation can be calculated. We adopted a value of 336 (almost 40 cells along the flow path) for flow accumulation in order to level and color the flow accumulation dataset. As the resolution of DEM dataset was 30 m × 30 m, the area of 40 cells represented 36,000 m², as shown in Equation (4),

$$Area=cells\times resolutio{n}_x\times resolutio{n}_y,$$

(4)

where $Area$ (its unit is m²) represents the final area, $cells$ represents the amount of cells along the flow path, and $resolutio{n}_x$ and $resolutio{n}_y$ represent the resolution of X and the resolution of Y.

Therefore, as shown in Figure 7, the flow accumulation was mapped into five kinds of colors (white, light green, dark green, naval blue, and dark blue): the white color represents flow accumulation values in the range of (0–336), where the flow area is 0.036 km²; the light green color represents flow accumulation values in the range of (336–3360), where the flow area is 0.36 km²; the dark green color represents flow accumulation values in the range of (3360–33,600), where the flow area is 3.6 km²; the naval blue color represents flow accumulation values in the range of (33,600–336,000), where the flow area is 36 km²; and the dark blue color represents flow accumulation values greater than 336,000, where the flow area is greater than 36 km².

From Figure 7, we can see that all of the natural water systems, such as lakes and rivers, had a flow accumulation across about 3.6 km². The accumulative areas of most parts of the Yangtze River, Lake Dong, Lake Yandong, Lake Yanxi, Lake Nan, and Lake Qingling were about 36 km², while the accumulative area of Lake Tangxun was larger than 36 km². We can see that Lake Huangjia, Lake Sha, and part of Lake Nan, part of Lake Yandong, and part of Lake Bei were colored dark green, which represented a flow area of 3.6 km². If we only consider the natural drainage by elevation drop, the lake areas are the main places to store rainwater during certain rainfall conditions. Furthermore, not all the lakes—even different parts of the same lake—can be concluded to converge rainwater at the same level. Among these lakes, Lake Tangxun may receive more runoff water than the other lakes in Wuhan. In Figure 7, the top buffer layer is the analysis result of Lake Tangxun derived by 6km from outside. This means the buffer area can represent the flow-in area of the Lake Tangxun basic unit, approximately. We can observe that the Lake Nan, Lake Huangjia, and even a little part of the Lake Dong southern area, are located in this buffer area. This phenomenon indicates that the Tangxun basic unit may receive the overflow water from Lake Nan, Lake Huangja, and part of Lake Dong, in theory. This means that the waterlogging disaster risk in the Lake Tangxun basic unit may be higher than the other three basic units (Lake Nan, Lake Huangja, and Lake Dong). We also observe that this buffer area almost contains all of waterlogging points distributed between Lake Dong and Lake Tangxun. This may confirm the fact that the waterlogging disasters which took place in this area were more serious than other places in the flooding and waterlogging event during early July in 2016.

3.3. Watershed Delimitation in Term of Precipitation

According to the datasets of merged basic units and flow direction, the watershed delimitation method can trace the flow paths. Those basic areas that can flow out can be considered as special big cells converging with the outside by their outlets, and those big special cells which cannot flow out might be treated as isolated basic sub-watersheds that only receive water from the outside normal cells or special big cells.

According to the calculated design precipitation of a 12 h intensity, we adopted four kinds of precipitation (in return periods of 1 year, 5 years, 20 years, and 100 years) to delimitate sub-watersheds, as shown in Figure 8. It can be seen that, as the amount of large sub-watersheds increased, the small sub-watersheds (especially those located inside the range of the large ones) seemed to diminish. Furthermore, the distribution of isolated sub-watersheds tended to decrease as the precipitation increased, as can be seen from Figure 8a–d.

From Figure 8, it can be concluded that, as precipitation increases, the average sub-watershed area will increase and the amount of sub-watersheds will reduce, as shown in

Table 4. Statistical information of four kinds of watershed delimitations derived by the corresponding design rainfall precipitations.

Index	Statistical Items		Precipitation (mm)
			73.6 (1 year)	138.1 (5 years)	205.6 (20 years)	336.2 (100 years)
1	Amount of sub-watersheds		16,106	7163	3747	1479
2	Sub-watershed area (cells)	Average	187.07	420.67	804. 06	2038.45
3		Minimum	6	6	6	6
4		First Quartile	24	27	30	33
5		Median quartile	49	60	68	75
6		Third quartile	111	153	178	201.5
7		Maximum	74,721	224,714	283,163	307,768

. There were more than 75% of sub-watersheds having smaller areas than the average value. The amount of sub-watersheds in the four schemes decreased as the precipitation increased from 73.6 to 336.2. For example, in the result for a precipitation of 73.6 mm, the average area of a sub-watershed was 187.07 cells; in contrast, for the result with a precipitation of 336.3 mm, the average area was 2038.45 cells. Comparing the minimum and maximum values, all of the minimum areas were 1 cell, and the maximum values for 73.6 mm and 336.2 mm precipitation were 74,721 cells and 307,768 cells, respectively.

We analyzed the quality of the sub-watersheds in the four schemes, as shown in Figure 9. All of them adopted the count value of the area to draw the histogram, mapped by the abscissa axis and the main vertical axis, and used the cumulative ratio of the area to draw the cumulative line graph mapped by the main abscissa axis and the secondary vertical axis.

From Figure 9a–d, we can see that the amount of the unique area value of sub-watersheds decreased as the precipitation increased. It is clear that the distribution of the histogram was sparser when the precipitation was 336.2 mm, compared to when the precipitation was 73.6 mm. This indicates that the amount of sub-watershed with 336.2 mm precipitation was less than that for 73.6 mm. We also found that the cumulative curve changed more gently from Figure 9a–d as the precipitation increased. This indicates that the sub-watershed area changed from small to large, becoming more and more evenly distributed, rather than almost all of the basic units being in a small area. This result shows that the proposed method seems to be consistent, compared with the case where the average area of the four schemes continued to increase, as in Figure 9a–d.

From Figure 9a–d, we can observe that the cumulative ratio of the sub-watersheds increases rapidly at the beginning, before leveling off for the remainder of the cumulative ratio curve. However, compared to the cumulative ratio of the basic units, the inflection point of the sub-watersheds’ cumulative ratio curves were located before that of the basic units’ cumulative curves, and the rates of change of the sub-watersheds were more steady than those of the basic units. In

Table 5. The cumulative ratio of the areas in the four schemes of precipitation to divide the sub-watersheds.

Scheme	Cumulative Ratio	Area (Cells)	Count	Ratio	Scheme	Cumulative Ratio	Area (Cells)	Count	Ratio
73.6	3%	9	381	2.37%	138.1	4%	11	86	1.20%
	6%	11	328	2.04%		9%	14	128	1.79%
	12%	14	167	1.04%		13%	17	84	1.17%
	25%	23	243	1.51%		24%	26	92	1.28%
	48%	46	109	0.68%		48%	56	41	0.57%
	72%	98	44	0.27%		72%	132	12	0.17%
	80%	140	27	0.17%		80%	199	8	0.11%
	90%	262	9	0.06%		90%	499	3	0.04%
	99%	2077	1	0.01%		99%	6176	1	0.01%
	100%	74,721	1	0.01%		100%	224,714	1	0.01%
205.6	4%	11	29	0.77%	336.2	3%	11	10	0.68%
	7%	13	63	1.68%		6%	13	23	1.56%
	16%	18	43	1.15%		12%	19	13	0.88%
	24%	29	38	1.01%		24%	31	14	0.95%
	48%	63	21	0.56%		48%	69	5	0.34%
	72%	152	6	0.16%		72%	169	2	0.14%
	80%	234	2	0.05%		80%	277	1	0.07%
	90%	557	1	0.03%		90%	767	1	0.07%
	99%	15,099	1	0.03%		99%	55,420	1	0.07%
	100%	283,163	1	0.03%		100%	307,768	1	0.07%

, we list parts of the area distribution of sub-watersheds in the four schemes.

From Table 5, we can observe that, as the cumulative ratio of the sub-watersheds increased from 73.6 mm to 336.2 mm, the size of the sub-watersheds evidently increased. For example, when the cumulative ratio is 80%, the areas with 73.6 mm, 138.1 mm, 205.6 mm, and 336.2 mm of precipitation are 140, 199, 234, and 277, respectively. In the four schemes, the areas of 90% of the basic units are less than 770. However, when the cumulative ratio is 99%, the areas with 73.6 mm, 138.1 mm, 205.6 mm, and 336.2 mm of precipitation are 2077, 6176, 15,099, and 55,420, respectively. This indicates that the precipitation levels directly affect the areas of the sub-watersheds and may result in large areas of sub-watersheds, as shown in Figure 10.

As shown in Figure 10, we can find that, as the precipitation increased from 73.6 mm to 336.2 mm, the area of the certain cumulative ratio increased. As the precipitation increased, this resulted in the occurrence of a large area of sub-watersheds.

3.4. Comparing the Sub-Watershed Delimitation with the Distribution of Waterlogging Points

In order to verify whether the proposed method was suitable for flood control practice, we analyzed the correspondence between the distribution of real waterlogging points and the watershed delimitation scheme. As the D8 algorithm is widely used in watershed delimitation, we adopted the watershed schemes produced by the ArcGIS 10.2 hydrology tool, which is a concept of the D8 algorithm, using the proposed method, as shown in Figure 11. The 40 locations of waterlogging points inside the third ring road of Wuhan in which flooding and waterlogging events occurred from 30 June 2016 to 6 July 2016, were identified.

The D8 algorithm, which uses a filled DEM, is the traditional, commonly used method to delimitate watersheds in hydrology research. As the hydrology tool can identify the sub-basin areas and sub-watersheds, and the sub-basins are the areas divided by the ridge lines of the terrain, which is similar to the sub-watershed delimitation provided in this paper, we adopted the basin delimitation overlapped with the waterlogging points shown in Figure 11a. We can see that the sub-basin in blue (followed by the name of Basin No. 2) contained most of the Yangtze River in the area (along the downstream at the junction of the Yangtze and Han Rivers), Lake Dong, Lake Yanxi, and Lake Yandong. The sub-basin in green (followed by the name of Basin No. 1) was at the north-west shore of the Yangtze River and the Han River. The sub-basin in yellow (followed by the name of Basin No. 3) contained the west part of the Yangtze River, the downstream of the Han River (pouring into the Yangtze River), Lake Nan, Lake Huangjia, and Lake Tangxun. The sub-basin in purple (followed by the name of Basin No. 4) was located at the south-west shore of the Yangtze River and south of the Han River, as well as Lake Moshui and Lake Sancha.

During the flooding and waterlogging event from 30 June 2016 to 6 July 2016, in Wuhan, the average cumulative precipitation was 532.8 mm. Thus, based on the concept of sink area outflow with the condition of “water level higher than outlet”, we defined the precipitation as 532.8 mm to delimitate the sub-watershed, as shown in Figure 11b. We can see that the Yangtze River was divided into two sub-watersheds (indicated as sub-watersheds No. 3 and No. 9). Although the waterlogging points were distributed in the Qiaokou and Jianghan Districts near the river, most of them were not in sub-watersheds containing the Han or Yangtze Rivers. We can also see that the Hongshan District was mainly divided into three sub-watersheds (indicated as sub-watersheds No. 5, No. 6, and No. 8). Sub-watershed No. 8 contains the Lake Nan and the Lake Tangxun, sub-watershed No. 6 is located at the edge of No. 8 and No. 5, and sub-watershed No. 5 contains the Lake Dong. We can see that most of the waterlogging points were located in sub-watersheds No. 8 and No. 6, and only three waterlogging points are in sub-watershed No. 5. According to the terrain characteristics, the lowest elevation on the shore of Lake Dong was higher than that of Lake Tangxun and Lake Nan. Therefore, this phenomenon indicated that the waterlogging points in Hongshan District near Lake Nan and Lake Tangxun were more serious than those at Lake Dong, which may be more acceptable and reasonable.

We list the distribution of the waterlogging points grouped by the two sub-watershed delimitations, as shown in

Table 6. Statistical information of four kinds of watershed delimitations derived by the corresponding design rainfall precipitations.

Method	Contain Range of River	Watershed	All Cars Cannot Pass			Small Cars Cannot Pass
			Amount	Subtotal	Ratio	Amount	Subtotal	Ratio
D8 algorithm	Not contain	No.1	9	14	40%	3	3	60%
		No.4	5			/
	Contain	No.2	8	21	60%	2	2	40%
		No.3	13			/
Proposed method	Not contain	No.1	10	32	91%	3	5	100%
		No.4	2			2
		No.5	4			/
		No.6	3			/
		No.7	5			/
		No.8	8			/
	Contain	No.2	2	3	9%	/	0	0%

. In this paper, we consider that the sub-watershed that does not connect with the main outlet may need to be flown outside to assist drainage. While the drainage ability will be restricted by the pumping work conditions, the sub-watersheds connected to the main outlet can drain rainwater according to the flow path independently. Therefore, we assume that these sub-watersheds do not connect to the main outlet, which may increase the risk of waterlogging events compared to those sub-watersheds connected to the main outlet directly. As the main outlet of Wuhan is the Yangtze River and Han River, we group the sub-watersheds as “contain range of river” and “not contain range of river”.

From Table 6, we can observe that in the scheme of the D8 algorithm—the situation of all cars cannot pass—there are 14 waterlogging points located in “not contain range of river” sub-watersheds, and there are 21 waterlogging points located in “contain range of river” sub-watersheds. In the situation whereby small cars cannot pass, there are three waterlogging points located in “not contain range of water” sub-watersheds and two waterlogging points located in “contain range of river”. This means that 60% of “all cars cannot pass” and 40% of “small cars cannot pass” waterlogging events take place in the “contain range of river” sub-watersheds. As in the scheme of the proposed method, we can observe that in the situation of “all cars cannot pass”, 32 waterlogging points were located in “not contain range of river” sub-watersheds and only three waterlogging points were located in “contain range of river” sub-watersheds; in the situation “small cars cannot pass”, all of the waterlogging points occur in “not contain range of river” watersheds. Therefore, 91% of “all cars cannot pass” and 100% of “small cars cannot pass” waterlogging events took place in “not contain range of water” watershed.

Therefore, the scheme of sub-watershed delimitated by the D8 algorithm performed less consistently with the distribution of the real waterlogging points, compared with the sub-watershed delimited by the proposed method. This indicated that the proposed method may be useful in the practice of analyzing the possible waterlogging points that would occur under certain intensities of rainfall.

4. Conclusions

In this paper, a dynamic watershed division method by identifying basic units (basic units contain the cells of certain sink areas and those cells that flow into the sink areas) and merging them in terms of the convergence of terrain and precipitation was proposed. This method adopted the original DEM without filling or removing the sink areas, such that it could retain the real terrain characteristics of low-lying and flat lands. This character indicates that the proposed method can be used in the rainwater and sedimentary carrying capacity analysis, which is sensitive in the original storage ability of sub-watersheds.

In this paper, we adopted Wuhan City as a study area, using the precipitation of design return periods of 1, 5, 20, and 100 years to derive four sub-watershed schemes. We verified that the watershed delimitated by our method is reasonable from the view of the consistency of the basic units and sub-watersheds delimitated by the proposed algorithm. According to the delimitation overlapped with the real waterlogging points from the flooding and waterlogging event during early July of 2016 in Wuhan, we found that sub-watershed delimitation based on the average accumulative precipitation by the proposed method can well explain the distribution of waterlogging points, while the sub-basin delimitated by the traditional D8 algorithm performed worse in the corresponding analysis with the waterlogging distribution. Flow accumulation analysis also showed that the provided method also showed the corresponding with the distribution of waterlogging points.

In this work, we adopted the AW3D30 DEM dataset, whose horizontal and vertical resolution were 30m and 1m, respectively. In theory, the proposed method could be used in different types of DEMs datasets, such as ASTER DEM or TanDEM-X. The resolution was coarse, as urban catchments usually have large areas of low-lying and flat areas, which need more refined DEM to describe the terrain features. This might be the reason why there were a large number of basic units in small areas, which would influence the derivation of isolated and small-scale sub-watersheds in the final delimitation. Furthermore, we did not distinguish the runoff production influenced by different land cover types and did not add the effects on flow direction changed by the drainage system. In the future, we will pay more attention to these problems and conduct research in resolving small-scale sub-watersheds and increasing quantitative effectiveness in our watershed divisions.