Identifying Cost-Effective Low-Impact Development (LID) under Climate Change: A Multi-Objective Optimization Approach

Abstract

Low-impact development (LID) is increasingly used to reduce stormwater’s quality and quantity impacts associated with climate change and increased urbanization. However, due to the significant variations in their efficiencies and site-specific requirements, an optimal combination of different LIDs is required to benefit from their full potential. In this article, the multi-objective genetic algorithm (MOGA) was coupled with the stormwater management model (SWMM) to identify both hydrological and cost-effective LIDs combinations within a large urban watershed. MOGA iteratively optimizes the types, sizes, and locations of different LIDs using a combined cost- and runoff-related objective function under both past and future stormwater conditions. The infiltration trench (IT), rain barrel (RB), rain gardens (RG), bioretention (BR), and permeable pavement were used as potential LIDs since they are common in our study area—the city of Renton, WA, USA. The city is currently adapting different LIDs to mitigate the recent increase in stormwater system failures and flooding. The results from our study showed that the optimum combination of LIDs in the city could reduce the peak flow and total runoff volume by up to 62.25% and 80% for past storms and by13% and 29% for future storms, respectively. The findings and methodologies presented in this study are expected to contribute to the ongoing efforts to improve the performance of large-scale implementations of LIDs.

1. Introduction

The expansion of impervious surfaces in urban areas, combined with the increased frequency and intensity of severe storms caused by climate change, often leads to an increase in stormwater runoff, deterioration of water quality, and a decrease in the arrival time of peak flows [1,2,3,4,5,6,7]. Several studies showed the growing vulnerability of urban drainage systems to flooding [8,9,10,11]. The damage to property and ecosystems from such flooding often exceeds the cost of improving the stormwater system [10]. Thus, upgrading the urban stormwater system is considered an effective way of mitigating major floods in urban areas [12].

Several emerging runoff management technologies, including the best management practices (BMPs), low-impact development (LID), and sustainable urban drainage system (SUDS), were recently introduced to preserve the predevelopment runoff [13] and improve stormwater systems [14]. Particularly, the LID has gained considerable recognition as a promising stormwater management practice that can effectively reduce runoff and improve water quality by handling the stormwater close to its sources [15,16,17]. It uses distributed, small-scale, onsite, and integrated land-management practices to reestablish the natural hydrological conditions that mimic the predevelopment conditions [18]. Several studies showed the effectiveness of LIDs in mitigating urban flooding [6,19,20,21] and their ability to provide multiple benefits to ecosystems in various climatic and geographic regions [21,22]. They are more effective in adapting to climate change than the traditional flood control measures [23] and can be utilized to replace conventional stormwater management practices for cities’ long-term sustainability [24,25,26,27]. Consequently, LIDs have gained acceptance worldwide to improve infiltration within a watershed, mitigate floods, control erosion, and maintain surface water quality [28]. In the United States, cities are required to include green infrastructures in their combined sewer overflow control plans to minimize stormwater pollutions. Many cities, including Seattle, Portland, Kansas, Washington, New York, and Philadelphia, have adopted this strategy to some level [29]. Seattle and Portland lead the efforts because of their relatively strict stormwater regulations [30]. In Melbourne, Australia, the stormwater management systems are mostly retrofitted with LIDs that include widespread uses of green roofs to meet runoff and water quality control targets.

The infiltration and evapotranspiration processes play essential roles in the performance and selection of sites for LIDs. Thus, factors affecting these processes, including soil type, permeability, roughness coefficient, land cover/use, precipitation, and slopes, are used to determine the appropriate LIDs. For example, infiltration trenches (IT) strongly depend on the underlying soil types and permeability, and are generally built in the areas where they are intended to intercept surface flow from roadside and parking lots [18]. The overall effectiveness and ability to manage the different parts of the hydrograph also depends on the types of LIDs, watersheds, and storm conditions. Reference [31] showed that a green roof outperforms permeable pavement in peak-flow reduction, while [32] demonstrated that bioretention is better than the green roof and permeable pavement in reducing the surface runoff from a small drainage area. Eckart et al. [33] showed that bioswale cells in combination with bioretention (BR) could provide higher protection against floods compared to the combined rain barrel (RB) and permeable pavement. The reduction in peak flow by green roofs can reach 22–70% depending on the intensity and duration of rainfall [34]. Oberascher et al. [35] evaluated the performance of a rain barrel (RB) in Austria by optimizing its installation site and found an 18–40% reduction in flood volume depending on where they are installed. Overall, green roofs, permeable pavements, and bioretention cells are found to be the most effective approaches to mitigating urban runoff [16,36,37]. Additional LID practices, such as rain gardens (RG), detention ponds, and infiltration systems, are also used to mitigate the issue of water quality involving non-point pollutants [18,38,39].

In addition to the hydrological performance, the cost of implementation is an important factor in deciding the type, size, and location of LIDs [40]. Thus, hydrological modeling and optimizations, considering the different factors, are often required to implement LIDs to achieve maximum reduction in peak flow and total runoff at the lowest possible cost [32,41]. Several models have been used to simulate and evaluate the effects of LIDs on hydrology [6,42,43,44,45]. In the United States, the stormwater management model (SWMM) is commonly used because of its reasonable amount of data and computational needs [46,47,48]. The model is particularly suitable for evaluating the performances of various stormwater systems under current and future storm conditions in urban areas [33,49,50].

For optimizing the LID types and sizes, multi-objective genetic algorithms (MOGA) are commonly used because of their advantages of representing different LIDs and handling multi-objective functions [51,52,53,54]. For example, Wu et al. [55] coupled the SWMM model with a nondominated sorting genetic algorithm (NSGA) to find the optimal combinations of the infiltration trench (IT), permeable pavement, and bioretention (BR) at the watershed scale that maximize the reduction in the polychlorinated biphenyl (PCB) load at minimum implementation cost. A similar approach was employed by [7,32,56] to improve stormwater runoff quantity and quality while minimizing the implementation costs.

In this study, we coupled MOGA with SWMM to optimize the implementation of different LIDs under current and future storm conditions. Unlike earlier studies, we incorporated different LIDs into existing stormwater systems to handle floods from both past and future storm events and considered several decision variables, including the LID types, numbers, sizes, and locations, to determine the optimal combinations and performance–cost tradeoff of LIDs and traditional stormwater systems. The paper is also among the few to demonstrate the applicability of the optimization technique to a relatively large urban watershed. Thus, the study not only extended the application of the optimization techniques to include more decision variables and the potential change in future rainfall but also contributed to engineering practices for addressing the stormwater challenges in larger cities.

The methodology was tested for the city of Renton, WA, USA, under both past and future climate change scenarios. The region experienced increasing trends in floods, recording 16 major flooding events since 1990. The extreme precipitations and floods in the region are expected to continue, increasing in the future because of climate change [57]. These flood trends and the frequent failures in the stormwater system show the challenges the city faced in adequately managing stormwater runoff in the current climatic conditions that will likely become severe in the future. Developing and implementing optimum stormwater management strategies by incorporating LIDs are critical to reducing the flood risk in the area. We first assess the performance of single LIDs and combinations of different LIDs in reducing the increased stormwater using the SWMM model. The LIDs are then further optimized for their types, locations, and sizes using a coupled SWMM and MOGA to identify hydrologically and cost-effective approaches. This paper aimed to assess the suitability and benefits of integrating LIDs into the existing stormwater management system to manage the stormwater challenges arising from increased urbanization and extreme precipitation because of climate change. Specifically, it aimed to identify an appropriate climate model and downscaling method for the city of Renton, adapt a suitable hydrologic model and its integration with various LIDs, assess the impacts of LIDs on runoff reduction during flood events under a changing climate, optimize LID designs to obtain maximum stormwater reduction benefits, and conduct a cost–benefit analysis of the LID-based stormwater management system.

2. Materials and Methods

2.1. Study Area

The study area is located in the city of Renton, WA, USA, which is a suburb of Seattle and is part of the lower Cedar River and Puget River Basin (Figure 1). It occupies nearly 53.82 km² of areas comprising 606 sub-catchments of undeveloped land (52%) and paved developed areas (48%). The eastern part of the study area mainly contains commercial buildings, while the western part is primarily commercial, multifamily, high-density residential, medium-density residential, low-density residential, grass, forest, and wetland areas. The boundary to the northwest is marked by a ridgetop connecting the city to Lake Washington. The Cedar River flows through the city and into Lake Washington, collecting most of the stormwater. The regional climate is classified as mild Pacific maritime and is greatly influenced by the Pacific Ocean. It receives approximately 46 inches (1168.4 mm) of annual precipitation, of which nearly two-thirds occurs during the winter season from October to March. During the winter months, the region is also affected by storm systems associated with mid-latitude cooled jet streams and low-level jet streams ahead of the cold front (also known as the atmospheric river), resulting in snow and severe storms.

Based on the Federal Emergency Management Agency (FEMA), the floodplain zones in the city are shown in Figure 1. The flood zone designation in the area is primarily Zone X (moderate flood hazard areas) and Zone AE (special flood hazard area), which were established based on 500-year and 100-year floodplain elevations, respectively [58]. The flood risks in the region are rising, with the Puget Sound region experiencing at least 16 federally declared flood-related disasters with significant flooding and closing of highways since 1990 [59,60] showed considerable change in the atmospheric rivers and increasingly severe storms in the region over the past 30 years (1980 to 2009). The recent major flood events include the November 2022 flood caused by a record-breaking 203.2 mm (8 in) of rainfall over 24 h and the February 2020 phase four (severe) flood caused by the fifth-highest daily rainfall (82 mm) and snowmelt on the mountain [59]. Many areas in the city are experiencing challenges in adequately and effectively managing the increasing stormwater and floods. The cost of the damages from these floods were in the tens of millions of dollars, with the 2020 flood alone costing the city over USD 27 million in damages and associated costs [61]. The region is expected to experience severe floods more frequently in the future because of climate change [62]. The sea-level rise in Puget Sound is also expected to increase the flooding risk in the area. According to [49], the region may witness a maximum sea-level rise of 0.4 feet (12.19 cm) by 2030 and 2.2 feet (67 cm) by 2100 under the SSP-585 high-emissions scenario.

The existing drainage systems mainly comprise underground pipes, catch basins, curbs, and gutters, which collect runoff throughout the stormshed and direct it to a detention wetland for retention and natural treatment. Failure in and overflow from the stormwater systems often cause the flooding of streets and private property. The city spends about USD 3.5 million each year to improve stormwater infrastructure and reduce flooding [63]. In 2007, the city of Renton adopted a FEMA-based comprehensive stormwater system improvement project to address current and future basement and street flooding by increasing the size of the culverts within the watershed. However, the upgraded culverts failed during the recent severe flooding. Moreover, there will be an approximately 26% increase in the total discharge in the future climate change scenario compared to the historical period in the area, which will further aggravate the flooding problems in the future [64].

2.2. Datasets

The base input data for the model and optimization include precipitation, watershed, and stormwater drainage data. The DEM and land-cover data with 0.5 m and 10 m resolutions, respectively, were obtained from the city of Renton. The topography in the study area is relatively flat, with an elevation range from 6 to 9 m. The catchment boundaries and sub-catchments were delineated using DEM. The land covers are predominantly impervious surfaces (48% with commercial buildings ranging from 2030 m² to 12,140 m²) and woodland (12%). The rest are grasses, small deciduous trees, and bushes. The high-resolution Soil Survey Geographical Database (SSURGO 2.2), obtained from the National Resources Conservation Service (NRCS) Geospatial Database, was used. The dominant soil type in the study area is Woodinville silty loam (Wo), with traces of Newberg silty loam (Ng) and Puget silty clay loam (Pu). All three types of soils have slow runoff potential and have high water holding capacity. The sub-catchment parameters, including area, slope, and infiltration characteristics, depth of depression storage for impervious and pervious surfaces, percent zero (percentage of impervious area that has no depression storage), Manning’s n values for impervious and pervious surfaces, and internal routing variables were used to determine the curve numbers for rainfall–runoff modeling. The 15-min precipitation data from 1 January 1995, through 1 May 2014, was obtained from the city of Renton database (https://green2.kingcounty.gov/hydrology/GaugeMap.aspx). The daily precipitation projection for 2020–2050 under the shared socioeconomic pathway (SSP585) emission scenario was downscaled from two CMIP6 climate models, namely MIROC6 and CMCC-ESM2. The projected daily precipitation data were biased-corrected and converted to 15 min data using the SCS Type IA distribution. The 15 min historical and future precipitation data were then used to develop the intensity–duration–frequency (IDF) curves of the storms in our study area.

2.3. Optimization Framework

A simulation–optimization modeling framework (Figure 2) was developed to reduce total runoff and peak flow across the watershed at a minimum implementation cost of LIDs. The framework couples the SWMM model with a genetic algorithm (GA) to find Pareto optimal solutions for different LIDs and design parameters (i.e., type, number, location, and size). The Borg multi-objective evolutionary algorithm [32] was used because of its effectiveness in handling highly nonlinear problems and many decision variables.

The optimization framework was developed using MathWorks produces mathematical computing software MATLAB R2022b which is established in Las Vegas, USA to coupling the SWMM model and MOGA. The framework starts with a randomly generated initial population of LID type, size, location, and number, which are represented by binary strings of numbers or chromosomes in the MOGA. Each individual in the population is typically referred to as a chromosome and represents the size of a specific solution for the investigated problem. There are no lower or upper limits for population size. The SWMM then evaluates the performance of the individuals in the population based on the resultant runoff and cost. The simulation results from SWMM were sent back to the MOGA for constructing the Pareto-optimal front and generating the next generation. This process continues till optimized solutions are found, or user-specified limits are reached.

2.4. SWMM Model Development

The recent version of SWMM (5.1), which includes a toolbox to incorporate several LIDs, was used to simulate the stormwater runoff and evaluate the LIDs’ performances in reducing runoff at the basin outlet. The model effectively simulates the rainfall-runoff processes and the hydrodynamics associated with the stormwater system for both event-based and continuous rainfall events [48]. SWMM was also applied previously to optimize the LID implementations [45,65].

The study area was divided into 606 sub-catchments based on the topography (Figure 3). The sub-catchment area and average width were determined using the auto-length feature of SWMM as the backdrop technique and then scaling the dimensions. The sub-catchment width is important to calculate the overland flow distances that the flow travels before becoming channelized. The accuracy of the distances was confirmed by comparing with the topographical details of the stormshed. The slope and impervious areas were configured using GIS and the city of Renton’s geodatabase (CORMAP). The curve number method coupled with dynamic wave surface runoff routing and Manning channel routing was used to simulate the watershed runoff because of their computational efficiencies and previous applications for stormwater system optimization [7,66]. The dynamic wave routing method has the advantage of considering the backwater effects, channel storage, entrance/exit losses, pressurized flow, and flow reversal.

The existing stormwater system, including the LIDs, as well as proposed new LIDs, were incorporated into the model. The infiltration at the LIDs was computed using the Green-Ampt method that uses the initial moisture deficit, ${\theta }_0$, parameter to represent the initial soil moisture condition before the storm. The initial values of ${\mathit{\theta }}_{\mathbf{0}}$ were determined based on the soil types and the depth of the groundwater water table. The King County database (https://kingcounty.gov/services/gis/Maps/imap.aspx, during 20 June 2019 and 7 October 2019 was used to determine the groundwater level, which varies in depth from 0 to 16 feet. Given that the study area has predominately sandy loam soil type, we have considered initial moisture deficit values ranging from 0.35 for sand to 0.25 for sandy loam [66].

2.4.1. Sensitivity Analysis and Model Calibration

The SWMM model outputs can be sensitive to several input data and parameters, making sensitivity analysis an essential step in the model development [67]. Earlier studies found that time to peak and peak flows are significantly affected by the width of the sub-catchment and the soil moisture content [68]. In addition, the percent of imperviousness areas and slope of the catchment significantly influences flow hydrographs, particularly the peak flows. In this study, we considered different hydrological and catchment-related parameters such as percent zero imperv (or percent of the impervious area with no depression storage), percent of imperviousness, curve number, manning coefficient, and channel width for sensitivity analysis. The study area has negligible slopes; hence, the slope is not considered in the sensitivity analysis.

Table 1. SWMM parameters used for the sensitivity analysis.

Parameters	Values	Initial Values Source
Sub-catchment width	Variable	Geometry
Manning’s n for impervious surface	0.05–0.15	SWMM manual
Impervious percentage with no depression (%)	5–15	SWMM manual
Curve number	Variable	DEM, land use, and soil data

provides lists of parameters considered in the sensitivity test.

The SWMM model was manually calibrated using the sensitive sub-catchments and routing parameters. The initial values of the curve numbers and other variables (e.g., Manning’s n, storage depth) were identified from the sub-catchment characteristics and the SWMM manual [66]. The storm events from 18–20 June 2019 and 7 October 2019 were used to calibrate the SWMM model, while the storm events from 8 August 2019 and 8–9 October 2019 were used to validate the model. After the validation, the model was used to simulate the 50-year and 100-year past and future climate change scenario design storms with a 15 min time step. The high-emission scenario (SSP585) was used for a conservative design of LIDs and to assess their performance under a potential large increase in future storm intensity considered to generate the future storms.

2.4.2. Existing LID Types in Our Study

Five LIDs, namely the infiltration trench (IT), rain barrel (RB), rain garden (RG), bioretention (BR), and permeable pavement (PP), were used because of their current popularity in the study area. They were identified based on a literature review, locally available data, and consultation with an engineering company that has worked on municipal infrastructure projects in the area. The LIDs are often hydrologically designed as a series of vertical layers, and their properties are defined using per unit area. SWMM 5.1 performs water balance computations during the simulation to keep track of how much water is stored or moved within each LID layer. The LIDs parameters, including available area, land use, soil infiltration properties, water-table depth, and slope, were determined based on watershed characteristics and LID design standards.

Table 2. The five existing LIDs considered in our study, their associated costs, sizes, and implementation strategies.

LIDs	Cost in WA	No. LIDs	Implementation Strategies
Bioretention Cell (BR)	4.8 USD/m2	78	Bioretention (BR) captures runoff from sidewalk and 1/4 of driveway and parking lots in commercial development.
Rain Barrel (RB)	150 USD/unit	130	Rain barrels (RB) collect water from the entire roof of a house.
Permeable Pavement (PP)	30.16 USD/m2	43	Permeable pavement captures runoff from the entire driveway and parking lot in commercial development.
Infiltration Trench (IT)	3.36 USD/m2	90	The infiltration trench (IT) captures runoff from half of the walkway along with the front houses and half of the driveway.
Rain Garden (RG)	3 USD/m2	105	Rain garden (RG) captures water from half of the house roof and half of the driveway.

shows the numbers of existing IT, RB, RG, BR, and PP in the study area and the estimated construction costs obtained from the city of Renton and online vendors.

2.4.3. Future Storm Considering Climate Change

The ability of LIDs to manage the expected increase in future climate-change scenario stormwater runoffs was assessed using the projected precipitations from the CMIP6 under SSP 585. The precipitation data from the CMIP6 have a coarse spatial resolution (ranging from 10 to 300 km) and are not suitable for direct use in the SWMM simulations. Therefore, the projected precipitation is commonly downscaled and bias corrected using statistical relationships with the observed precipitation [69]. Many statistical techniques are available to downscale and bias correct the projected precipitation, with the delta change (DC), quantile mapping (QM), and equidistance quantile mapping (EQM) methods being the most commonly used [70].

The DC method (Equation (1)) downscales and bias corrects the projected precipitation (2020–2050) by multiplying it with the ratio of observed and computed mean precipitations from the training period (1994–2014) [71].

$${\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}^{\mathit{D}\mathit{e}\mathit{l}\mathit{t}\mathit{a}}={\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{S}\mathit{S}\mathit{P}}\times \frac{{\overline{\mathit{P}}}_{\mathit{O}\mathit{b}\mathit{s}}}{{\overline{\mathit{P}}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{h}\mathit{i}\mathit{s}\mathit{t}.}}$$

(1)

where ${\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}^{\mathit{D}\mathit{e}\mathit{l}\mathit{t}\mathit{a}}$ is the bias-corrected precipitation projection, ${\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{S}\mathit{S}\mathit{P}}$ is the precipitation projection from the climate models, and ${\overline{\mathit{P}}}_{\mathit{O}\mathit{b}\mathit{s}}$ and ${\overline{\mathit{P}}}_{\mathit{O}\mathit{b}\mathit{s}}$ are the means of observed and estimated precipitation from the climate models for the training period, respectively.

The DC method is easy to implement but is often limited in reproducing extreme precipitation because of its emphasis on the mean precipitation. In contrast, QM is capable of reproducing the full spectrum of the precipitation distribution as it corrects precipitations of different quantiles. However, the QM method is also limited because of its assumptions of constant variance and skewness of the climatic variables over the projection and training periods. This led to an extension of the QM method to EQM. The adjustment function (i.e., differences between observed and model values for any percentile during the training period) remains the same for the training and future periods in the QM method, whereas the EQM method also considers the differences between the cumulative distribution functions (CDFs) for the training and future periods. Mathematically, QM and EQM methods can be written as shown in Equations (2) and (3), respectively.

$${\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}^{\mathit{Q}\mathit{M}}={\mathit{F}}_{\mathit{O}\mathit{b}\mathit{s}-\mathit{h}\mathit{i}\mathit{s}\mathit{t}}^{-\mathbf{1}}\left[{\mathit{F}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{h}\mathit{i}\mathit{s}\mathit{t}}\left({\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}\right)\right]$$

(2)

$$\begin{gathered} {\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}^{\mathit{E}\mathit{Q}\mathit{M}}={\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}+{\mathit{F}}_{\mathit{O}\mathit{b}\mathit{s}-\mathit{h}\mathit{i}\mathit{s}\mathit{t}}^{-\mathbf{1}}\left[{\mathit{F}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}\left({\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}\right)\right]- \\ {\mathit{F}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{h}\mathit{i}\mathit{s}\mathit{t}}^{-\mathbf{1}}\left[{\mathit{F}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}\left({\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}\right)\right] \end{gathered}$$

(3)

where ${\mathit{F}}_{\mathit{O}\mathit{b}\mathit{s}-\mathit{h}\mathit{i}\mathit{s}\mathit{t}}^{-\mathbf{1}}$ and ${\mathit{F}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{h}\mathit{i}\mathit{s}\mathit{t}}$ are the CDFs for the observed and predicted precipitation for the training period (1994–2014), respectively, ${\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}^{\mathit{Q}\mathit{M}}$ and ${\mathit{P}}_{\mathit{G}\mathit{C}\mathit{M}-\mathit{f}\mathit{u}\mathit{t}}^{\mathit{E}\mathit{Q}\mathit{M}}$ the biased corrected precipitation projections for 2020–2050.

We tested and compared the DC, QM, and EQM methods for downscaling the maximum daily precipitation. Two climate models (MIROC6 and CMCC-ESM2) from the CMIP6 were used because of their relatively good performances for the region. The Nash–Sutcliffe (NSE) and correlation coefficient (R2) are used to quantify the degree of agreement between observed and downscaled precipitation, with the value ranging between 0 and 1.

The downscaling approaches were further evaluated based on their performances to reproduce not only the mean precipitations but also the observed extreme precipitations. Frequency analysis of annual maximum daily precipitation was performed for the observed and downscaled precipitations. The results were compared to identify the best downscaling method, and the corresponding downscaled precipitation was used for developing the future IDF curves and design storms used for the SWMM simulations. For a particular return period, a typical IDF relationship can be described by a general mathematical statement provided in Equations (4)–(11).

$$\mathit{i}=\frac{\mathsf{\omega }}{{\left({\mathit{d}}^{\mathit{v}}+\mathsf{\theta }\right)}^{\mathsf{\eta }}}$$

(4)

where $\mathit{\omega }, \mathit{v}, \mathit{\theta }, \mathrm{and}$η are non-negative coefficients with $\mathit{v}\mathit{\eta }$ ≤ 1. The equation was derived from the Sherman (1931) IDF equation:

$$\mathsf{i}=\frac{\mathit{a}\left(\mathit{T}\right)}{\mathit{b}\left(\mathit{d}\right)}$$

(5)

where a is a function of the return period (T) and represents the quartile for the cumulative probability distribution associated with T. The denominator function b(d) can be expressed as:

$$\mathit{b}\left(\mathit{d}\right)={\left(\mathit{d}+\mathsf{\theta }\right)}^{\mathsf{\eta }}$$

(6)

where d is is the time duration (usually minutes), and θ and η are model parameters (θ > 0, 0 < η < 1).

a (T), which is ω in Equations (3)–(10), can be calculated using the following two empirical relations [72]:

$$\mathit{a}\left(\mathit{T}\right)=\mathsf{\lambda }{\mathsf{T}}^{\mathit{K}}$$

(7)

$$\mathit{a}\left(\mathit{T}\right)=\mathit{c}+\mathsf{\lambda }\mathbf{ln}\mathit{T}$$

(8)

where parameters λ and K depend on the return period T. a(T) is determined from the probability distribution function of the maximum rainfall intensity I(d) of a given duration d. The intensity of rainfall I(d) has a cumulative distribution of ${\mathit{F}}_{\mathit{I}\left(\mathit{d}\right)}$ that is the same as the cumulative distribution of variable $\mathit{Y}=\mathit{I}\left(\mathit{d}\right)\mathit{b}\left(\mathit{d}\right)$, which represents the rainfall intensity rescaled by $\mathit{b}\left(\mathit{d}\right)$. Mathematically, it can be expressed by:

$$\mathit{a}\left(\mathit{T}\right)=\mathit{P}\left\{\mathit{I}\left(\mathit{d}\right)\le \mathit{i}\right\}=\mathit{P}\left\{\mathit{I}\left(\mathit{d}\right)\mathit{b}\left(\mathit{d}\right)\le \mathit{i}\mathsf{x}\mathit{b}\left(\mathit{d}\right)\right\}=\mathit{P}\left\{\mathit{Y}\le {\mathit{y}}_{\mathit{T}}\right\}$$

(9)

where P is the cumulative probability calculated using:

$${\mathit{F}}_{\mathit{I}\left(\mathit{d}\right)}\left(\mathit{i},\mathit{d}\right)={\mathit{F}}_{\mathit{Y}\left({\mathit{y}}_{\mathit{T}}\right)}=\mathbf{1}-\frac{\mathbf{1}}{\mathit{T}}$$

(10)

Hence, ${\mathit{y}}_{\mathit{T}}$ or $\mathit{a}\left(\mathit{T}\right)$ can be computed using:

$${\mathit{y}}_{\mathit{T}}=\mathit{a}\left(\mathit{T}\right)={\mathit{F}}_{\mathit{Y}}^{-\mathbf{1}}\left(\mathbf{1}-\frac{\mathbf{1}}{\mathit{T}}\right)$$

(11)

Different cumulative distribution models are used to determine the $\mathit{a}\left(\mathit{T}\right)$.

Both the past and future data were used to construct the IDF curves. The future scenario for the annual maximum rainfalls for the years 2019–2049 with 1 h, 3 h, 6 h, 12 h, and 24 h durations were developed using the cumulative distribution functions (CDFs).

2.5. MOGA Optimization

The MOGA algorithm was used to optimize the implementation of new LIDs under the current and future storm conditions. The MOGA parameters are given in

Table 3. MOGA parameters used for the LIDs optimization.

Parameter	Value
Initial population size	100
Maximum population size	200
Probability for crossover	1
Probability for mutation	0.03
Maximum generations per run	25
Maximum number of function evaluations	10,000

and the description of the terms used is given in Supplementary Material Section S1. The optimization starts with an initial set of randomly generated N population (or chromosomes) made of genes representing the different values for the LIDs’ type, number, size, and location throughout the watershed. The decision variables were represented using binary encoding and decoding to facilitate data exchange between SWMM and MOGA. We used 10-bit chromosomes with 1024 (210) possible solutions as this avoids the duplicity of the solutions while providing sufficient variations to represent the decision variables. After the SWMM’s performance evaluation, crossover and mutation operators were applied to the initial population and generated the next population. First, better-performing individuals were identified as potential parents to produce offspring for the next generation. Then, some of the bites or genes in those individuals were randomly changed (mutated) before mixing them or applying crossover. Such selection and random gene changes ensure the survival of the fittest individuals and diversity that could lead to globally optimal solutions. Each chromosome bit was assigned a random real number ranging from 0.0 to 1.0. The mutation alters one more bit corresponding to the random numbers less than the specified mutation probability (0.001 to 0.05). In general, the mutation probability should be set to a low value, typically less than 0.5, to avoid the algorithm from switching to a random search. The optimization was done for four different design storms, representing the past (HIST 50 and HIST 100) and future (CC 50 and CC 100) design storms corresponding to the 50-year and 100-year return periods. The MOGA provides both the dominated and non-dominated solutions, which are used to develop the Pareto tradeoff between flow reduction and cost. The best solution from the Pareto front is identified as the point where there is a sudden change in the front, which is often known as the knee-of-curve point. This was done by estimating and comparing the slopes at each point on the Pareto front.

2.5.1. Implementation Strategies for New LIDs

The implementation of new LIDs was determined based on the availability of spaces, which is expressed as the LID area ratio (or ratio of LID area and available open spaces in each sub-catchment). Current land cover categories (Figure 4) obtained from the city of Renton and King County include mostly commercial, multifamily, high density residential, medium density residential, low density residential, grass, forest, and wetland. Multifamily areas are defined by attached housing with greater than seven dwelling units per acre, high-density areas have four to seven dwelling units per acre, medium-density have one to three dwelling units per acre, and low-density have less than one dwelling unit per acre. The blocks in the sub-catchments were categorized into commercial and residential zones (Figure 5) to facilitate the LID implementations. The residential zones contain single-family areas, detached, and attached multi-units, covering 470 sub-catchments with a recommended runoff coefficient of 0.25–0.35. The residential zone was further categorized based on the soil type (Figure 4), as sandy and heavy soil with no crop, crop, pasture, and woodlands. The commercial zone, which covers 136 sub-catchments, consists of land use with a high percentage of impervious areas and runoff coefficient ranging from 0.70–0.95. The land cover and availability were used to constrain the choice of LIDs. For example, where there is not enough available land for LID expansion, RGs were chosen because of their minimum space requirements and low installation and maintenance costs. RBs were chosen for the private area as they take relatively smaller space but have higher costs of implementation. On the other hand, IT and BR are most suitable for sub-basins with parks and schools. The PP is a candidate to replace roads, sidewalks, and parking spaces in the industrial zone. In addition, the soil types were also used as constraints, with BR and IT recommended for area with sand, clay, and clay loam, while RG were considered for sites with 30% sandy soil and mixed with loamy topsoil to obtain proper drainage conditions.

The IT area ranges from 10 m² to 300 m² depending on the percentage of impervious area within a sub-watershed. They are typically implemented within 20% to 50% impervious area and can take up to 50% of the front yard and 50% of the driveway of private properties, while the width was considered to be between 1.0 and 1.5 m. RG was applied in areas with land slope below 4%, and its depth ranges from 75 mm to 150 mm based on the surface topography. The soil conditions also influence the RG sizes, with a poorly drained soil requiring a larger surface area with shallow depth, while a well-drained soil requires a smaller surface area with greater depth. The RG is considered to take 50% of the impervious runoff from driveways and 25% of the runoff from rooftops. The other 25% runoff from rooftops will be routed to RB, where there are 0 to 4 RB per house. Generally, the BR is more effective in areas with underlying clay or sandy soil.

The majority of the study area was considered moisture sensitive because of the dominant silt–loam soil and thicker topsoil in the area. During the wet season (October to March), an increase in soil moisture content can cause significant reduction in the infiltration capabilities. In addition, the wet soil dries at a slow rate, affecting the soil permeability. Thus, BR in our study area needs an underdrain at the bottom to allow more infiltration and drainage of the soil layers. PP (porous concrete and asphalt) were used for roads, private driveways, and parking lots. A commercial parking lot requires 7 cm thick permeable concrete pavers or 7–18 cm thick porous asphalt depending on vehicle types and runoff volume [73]. The storage and infiltration capability of PP generally depends on the number of LIDs per house and on the type of pavement. In general, 0 to 1 units of RG, BR, and PP were used per house. The different LIDs that currently exist in the study area and are used in the SWMM model are shown in Table S1. The potential locations of the different LIDs based on the above criteria (i.e., prior to the optimization) are shown in Figure 5.

Before optimizing the LID implementation, we tested the potential effectiveness of adding LIDs into the existing stormwater system and LIDs to reduce the historical and future climate-change scenario peak and total flow in the study area. We have conducted a preliminary study using six hypothetical LID scenarios, consisting of individual LIDs and combinations of different LIDs (

Table 4. Hypothetical LID scenarios to evaluate the effectiveness of redesigning the existing LIDs.

Scenario		Description
S1	RG only	Use only rain garden (RG) to collect runoff from the houses.
S2	IT + RB	Use infiltration trench (IT) and rain barrel (RB) to collect the runoff from roof and surrounding areas of houses.
S3	PP	Use permeable pavement to collect runoff from commercial regions.
S4	PP + BR	Use permeable pavement and bioretention (BR) to collect runoff from commercial regions.
S5	BR only	Use only bioretention to collect runoff from small commercial areas.
S6	RB + BR + IT	Use combination of rain barrel (RB), bioretention (BR), and infiltration (IT) trench to collect runoff from both private and commercial areas.

).

2.5.2. Objective Functions

In general, the multi-objective function can be represented as follows:

$$F\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right),\dots ..f_n\left(x\right)\right)$$

(12)

where $\left(x_1,x_2,\dots .x_n\right)$ are the design variables and $f\left(x\right)$ is the different objective functions. The multi-objective optimization problems can be minimizing or maximizing the above objective function and are solved iteratively by identifying sets of solutions comprised of the trade-off between the various objectives. In this study, we used three objective functions, including total flood volume and peak-flow discharge at the basin outlet, and construction costs of LIDs. The primary goal of the optimization process is to identify combinations of LIDs with variable sizes and locations that optimally decrease the total flow and peak flow at a minimal cost of construction. Mathematically, the three objective functions are given in Equations (13)–(15).

$$\mathrm{Objective} 1 (\cos \mathrm{t} \mathrm{minimization}): min{\sum }_{J=1}^m{\sum }_{k=1}^n{x}_j^k\left(S,N\right)$$

(13)

where $x_j^k$ represents the cost of LIDs type k per lot within sub-catchment j. This cost varies with the LID volume S and number N in the sub-catchment; n is total number of the different LIDs within the sub-catchment and m is the total number of sub-catchments. The average cost of each unit was estimated based on the standard cost presented in Table 2.

$$\mathrm{Objective} 2 (\mathrm{peak} \mathrm{flow} \mathrm{minimization}): min\left(Q_{p_i}\right)$$

(14)

$$\mathrm{Objective} 3 (\mathrm{total} \mathrm{flow} \mathrm{minimization}): min{\sum }_{t=0}^k{\sum }_{i=1}^n{R}_i$$

(15)

where $Q_{p_i}$ and $R_i$ are the peak flow rate and flow volume, respectively, generated from the design storm from sub-catchment i. The total runoff volume is the sum of the volumes from each sub-catchment and over time.

2.5.3. Decision Variables

The decision variables include type, size, location, and number of LIDs. As described in

Table 5. Sensitivity analysis of SWMM parameters on peak flow.

Parameter Test	Percentage Change (%)	Reduction in Peak Flow (%)	Percentage Change (%)	Reduction Peak Flow	Percetage Change (%)	Reduction in Peak Flow (%)
Zero Impervious	5	−13.89%	10	−16.9%	15	−22.44%
Width	5	−10.75%	10	−14.60%	15	−18.31%
Manning n	5	−15.20%	10	−19.86%	15	−23.51%
CN	5	−17.29%	10	−20.69%	15	−24.74%

, the implementation of these decision variables was constrained by soil type, land cover, available space, and other conditions. The list of decision variables, constraints, and feasible regions for number and sizes of LIDs are provided in Table S2. The feasible regions for the decision variables are defined as the sets of solutions that fulfill any constraints imposed on the variables [74].

3. Results and Discussions

3.1. Climate Change Impact on Future Storms

The observed and downscaled maximum daily precipitation from the MIROC6 and CMCC-ESM2 climate models are shown in Figure S1. Overall, the downscaled precipitation from the MIROC6 shows better performance compared to the CMCC-ESM2 in reproducing the observed precipitation. The MIROC6 has R2 and NSE values of 0.96 and 0.91, while CMCC-ESM2 has 0.91 and 0.77, respectively.

Additionally, Figure S2 compares the frequency analysis results for observed and downscaled annual maximum daily precipitation using the three downscaling methods for the MIROC6 and CMCC-ESM2 projections. The results show that in both cases the EQM method outperforms the other methods because of its ability to consider the different quantiles of the precipitation distribution separately during the bias correction. The delta method only considered the mean of the distribution and often underestimated the extreme precipitation. The MIROC6 also performs better than CMCC-ESM2, capturing the observed extreme precipitation.

The design storms, estimated from the IDF curves, are commonly used to design and manage stormwater systems. The 15-min precipitation data, collected from the city of Renton database (https://green2.kingcounty.gov/hydrology/GaugeMap.aspx, 1 January 1995 through 1 May 2014), was used to develop the historical IDF curves. The bias-corrected daily maximum precipitation from the MIROC6 model was converted to 15 min precipitation using the SCS Type IA distribution and used to develop the future IDF curves. Using these two datasets, the IDF curves for 2-, 5-, 10-, 20-, 50-, and 100-year return periods are constructed for both historical and future-scenario periods (Figure S3). The results show an increase (2.5% to 30%) in the storm intensities in the future scenario for all the return periods, with the increase being relatively more significant for higher return periods like the 50-year and 100-year storms.

3.2. Sensitivity Analysis

The detailed results of sensitivity analysis of SWMM are given in Table 5. In the present study, 5–15% perturbations in the sub-catchment’s width reduced the peak flow by 10.75–22.44%, whereas identical changes in CN reduced the peak flow from 17.29–24.74%. Sub-basins under consideration in this study have slopes within the shallow range. However, the slope parameter was not considered for the model calibration, but an increase in the depth of depression storage leads to a delayed flow after a storm event, minimizing total runoff and peak flows [69,71,73]. Zero percentage (impervious area with no depression storage) minimized the peak flow from 13.89–22.44% on 5–15% perturbance while the percent imperviousness decreased the peak flow by 15.25 to 23.51%. In this study, perturbations in CN indicated a significant impact on the peak flow while percent was the least sensitive parameter. Factoring other parameters under consideration, decreasing the catchment’s width leads to a reduction in peak flows, makes peaks less abrupt, and intensifies the hydrograph trailing arm. High variations in curve number magnitudes can be considered significant, similar to the effect of the percent imperviousness, while reduction in the percent width parameter effectively increases overland flow length and overland flow travel time.

3.3. Calibration and Validation of the SWMM Model

The SWMM model was calibrated and validated manually using the observed flow from four different rainfall events. Two rainfall events were used for manual calibration (Figure S4a,b), and the other two were used for validation of the model (Figure S4c,d). The calibration and validation resulted in good agreements between simulated and observed hydrographs, with NSE statistics ranging from 0.79 to 0.84 and R2 ranging from 0.81 to 0.87 for the four storms.

3.4. Potential of LID to Reduce Floods

Prior to optimizing the LIDs, we tested their potential to reduce floods in our study area using six hypothetical LID implementation scenarios. Figure 6 shows the resulting hydrographs from the 50-year and 100-year design storms for the historical period (HIST 50 and HIST 100) and future-scenario period (CC 50 and CC100) compared with the baseline results that considered only the existing stormwater system. All the scenarios reduced the peak flows from historical and future-scenario storms. Overall, the peak flow and total runoff volume can be reduced by up to 13% and 29% under future climate-change scenarios and by up to 62% and 80% under the current storm conditions. Comparing the five individual LIDs, the infiltration trenches (IT) performed better in reducing the peak flow and runoff volume per construction cost, while the rain garden (RG) and bioretention (BR), respectively, are less effective in reducing the peak flows and total flow. The maximum reduction of runoff volume (80%) and peak flow (65%) were observed when a combination of rain barrel (RB), bioretention (BR), and infiltration trench (IT) was used (S6). The moderate performance of the bioretention (S5) was due to its small surface areas, causing overtop of stormwater when receiving a large portion of the surface water from the impervious area. The combination of infiltration trench and rain barrel (S2) provides a peak-flow hydrograph slightly lower than the individual LID under 50- and 100-year historical events. The peak discharge time for rain garden (S1), bioretention (S5), and mix of infiltration trench and rain barrel (S2) is nearly identical to the baseline scenario, while the other scenarios have delayed the peak discharge time.

3.5. Cost–Benefit Curves for Peak Flow

The multi-objective optimization provides optimal solutions and their potential trade-off for different objectives. In this study, the multi-objective consists of maximizing the peak flow and volume reduction and minimizing the construction cost of different LIDs. The performance of a given LID can vary depending on its size, number, and location. Figure 7 presents the optimal trade-offs (Pareto curves) between the percentage of peak-flow reduction and the construction costs under the historical and future climate scenarios. Each point in the figure represents a solution associated with an LID with a given size, number, and location. The red points represent the most cost-effective solutions. It has been found that the cost-effective curves generally show a significant increase in flow reduction with added investment on the LIDs up to a certain point, also known as knee of curve (blue dot), after which the LIDs’ performance does not improve considerably with the added investment. Therefore, the knee of curve is considered as the optimum solution. With an approximately USD 5 million investment on LIDs, the peak flows from past 50- and 100-year storms can be reduced by about 70% and 45%, while the peak flows from future-scenario 50- and 100-year storms can be reduced by about 42% and 35%, respectively. The additional investment beyond the USD 5 million will not significantly reduce the peak flows. Declining investment returns occur for investing beyond the effective solutions with the increase in the peak-flow reduction can only be obtained by investing in less efficient LIDs. The level of peak-flow reduction that can be achieved before a significant point of diminishing returns is reached depends on the limitations of the most efficient LID types in each scenario.

The results also show that the cost-effective solution varies with the storm sizes. For the HIST 50 storm, a 62.5% maximum reduction of peak flow was achieved with LID costing USD 3.9 million (Figure 7a). Whereas, for the HIST 100 storm, the maximum peak-flow reduction achieved was about 42% with a cost of USD 4.1 million (Figure 7b). These indicate that the LID practices performed differently for different design storms, and, as expected, flow reduction is higher for design storms with lower recurrence intervals (or moderate storms). For the future-scenario storms (CC100 and CC50), the peak flows are reduced by up to 33% and 42% at a total cost of USD 5.20 and USD 4.78 million, respectively (Figure 7c,d).

Figure 8 is breakdown to some of the LID mixes and decision variables contained in the cost-effective curves (red dots in Figure 7). The infiltration trench (IT) and its mix with other LIDs, such as rain gardens (RG), dominate the optimum solutions. The infiltration trenches (IT), with sizes ranging from 10 m² to 300 m² depending on the percentage of impervious land in the sub-catchment, are the most efficient and effective that can allow for further reduction of the peak flow by adding a rain garden (RG). Most trenches need continuously perforated pipes and relatively flat topography (slope < 1%). Higher slope reduces the infiltration and requires longer trenches at higher cost. Rain barrels (RB) are also present in some of the solutions but they are less effective at reducing peak flows than they are at reducing total runoff because they fill up before the most intense part of the precipitation event. Overall, a reduction of up to 50% peak flow from the past storm (HIST 50) was achieved with the use of IT and few RB that cost below USD 96,000. To further reduce the peak flow by up to 76%, the IT needs to be mixed with RG that costs about USD 269,000. The RB is relatively less effective in reducing peak flow from the future-scenario storm (CC-100) than HIST-50, achieving a maximum reduction of 30% with the related cost of USD 131,000. Although the combined uses of IT, RG, and RB provide high-flow reduction, they incur significantly high costs.

3.6. Cost–Benefit Curves for Total Runoff

The multi-objective solutions and the cost-effective curves for total runoff reduction are shown in Figure 9. The solutions are different from the peak-flow reduction solutions. Compared to the existing system, the optimized solutions for the HIST 50 and HIST 100 storms show a total runoff reduction of 70% and 58.2% with corresponding costs of USD 2.67 million and USD 3.5 million, respectively (Figure 9b,c). The reduction of total runoff for a CC 100 storm is around 6% lower than that of the CC 50, while the total cost is significantly increased because of the expected higher storm intensity for a CC 100 storm compared to the CC 50 storm.

The LIDs’ performances in reducing the total runoff were further evaluated using selected decision variables in Figure 10. The result shows that the increasing rate in cost-effective curves for the total flow reduction is steeper than that of the peak-flow reduction. This indicates that achieving a peak-flow reduction above the knee-of-curve point generally requires a higher cost than the total runoff reduction for all the four storm conditions we considered. In addition, the percentage reduction for total runoff for the identified LID practices is better overall than that of the peak runoff. Similar to the peak-flow result, the IT provides the most optimum cost–benefit solution for the total runoff reduction. Additionally, RB and RG are frequently found in a mix with IT to reduce the total runoff effectively. A higher reduction in flood volume was observed for infiltration trenches and permeable pavement (IT + PP) under the 100-year future-scenario storm, but they incur significantly increased costs (Figure 10).

In summary, based on the knee-of-curve results, it was observed that:

• Nearly 70% of the runoff volume from the 50-year past storm can be reduced using a combination of 30 ITs with a surface-area range of 4–80 m², 70 RGs with a surface-area range of 4–10 m², 18 BRs with a surface-area range of 4–10 m², 9 PPs for driveways with a surface-area range of 20–80 m², 100 RB with commercial size (200 L) at a total cost of USD 2.67 million;
• Approximately 58% of runoff volume from a 100-year past storm can be reduced using 79 ITs with a surface-area range of 10–80 m², 77 RGs with a surface-area range of 8–10 m², 12 PPs for driveways with a surface-area range of 20–80 m², and 150 RBs with the commercial size 200–300 L at the cost of USD 3.5 million.
• Approximately 50% of runoff volume from a 50-year past storm can be reduced using 86 ITs with a surface-area range of 25–80 m², 120 RGs with a surface-area range of 8–10 m², 80 BRs with a surface-area range of 4–10 m², 17 PPs for driveways with a surface-area range of 20–80 m², 168 RBs with commercial size 200 L at the cost of USD 3.5 million dollars.
• The peak flows from the past 50-year and 100-year storms were decreased by 63% and 61%, with LIDs costing about USD 3.9 million and USD 4.3 million, respectively. Meanwhile, for the future-scenario 50-year and 100-year storms, the peak flow was reduced by 42% and 31% with USD 4.8 million and 5.2 million total costs, respectively.

4. Conclusions and Recommendations

The expansion of urban areas and increased extreme storms caused by climate change are causing severe floods to become more frequent in most regions. Many cities around the globe recognize the increasing trends in severe floods and have started upgrading the existing stormwater systems with relatively effective and sustainable LIDs. However, the effectiveness and cost of LIDs vary significantly depending on the characteristics of watersheds and design storms. For example, the widely used infiltration trench is only effective, both hydrologically and in cost, in areas where the topography has a gentle slope and the soil is deep with higher hydraulic conductivity. Thus, in most regions, mixing different LIDs and existing conventional stormwater systems and optimizing the size, location, and number of LIDS are essential to improve their performance and minimize cost.

The integrated SWMM–MOGA framework was first used to assess the performance of commonly used LIDs and then identify the hydrologically and cost-effective mixes of LIDS under both past- and future-scenario storm conditions in Renton City, WA. After confirming the benefit of adding LIDs to the existing stormwater system, we have optimized their applications by considering their sizes, locations, numbers, and watershed conditions. The reduction in peak flow, total runoff volume, and construction costs were used as objective functions for the optimization.

The results indicated that incorporating LIDs can significantly reduce the peak flow and total runoff in our study area. Overall, they appeared to be more effective in reducing total runoff than the peak flows, with the total runoff and peak-flow reductions ranging from 42% to 70% and 33% to 62.5% for the future 100-year and past 50-year storms, respectively. Combined uses of rain barrels, bioretention, and infiltration trenches represent the most effective LIDs to reduce the peak flow and volume in our study area. The effectiveness of a given type of LID is nonlinearly related to its design parameters, the combination of different LID practices, and storm magnitudes. For example, the optimized solution for the 100-year future storm indicated that a 50% runoff reduction (or 1200 m³) can be achieved with an implementation cost of USD 91,500. It costs three times more (USD 341,000) to further reduce the runoff by 2400 m³ from the same storm event.

The simulation–optimization framework from the study can be a useful decision-making tool to quantify the hydrological- and cost-effectiveness of different LIDs. It can assess the effect of different LIDs’ design parameters (e.g., placement and sizing of LIDs) on their performance and optimize the integrated LID and traditional stormwater system under different storm and urban development scenarios. Such a modeling framework allows the identification of an optimal LID implementation strategy for adapting to future climate change and providing quantitative information for urban planners, stormwater engineers, decision-makers, and stakeholders to implement an efficient and sustainable stormwater management system.

5. Limitations and Future Aspects

The presented study considered only a few popular LIDs in our study area. Including the other less commonly used LIDs (detention vaults and ponds) from the current stormwater system may further improve the overall performance. The water-quality benefits of the LIDs were not investigated in this study. Including water quality in the future-scenario study will give the decision-makers complete cost–benefit information on their investment in LIDs. In addition, when groundwater interferes with the LIDs because of a shallow water table, the model (which assumes conditions favorable to stormwater infiltration at LID sites) provides a poor simulation of the actual field conditions. This issue can be addressed to some level by adding an impervious layer underneath the LIDs to limit the infiltration. In addition, improving the modeling framework to account for the different sources of uncertainty will enhance its use as a decision-making tool. In modeling the future scenario, we only considered changes in storm events while considering the same land cover and static performance by LIDs. However, the continuously growing urban areas will impact the land cover by expanding the impervious surface. The current version of the framework does not explicitly account for the potential uncertainty in the climate projections. However, the framework can be run under different climate projections, representing the uncertainty range, to generate alternative implementations of LIDs. Finally, the study did not consider the flood risk caused by a sea-level rise since such a flood risk is relatively low in the city.