The Effects of Low-Impact Development Best Management Practices on Reducing Stormwater Caused by Land Use Changes in Urban Areas: A Case Study of Tehran City, Iran

Abstract

Urbanization growth and climate change have increased the frequency and severity of floods in urban areas. One of the effective methods for reducing stormwater volume and managing urban floods is the low-impact development best management practice (LID-BMP). This study aims to mitigate flood volume and peak discharge caused by land use changes in the Darabad basin located in Tehran, Iran, using LID-BMPs. For this purpose, land use maps were extracted for a period of 23 years from 2000 to 2022 using Landsat satellite images. Then, by using a combination of geographic information system-based multi-criteria decision analysis (GIS-MCDA) method and spatial criteria, four types of LID-BMPs, including bioretention basin, green roof, grass swale, and porous pavement, were located in the study area. Next, rainfall–runoff modeling was applied to calculate the changes in the mentioned criteria due to land use changes and the application of LID-BMPs in the area using soil conservation service curve number (SCS-CN) method. The simulation results showed that the rise in built-up land use from 43.49 to 56.51 percent between the period has increased the flood volume and peak discharge of 25-year return period by approximately 60 percent. The simulation results also indicated that the combined use of the four selected types of LID-BMPs will lead to a greater decrease in stormwater volume and peak discharge. According to the results, LID-BMPs perform better in shorter return periods in a way that the average percentage of flood volume and peak discharge reduction in a 2-year return period were 36.75 and 34.96 percent, while they were 31.37 and 26.5 percent in a 100-year return period.

1. Introduction

The ever-increasing growth in the world population, especially in cities, and global warming, on the one hand, have caused the growth in urbanization, increasing the impermeability of urban surfaces, significantly reducing the ratio of infiltration capacity, raising the intensity of surface stormwater, decreasing the time of concentration, and declining rainfall [1], and on the other hand, it has caused an increase in atmospheric humidity, evaporation, and transpiration and related changes in rainfall patterns [2]. Climate change, even on a small scale, can increase the intensity of rainfall in urban areas [3]. The physical growth of cities (especially metropolitan cities), rapidly shifting patterns of land use, and demographic increase cause the natural vegetation cover and the capacity of rainfall infiltration to decrease on the surface of the earth; therefore, the risk of flooding in megacities has raised [4,5]. Although most metropolitan cities now use pumping stations and sewage systems to facilitate flood evacuation from urban areas, when the duration and intensity of a flood event exceed the maximum loading capacity of the urban sewage system, the residents of the surrounding areas may suffer from financial and life losses caused by the flood [6,7].

One of the effective ways to reduce the damage caused by floods in urban areas is to upgrade and expand the capacity of drainage systems; however, it has been proven in various studies that increasing the capacity of drainage systems can be costly, unsustainable, and even impractical, especially in dense urban areas [8,9,10]. In addition, drainage systems that have already been designed in urban areas are usually based on old climatic data and do not consider changes in the nature of rainfall over time [11,12]. Therefore, new techniques such as bioretention basin, green roof, grass swale, and porous pavement have been developed to deal with the problem of urban stormwater during floods. All of these techniques are called low-impact development best management practices (LID-BMPs) [13]. LID-BMPs, as innovative approaches in environmental management and urban improvement, play a vital role in strategic planning and urban sustainability. These strategies, including permeable pavements, green roofs, rain gardens, and other similar solutions, significantly enhance groundwater recharge, reduce runoff [14], and mitigate urban heat islands while integrating ecological [15] and social equity into urban design [16]. They strengthen cities’ resilience to climate change by effectively managing stormwater, reducing flooding risks [17], and benefiting underserved areas through equitable urban development [18]. Additionally, these practices improve public health by creating green spaces, lower long-term infrastructure costs by reducing reliance on traditional systems, and foster sustainable economic growth aligned with long-term development goals [19]. Although LID-BMPs offer numerous advantages, their implementation is not without challenges. In some cases, the costs of extending stormwater systems with devices aimed at restoring natural water conditions can exceed those of traditional drainage systems [20]. Additionally, space limitations, especially in densely populated urban areas, may hinder their development [21]. Factors such as soil permeability and potential contamination risks also significantly influence the feasibility of these systems. Addressing these considerations ensures the effective and sustainable integration of LID-BMPs into urban environments.

Overall, LID-BMP serves as a multifaceted tool within strategic urban planning frameworks, providing environmental, social, and economic benefits that are essential for sustainable development. LID-BMP is one of the stormwater management strategies to increase permeability, reduce surface stormwater, and reuse water [22], which has been implemented in recent years to boost urban resilience and diminish flood risk [23,24]. In various studies, the efficiency of LID-BMPs has been reported as an effective method to lessen the effects of urban floods [2,18,25,26,27]. In other words, these studies show that different LID-BMP strategies can have different benefits for urban ecosystems in different geographical locations. In addition to decreasing the effects of urbanization on the basin’s hydrology, LID-BMPs increase the level of groundwater, reduce stormwater, and preserve the natural environment. For instance, a study concluded that a green roof can reduce 60 to 70% of the volume of stormwater compared to conventional roofs [23]. Similarly, an experiment conducted found that green roofs can reduce stormwater volume by 30–78% in comparison with conventional roofs, in addition to reducing discharge time [28]. However, the optimal combination of LID-BMP measures varies in different severity zones, highlighting the need for a balanced approach to achieve sustainable development in sponge cities [29]. LID-BMPs can also reduce or eliminate the need to make ponds that are created to store rainwater in urban areas and usually occupy a large area from the occupancy level [2]. They provide the possibility of high-density urbanization with the least negative hydrological effects, and as a result, they can be considered as a positive measure in urban development in the direction of infrastructure development [30]. Another important feature of LID-BMPs is their affordability. For example, the economic importance of LID-BMPs was investigated in Hong Kong by presenting a cost–benefit model and concluded that USD 1.2 billion and 5.3 billion of environmental and economic benefits have been saved in a 30-year period [31]. According to the mentioned studies, LID-BMPs can be considered as an alternative to the old and traditional ways of managing stormwater at the city level and provide suitable solutions for the coping capacity against climate change and the sustainable development of cities.

Tehran, receiving considerable rainfall, is located at the foot of the Alborz Mountain range. Many rivers drain the water from the rain and transfer it to the downstream plains. The general direction of the rivers and canals is mainly towards the urban area, and this has caused the residents to witness waterlogging of the streets due to torrential rains and floods in case of heavy and continuous rains. At the same time, the urban area of Tehran is the passageway of many channels and rivers that are responsible for the water drainage of the upstream basins and also the collection of rainwater in the urban area. In some parts of the city where the cost of land is high, beds of rivers and canals inside the city have been encroached upon and their cross-section has been reduced and limited. This issue can affect the natural regime and flow of the river and cause irreparable damages in the time of heavy rainfalls. Therefore, the purpose of this study is to use LID-BMPs to compensate for the decline in permeability of the surfaces due to the rise in urban land use and the escalation in stormwater due to rainfall in a basin located in the northeast of Tehran. Among the innovations of this study is an integrated investigation of the cause (change in land use) to the solution (use of LID-BMPs) for urban stormwater management all in one study. Another novelty of this study is the use of a combination of the GIS-MCDA method and spatial criteria for the site selection of LID-BMPs. In addition, eight criteria, including flooding potential, were used to select the optimal locations for installing the LID-BMPs. In this study, four proposed scenarios of the use of bioretention basin, green roof, grass swale, and porous pavement are investigated.

2. Materials and Methods

2.1. Study Area

The Darabad basin, situated in the northeast of Tehran (35°49′31″ N, 51°30′12″ E), covers an area of 3987.39 hectares and receives an average annual rainfall of 380 mm. Tehran, the capital of Iran, is a highly populated city, with a population of 15.7 million, accounting for one-fifth of the country’s total population, based on the 2016 census [32]. The northern area of Tehran city is connected within the Alborz Mountain ranges and has a length of about 40 km. The Darabad river originates from the Alborz Mountain ranges at an altitude of 3160 m above the mean sea level (amsl), and then, it is directed toward the inner city channels, and after the city sewage enters it, it finally terminates in the Kavir plain at an elevation of about 1700 m amsl. The lack of proper land use planning in the region can be deduced from the following: i. soil erosion, ii. inadequate urban drainage system, iii. building urban settlements up to the boundary of the protected area, iv. The poor management of public green spaces, and v. the absence of LID-BMPs and BMPs for managing stormwater [33]. Human activities along with natural influences (intense extreme events) have aggravated the vulnerability of the basin [34]. In recent years, as urbanization progresses, the basin has witnessed the transformation of large areas of gardens, grass land, and bare land into construction land, which has increased the amount of impervious area. A display of the study area is shown in Figure 1.

2.2. Data

The materials used in this study mainly consist of land use data, meteorological data, and basic geographic information. In terms of land use data, two sets of temporal multispectral satellite images of 30 m resolution for the years 2000 and 2022 were acquired from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/, accessed on 24 December 2024). The average annual rainfall and the average maximum rainfall of 24 h of the Darabad basin were obtained from I.R. of Iran Meteorological Organization (https://www.tehranmet.ir, accessed on 24 December 2024). In addition, the maximum instantaneous discharge statistics of Qalak rain gauge, which is located inside the basin, were acquired from Regional Water Company of Tehran (https://thrw.ir/, accessed on 24 December 2024). The digitized topographic map of the basin with a spatial resolution of 5 m was taken from the Municipality of District 1 of Tehran city.

2.3. Methodology

The general procedure carried out in this research is depicted in a flowchart in Figure 2.

2.3.1. Land Use Land Cover (LULC) Change

Recognizing LULC is a crucial component of environmental analysis, planning, and management [35], since these changes are regarded as one of the most significant environmental problems facing the world today [36]. Meanwhile, satellite imagery has been widely utilized to track environmental degradation and changes in LULC over time [37]. The time range of the Landsat series of satellites, from 1975 to the present time, has made it the best option for studying urban physical development [38]. The details corresponding to the satellite data are given in

Table 1. Details of the satellite data.

Satellite	Sensor	Spectral Bands	Month/Year	Path/Row
Landsat 7	ETM+	8	July 2000	164/35
Landsat 8	OLI/TIRS	11	July 2022	164/35

. The screening and processing of enormous amounts of remote sensing data are essential in remote sensing technology-based research and remote sensing image analysis [39]. After applying radiometric and atmospheric correction to the images, a total of five LULC classes, i.e., built-up (residential, commercial, industrial, and traffic land), barren land (bare soil without any cover), vegetation (parks and gardens), rangeland (shrubs and grasses), and ridge (mountainous lands), were identified. The LULC classification and mapping process includes the post-classification accuracy assessment, which is used to evaluate the accuracy of the classified maps. The accuracy of a classified map should not be less than 80% for accurate interpretability and identification [40]. Utilizing 400 randomly chosen points, the Kappa coefficient technique is employed in this study to assess the precision of LULC maps of the Darabad basin. In the selection of the points, it was performed in a way that covers all the study area and all the classes. Since field data were not available for 2000, the ground observations for this year were acquired from the Google Earth Pro platform. Additionally, for 2022, the ground observations were obtained from a combination of field visits and utilizing the Google Earth Pro platform. The accuracy assessment results indicate that the overall accuracy level was 93.56 percent and 96.80 percent for the years 2000 and 2022, respectively. The corresponding Kappa statistics were 0.8884 and 0.9441, respectively.

2.3.2. Hydrology

Hydrological Divisions of the Area

In order to extract more accurate and uniform information, the study area was divided into 7 hydrological units and 1 non-hydrological unit. A non-hydrological unit is a basin that does not have an independent stream, and water enters it from adjacent basins. The division of sub-basins along with the digital elevation model (DEM) of the Darabad basin is shown in Figure 3.

Longitudinal Profile of the Main Stream and Slope Calculation

Knowing the slope of different points of a river’s path can give a clear picture of the river’s destructive power in different places. In addition, in order to calculate the time of concentration (TOC) of streams in the sub-basins, it is necessary to calculate their slope first. For this purpose, the longitudinal profile of streams was calculated and drawn separately for each sub-basin. The longitudinal profile of the Darabad river is obtained by transferring the length of the river on the X-axis and the height of its different points on the Y-axis and with the help of contour lines from topographic maps. In drawing the longitudinal profile of the main river of the study area and its hydrological units, the digitized topographic map (with a spatial resolution of 5 m) of the area was used. According to the profile, sub-basin A has the highest and sub-basin G has the lowest slope with 1.09 and 36 percent, respectively.

Calculation of Stormwater Using the Soil Conservation Service Curve Number (SCS-CN) Method

One approach for calculating the direct stormwater depth for a rainfall event is SCS-CN. Equations (1) and (2) are used in surface stormwater prediction equations based on the SCS-CN approach [41]. CN is determined based on the soil type, land use and land cover, hydrological conditions, and antecedent moisture conditions of the soil acquired from the SCS [42]. In this way, after preparing the map of the hydrological groups of the basin and integrating it with the land use map, the value of CN and then its weighted average were calculated using Equations (1) and (2).

$${\mathcal{Q}}_{SCC-CN}={\displaystyle \frac{{(P-0.2S)}^2}{P+0.8S}}$$

(1)

$$S={\displaystyle \frac{25400}{C_N}}-254$$

(2)

where $\mathcal{Q}$—direct stormwater (mm); $P$—total is precipitation (mm); $S$—maximum retention capacity obtained as per the SCS-CN method (mm); and $C_N$—curve number that depends on the soil type [43].

TOC

One of the most pivotal and widely used physical parameters of a basin, especially in hydrological equations, is the TOC. TOC is the interval of time between the centroid of rainfall excess and the recession’s inflection point for direct stormwater [44]. The time needed for a stormwater particle to move from the watershed boundary that is hydraulically farthest away along the longest watercourse to the basin outlet is the physical definition of TOC [45]. The TOC in each basin depends on several factors, including the slope of the basin, the slope of the waterway, the geological structure, the type and amount of vegetation, the intensity of rainfall, and the amount and distribution of rainfall in time and place. This study employed the Kirpich method (Equation (3)) for calculating the TOC.

$$T=0.0663L^{0.77}{S}_o^{-0.305}$$

(3)

where $T_c$ is time of concentration by the Kirpich method in hours, $L$ is the main channel length in kilometers, and $S_o$ is the main channel slope in meters per meters [46]. In this study, the watershed was divided into eight distinct sub-basins to enhance the precision of the hydrological analysis. For each sub-basin, both the basin slope and the main channel slope were computed. The results revealed minimal differences between the two parameters within most sub-basins, with the exception of sub-basin A (

Table 2. Comparison of the average channel slope and basin slope for sub-basins.

Sub-Basin	A	B	C	Dint	E	F	G	H
Basin slope (%)	47.79	41.11	45.03	20.21	24.14	6.30	4.31	13.16
Channel slope (%)	44.93	39.87	43.67	19.30	23.34	6	3.93	12.38
Difference (%)	2.86	1.24	1.36	0.91	0.80	0.30	0.38	0.78

). As a result, the main channel slope was chosen as the primary parameter for calculating the specific TOC for each sub-basin. This approach was deemed both practical and representative of the flow dynamics within the sub-basins, reflecting the key characteristics driving water movement.

2.3.3. Meteorology

In order to implement the SCS-CN rainfall–runoff model, the maximum rainfall at the TOC is required, and its calculation requires the average annual rainfall of the basin and the average maximum rainfall of 24 h. The minimum, average, and maximum annual rainfall statistics of the basin were calculated from three meteorological stations around the basin from 1990 to 2022 and can be seen in Figure 4.

Rainfall Gradient in the Area

In a general sense, the gradient represents the change in a phenomenon. In the study of rainfall gradient, a logical relationship is established between rainfall as the dependent variable and altitude or longitude and latitude as independent variables, to estimate the amount of rainfall in different parts of the basin. To determine the rainfall gradient in the area using meteorological station data, the relationship between rainfall and altitude in the basin and sub-basins was investigated. Due to the absence of stations within the study area, data from the nearest available stations were used. These stations, as shown in Figure 5, are located at varying elevations. The presence of meteorological stations at different elevations helps improve the accuracy of analyses, covers diverse areas, and enhances hydrological modeling. Although the trend line is based on three points, these points were selected as the most representative data available for the study area to ensure the reliability of the analysis within the given constraints. The results are presented separately below. Using the 33-year average rainfall statistics, the relationship between rainfall and altitude in the study area was analyzed, and the rainfall gradient was finally obtained. In Figure 5, the relationship and correlation diagram between rainfall amount and the elevation of meteorological stations are shown. The coefficient of determination obtained for the rainfall gradient equation in the study area is 0.95, which is acceptable.

Average Maximum 24-Hour Rainfall in Different Return Periods

This parameter, which is very important in the design of structures and calculations related to the maximum instantaneous discharge, can be analyzed using the data recorded in nearby meteorological stations. Kolakchall (28 years), Niavaran (30 years), and Shahid Abbaspur (24 years) stations were the stations with the data of average maximum 24-hour rainfall around the Darabad basin, and all their data were used in this research. Using the EasyFit 5.6 software, the average maximum 24-hour rainfall was calculated in different return periods. EasyFit is a software tool used for fitting probability distributions to data. It provides a wide range of statistical distributions and helps in selecting the best-fitting distribution based on the data. In

Table 3. Selected statistical distributions used for calculating average maximum 24-hour rainfall (along with distribution parameters) for index stations.

Meteorological Station	Average Annual Rainfall	Statistical Distribution
Kolakchall	511.7	Wakeby (ξ: 189.15; α: 11.046; β: 8.697; γ: 0.19457; δ: 16.588)
Niavaran	243.3	LogLogistic (α: 3.7741; β: 19.409; γ: 7.1402)
Abbaspour	334	LogLogistic (α: 2.558; β: 15.643; γ: 17.546)

, the selected statistical distributions of average maximum 24-hour rainfall for the meteorological stations are given. The Wakeby and LogLogistic distributions are statistical models often used in hydrology and other fields to analyze extreme events like floods or rainfall. The Wakeby distribution is a five-parameter model defined by ξ (location), α (shape), β (shape), γ (scale), and δ (shape). The LogLogistic distribution, also used for skewed data such as rainfall or survival times, has parameters such as α (scale), β (shape), and γ (location). In EasyFit software, these parameters are derived from dataset and displayed with their values to help model the data accurately. In

Table 4. The average maximum 24-hour rainfall (mm) in different return periods for the selected stations along with the gradient of annual rainfall frequency.

Meteorological Station		Average Maximum 24-Hour Rainfall with Return Period (Year)
		2	5	10	25	50	100
Kolakchall		40.2	50.2	59.0	72.6	84.7	98.5
Niavaran		26.6	35.2	41.9	52.2	61.6	72.7
Abbaspour		33.2	44.4	54.5	71.7	89.2	111.8
Rainfall gradient	R2	0.99	0.95	0.90	0.78	0.65	0.50
	Gradient of a line	0.0491	0.0526	0.0584	0.0663	0.0708	0.0724
	Intercept	15.47	24.17	30.58	41.44	52.78	68.06

, the average maximum 24-hour rainfall (mm) in different return periods for the selected stations along with the gradient of annual rainfall frequency has been calculated.

Analysis of Short-Term Rainfall and the Intensity, Duration, and Frequency (IDF) Curve

The intensity of rain in a rain storm is the amount of rain per unit of time, which is usually expressed in millimeters per hour or inches per hour. In general, the shorter the duration of the rainfall, the greater its intensity and vice versa. Rainfall intensity–duration equations change with frequency or return period. The longer the return period, the more intense rains are expected. Due to the lack of suitable stations containing information near the study area, IDF diagrams in this report were drawn using experimental and calibrated formulas. Risk management requires a precise and trustworthy estimation of the amount of rainfall indicated by the points on the IDF curves [47]. By analyzing the relationship between the intensity, duration, and frequency of observed rainfall data, it is possible to estimate the intensity of extreme rainfall events with a particular duration and return period [48]. There are diverse equations to calculate the intensity of short-term rainfall with different return periods, and in this report, Ghahraman and Abkhezr method [49], which is a generalized Bell method for Iran, has been used. In this method, Iran is divided into six regions, and for each of these regions, specific equations are provided to calculate one-hour rainfall with a return period of 10 years. Based on this zoning, the Darabad basin is located in District 2. The relation between one-hour rainfall calculation and 10-year return period in this area is given in Equation (4):

$$P_{60}^{10}=9.99+0.212P_{1440}^2$$

(4)

where $P_{60}^{10}$ is 60-minute rainfall with a return period of 10 years, and another parameter of the equation is the average maximum daily rainfall in the study area. To determine the amount of rainfall in the return period and the duration of different rainfall, Equations (5) and (6) are used for the duration of rainfall of less and more than 2 h:

$$P_t^T={0.1372t}^{0.4778}\left[0.4608+0.2349Ln(T-0.62)\right]P_{60}^{10}$$

(5)

$$P_t^T={0.2009t}^{0.3937}\left[0.5565+0.1948Ln(T-0.8)\right]P_{60}^{10}$$

(6)

Equation (5) is used for rainfalls of up to 2 h, and Equation (6) is used for rainfalls of more than 2 h. In the above equations, $t$ is the duration of rainfall in minutes, $T$ is the return period in years, and $t$ is minutes of rainfall with a return period of $T$ years. After calculating the TOC using the Kirpich method and placing it in the above equations, the amount of rainfall and its intensity in different return periods were obtained separately for each sub-basin. The diagram related to the IDF of the basin is given in Figure 6.

Unit Hydrograph of Sub-Basins by SCS-CN Method

One of the most important steps in hydrological analysis is the calculation of the unit hydrograph. In this research, the SCS-CN empirical method and the following relations are used to calculate the unit hydrograph.

$$\left\{\begin{array}{c}{Q}_p={\displaystyle \frac{0.2083AQ}{T_p}}\\ T_p={\displaystyle \frac{D}{2}}+T_i\\ D=0.1333T_c\\ T_i=0.6T_c\\ T_b=2.67T_p\end{array}\right.$$

(7)

In above equations, $Q_p$ is the peak discharge of the unit hydrograph (m³/s), $A$ is the area of the basin (km²), $Q$ is the surface stormwater height (mm), $T_p$ is the time to reach the highest point of the unit hydrograph (h), $D$ is the duration of effective rainfall of unit hydrograph (h), $T_l$ is the lag time of the basin from the center of gravity of rainfall to the peak of the hydrograph (h), $T_b$ is the base time of the unit hydrograph (h), and $T_c$ is the time of concentration of the basin (h) [50].

To obtain unit hydrograph dimensions, the value of TP and QP must be multiplied in the SCS-CN dimensionless table; in this way, the dimensions of the hydrograph corresponding to each hydrological unit can be obtained. The dimensions of the unit hydrograph of each sub-basin, as well as the entire basin, are given in

Table 5. Unit hydrograph dimensions of sub-basins for one centimeter of stormwater by SCS-CN method.

Sub-Basin	Area (km2)	Tc (h)	D (h)	Ti (h)	Tp (h)	Qp (m3/s)
A	19.81	0.66	0.09	0.40	0.44	9.31
B	0.96	0.15	0.02	0.09	0.10	1.98
C	5.13	0.25	0.03	0.15	0.17	6.37
Dint	3.24	0.38	0.05	0.23	0.26	2.63
E	3.99	0.45	0.06	0.27	0.30	2.75
F	3.78	0.62	0.08	0.37	0.41	1.91
G	1.40	0.54	0.07	0.32	0.36	0.81
H	1.57	0.19	0.03	0.11	0.13	2.60
Darabad basin	39.87	0.93	0.12	0.56	0.62	13.39

. Using the dimensions of the unit hydrograph obtained by the SCS-CN method and also the amounts of rainfall excess for different return periods, the dimensions of the flood hydrograph and the maximum values of the flood discharge were estimated.

2.3.4. Selection of Suitable LID-BMPs

There are many types of LID-BMPs and possible sites for their installation, and the process of selecting LID-BMPs is a complex one, given the diversity in site conditions, performance, and cost of LID-BMP [51]. In this research, the multi-criteria index system (MCIS) method [52] was used to select LID-BMPs. The evaluation process was conducted based on three main criteria, i.e., site suitability, stormwater control benefits, and cost and maintenance, along with a set of sub-criteria. The sub-criteria under site suitability included land availability, groundwater characteristics, site conditions, slope, and soil type, while stormwater control benefits were evaluated based on runoff reduction efficiency, water quality improvement, and flood risk mitigation. The cost and maintenance criterion included construction costs, maintenance requirements, and lifespan.

To ensure a robust evaluation, a panel of experts, including water resource engineers, urban planners, environmental civil engineers, GIS experts, and environmental engineers, was formed to score each LID-BMPs against the defined criteria and sub-criteria. A range of potential LID-BMPs was initially considered, including constructed wetlands, infiltration basins, wet detention ponds, dry detention ponds, bioretention basins, green roofs, grass swales, and porous pavements. However, certain options, such as constructed wetlands, infiltration basins, wet detention ponds, and dry detention ponds, were excluded from the final list due to the need for a large area, which is practically not feasible in the Darabad basin due to the high price of land.

The scoring process was performed using a weighted average approach, where weights were assigned to each criterion based on their importance to the study objectives. As a result, four LID-BMPs with the highest scores were selected as suitable options for the study area: 1—bioretention basin; 2—green roof; 3—grass swale; and 4—porous pavement.

2.3.5. Optimal Location Identification for LID-BMPs

In order to determine the optimal locations for LID-BMPs, 8 criteria of slope, land use, rainfall, flooding potential, groundwater level, distance from the waterway, distance from the street, and distance from the fault line were used. Since the unit of measurement of each of these criteria and their valuation range is different, it is necessary to compare all the criteria with the same unit of measurement; fuzzy logic was used for this purpose. A fuzzy logic system is distinctive in that it can manage linguistic information and numerical data at the same time. To use fuzzy logic in ArcGIS, it is necessary to assign a weight to each criterion (layer), and the greater weight indicates the greater importance of that criterion. The assigned weights were extracted by the analytical hierarchy process (AHP) using Expert Choice 11 software. To derive these weights, a structured questionnaire was designed, which was completed by 25 experts specializing in water resource engineering, urban planning, environmental civil engineering, GIS, and environmental engineering. The experts were selected based on their professional experience and familiarity with LID-BMP implementation in urban basins.

Each expert provided pairwise comparisons of criteria and sub-criteria based on their judgment, and the final results were obtained by aggregating these individual judgments using the geometric mean method. This approach ensures a balanced and representative evaluation, reflecting the collective expertise of the panel [53]. The consistency ratio (CR) was also calculated to confirm the reliability of the judgments, with all values falling below the acceptable threshold of 0.1, indicating consistent responses [54,55]. To assign optimal locations to the LID-BMPs using the fuzzy-AHP method, it is necessary to fuzzify each criterion according to its characteristics, with fuzzy membership functions. It should be noted that the type of membership function for each criterion has been determined using a questionnaire and experts’ opinions.

Table 6. Fuzzy membership functions and weights of criteria.

Criterion	Fuzzy Function		Weight
			Bioretention Basin	Green Roof	Grass Swale	Porous Pavement
Slope	Linear (decrease)	Min: 20 Max: 0	0.217	0.097	0.243	0.297
Groundwater level	Linear (increase)	Min: 1.22 Max: 97.3	0.062	0.04	0.138	0.115
Distance from waterway	Linear (increase)	Min: 0 Max: 2220	0.04	0.041	0.057	0.035
Distance from street	Linear (decrease)	Min: 6556.7 Max: 0	0.122	0.045	0.303	0.224
Distance from fault line	Linear (increase)	Min: 0 Max: 5.790	0.083	0.068	0.064	0.057
Flooding potential	Linear (increase)	Min: 0.66 Max: 1.28	0.073	0.14	0.053	0.078
Rainfall	Linear (increase)	Min: 432 Max: 555.12	0.032	0.231	0.034	0.056
Land use	User defined		0.371	0.338	0.108	0.138

shows the fuzzy functions and weight of each criterion.

For the land use criterion, since 4 different types of LID-BMPs have been examined in this research and the priority of use for each of them is different, the land use scoring map for each LID-BMP was given separately. It should be noted that a land use type with a lower priority number in

Table 7. Locative priority for each type of LID-BMP.

Land Use	Location Priority
	Bioretention Basin	Green Roof	Grass Swale	Porous Pavement
Public services	5	2	4	4
Pertaining to the police	10	10	10	10
Residential garden	3	5	5	9
Residential	4	1	6	8
Green space	1	-	2	1
Recreational	2	-	3	2
Religulous	5	4	7	4
Parking	6	6	7	3
Bare land	6	-	1	-
Rocky land	5	-	5	-
Commercial–industrial–therapeutic–infrastructure	8	3	9	5

has a higher number of points in the land use scoring map. To prepare Table 7, a comprehensive questionnaire was designed and distributed among experts in water resource engineering, urban planning, environmental civil engineering, GIS, and environmental engineering. The responses were analyzed using statistical methods to extract the priority of different land use types for each LID-BMP. The resulting priorities were then standardized on a scale from 1 to 10, where a higher score indicates better compatibility of a specific land use type with the corresponding LID-BMP. This table serves as the basis for fuzzifying the land use layer, with each priority level reflecting its contribution to the overall suitability analysis.

After preparing the fuzzy layers of criteria, it is necessary to multiply them by their obtained weights using the AHP method to obtain weighted fuzzy layers. According to the different functions of LID-BMPs, the relative weight of each criterion is also different in these systems. For example, for grass swale, the distance from the street is the first priority, and then, the slope criterion has the highest relative weight, while for the bioretention basin, due to the possibility of using the stored stormwater in the irrigation of green spaces or the car wash, the land use criterion holds the highest importance [56]. Table 6 shows the fuzzy membership functions and the resulting weight for each criterion separately for 4 types of LID-BMPs.

To fuzzify the land use layer, since the priority of each kind of land use is different for each type of LID-BMP, the User Defined function is used in the TerrSet 19.0 software. TerrSet is an integrated geospatial software designed for land planning, monitoring, and analysis, widely used for environmental management and urban planning applications [57]. The User Defined function in TerrSet allows for customized fuzzy membership functions to be applied, ensuring the accurate representation of expert-derived priorities.

In Table 7, the priority of each kind of land use for each LID-BMP is shown, which was extracted using the aforementioned questionnaire and expert opinions. The values in this table range from 1 to 10, where a higher number indicates a better compatibility between the land use and the corresponding LID-BMP.

In the next stage of locating the LID-BMPs, it is necessary to overlay the weighted fuzzy layers of the criteria using one of the fuzzy overlay operators. The common areas of the 8 criteria will indicate the appropriate place for placing the LID-BMP. In this research, gamma operator is used to overlay the layers. The gamma operator, applied in this study, is a powerful fuzzy aggregation operator used to combine multiple criteria layers. By adjusting the gamma value, it balances the trade-off between the minimum (conservative) and maximum (optimistic) aggregation methods, providing a more flexible approach to decision making [58,59]. This operator was instrumental in generating the final suitability maps, as it ensured that the specific characteristics and priorities of each LID-BMP type were appropriately integrated.

2.3.6. Modeling of the Selected LID-BMPs

Since the SCS-CN rainfall–runoff model calculates the infiltration loss using CN, LID-BMP can be defined as a permeable sub-basin (with a change in the infiltration rate) such that the rainfall excess in this unit is less compared to other units and its impact is observed in the outflow. In this method, each modeled sub-basin consists of several different parts of the CN (places with different permeability), and by calculating the new CN obtained cumulatively, the amount of new stormwater is calculated. The estimation of the new CN for LID-BMP is also based on its type as well as the review of previous research. A study by Jia, Yao, Tang, Yu, Zhen, and Lu [51] was conducted on the development of an MCIS to select the best urban stormwater management practices and concluded that the LID-BMP was effective in controlling the quantity of stormwater as shown in

Table 8. LID-BMP’s stormwater quantity control effectiveness.

LID-BMP	Stormwater Quantity Control
	Stormwater Volume	Peak Flow	Flow Rate	Total
Bioretention basin	16.25	18.3	14.3	48.75
Green roof	14.75	16.8	14	45.5
Grass swale	15.75	17.8	13.8	47.25
Porous pavement	16	18.5	13.3	47.75

. The LID-BMPs used in this research, which include bioretention basin, green roof, grass swale, and porous pavement, are all effective in reducing the quantity of rainfall.

3. Results

3.1. Spatio-Temporal Land Use Change

The analysis of land use changes using remotely sensed data helps us gain a deeper understanding of the historical interactions between human activities and the environment. It also enables us to predict potential landscape conversions, which in turn can aid in the creation of realistic and multi-dimensional scenarios for developing sustainable environmental policies and management practices. As it is shown in

Table 9. Changes in all land use classes during a period of twenty-two years in hectare (ha).

Year	Built-Up		Barren Land		Vegetation		Rangeland		Ridge
	ha	%	ha	%	ha	%	ha	%	ha	%
2000	662	43.49%	262	67.17%	357	52.19%	1332	52.44%	1380	48.43%
2022	860	56.51%	128	32.83%	327	47.81%	1208	47.56%	1469	51.57%
Total change	+198	+13.02%	−134	−34.34%	−30	−4.38%	−124	−4.88	−89	+3.14

, the most considerable change pertains to the transformation of barren land to built-up areas. Thus, during the time period of the study, 134 ha of barren land were converted into urban land use. The second significant change relates to 128 ha increase in built-up area, which was made possible by the conversion of barren land and vegetation areas. The considerable decrease in vegetation area in the urban part of the basin was moderated to some extent by the creation of gardens in its mountainous part. Gardens were also responsible for the decline in rangeland area; however, it mostly decreased naturally, which was also the reason for the increase in mountainous areas. The overall spatio-temporal LULC changes in the Darabad basin are shown in Figure 7.

3.2. Annual Rainfall Values and Average Maximum 24-Hour Rainfall

After calculating the rainfall gradient in the study area, the average annual rainfall in the Darabad basin and its sub-basins was obtained, as shown in Figure 8. As observed in Figure 8, the Darabad basin has the highest average height (2617.9 m) and annual rainfall (527.2 mm) among all sub-basins. In contrast, sub-basin E exhibits the lowest values, with an average height of 1704.5 m and annual rainfall of 354.9 mm. Significant differences in annual rainfall are evident between the sub-basins. For example, sub-basin C, with an average height of 2725.6 m, receives an annual rainfall of 463.4 mm, highlighting the direct relationship between elevation and rainfall in the region. Additionally, sub-basins with lower elevations, such as E and F, experience annual rainfall below 400 mm.

To calculate the maximum 24-hour rainfall in different return periods for the sub-basins, linear relationships were established between the maximum 24-hour rainfall of the stations in each return period and their average annual rainfall. Figure 9 illustrates the maximum 24-hour rainfall values of the sub-basins for return periods of 2, 5, 10, 25, 50, and 100 years.

As shown in Figure 9, the Darabad basin exhibits the highest maximum 24-hour rainfall values across all return periods, with 41.36 mm for the 2-year return period, 51.88 mm for the 5-year return period, and 106.27 mm for the 100-year return period. In contrast, sub-basin F demonstrates the lowest rainfall values, with 31.72 mm for the 2-year return period, 41.57 mm for the 5-year return period, and 92.03 mm for the 100-year return period.

Notably, the differences between sub-basins become more pronounced as the return period increases. For example, the rainfall for the 100-year return period in sub-basin C reaches 101.61 mm, while in sub-basin E, it is only 93.84 mm. These variations highlight the spatial heterogeneity of extreme rainfall events across the study area, which can be attributed to differences in elevation, topography, and climatic conditions.

This analysis provides valuable insights into the potential impacts of extreme rainfall events in different parts of the Darabad basin, which is crucial for designing appropriate flood management strategies and hydraulic infrastructure.

3.3. Investigating Flooding Potential in the Basin

After the initial estimation of stormwater volume and peak discharge with the created model, the SCS-CN model was calibrated, and the final model was prepared. The maximum instantaneous discharge statistics of Qalak station were used to calibrate the model. The parameters of soil maintenance coefficient and initial loss were considered as parameters of calibration or optimization of the model with 20% variance. The initial storage amount and the maximum storage capacity depend on various factors such as soil type, vegetation type, permeability, pit storage, and initial soil moisture. After calibrating the model, the peak discharge and flood volume in the sub-basins were calculated.

After preparing the SCS-CN rainfall–runoff model and calibrating it in the study area, the values of the maximum specific instantaneous discharge with different return periods were calculated for each sub-basin. The results obtained were used to classify flooding potential and prepare the corresponding layer, as one of the criteria for site selection of LID-BMPs. The maximum specific instantaneous discharge of each sub-basin is shown in Figure 10. After preparing the SCS-CN rainfall–runoff model and calibrating it for the study area, the values of the maximum specific instantaneous discharge for different return periods were calculated for each sub-basin. These results were utilized to classify flooding potential and prepare the corresponding layer as a criterion for the site selection of LID-BMPs. The maximum specific instantaneous discharge values for each sub-basin are shown in Figure 10.

As shown in Figure 10, sub-basin E exhibits the highest maximum specific discharge values across all return periods, with 0.23 m³/s/km² for the 2-year return period, 0.6 m³/s/km² for the 5-year return period, and 2 m³/s/km² for the 100-year return period. On the other hand, sub-basins F and B show the lowest values, with sub-basin F recording only 0.11 m³/s/km² for the 2-year return period and 1.06 m³/s/km² for the 100-year return period. The differences in maximum specific discharge between the sub-basins are influenced by factors such as land use, soil type, and topographic conditions. For instance, sub-basin A, with a maximum specific discharge of 1.8 m³/s/km² for the 100-year return period, reflects moderate flooding potential, while sub-basin G, with a value of 1.33 m³/s/km² for the same return period, shows relatively lower flooding potential. These findings indicate significant spatial variations in flooding potential within the study area. The classification of flooding potential based on these results plays a critical role in the strategic placement of LID-BMPs to mitigate flood risks effectively.

3.4. The Relationship Between Land Use/Land Cover Change and the Increase in Peak Discharge and Flood Volume

For investigating the relationship between land use change and the increase in peak discharge and flood volume in the study area, the years of 2000 and 2022 were selected as index years. The curve number values were calculated for each of the years, and according to the rainfall values at the time of concentration and the SCS-CN calibrated model, the peak discharge and flood volume values were estimated for the years. As it can be seen in Figure 11 and Figure 12, the amounts of changes in peak discharge and flood volume are the highest in the 2-year return period, and as the return period increases, the amounts of changes are reduced. The reason is that, as the return period becomes longer, the peak discharge and, naturally, the destructive power of the flood increase, for instance, in a rainfall with a return period of 100 years, which is extremely powerful, and the land use parameter is less effective on the peak discharge and flood volume. Only urban sub-basins that were affected by land use changes are shown in this graph. According to the graph, it can be seen that the lowest amount of change in peak discharge and flood volume per land use change belongs to sub-basin B, because there was a slight change in land use. The highest amounts of changes in stormwater volume and peak discharge per land use change occurred in sub-basin F because it has been affected more by transformations caused by urbanization than other sub-basins.

3.5. Site Selection of Optimal Low-Impact Developments

The fuzzy maps of site selection criteria are shown in Figure 13. Map values are between 0 and 1, where 0 represents the lowest value and 1 represents the highest value. The criteria map shows that the high values are different in the study area. For example, in the slope map, the areas with high values are mostly in the central and southern regions of the studied region. In the standard map of rainfall, the northern regions have high values compared to other regions. By combining these maps, it becomes evident that different criteria emphasize distinct areas of the study region, underlining the importance of multi-criteria decision-making approaches in site selection processes. This integrated analysis allows for the identification of optimal locations that balance the varying influences of slope, rainfall, and other relevant factors.

After determining the appropriate areas for each LID-BMP, their location raster map was divided into five classes to ascertain the best class that represents the best areas for each LID-BMP. After extracting the best areas for each LID-BMP, these areas were put together. In the areas that were found to be suitable for two or three LID-BMP types, according to the characteristics of the LID-BMPs and the priority of the land use for each LID-BMP, these area was assigned to one LID-BMP model. In Figure 14, the final map of the suitable locations for all four LID-BMPs of bioretention basin, green roof, grass swale, and porous pavement is shown together in the urban part of the Darabad basin. Finally, 119.85 ha for grass swale, 67.77 ha for bioretention basin, 37.31 ha for green roof, and 18.84 ha for porous pavement are proposed in the urban part of Darabad basin. The ratio of the area of LID-BMPs to the area of each sub-basin can be seen in

Table 10. The ratio of the area of LID-BMPs to the total area of each sub-basin.

Sub-Basin	Area of LID-BMP (ha)	Total Area (ha)	Ratio
Dint	27.60	323.98	8.51%
E	49.72	398.53	12.47%
F	63.23	377.67	16.74%
G	56.43	140.03	40.29%
H	15.04	118.82	12.65

. In sub-basins A, B, and C, which are non-urban, there is no effect on the pick discharge and stormwater volume, so in the continuation of the research, the numbers related to these three basins have been avoided.

After determining appropriate sites for each LID-BMP, it is time to model the SCS-CN model assuming the use of LID-BMPs in the sub-basins. The effectiveness of each LID-BMP model is calculated according to Table 2. As mentioned earlier, after applying LID-BMPs in each sub-basin, the CN values change for them. In

Table 11. Values of the CN in the sub-basins after applying LID-BMPs.

Sub-Basin	A	B	C	Dint	E	F	G	H	Darabad Basin
CN values	94.54	94.93	94.46	95.81	95.66	87.62	90.33	93.11	93.72

, the CN values of the sub-basins after the use of LID-BMPs are discussed.

3.6. Simulation Results

The results show that the greatest impact of LID-BMPs on the peak discharge in all return periods occurred in sub-basin G because the ratio of the area where LID-BMPs were applied to the total area of the sub-basin was higher in sub-basin G. The simulation results of the rainfall–runoff model showed that, if LID-BMPs are applied, the average percentage of peak discharge reduction in urban sub-basins in a return period of 25 years will be 28.82%. Figure 15 show the percentage of change in peak discharge with the application of LID-BMPs. It is also evident from the graphs that, under the condition of using LID-BMPs, the percentage of peak discharge changes in shorter return periods is higher and better than that in longer return periods. In other words, LID-BMPs perform better in shorter return periods.

As shown in Figure 16, if LID-BMPs are used, the average percentage of flood volume reduction in urban basins in the 25-year return period will be 32.72%. Again, a better performance of LID-BMPs in shorter return periods can be seen here.

4. Discussion

The methodological approach has provided sufficient evidence that there is a relationship between temporal LULC change and changes in important hydrological parameters of stormwater volume and peak discharge. According to a study in Indiana, a 46% increase in urban land use can cause a 75% growth in stormwater volume [60]. Another survey conducted in Florida from 1920 to 1990 showed a 22–55% rise in stormwater volume due to land use changes from green cover to built-up [61]. A report from China illustrated that a 10% increase in built-up area led to an 11% increase in surface stormwater due to urbanization effects [62]. It is concluded that there is a connection between flow rate, as a measure of flood hazard intensity, the type of LULC, the intensity of rainfall, and return periods.

As a natural phenomenon, flooding cannot totally be avoided and leaves a significant impact on local livelihoods and socioeconomic development. In this research, after determining the level of vulnerability of each of the sub-basins against floods as a result of LULC change, their vulnerability in higher return periods was also estimated. The use of LID-BMPs was suggested as a management solution to control stormwater volume and peak discharge, and the selected models were located according to their functional characteristics and the spatial characteristics of the sub-basins. In addition, the impact of LID-BMPs in reducing flood potential was calculated by reducing the two parameters of stormwater volume and peak discharge for different return periods.

4.1. Urban LULC Change

Our results show a significant increase in urbanization in Tehran, particularly in the peri-urban areas. This trend has led to substantial changes in LULC, similar to the findings of other global studies. For instance, a study by Olivera and DeFee [63] observed a 181% increase in impervious areas in the Whiteoak Bayou watershed in Texas, USA, which led to a 77% increase in annual runoff. Similarly, Huq and Abdul-Aziz [64] found that urban expansion led to a marked increase in flood volume, supporting the conclusion that urbanization significantly impacts stormwater management.

In comparison to previous studies, such as the one by Moniruzzaman et al. [65] which reported a 32.92% increase in stormwater volume due to urbanization, our results for Tehran’s Darabad basin show a notable shift from vegetative and barren lands to built-up areas, which has directly influenced both stormwater volume and peak discharge. This highlights the accelerated urbanization pattern in Tehran, which is more rapid than observed in many other cities. In particular, sub-basin F, being the most urbanized area in our study, exhibited the highest changes in flood volume, aligning with similar patterns reported globally where urban areas have been most affected by changes in LULC.

4.2. Effect of LULC Change on Urban Stormwater

The regional hydrological cycle and flood control have been significantly impacted by rapid urban development [66]. According to a prediction, if the expansion of built-up areas is not controlled, the risk of flooding is likely to rise in the future because of the increase in stormwater [67]. Tehran’s geomorphological features include mountains in the north, a Piedmont zone in the center, and a desert in the south. These features result in several rivers and watercourses flowing through the city, which could potentially lead to flooding. Due to the lack of proper drainage systems, flash floods caused by overflowing rivers or extreme rainfall occur frequently in Tehran. The analysis of satellite images of the study area indicates the reduction in vegetation cover and bare land and their conversion to built-up area. The results of rainfall–runoff modeling showed that the lowest and highest changes in peak discharge and flood volume occurred in sub-basin B and F, respectively. This is because almost half of sub-basin B is located in the mountainous part of the area, while sub-basin F is located in the most urban part of the study area and, consequently, is more affected by developmental land use changes.

4.3. Site Selection of Low-Impact Developments

Determining the location of LID-BMPs is inherently a spatial decision and depends on various criteria such as topography, access, and environment. Therefore, in this study, a GIS-MCDA model was used to locate different LID-BMPs. The integration of MCDA and GIS creates synergy by combining spatial data analysis with user preferences, enabling structured and creative decision making [68,69,70]. The accuracy and efficiency of GIS-MCDA models depend on the considered effective criteria, the weight determined for the criteria, and the efficiency of aggregation models [71,72]. By reviewing past studies and considering experts’ opinion, effective site selection criteria and the weight of each criterion were calculated using the AHP method. Among the advantages of this method is the simultaneous consideration of quantitative and qualitative criteria, the possibility of checking the consistency of the results, and easy implementation.

4.4. Limitation and Uncertainties

One of the limitations of this study is the lack of access to free satellite images with high spatial resolution, which can cause the impact of urban infill on flooding to be underestimated or overestimated [73]. Another argument inferred that, in order to accurately represent the unique sub-catchment responses to storm events, detailed topological and geographic information is required rather than simply classifying the impervious land cover [74].

The lack of knowledge about the spatial vulnerability of floods can lead to a poor understanding of its management. Also, the lack of meteorological and rain gauge stations in the study area were other limitations of this study in preparing relevant data, and therefore, rainfall–runoff modeling was performed using the statistics of the nearest stations, which increased the uncertainty.

5. Conclusions

The results of this study assessed the impact of land use changes and the application of low-impact development best management practices (LID-BMPs) on stormwater management in an urban basin located in Northern Tehran. The findings showed that, between 2000 and 2022, green and barren land uses decreased from 357 and 262 hectares to 327 and 128 hectares, respectively, while urban land use increased by 13.02%. These changes highlight the urgent need for effective management strategies to address the challenges of urbanization. The study evaluated four types of LID-BMPs, including bioretention basins, green roofs, grass swales, and porous pavements. Site selection analysis revealed that land use is the most important criterion for choosing bioretention basins and green roofs, while slope and distance from streets are the most decisive factors for grass swales and porous pavements. The results indicated that the study area has the highest potential for grass swales, with approximately 120 hectares of urban land identified as suitable for their implementation. Changes in peak discharge and flood volume before and after the application of LID-BMPs demonstrated that these practices significantly reduce peak discharge and flood volume, particularly for shorter return periods, highlighting their effectiveness in managing frequent, low-intensity storm events. Based on these findings, it is recommended that urban planners prioritize the implementation of grass swales in the identified areas and integrate different types of LID-BMPs into urban designs based on the specific characteristics of each area. Additionally, due to the positive impact of these practices on reducing short-term floods, they should be used to manage frequent storm events. Finally, future research should focus on integrating high-resolution satellite imagery and detailed topographic data for the more accurate modeling of LULC changes and hydrological impacts. Expanding meteorological monitoring networks, examining the combined effects of urbanization and climate change, and employing advanced AI-based modeling techniques could significantly enhance flood risk assessments. Additionally, incorporating community participation and policy integration will help develop more effective and sustainable urban flood mitigation solutions.