Highly Resolved Runoff Path Simulation Based on Urban Surface Landscape Layout for Sub-Catchment Scale

: The present study explored the regularities of the path and network structure of surface runoff formed under the inﬂuence of urban surface landscapes. We used unmanned aerial vehicle sensors to examine terrain and land use/cover change. The sub-catchments of a typical city, Luohe, China, were evaluated for the effect of landscape on surface runoff. Landscape and topographic parameters from 166 urban sub-catchments in Luohe were obtained by measuring digital surface models and orthophoto maps. The minimum cumulative resistance model was used to simulate potential runoff and 491,820 potential runoff paths, connected upstream and downstream, were obtained in 166 sub-catchments. The chi-square test was used to compare simulation runoff paths and actual runoff depth, with the results showing that they led to the same distribution trend. When the gravity coefﬁcient was greater than 18.93, path disconnection occurred among 166 subcatchments, with a decrease in channels. The potential runoff distribution appeared in aggregation; as the gravity coefﬁcient increased from low to high, aggregation showed a trend of increasing initially but subsequently decreasing. The initial runoff formed sub-catchments with high gravity coefﬁcients, then accumulated and spread to the others. It is important that proper measures are taken to establish a uniﬁed planning of the city’s surface landscape in order to produce suitable surface runoff distribution.


Introduction
Rapid urbanization has caused huge changes in urban sub-catchments and profoundly affected surface runoff. It not only leads to a decrease in infiltration, but also increases runoff and ultimately affects the hydrological cycle and urban environment [1,2]. Studies have confirmed that land use/cover change (LUCC) induced by urban expansion and the subsequent increase in the areas of impermeable ground surface are considered the main causes of changes in urban hydrology. In addition, the expansion of impermeable areas is linked to urban disasters, such as flash floods, non-point source pollution, water shortages, and river deterioration, as well as climatic change [3,4]. Studies have shown that without factoring in the potential impacts of climate change, cities around the world will be 2.7 times more likely to suffer flood hazards by 2030. For example, a devastating flash flood occurred in Beijing, the capital of China, on 21 July 2012. This flooding event resulted in a loss of 79 lives and USD 1.67 × 10 12 in loss. The connection between surface runoff and landscape has been steadily established through further research. Impervious landscapes contribute to increased accumulation of stormwater and urban floods [5], and reducing impervious surfaces and planting scattered vegetation that can maintain and purify water quality and quantity [6]. Impervious areas amplify the spread of surface runoff into urban areas. The common conventional solution is often to discharge water directly through channels and ment is based on the application of UAVs to acquire images of high spatial resolution. UAV sensors can acquire surface landscape data which is needed to model the urban surface runoff within a city for analysis. The strength of our analysis is that we acquired a 0.09 m resolution orthogonal projection and 0.09 m vertical resolution digital surface model (DSM) of the entire urban study area through the UAV [15].

Study Area
Luohe is located in middle Henan Province, China ( Figure 1). The administrative division has an area of 2617 km 2 , with 78 km 2 dedicated to the city with a population of 2.84 × 10 7 , where approximately half of the people lived in the central district at the end of 2018. Our study focused on this central district of 169.38 km 2 , where the urbanization rate is 52.47%. The construction area of Luohe increased from 24 km 2 in 1999 to 83 km 2 in 2018, with an annual expansion rate of 3.11 km 2 /year. Green space accounts for 39.8%, with 12.8 m 2 of green space per capita of the city. While the green space increased from 37.4% in 1999 to 41.3% in 2019, this increase was due to the expansion of the city area rather than increased green space within the city itself, and it distribution still remains imbalanced [16]. The city is an alluvial plain by construction of embankments along the confluence of the Sha-li River, increasing flood discharge capacity to protect the city from most rivercresting events. However, the decrease of permeable surfaces has led to an increment of urban waterlogging events, in terms of increased frequency and severity inside the city. Landscape layout is one of the important factors that affect urban waterlogging.

Land Use and Cover Classification
We acquired a digital land use/land cover map using an aerial photography image from a Trimble UX5, fixed-wing, unmanned aerial vehicle (UAV), with a 48 megapixels, silent frame and 0.09 m ground resolution, using a 25 mm focal length lens. The total extent of the image was 169.38 km 2 . Details of the UAV flight included: flight altitude 300 m, a single flight covering an area of 3.7 km 2 flight mode for the 'Z'-type lane, 5:1 aspect ratio of a rectangle, and 80% image overlap rate. The post-processed kinematic and global navigation satellite technology in the UAV system established pinpoint image locations (accuracy > 0.15 m). To obtain a clear land classification and patch boundary identification, aerophotography was conducted during the defoliating period. We converted raster images into vector layers for classification, which were interpreted and verified manually using updated planning maps. This work was based on ArcGIS 10.4. Surface landscape factors were allocated into two groups-permeable (group 1) and impermeable (group 2) -which were categorized into seven classes and classified by dominant human land use based on their influence on surface runoff.

Generation of a High-Resolution Digital Surface Model
In urban GIS and management, digital surface model (DSMs) are used for generating object models. A DSM was generated to obtain the inverted elevations and slopes of the stormwater structures. The preliminary step was to segment the interesting areas. We acquired a 0.09 m vertical resolution DSM of the entire urban study area through the UAV. Then, the terrains could be extracted from DSM ( Figure 3). Since urban appearance is complicated with buildings, green space, water, etc., terrain extraction was implemented through several stages. In the first stage, building roofs were extracted from the DSM. In the second stage, marker-based sub-catchments were implemented to get the boundaries of the objects above the ground [17]. The result of marker-based buildings and sub-catchments were merged to improve the urban terrain extraction accuracy.

Sub-Catchment Division
Sub-catchments were defined according to physical parameters, forming units of urban surface runoff [18,19]. Usually, sub-catchments are divided based on outlets to the sewer system, DSMs, and the urban pipe network in the urban central district within the clear flow direction [20]. Therefore, we acquired 166 sub-catchments. There was a high density of buildings and a network of roads, so the sub-catchments were small and in large numbers. However, the areas at the edge of the city were larger. Among the 166 sub-catchments we analyzed, the largest was 11.02 km 2 , the smallest was 0.072 km 2 , and the average area was 1.03 km 2 ( Figure 4).

Minimum Cumulative Resistance (MCR) Model
With respect to the distribution of runoff process, an urban landscape layout could be divided into a "source" or "sink" [21]. The source consists of those processes in the upstream area from which surface runoff is produced, and the sink consists of the downstream area that receives stormwater runoff [22]. Although source and sink of landscapes have opposite qualities, a sink in one process may become a source in another. The main purpose of the proposed the source and sink landscape theory was to explore how the status quo of the urban landscape would affect the distribution of the runoff process. In turn, the runoff process would help in the identification of a suitable spatial pattern and overall urban landscape in every region [20]. The transformation between source and sink occurs as a mutual conversion process for the path of surface runoff, which was achieved by overcoming various forces of resistance. They were then integrated to reflect the resistance surface [23].
The MCR model essentially reflects the degree of resistance in landscape layout. The cost distance in the MCR model is different from the actual distance. Instead, it emphasizes the relative spatial relationship between points. [24]. This cost distance is based on calculating the resistance coefficient when the target path passes through different landscape units [25,26]. Based on the resistance value of MCR when a target path passes through a specific unit, we can determine the connectivity between the two points. Usually, the "source" is the initiator of runoff distribution and flows into the "sink" through different landscape units. To determine the minimum cumulative resistance value, MCR, the following modified formula is: where f is a function of the positive correlation that reflects the relation of the least resistance from any point to any point in space and the characteristics of the landscape base surface, min is the minimum value of cumulative resistance produced in any processes of source unit i transforming into a different sink unit j, D ij is the spatial distance between unit i and source unit j, and R i denotes the resistance coefficient that exists in transition from source unit, i, to sink unit, j. R i can be determined by the value of i, and the path from unit i to unit j produces different resistance values. When i has been determined, calculating MCR requires selecting the path with the least resistance value of the spatial distance.
In this paper, source refers to the landscape units that initiate the process by which surface runoff occurs. Sink refers to the landscape units that receive inflow from the source runoff. During MCR calculation, the input of the source and sink can be a patch or a combination of patches, and sources could be connected or unconnected in space. We assumed that the divided sub-catchments were the sources and sinks of each other. Surface runoff has the common characteristics of liquid flow and can flow between different sub-catchments [27].

Chi-Square Test of Simulation Results
The chi-square test is a widely used hypothesis testing method, which is used to calculate the degree of fitting between the observed sample and the simulated value [28]. The chi-square test was used to measure the relationship between simulation value and actual runoff depth. The runoff depth values were, respectively, the runoff depth on 27 June, 22 July, and 30 August 2020 at the 36 flooding points in Luohe [16]. Simulation values were potential runoff paths from the MCR model.

Gravity Model
Gravity has been one of the most successful empirical models in spatial analysis. Interactions between spaces, such as information flow and material flow have been related to size and distance. In urban network structures, interactions between objects such as cities or areas lead to models remembering Newton's gravity law, where the sizes of the objects play the role of mass [29]. Incorporating deeper theoretical foundations of gravity into recent practice lead to a richer and more accurate estimation and interpretation of the spatial relations described by gravity. Wider acceptance has followed, and suggestions are made for promising future research [30].
where G ab is the interaction between surface runoff paths a and b of the sub-catchments, N a and N b are the weights between a and b, D ab is the standardized value of potential runoff path resistance in a and b, P a is the resistance of sub-catchment a, P b is the resistance of sub-catchment b, S a is the area of sub-catchment a, S b is the area of sub-catchment b, L ab is the accumulated resistance value of the surface runoff path from a to b in the sub-catchment, and L max is the maximum resistance value of the potential surface runoff path in the entire study area [30].

Spatial Autocorrelation
When the values of variables observed at spatial positions resemble each other more than expected for a randomness model, the variable is said to be spatially autocorrelated [31]. It is a measure of the degree to which objects, situated in close proximity, have a tendency for similar values of a given index. In studies of surface runoff, spatial autocorrelation may arise from the underlying causal variables of landscape layout and various other factors [32]. The spatial distribution of runoff formation can be revealed by analyzing runoff density in each sub-catchment. The potential runoff density is d = l/a, where d is the potential runoff density in each sub-catchment, l is the length of the potential runoff, and a is the area of each sub-catchment.
Moran's I is one of the most popular indices for assessing spatial autocorrelation in spatial data, and is defined as follows: In expression (3), where pi is the runoff density in the ith units, for i = 1..., m, and W ij is a weight assigned to every pair of units. Conventionally, w ii = 0. This weight indicates the importance of the locality pair in measuring the spatial autocorrelation. Usually, units that are adjacent or closer together than some distance thresholds have weights assigned to unity, with the other weights being zero. Other weight structures are also possible. Moran's I usually ranges between 1 and −1, with positive values indicating similarity between geographic neighbors, negative values indicating dissimilarity, and values close to 0 indicating a random pattern. In recent years, spatial autocorrelation has come into its own in studies on urban rainwater distribution.
Industrial zones, business zones and high-density residential zones tend to have a greater influence on runoff [33,34]. Chang et al. [35] examined the spatial patterns of annual runoff ratios and their variability and identified the determinants of runoff indices for 238 reference basins with low levels of anthropogenic influence and 1352 non-reference basins with substantial levels of anthropogenic influence variability across the contiguous U.S. The research results showed local specific relationships between runoff indices and landscape factors [35]. The geo-statistical properties of soil moisture patterns from five different locations in Australia and New Zealand have been researched. The results showed that the processes controlling spatial patterns can change between temporal and spatial scale with catchment moisture status; however, when similar general phenomenon reoccur in a catchment, similar spatial patterns result [36].

Identification and Analysis of Urban Surfaces
In the course of rainfall, surface runoff mainly spreads in the horizontal direction, and the horizontal movement mechanism is the method of redistributing surface runoff in the process of rainfall. It is also the basic form of the catchment confluence process.
In the process of horizontal confluence, surface runoff will be influenced by the surface landscape, including terrain factors (altitude, slope, relief amplitude, and roughness) and LUCC, which produce resistance along the way and local resistance, making the runoff continually change direction to form a network. The network path of surface runoff is determined by the quality of the source and sink, the resistance between the source and sink, and the distance between the source and sink. Among them, the surface landscape plays an important role in the flow process of surface runoff. Combined with literature materials, this study evaluates the resistance of different surfaces generated along the path and local resistance and obtains the confluence resistance consumption surface [27].
First, we extracted the single factor cost surface through satellite imaging and the DSM. Then, values of 120-1000 were assigned to different indices according to different factors (Table 1) and the AHP method was used to calculate the weight base on references ( Table 2) [37][38][39]. Finally, different single factor cost surfaces were superimposed by raster calculation to form the comprehensive cost surface [40,41] (Figure 5).

Distribution Characteristics of Potential Surface Runoff Paths
The potential runoff paths constructed by the MCR model revealed 491,820 runoff connections in total, distributed among 166 sub-catchments. The potential runoff paths were linear, and not only connected the upstream and downstream sub-catchments but also showed clear directions, pointing from sources to sinks ( Figure 6). The simulation results showed that 1.45% of the runoff intersected with the buildings but had an error rate of less than 5%. Causes of conjecture errors include the following: (1) only features and landscape elements were used as resistance characteristics for simulation in this study, (2) aerial imagery accuracy, and (3) the calculation of frictional resistance and local resistance to runoff is more applicable to the calculation of closed pipes or fixed transmission channels. The cost ratio index is used to quantify the average consumer cost of a network, and it mainly reflects the complexity of the network. The closer the cost ratio index is to 1, the more complex the network structure. The formula is as follows: cost ratio = 1 −L/d, (L is the quantity of runoff; D is the runoff length). The resulting cost ratio index was 0.71, and the network costs were relatively high.
Assumption was that potential paths were associated with depth of runoff. The depth of flooding points (n = 36) was measured ( Figure 6) and Chi-square tests show that the signification (p) of repeat the experiment three times was greater than 0.05 (Table 3), respectively. That was to accept the null hypothesis, indicating that the overall efficiency of the simulated value of the null hypothesis was equal to the measured value. The difference between the potential paths and depth of runoff was not statistically significant, and there was no significant difference between the simulated value and the observed value, so the hypothesis of correlation was valid. Meanwhile, the verification result was consistent with my previous research tendency by Storm Water Management Model (SWMM) [16].

Distribution Characteristics of Potential Surface Runoff Based on the Gravity Model
The results of the gravity model showed that there was a great difference in the interaction intensity among sub-catchments. A gravity matrix was established to analyze the relationships among them (Table 4).  After screening, it was found that when the gravity coefficient was 18.93, the potential surface runoff network in 166 sub-catchments broke. Therefore, when the gravity coefficient was less than or equal to 18.93, all sub-catchments were connected. This was the maximum threshold of the connection of all the sub-catchments. Then, we selected gravity coefficient scenarios of 0, 18.93, 30, 40, 50, 100, 300, and 500 and calculated runoff under different gravity coefficient orders, namely, the gravity coefficient. The results showed that the greater the sub-catchment connection strength between catchments, the more prone they were to surface runoff communication. From this outcome, we can infer the actual situations. When rainfall begins, the sub-catchments with the highest gravity coefficients generate runoff first, and with the increase in rainfall intensity, other sub-catchments generate runoff according to this parameter.

Spatial Autocorrelation Analysis of Sub-Catchments
The potential runoff network density (d) in eight scenarios (0, 18.93, 30, 40, 50, 100, 300, and 500 of the gravity coefficient) was calculated (Figure 7). On this basis, a spatial autocorrelation analysis (Moran's I coefficient) of runoff density in sub-catchments was carried out. It was found that (1) in eight scenarios, Moran's I coefficient was greater than 0, and Z > 2.58; p < 0.01, indicating that the distribution of potential runoff density was significantly positively correlated with the distribution of catchments and showed a significant aggregation, with a gradually decreasing trend from the inner central district to the outer city edge. (2) With the increase in the coefficient of gravity, Moran's I presented a trend of first increasing and then decreasing. First, the spatial distribution of runoff density was positively correlated with the gravity coefficient. The sub-catchment areas with high gravity coefficients first generated runoff and showed high aggregation. Second, with the increase in runoff density in the region with a low gravity coefficient, the spatial distribution of runoff density in the whole study area was averaged, and the spatial aggregation decreased. This result not only reflected the spatial aggregation characteristics of runoff but also suggested that the potential surface runoff is mainly distributed in the urban center at the initial stage (Table 5).
This outcome indicated that the impermeable surface area, i.e., the number of buildings and streets in the city center, was higher than that in the surrounding areas, resulting in a high spatial dependence among sub-catchments, and the potential surface runoff connections were more likely to interact with each other. In general, the potential surface runoff density in the central city under the condition of different gravity coefficients was significantly higher than that outside the city center. Further observations showed that the distribution of potential runoff density revealed a significant multicore and multistage sequence distribution pattern, which is consistent with the characteristics of land change in the central urban area.

Advantages of Potential Surface Runoff Simulation
For urban areas, improving the conditions of surface runoff distribution and constructing land expansion are coexisting processes. The resulting orthophoto maps and DSM model was satisfactory and provided the sufficient details to define runoff paths based on elevation. This model was validated by Chi-square analysis of runoff along the urban surface landscape. This has reinforced the conclusions reached about the accuracy of high precision elevation data from Hoffmann et al. [42], Chigbu et al. [43], Sami et al. [17], Hasan [44] and my own research by SWMM [16]. It is worth mentioning that SWMM only requires the permeability coefficient of landscape layouts [45]. Further, our modeling methods consisting of the comprehensive surface landscape information, the orthophoto maps, and DSM's resolution were effective enough to obtain its spatial layouts, i.e., the slope, relief amplitude, roughness, and LUCC of study area. Studies have shown that runoff increases are caused by urban LUCC changes, and high runoff aggregation is mainly concentrated in the Luohe central district. Industrial areas, commercial districts, and highdensity residential areas tend to have a greater impact on runoff. Some research results showed specific local relationships between runoff indices and landscape factors [46]. Controlling landscape pattern processes could change runoff distribution between temporal and spatial; when similar general cases reoccur between sub-catchment areas, similar spatial patterns result. Currently, the increase in impermeable land surfaces caused by urban expansion is an inevitable trend. How to construct urban landscape structure reasonably and effectively, including permeable and impermeable landscapes, will be the key point of urban rainwater resource allocation in the future.
The aggregation distribution of surface runoff in the Luohe central district was closely related to urbanization development. This phenomenon appeared to occur frequently in high-density areas [47]. The permeable landscape is the basis of GI and usually includes green space, farmland, woodland, and water, which can reduce infiltration excess (Hortonian) and the lower the rainfall-excess rates, the longer the surface runoff concentration time and the time of the peak, and the smaller the peak runoff value at any distance [48]. So, effectively improving the collection, transportation, utilization, and regeneration of stormwater requires not only the improvement of municipal pipe networks but also the reasonable and adequate utilization of landscape layouts. Studies of the overall surface landscape can reflect the functional strength of urban stormwater carrying capacity. Surface runoff is significantly influenced by its landscape layout-permeable and impermeable. These effects were treated through the landscape resistance [49]. The effects of the resistance on surface runoff could be negative or positive [50]. With the help of "source" to "sink" identification (from a sub-catchment to others) and resistance surface (MCR model), the runoff transmission paths could be analyzed [51,52]. Thus, the MCR model was employed to analyze the distribution of surface runoff. Through statistical analysis, it was found that the capacity of the Luohe green space to deal with runoff was obviously insufficient; green space accounts for 14.64% of the studied district, but only 5.65% of runoff goes through the green spaces. The length of runoff in green spaces accounts for 7.51% of the total runoff. Regarding the hydrologic response, accurate information regarding the position and layout of the landscape structures was needed to simulate its geometry of surface runoff [17]. The potential runoff spread in different scale of sub-catchments was investigated. Although such data could be easily gathered in small areas, this task could be time-consuming for larger study sites, remaining as a challenge if upscaling of the results is planned [53]. Nevertheless, these methods require a high-resolution surface topography data and satellite images [54]. In our case, it was shown to be accurate enough to explain the relation between landscape layout and runoff distribution. By simulating the distribution of potential runoff, the key areas to be transformed can be identified, and the diversion design can be carried out at the intersection of surface runoff [55]. It is important to take full advantage of the existing permeable landscape. The simulation results can accurately locate the inlet and outlet of green space, guide the transformation of microtopography or landscape, guide the surface runoff into the permeable surface landscape, and increase the utilization rate of permeable green space. Thus, simulating the local hydrological effect of landscape layout could be an important tool used by managers to improve their management efficiency of such measures and for local promotion of their benefit [2,56].

Strategies for the Protection of Suitable Distribution of Urban Stormwater
The organic combination of urban permeable landscapes and impermeable landscapes plays an important role in the stormwater security of urban areas [57]. Different landscapes have different resistances to runoff, and different sub-catchments have different gravity impacts on runoff. Therefore, the values and map of the potential runoff distribution of Luohe provides guidance for urban landscape planning. There are a number of measures that could be considered for disaster reduction [58]. Setting up protection green corridors along the urban trunk road and branch road would be beneficial for organic evacuation of surface runoff. Guided spaces and protected green belts around the urban area next to the river could help prevent the spread of urban flood pollution. Establishment of rain gardens, bio-retention, and permeable areas combined with existing natural resources and even gardens on roofs, which account for 20.09% of the total area of the Luohe central district would improve the quality of life for urban residents and protect the environment. In addition, the organic combination of urban surface landscapes plays an important role in urban security and maintenance of urban biodiversity. Luohe construction continues to grow, showing spatial agglomeration. The central urban area has the most direct impact on the distribution of surface runoff, with the most significant feedback effect. The distribution of surface runoff in Luohe was consistent with the expansion trend of urban construction. Currently, the permeable surface landscape of Luohe city, especially the green space, is almost entirely dependent on the new development of communities and buildings, so in the LUCC map, it is difficult to observe a relatively continuous greenway or green belt. The simulation results can accurately locate the inlet and outlet of green space, guide the transformation of microtopography or landscape, guide the surface runoff into the permeable surface landscape, and increase the utilization rate of permeable green space. Establishing greenbelts on both sides of roads and rivers within the region and strengthening the functional connection of permeable space through integrative green space could increase infiltration, retention, and purification processes of runoff. Building a hydrologic park suitable for public events would increase the recreational functions of permeable space. Removing unnecessary urban structures in central district with connected runoff channels could guide organic evacuation of runoff from urban spaces and avoid growth in submerged areas and low land-use efficiency [59].

Conclusions
The MCR model was used to identify the runoff path from the surface landscape, and we compared the potential runoff between different sub-catchments using the gravity model and spatial autocorrelation from landscape properties in Luohe.
Our results showed that total permeable land (including green space, farmland, and unused land) in the central urban area accounts for 60.98% of the city's total land area, and green space occupies only 14.64%. There is a large and growing area of construction land. To meet requirements for rapid urbanization and to protect suitable ecological areas and green spaces, a reasonable allocation of both areas and a suitable spatial pattern of the two LUCCs are important.
For suitable surface runoff distribution, time must be taken to establish a unified planning of the city's surface landscape. Permeable landscape in urban areas and protection is primary, but still needs scientific urban planning to ensure landscape connectivity and the spatial integrity of critical land. Suitable construction land and rational development are needed along with increasing the infiltration efficiency of land use.
Finally, HD remote sensing was used to obtain surface landscape data and is an efficient method. However, due to the error of LUCC accuracy, the image is still different from reality; selection factors and weight allocation of the resistance surface in the MCR model revealed some subjective elements, and selection factors were not comprehensive enough because of the limitations of information and data collection in Luohe. Therefore, in future studies, it is necessary to improve the experiment and investigate more parameters.
Author Contributions: T.B. and Y.W. conceived of the presented idea and wrote the manuscript, with support from K.B. and J.Z. helped supervise the manuscript. All authors discussed the results and contributed to the final manuscript. All authors have read and agreed to the published version of the manuscript.