Flood-Prone Area Delineation in Urban Subbasins Based on Stream Ordering: Culiacan Urban Basin as a Study Case

: Urban development decreases inﬁltration, increases the runoff velocity, and reduces the concentration times. This situation increases the ﬂood risk in urban watersheds, which represent a management challenge for urban communities and authorities. To increase the resilience of communities due to modiﬁcations of the hydrological cycle produced by climate change and urban development, a methodology is proposed to delineate ﬂood-prone areas in urban basins. This methodology is implemented in an urban subbasin of Culiacan, Mexico, and is based on stream order. A high-resolution digital elevation model was used, which was validated independently through a photogrammetric ﬂight with an unmanned aerial vehicle and ground control points obtained with GNSS (global navigation satellite systems) receivers. Morphometric parameters related to geometry, shape, relief, and drainage network aspects of the subbasin were determined and analyzed. Then, ﬂood-prone area zonation was carried out based on stream-order classiﬁcation and ﬂow direction. Fieldwork was also carried out for the inspection of the sewage network conditions. This methodology simpliﬁes the identiﬁcation of the ﬂood-prone areas in urban subbasins without carrying out complex hydraulic calculations. The present study proposes a simple methodology to identify and delineate the ﬂood-prone areas based on the topographical characteristics and stream order Author Contributions: Conceptualization, A.J.S.-G.; methodology, S.A.M.-A.; validation, A.J.S.-G. and S.A.M.-A.; formal analysis, S.A.R.-G.; investigation, S.A.M.-A. and A.J.S.-G.; data curation, Z.D.M.-F.; writing—original draft preparation, Z.D.M.-F. A.J.S.-G.; writing—review and editing,


Introduction
Hydrometeorological events arise as one of the most important natural hazards. According to Zúñiga et al. [1], more than 50% of the disasters registered in the last three decades worldwide are related to floods. These events produce intense short-duration rainfall or even prolonged-duration rainfall caused by tropical cyclones. They generate water flows that exceed the infiltration capacity of the soil [2].
The dynamic processes of land occupation increase the flood risk in urban areas since the infiltration capacity of soils is reduced. The transition from rural to urban areas is favored by population increase, which promotes the concentration of population in urban environments. However, this transition is carried out without considering the natural hydrological processes present in the original watershed. Hence, mainly the infiltration capacity of subsoil and the natural hydrography is modified [3].
Despite hydrometeorological events contributing positively to increasing the water levels of dams and lakes, they also cause serious problems in urban areas [4,5]. The occurrence of these events in urban areas can generate disruption to community affairs, danger to human life, and health and economic losses. The transformation of precipitation into the runoff is a process that involves factors such as the morphology of the watershed, vegetation covers, and land uses [6]. However, other aspects related to flooding problems in an urban area could be urban solid waste, urban drainage infrastructure, urban vegetation, and housing development [7].
The flood risk assessment has been widely used to support the decision-making process to develop better management and control of urban floods and for the development of plans for urban surface water management [8,9]. However, this experience has been focused on river modeling, while urban flood studies are comparatively scarce [10]. The urban flood risk assessment must include the effect of sewer systems, seepage, infiltration devices, sources and sinks, and the role of rivers in the studied area. Flood risk analysis must reflect various rainfall events with different rainfall intensities and periodicities, where flow velocities and water depth are simulated for selected scenarios of flood with typical recurrence periods [11,12].
From flood risk assessment, flood hazard maps could be generated. Unfortunately, the generation of flood hazard maps is a multi-criterial and complex issue [13,14]. Different approaches based on GIS hazard maps [15], flood frequency [16], flood depth and velocity [17], and remote sensing techniques [18], among others, have been proposed worldwide. Four approaches are recommended to flow modeling in urban areas [19]: one-dimensional (1D) hydrodynamic models, dual one-dimensional (1D/1D) models, and 1D/2D-coupled models simulate flow in minor systems (sewer network) and major systems (a network of surface flow paths and ponds), while Geographic Information Systems (GIS) models provide insight in flow paths and depressions for a certain rainfall amount in the major system [20].
In 1D models, such as the Storm Water Management Model (SWMM), the surcharged flow (major system flow) is kept atop of the manholes and only allowed to be drained back to the sewer once the sewer capacity is made available. This is a limited representation of flow behavior in the major systems since flow is not allowed to move across the surface. Dual one dimensional models, such as MIKE-Urban or Infoworks models, overcome this inconvenience by connecting minor and major systems through flow simulation of weirs, orifices, or sinks [1]. Moreover, 1D/2D models could be considered as the best representation of real urban flow conditions since they couple a 1D sewer model with a 2D overland model [21]. GIS models represent overland flow by flow paths and depths assuming that a portion of precipitation is drained through sewers.
Despite 1D/2D models being the most realistic simulation of urban flow, their calculations require a huge effort on both the data processing and computational demands, which constitutes an important downside that makes them unsuitable for operational management and quick predictions [19,20]. This situation is also observed in 1D/1D models, where a high data demand is required to interpret 2D flow into a 1D simulation. Alternatively, GIS models require less calculation and data resources but commonly are associated with low accuracy and fail to consider the interactions between sewer systems and surface runoffs [19]. However, GIS models only determine flow paths and water depth and it is not possible to determine flooding duration [20].
Despite the advances in flood hazard mapping tools and computational resources, high-quality flood mapping remains challenging for studies with a large analysis scale, with limited evaluation times and scarce data conditions [19]. GIS models can support strategies and activities for flood risk management in regions with limited resources [22]. In this sense, GIS models can make a quick scan of pluvial flooding with limited effort in a few minutes for medium-sized urban areas. These models use GIS capabilities to assess a primary hydrological analysis based on depression extraction (DE) and topographic wetness index (TWI) to predict the surface flow directions and spatial distributions of the accumulated flows.
Since these maps are crucial to develop social resilience and to cope with extreme weather, this study proposes a less complex methodology for the identification and delineation of flood-prone areas in urban basins. The methodology proposed in this paper corresponds to the GIS models. This methodology is based on the delimitation of an urban subbasin obtained from a DEM (digital elevation model) and relevant drainage network aspects, such as stream order and surface flows direction. The drainage network aspects were determined using Geographic Information Systems (GIS), cartography, and UAV photogrammetry. Essential morphometric parameters were calculated to obtain a perspective of the urban subbasin under study with respect to other subbasins. Thus, this methodology simplifies the generation of flood-prone maps without carrying out the hydraulic calculations needed in traditional flood risk studies. Given the simplicity of the methodology, the community could implement it with limited resources and calculation work. In this study, a demonstrative implementation of this methodology is presented in the Culiacan urban area as an urban subbasin case study.

Study Area
This study was carried out in the city of Culiacan, Mexico, which is located between the Universal Transversal de Mercator (UTM) coordinates X:242,000, Y:27,555,000, and X:269,500 Y:2,733,500 in zone 13N. Based on the Risk Atlas of Sinaloa [23], certain sectors of Culiacan city are vulnerable to intense rainfall according to the return period analyses of 2, 5, and 10 years. In these sectors, severe economic losses and human health damages have been registered due to the occurrence of intense rainfall. Figure 1 shows the geographic location of the Culiacan city in Mexico, where the Tamazula, Humaya, and Culiacan watersheds meet. The Humaya and Tamazula rivers flow through the Culiacan urban area and intersect in Culiacan city downtown. The river formed from this confluence takes the name of Culiacan River [24].

Dataset Acquisition
A Digital Elevation Model was used to describe the elevations in the study area [25]. This DEM has a spatial resolution of 5 m at a scale of 1:10,000. Likewise, a vector dataset of the terrain graphic entities provided by INEGI [26] was used. These data include the delimitation of the urban zone in the study area, and the neighborhoods', blocks', and municipalities' shapefiles.
Similarly, satellite images from Google (https:/earth.google.com) accessed on 16 September 2021 and ESRI (http://goto.arcgisonline.com/maps/World_Imagery) accessed on 16 September 2021 were used to discriminate flood-prone areas in the study area. The hydrological and geomorphological analysis was carried out using QGIS software (version 3.16, OPENGIS.ch: Laax, Switzerland).

DEM Validation
A DEM was used to characterize the surface terrain through the mathematical values of height and position of the surface. In this sense, DEM quality must meet the quantitative descriptors of relief, such as elevation, shape, and other topological characteristics such as slope and curvature, among other surface features [27]. Hence, a validation process was carried out to determine the DEM quality obtained from INEGI following the criteria established by the Geospatial Positioning Accuracy Standards [28]. For the validation process, the spatial data of the DEM were compared against another independent

Dataset Acquisition
A Digital Elevation Model was used to describe the elevations in the study area [25]. This DEM has a spatial resolution of 5 m at a scale of 1:10,000. Likewise, a vector dataset of the terrain graphic entities provided by INEGI [26] was used. These data include the delimitation of the urban zone in the study area, and the neighborhoods', blocks', and municipalities' shapefiles. Similarly, satellite images from Google (https:/earth.google.com) accessed on 16 September 2021 and ESRI (http://goto.arcgisonline.com/maps/World_ Imagery) accessed on 16 September 2021 were used to discriminate flood-prone areas in the study area. The hydrological and geomorphological analysis was carried out using QGIS software (version 3.16, OPENGIS.ch: Laax, Switzerland).

DEM Validation
A DEM was used to characterize the surface terrain through the mathematical values of height and position of the surface. In this sense, DEM quality must meet the quantitative descriptors of relief, such as elevation, shape, and other topological characteristics such as slope and curvature, among other surface features [27]. Hence, a validation process was carried out to determine the DEM quality obtained from INEGI following the criteria established by the Geospatial Positioning Accuracy Standards [28]. For the validation process, the spatial data of the DEM were compared against another independent dataset.
The dataset for the DEM validation was acquired independently by using unmanned aerial vehicle (UAV) photogrammetry and ground control points (GCP) set with GNSS receivers (CHC model x900). A DJI Phantom 4 RTK Quad Copter (SZ DJI Technology Co., Shenzhen, China) was used for image acquisition. A Sony RX100 camera with a CMOS sensor and a focal length of 35 mm was mounted on the UAV. This camera acquired images with an effective resolution of 20.2 megapixels. This UAV had a global positioning system (GPS/GLONASS) and a gimbal stabilizer to avoid producing images with distortions.
The flight planning and the image processing algorithms suggested by Mora-Félix et al. [29] were applied to obtain the orthophoto. This orthogonal projection was used to visualize the 29 points marked and distributed on the ground over the study area by accessing the geometric properties (graphic entities) of the cartographic plane.
The root mean square error (RMSE) was calculated to estimate the positional accuracy as suggested by Mora-Félix et al. [29] and Polidori and El Hage [27]. The RMSE is the square root of the average of squared errors between the values observed in the DEM and the values of the coordinates (x, y, z) obtained from GCP according to Equations (1) and (2). The RMSE is used to evaluate the accuracy of the DEM in predicting the coordinates (x, y, z) measured in the field.
where n refers to the total of control points that were measured; x o , y o , z o , correspond to the x, y, z coordinates of the points marked on the ground that were observed in the DEM; x GCP , y GCP , z GCP correspond to the x, y, z coordinates from the independent source measured with GNSS equipment.

Raster and Vector Processing Raster Pre-Processing
After the DEM accuracy was validated, a pre-processing of the DEM was carried out to delineate the subbasins within the study area. The main river course, as well as its tributaries, were also identified at this stage. The pre-processing of the DEM was carried out using a sequential calculation method with specialized tools for hydrological analysis of watersheds. The DEM was processed with the open-source software QGIS 3.16, using the mathematical algorithms of "watershed basin analysis" and the integrated GRASS tool (GRASS GIS 7.4.1, GRASS Development Team: Bonn, Germany).

Raster Processing
First, the study area was delimited using a mask which generated a polygon of the area under study. Subsequently, a filtering of the raster pixel values was applied to eliminate high variations in the different tones that appear in an image. The objective of this filtering is to eliminate errors that are present in the raster image due to the instrument used to acquire the information. These errors could be the presence of areas without information, flat areas, and/or depressions. In this study, a filtering process was carried out to eliminate the sinks in the DEM by filling. These artificial depressions (sinks) were removed using the r.fill.dir algorithm proposed by Wang and Liu [30]. The algorithm generated a depression-free elevation map.

Delineation of the Urban Subbasin
Once the validation and pre-processing of the DEM had been carried out, the urban subbasin of Culiacan was delineated. According to Walesh [31], the delineation of a subbasin modified by anthropogenic activities is determined by surface topography and by the configuration of swales, storm sewers, channels, and other drainage elements. In this study, the delineation of the urban subbasins was carried out based on topographic and hydrographic criteria, where the highest elevations were established as the boundaries of the subbasins, without considering the influence of the surface drainage network. The urban subbasins were delineated using the DEM and the "watershed" algorithm in QGIS, where a minimum cell size was selected to generate the smallest portions of the hydrological basin.

Morphometric Analysis
The physical properties of a subbasin and its streams influence the amount of surface runoff generated. These properties are the result of multiple natural and anthropogenic factors. According to Esper-Angillieri et al. [32], the morphometric analysis includes the study of a set of geometry, shape, and relief aspects of a basin. These physical characteristics are used to describe the hydrological functionality of a subbasin, establish comparisons between subbasins and implement strategies for adequate surface water management. In this study, the characterization process of the hydrological and morphological parameters of the urban subbasin under study was carried out following the traditional methodologies proposed by Horton [33,34], Miller [35], Schumm [36], Gravelius [37], Strahler [38], and Kirpich [39]. Table S1 shows the mathematical equations used to characterize the geometry, shape, relief, and drainage network morphometric parameters of the urban subbasin under study.
The basin geometry morphological parameters were perimeter, area, the maximum length of the basin, length of the mainstream, and basin width. The shape aspects of the basin were obtained, such as form factor, Gravelius compactness, and circularity ratio. The morphometric parameters related to the basin relief, such as maximum and minimum elevations, mean elevation, most frequent elevation, and mean slope of the basin, were also determined. Finally, the stream order and length, drainage density, stream frequency, coefficient of torrentiality, mean slope of the mainstream, and time of concentration were determined. These parameters are related to the drainage network.

Stream Ordering
Flood-prone areas were delineated according to the stream order of the urban watercourses. First, the streamflow directions were identified in the urban subbasin using the DEM. Subsequently, the stream orders of the urban watercourses were determined by using a QGIS tool that implements the numerical algorithm proposed by Garbrecht & Martz [40] based on topographical drainage patterns.

Sewage Network Inspection
Since the watercourses in an urban basin correspond to streets, fieldwork was necessary to contrast the information resulting from the DEM. The maximum and minimum elevations in the terrain were observed during the field visit and the possible water flow direction on the streets was validated. Likewise, the location of the main storm drainage infrastructure was identified by fieldwork.

Identification and Delineation of Flood-Prone Areas
The identification of flood-prone areas in urban subbasins was carried out based on the urban surface topography, the urban watercourses, and the drainage infrastructure located within the urban subbasin under study. Delineation of flood-prone areas was carried out using QGIS software according to stream-order criteria. Higher-order streams were delineated as flood-prone areas in the urban subbasin.  Figure 3 shows a graphical description of the methodology used to extract the UTM coordinates from the DEM using an orthophoto generated with a photogrammetric flight. The orthophoto was georeferenced in the DEM using the marks on the ground in such a way that the UTM coordinates of these marks were extracted from the DEM. Its coordinates were compared with the UTM coordinates obtained by using GNSS equipment at the points marked on the ground. The results of the comparison of the x, y, z coordinates from both sources are shown in Table S2. The accuracy of the DEM was estimated by using the root mean square error (RMSE). Since the GPS surveyed points that were assumed as the true terrain coordinate values of ground control points (GCP), the residuals were used to determine the altimetric (RMSE z ) and planimetric (RMSE xy ) accuracies of the DEM according to Equations (1) and (2).  Figure 3 shows a graphical description of the methodology used to extract the UTM coordinates from the DEM using an orthophoto generated with a photogrammetric flight. The orthophoto was georeferenced in the DEM using the marks on the ground in such a way that the UTM coordinates of these marks were extracted from the DEM. Its coordinates were compared with the UTM coordinates obtained by using GNSS equipment at the points marked on the ground. The results of the comparison of the x, y, z coordinates from both sources are shown in Table S2. The accuracy of the DEM was estimated by using the root mean square error (RMSE). Since the GPS surveyed points that were assumed as the true terrain coordinate values of ground control points (GCP), the residuals were used to determine the altimetric (RMSEz) and planimetric (RMSExy) accuracies of the DEM according to Equations (1) and (2).  Table S2 shows that altimetry revealed the greatest differences, while planimetry remained practically constant for most of the points. The DEM elevations (z) presented a difference variability between a range from −3.83 m to 3.63 m that correspond to GCP 3 and 9, respectively. They showed the greatest elevation difference (|Δz| > 3.6 m), while GCP 10, 22, and 28 presented a moderate difference (|Δz| < 3.0 m). This variation could be related to the measuring times between the coordinates, the presence of slopes, or the interference of the GPS signal. The accuracy obtained for altimetry (RMSEz) is below 3 m, which is within the range established for a representation at a working scale of 1/1000 (10 m in the terrain). Hawker et al. [41] showed that high-accuracy DEMs give better flood estimations. Annis et al. [42] also demonstrated that DEM accuracy obtained from UAVs has a crucial impact on flood dynamics. However, many financial and practical challenges  Table S2 shows that altimetry revealed the greatest differences, while planimetry remained practically constant for most of the points. The DEM elevations (z) presented a difference variability between a range from −3.83 m to 3.63 m that correspond to GCP 3 and 9, respectively. They showed the greatest elevation difference (|∆z| > 3.6 m), while GCP 10, 22, and 28 presented a moderate difference (|∆z| < 3.0 m). This variation could be related to the measuring times between the coordinates, the presence of slopes, or the interference of the GPS signal. The accuracy obtained for altimetry (RMSE z ) is below 3 m, which is within the range established for a representation at a working scale of 1/1000 (10 m in the terrain). Hawker et al. [41] showed that high-accuracy DEMs give better flood estimations. Annis et al. [42] also demonstrated that DEM accuracy obtained from UAVs has a crucial impact on flood dynamics. However, many financial and practical challenges in producing a wide area of high-resolution DEMs are recognized, especially in developing countries [22]. Muthusamy et al. [43] suggest that the highest resolution DEM is not necessary. They used a 30 m resolution DEM for flood modeling, and an overall RMSE range from 2.6 m to 0.9 m was reported in flood depths. Despite that no flood estimations were carried out in this study, the RMSE z obtained in this study is close to the RMSE previously reported in flood depth. In this sense, this could be a cost-effective solution for locations where higher resolution DEMs may not be available [43]. Therefore, the DEM used in this study is suitable for morphological and hydrological analysis in the studied urban subbasin, since the minimum unit to represent is much higher than the RMSE accuracy obtained. In planimetry, a root mean square error (RMSE xy ) of 0.104 m was obtained. The greatest position difference was found at control point 23 with a ∆x < 0.555 m, but ∆y was less than 0.013 m in all GCPs. Based on the above, the DEM provided by INEGI [21] has a good planimetric and altimetric precision since it is below the representation scale of the final project.

DEM Depression Elimination
The depression elimination process consisted of filling the sink in the DEM using the r.fill.dir algorithm of QGIS. The removal of low elevation areas in DEMs is critical since sometimes the flow direction could not follow the terrain slope because of the presence of pools. According to Garbrecht & Martz [40], this step is needed before the use of the flow direction routine. This process was repeated three times until the algorithm generates a depressionless elevation map ( Figure 4). This product is essential to carry out the subsequent hydrological analysis. in producing a wide area of high-resolution DEMs are recognized, especially in developing countries [22]. Muthusamy et al. [43] suggest that the highest resolution DEM is not necessary. They used a 30 m resolution DEM for flood modeling, and an overall RMSE range from 2.6 m to 0.9 m was reported in flood depths. Despite that no flood estimations were carried out in this study, the RMSEz obtained in this study is close to the RMSE previously reported in flood depth. In this sense, this could be a cost-effective solution for locations where higher resolution DEMs may not be available [43]. Therefore, the DEM used in this study is suitable for morphological and hydrological analysis in the studied urban subbasin, since the minimum unit to represent is much higher than the RMSE accuracy obtained. In planimetry, a root mean square error (RMSExy) of 0.104 m was obtained. The greatest position difference was found at control point 23 with a Δx < 0.555 m, but Δy was less than 0.013 m in all GCPs. Based on the above, the DEM provided by INEGI [21] has a good planimetric and altimetric precision since it is below the representation scale of the final project.

DEM Depression Elimination
The depression elimination process consisted of filling the sink in the DEM using the r.fill.dir algorithm of QGIS. The removal of low elevation areas in DEMs is critical since sometimes the flow direction could not follow the terrain slope because of the presence of pools. According to Garbrecht & Martz [40], this step is needed before the use of the flow direction routine. This process was repeated three times until the algorithm generates a depressionless elevation map ( Figure 4). This product is essential to carry out the subsequent hydrological analysis.

Delineation of the Urban Subbasin
One of the most important parameters to delineate a watershed for hydrological analysis is the minimum cell size. This parameter is used by the GIS algorithm to generate the subbasins. This size depends on certain factors such as the objective and scope of the study, the resolution of the digital elevation model, among others [44]. In this work, the subbasins were identified within the Culiacan city urban area. These subbasins were delimited semi-automatically using a surface size of 100 ha. In this sense, the algorithm generated subbasin delimitations between 1 and 100 ha, where the minimum cell size used was 40,000 pixels. Due to the cell size used, 38 subbasins were delimited automatically in the Culiacan urban area. Figure 5 shows some of the subbasins identified in the Culiacan urban zone. The urban subbasin selected as the case of study is contoured in Figure 5. This area partially coincides with the subbasins previously identified and delineated by local

Delineation of the Urban Subbasin
One of the most important parameters to delineate a watershed for hydrological analysis is the minimum cell size. This parameter is used by the GIS algorithm to generate the subbasins. This size depends on certain factors such as the objective and scope of the study, the resolution of the digital elevation model, among others [44]. In this work, the subbasins were identified within the Culiacan city urban area. These subbasins were delimited semi-automatically using a surface size of 100 ha. In this sense, the algorithm generated subbasin delimitations between 1 and 100 ha, where the minimum cell size used was 40,000 pixels. Due to the cell size used, 38 subbasins were delimited automatically in the Culiacan urban area. Figure 5 shows some of the subbasins identified in the Culiacan urban zone. The urban subbasin selected as the case of study is contoured in Figure 5. This area partially coincides with the subbasins previously identified and delineated by local civil protection authorities [23], as shown in Figure 6. This figure also presents the names (codes) of the subbasins identified by local authorities in the study area with their respective drainage area in km 2 . According to these authorities, this urban subbasin generates frequent flooding problems in areas with possible drainage deficiencies. These streams drain water towards a mainstream located at the lowest part of the subbasin.
Sustainability 2021, 13, x FOR PEER REVIEW 10 of 24 civil protection authorities [23], as shown in Figure 6. This figure also presents the names (codes) of the subbasins identified by local authorities in the study area with their respective drainage area in km 2 . According to these authorities, this urban subbasin generates frequent flooding problems in areas with possible drainage deficiencies. These streams drain water towards a mainstream located at the lowest part of the subbasin.    civil protection authorities [23], as shown in Figure 6. This figure also presents the names (codes) of the subbasins identified by local authorities in the study area with their respective drainage area in km 2 . According to these authorities, this urban subbasin generates frequent flooding problems in areas with possible drainage deficiencies. These streams drain water towards a mainstream located at the lowest part of the subbasin.

Analysis of Morphometric Parameters
3.4.1. Geometry and Shape of the Urban Subbasin

Geometry and Shape of the Urban Subbasin
The values obtained for the main morphometric parameters associated with the urban subbasin geometry and shape are shown in Figure 7 and Table 1. The total urban subbasin area is 2.517 km 2 . This value corresponds to the drainage area, which is a key parameter for hydrologic analysis because larger drainage areas lead to greater discharge. The urban subbasin length is 2.8 km, which is related to the maximum length of the subbasin measured parallel to the main drainage line. The values obtained for the main morphometric parameters associated with the urban subbasin geometry and shape are shown in Figure 7 and Table 1. The total urban subbasin area is 2.517 km 2 . This value corresponds to the drainage area, which is a key parameter for hydrologic analysis because larger drainage areas lead to greater discharge. The urban subbasin length is 2.8 km, which is related to the maximum length of the subbasin measured parallel to the main drainage line.  Watershed shape parameters reflect the runoff behavior at the outlet. According to the values obtained for the form factor, the compactness coefficient, and the circularity ratio, the urban subbasin has a slightly elongated shape. This subbasin presents a compactness coefficient (Gc) value of 2.39 which theoretically indicates that this subbasin is prone to generate peak flow rates of low magnitude in response to intense rains. However, multiple factors could modify the hydrologic behavior of the urban subbasin, such as the presence of paved roads, scarcely vegetated areas, blocked drainage systems, and improper drainage structures that obstruct water flows, among others. Due to the presence of these factors, it can be assumed that this urban subbasin behaves as a rounded-shape  Watershed shape parameters reflect the runoff behavior at the outlet. According to the values obtained for the form factor, the compactness coefficient, and the circularity ratio, the urban subbasin has a slightly elongated shape. This subbasin presents a compactness coefficient (Gc) value of 2.39 which theoretically indicates that this subbasin is prone to generate peak flow rates of low magnitude in response to intense rains. However, multiple factors could modify the hydrologic behavior of the urban subbasin, such as the presence of paved roads, scarcely vegetated areas, blocked drainage systems, and improper drainage structures that obstruct water flows, among others. Due to the presence of these factors, it can be assumed that this urban subbasin behaves as a rounded-shape rural subbasin. A rounded watershed shape has been associated with greater flooding risk. Runoff from rounder watersheds tends to reach the outlet more quickly and with greater erosive power and velocity [45,46].

Relief Aspects of the Urban Subbasin
The watershed relief has a great influence on the outflow velocity. Under the same rain conditions, a greater watershed slope generates flows with higher velocities and, in consequence, higher erosion and sediment transport capacities. Figure 8 shows the altitude ranges in the urban subbasin under study. The subbasin was divided by altitude ranges, and then the area was calculated for each of these elevation ranges to analyze the subbasin relief morphometric parameters. In this way, the hypsometric curve and the altimetric frequency histogram were constructed (Figure 9). The hypsometric curve describes the relationship between the altitude and the subbasin area [47]. rural subbasin. A rounded watershed shape has been associated with greater flooding risk. Runoff from rounder watersheds tends to reach the outlet more quickly and with greater erosive power and velocity [45,46].

Relief Aspects of the Urban Subbasin
The watershed relief has a great influence on the outflow velocity. Under the same rain conditions, a greater watershed slope generates flows with higher velocities and, in consequence, higher erosion and sediment transport capacities. Figure 8 shows the altitude ranges in the urban subbasin under study. The subbasin was divided by altitude ranges, and then the area was calculated for each of these elevation ranges to analyze the subbasin relief morphometric parameters. In this way, the hypsometric curve and the altimetric frequency histogram were constructed (Figure 9). The hypsometric curve describes the relationship between the altitude and the subbasin area [47].  Based on the hypsometric curve (Figure 9a), the urban subbasin under analysis could be classified as an early mature stage subbasin. This curve presents a concave shape towards the middle and lower parts, which is a typical characteristic of a highly eroded region. This classification is also evidenced in the altimetric frequency histogram (Figure 9b) since medium to low elevation ranges predominate in the study area. This histogram shows that 50% of the total subbasin area is located above 66 m.a.s.l. The elevation range between 51 and 59 m.a.s.l. concentrates the highest frequencies, but this elevation range only represents 18% of the total area of the subbasin.
The morphometric parameters related to the subbasin relief are shown in Table 2. In this study, a mean stream slope (Smed) value of 2.91% was obtained. This parameter describes the mean surface runoff velocity, which is highly related to the runoff erosion capacity. The value obtained in the present study is characteristic of urban subbasins, where paved areas cause a reduction in the infiltration rate and increase the risk of the existence of flood-prone areas [48]. Sustainability 2021, 13, x FOR PEER REVIEW 13 of 24 Based on the hypsometric curve (Figure 9a), the urban subbasin under analysis could be classified as an early mature stage subbasin. This curve presents a concave shape towards the middle and lower parts, which is a typical characteristic of a highly eroded region. This classification is also evidenced in the altimetric frequency histogram ( Figure  9b) since medium to low elevation ranges predominate in the study area. This histogram shows that 50% of the total subbasin area is located above 66 m.a.s.l. The elevation range between 51 and 59 m.a.s.l. concentrates the highest frequencies, but this elevation range only represents 18% of the total area of the subbasin.
The morphometric parameters related to the subbasin relief are shown in Table 2. In this study, a mean stream slope (Smed) value of 2.91% was obtained. This parameter describes the mean surface runoff velocity, which is highly related to the runoff erosion capacity. The value obtained in the present study is characteristic of urban subbasins, where paved areas cause a reduction in the infiltration rate and increase the risk of the existence of flood-prone areas [48].

Drainage Network Aspects of the Urban Subbasin
In this study, the method proposed by Horton [34] was applied to determine the stream order in the urban subbasin ( Figure 10). According to the stream order analysis,

Drainage Network Aspects of the Urban Subbasin
In this study, the method proposed by Horton [34] was applied to determine the stream order in the urban subbasin ( Figure 10). According to the stream order analysis, the subbasin has 53 first-order streams, 20 second-order streams, 11 third-order streams, and 30 fourth-order streams, with a total length of 17.6 km. the subbasin has 53 first-order streams, 20 second-order streams, 11 third-order streams, and 30 fourth-order streams, with a total length of 17.6 km. The results obtained for the morphometric parameters related to the subbasin drainage network are presented in Table 2. Based on the results and considering the drainage density criteria specified by Sujatha et al. [49], the subbasin under study could be considered as a natural well-drained subbasin. Regarding the stream frequency, the value obtained (Fs = 23.04 streams/km²) is considered high and is related to the natural hydrographic configuration of the watershed since a considerable number of the original streams are currently used as urban streams, which were transformed into streets; however, the total number of streams is finally determined by the urban layout. According to this, the subbasin under study theoretically had a low probability of presenting flood problems due to intense rainfall. However, it was observed that most of the morphometric parameters related to drainage network are meaningless when applied to an urban subbasin, which is modified by anthropogenic activities. This is the case of the urban subbasin under study; this subbasin is well-drained up to the flooding zone throughout the hydraulic network, but it is poorly drained by infiltration. This is because of the predominant presence of impermeable soils mainly constituted by paved streets and urban buildings.  The results obtained for the morphometric parameters related to the subbasin drainage network are presented in Table 2. Based on the results and considering the drainage density criteria specified by Sujatha et al. [49], the subbasin under study could be considered as a natural well-drained subbasin. Regarding the stream frequency, the value obtained (F s = 23.04 streams/km 2 ) is considered high and is related to the natural hydrographic configuration of the watershed since a considerable number of the original streams are currently used as urban streams, which were transformed into streets; however, the total number of streams is finally determined by the urban layout. According to this, the subbasin under study theoretically had a low probability of presenting flood problems due to intense rainfall. However, it was observed that most of the morphometric parameters related to drainage network are meaningless when applied to an urban subbasin, which is modified by anthropogenic activities. This is the case of the urban subbasin under study; this subbasin is well-drained up to the flooding zone throughout the hydraulic network, but it is poorly drained by infiltration. This is because of the predominant presence of impermeable soils mainly constituted by paved streets and urban buildings.
The time of concentration (T c ) can be defined as the time required for a "water particle" to travel from the watershed boundary along the longest watercourse to the watershed outlet. A T c value of 0.72 h was obtained in this urban subbasin which suggests that surface runoff has a short time to contribute to the peak discharge at the subbasin outlet. In addition, the coefficient of torrentiality obtained in this study was high (C t = 21.05 streams/km 2 ). This coefficient indicates that a high number of first-order streams are present in a basin, therefore shorter surface runoff generation times are observed. This situation could explain the conditions for the generation of flood-prone areas in the urban subbasin under study. Figure 11 shows the flow direction analysis within the urban subbasin. The flow direction was determined by using the DEM. Maximum elevations were identified in some streets by analyzing the longitudinal profile. This figure shows that several streets have two opposite surface water flow directions starting in sites of relatively high elevations (ovals marked in white). This morphological analysis helps to identify the first-order streams in the urban subbasin. The flow direction analysis demonstrated that this subbasin maintains its original main river course that drained its territory in natural conditions.

Identification and Delineation of Flood-Prone Areas in the Urban Subbasin
3.5.1. Flow Direction Figure 11 shows the flow direction analysis within the urban subbasin. The flow direction was determined by using the DEM. Maximum elevations were identified in some streets by analyzing the longitudinal profile. This figure shows that several streets have two opposite surface water flow directions starting in sites of relatively high elevations (ovals marked in white). This morphological analysis helps to identify the first-order streams in the urban subbasin. The flow direction analysis demonstrated that this subbasin maintains its original main river course that drained its territory in natural conditions. Figure 11. Identification of flow direction within urban subbasin. Figure 12 shows the stream order classification within the urban subbasin under analysis. In the study area, the mainstream accumulates water from lower-order streams. Therefore, the mainstream course (main channel) meets with third-order and fourth-order streams. This situation is shown in Figure 12, where the end of the third-and fourth-order streams are indicated with orange circles and red diamonds, respectively, and most of them are located in the mainstream course. This relationship is fundamental since the stream order is the basis for the identification of flood-prone areas in the subbasin. However, this study also demonstrated that the stream order in an urban basin can decrease  Figure 12 shows the stream order classification within the urban subbasin under analysis. In the study area, the mainstream accumulates water from lower-order streams. Therefore, the mainstream course (main channel) meets with third-order and fourth-order streams. This situation is shown in Figure 12, where the end of the third-and fourth-order streams are indicated with orange circles and red diamonds, respectively, and most of them are located in the mainstream course. This relationship is fundamental since the stream order is the basis for the identification of flood-prone areas in the subbasin. However, this study also demonstrated that the stream order in an urban basin can decrease rather than increase. This situation can be explained since the streams in an urban subbasin are related to the street layout, which can divide the flow instead of merging it. This behavior differs from what is generally observed in natural basins. rather than increase. This situation can be explained since the streams in an urban subbasin are related to the street layout, which can divide the flow instead of merging it. This behavior differs from what is generally observed in natural basins.  Figure 13 shows the zonation of urban streams based on the stream-order criteria. Four flood risk areas were identified. The first-order areas correspond to the highest elevations of the urban area, where a low possibility of flooding is expected. The secondorder stream areas are considered intermediate or transition zone, and these zones do not represent a high flood risk. The third and fourth stream-order areas receive water from the other lower-order streams. These areas are represented with orange and yellow colors in Figure 13, respectively. The higher-order streams are low elevation zones where a large volume of water is accumulated. Therefore, these areas are more susceptible to floods. According to this zonation, the fourth-order zones are more vulnerable to flood since flooding first occurs in these areas.  Figure 13 shows the zonation of urban streams based on the stream-order criteria. Four flood risk areas were identified. The first-order areas correspond to the highest elevations of the urban area, where a low possibility of flooding is expected. The second-order stream areas are considered intermediate or transition zone, and these zones do not represent a high flood risk. The third and fourth stream-order areas receive water from the other lowerorder streams. These areas are represented with orange and yellow colors in Figure 13, respectively. The higher-order streams are low elevation zones where a large volume of water is accumulated. Therefore, these areas are more susceptible to floods. According to this zonation, the fourth-order zones are more vulnerable to flood since flooding first occurs in these areas. The flood-prone areas identified in this study coincide with those reported by local Figure 13. Stream zonation in the urban subbasin.

Field Inspection of the Sewage Network Conditions
The flood-prone areas identified in this study coincide with those reported by local civil protection authorities, where flooding problems have been registered caused by intense rainfall. According to the Risk Atlas of Culiacan city [23], storm sewer structures installed in the study area do not meet the hydraulic recommendations to efficiently drain it. The fieldwork confirmed the presence of sewers in these zones. Culverts are located at strategic points on sidewalks and the ground level of streets and avenues. However, the reduced dimensions or improper construction of the drainage network do not allow draining the large volume of water accumulated by these streams.
This situation is evidenced in Figure 14, where an open main channel with a trapezoidal cross-section is reduced to a pipe culvert at point "A" (first transition). The trapezoidal cross-section dimensions are as follows: top width 6.50 m, bottom width 2.05 m; flow depth 1.20, side slope 1:1.85, and a freeboard height of 0.25 m. The main channel undergoes a channel transition toward a rectangular section with a width of 2.00 m, in which three circular section ducts are arranged. This is a pipe culvert under a vehicular avenue. The maximum flow area of the main channel with a trapezoidal section is 5.13 m 2 . This area is reduced to three pipes with a diameter of 1.00 m, equivalent to a flow area of 0.785 m 2 for each pipe with a total flow area of 2.36 m 2 . The pipe culvert area represents 46% of the main channel area. The transition to a culvert reduces the channel conveyance, which causes the main channel overflows that in turn generate floods when intense rainfall occurs. The pipe culver discharges in point B to another trapezoidal channel of smaller dimensions (second transition) than the above described for the main channel: top width 5.60 m, bottom width 2.00 m; flow depth 1.15, side slope 1:1.565, with no freeboard height and a flow area of 4.37 m 2 . However, at this point, a bottom rack located in the avenue discharges surface water to this channel. This channel conducts the water flow to point "C", located at the exit of a local university campus. At point "C", a single pipe with a diameter of 1.55 m is located. This is another channel transition (third transition) with a smaller flow area (1.89 m 2 ) due to an incorrect sewer design. This situation causes flooding within the University campus and constitutes a problem that occurs regularly in this area.

Discussion
This manuscript proposes a methodology to identify and delineate the flood-prone areas in an urban subbasin based on stream order. The morphometric characteristics related to the urban subbasin relief obtained in this study coincide with those reported by Ali Shaikh et al. [52]. Anthropogenic characteristics of the urban subbasin, such as the presence of paved soil, reduce the infiltration rate and originate flow currents that are difficult to capture by the storm drainage systems. Therefore, water volume is accumulated in the lower parts of the urban subbasin under study and causes floods. This situation also coincides with Talchabhadel et al. [53], who suggest that the land-use change has a significant impact on the occurrence of floods.
The influence of anthropogenic activities on drainage network aspects is also evidenced by Ercoli et al. [54]. This study found significant changes in the erosive processes in the municipality of Nova Lima, Brazil, due to a combination of factors such as the removal of vegetation cover, the insufficient drainage system for rainwater, the development of luxury housing, among others. This situation is similar to the one reported in the

Discussion
This manuscript proposes a methodology to identify and delineate the flood-prone areas in an urban subbasin based on stream order. The morphometric characteristics related to the urban subbasin relief obtained in this study coincide with those reported by Ali Shaikh et al. [52]. Anthropogenic characteristics of the urban subbasin, such as the presence of paved soil, reduce the infiltration rate and originate flow currents that are difficult to capture by the storm drainage systems. Therefore, water volume is accumulated in the lower parts of the urban subbasin under study and causes floods. This situation also coincides with Talchabhadel et al. [53], who suggest that the land-use change has a significant impact on the occurrence of floods.
The influence of anthropogenic activities on drainage network aspects is also evidenced by Ercoli et al. [54]. This study found significant changes in the erosive processes in the municipality of Nova Lima, Brazil, due to a combination of factors such as the removal of vegetation cover, the insufficient drainage system for rainwater, the development of luxury housing, among others. This situation is similar to the one reported in the present study. The urban subbasin under study is located in the city of Culiacan in Mexico and presents impermeable soils due to urbanization [55]. This situation explains that the urban subbasin under study is well-drained by its streams but poorly drained by infiltration, with high runoff velocity, and is susceptible to flooding.
The identification of flood-prone areas was carried out based on the flow direction and the stream order criteria at the street level. The stream order classification of the urban subbasin was based on Horton's topological stream order where the stream order number increases by one at every confluence [56]. In the study case, the stream order classification was carried out at street level since the original streams of the natural watershed have been modified and mainly transformed into streets. The flow direction and stream order analysis in this work reveals that the Horton criterion is modified in an urban subbasin since the stream order can decrease downstream because of flow division among streets.
The stream order classification identified that the first-order streams tend to be less prone to flooding since they are in higher elevation areas. In this sense, the higher stream order areas (third and fourth stream order) receive water from the lower order streams. These areas are in the lowest parts of the urban subbasins, where, theoretically, sewerage works should be located to drain the surface water runoff. However, fieldwork and local civil protection authorities [23] suggest that there are drainage structures that do not meet the hydraulic recommendations to discharge the water accumulated at lower parts of the urban subbasin. Therefore, this drainage structure eventually collapses due to its incorrect dimensioning and can cause flooding. In addition, the characteristics associated with the size and shape of the subbasin determine the response and velocity of runoff.
Mokarram & Sathyamoorthy [57] used the methodology proposed by Horton [34] to carry out the morphometric and hydrological characterization of an urban subbasin. Conventional morphometric parameters, such as the shape factor, compactness coefficient, and circularity coefficient, were used to indicate that the urban subbasin is not prone to intense rains and is subject to small magnitude floods. However, this study did not consider the urban conditions. In the present study, this perspective was considered to carry out an adequate morphometric characterization in the urban subbasin of Culiacan. The presence of multiple factors such as paved streets, scarcely vegetated areas, and surface runoff obstructions, among others, have a great influence on the runoff generated. This situation completely modifies the behavior of the original watershed in such a way that the urban subbasin must be analyzed differently from that of natural watersheds. However, this study demonstrated that a full basin characterization through the traditional morphologic parameters is not essential to recognize flood-prone areas in an urban basin. Instead, watershed delineation, drainage network layout, flow direction, and stream ordering can be used to recognize flood-vulnerable zones according to their predominant stream order. However, a morphological characterization gives a hydrologic perspective of the urban subbasin under study with regard to other subbasins.
Several studies have been carried out to delineate flood-prone areas. Comprehensive flood risk assessment studies are focused on quantifying the uncertainty in flood hazard estimation and mapping inundation [58,59]. The flood-affected areas have been also mapped using remote sensing satellite data [60]. Despite the importance of flood hazard mapping being well recognized, several problems in elaborating useful flood hazard maps in developing countries are also acknowledged [61]. The classical approaches to flood hazard estimation combine the use of hydrological and hydrodynamic models [62]. However, these studies are limited by high computational costs with massive high-resolution datasets and over large scales. The present study proposes a simple methodology to identify and delineate the flood-prone areas based on the topographical characteristics and stream order classification in an urban subbasin. This methodology could also be used to identify the sewage locations in an urban basin to mitigate the possible impacts of weather events that could cause flooding problems.
Since rainfall is not an input of this method, it is not possible to establish scenarios depending on return periods of precipitation. Hence, it is assumed that the method enlights a scenario of an intense rainfall that surpasses the sewer capacity (minor system) and needs the full capacity of the surface network paths and ponds (major system) to drain the rainfall excess. In this sense, the flood scenario is timeless but space-dependent, and the expected flood extent is governed by the stream order of the urban zone.
Flood-prone areas mapped by 1D, 1D/1D, 1D/2D, and GIS models are different from the GIS model proposed in this investigation since they are essentially delineated according to the flow depths reached on streets and other urban low-elevation facilities. Flow depths vary in a continuous pattern over the flooded surfaces and therefore flood-prone areas are delineated by rather sinuous contours. Conversely, in the method proposed in this work, flood-prone areas do not depend on flow depth. Instead, they are identified by the order of the streams which commonly corresponds to streets. Hence, flood-prone areas are adjusted to the street layout and, in consequence, their contours are primarily angulate.

Conclusions
The present study proposes a methodology to identificate and delineate the floodprone areas in an urban subbasin of the Culiacan city in Mexico based on stream order analysis. Traditional morphological parameters were implemented for morphometric characterization of the urban subbasin. However, some morphometric parameters that are normally applied to natural watersheds could not be used to characterize urban subbasins, since they are anthropogenically modified. According to the proposed methodology, the flood-prone areas in an urban subbasin are highly related to high-order streams. Third and fourth-order streams receive the most water flows from lower-order streams, and are responsible for generating floods that often cause great human and economic damage in the urban area. These flood-prone areas coincide with flood risk zones previously identified by local civil protection authorities.
Accuracy, data availability, time, and calculation complexity are current issues to choose the most appropriate modeling technique for flood mapping. Accordingly, this methodology is a simple tool that contributes to the identification of flood-prone zones and, in consequence, to recognize suitable locations for the development of redesign drainage systems projects. An improved channel conveyance system could be designed to avoid water accumulation in the streets and to reduce flooding problems in the urban subbasin. The proposed methodology can be replicated to the other urban subbasins identified in the Culiacan city or to other urban development projects to prevent flood problems due to the modification of hydrological basins. This methodology could also minimize the time and costs of field surveys when obtaining distances and elevations. In combination with flood simulations, this study could also propose hydraulic works to solve existing flood problems in the study area and to prevent the development of new flood problems. Likewise, this methodology could contribute to the decision-making process for the development and application of structural and non-structural measures to reduce the potential impacts of intense rains and the development of emergency action plans to safeguard human lives and increase resilience to intense rains.