Integrating the Spatial Configurations of Green and Gray Infrastructure in Urban Stormwater Networks

Green infrastructure (GI) practices improve stormwater quality and reduce urban flooding, but as urban hydrology is highly controlled by its associated gray infrastructure (e.g., stormwater pipe network), GI's watershed‐scale performance depends on its siting within its associated watershed. Although many stormwater practitioners have begun considering GI's spatial configuration within a larger watershed, few approaches allow for flexible scenario exploration, which can untangle GI's interaction with gray infrastructure network and assess its effects on watershed hydrology. To address the gap in integrated gray‐green infrastructure planning, we used an exploratory model to examine gray‐green infrastructure performance using synthetic stormwater networks with varying degrees of flow path meandering, informed by analysis on stormwater networks from the Minneapolis‐St. Paul Metropolitan Area, MN, USA. Superimposed with different coverage and placements of GI (e.g., bioretention cells), these gray‐green stormwater networks are then subjected to different rainfall intensities within Environmental Protection Agency's Storm Water Management Model to simulate their hydrological benefits (e.g., peak flow reduction, flood reduction). Although only limited choices of green and gray infrastructure were explored, the results show that the gray infrastructure's spatial configuration can introduce tradeoffs between increased peak flow and increased flooding, and further interacts with GI coverage and placement to reduce peak flow and flooding at low rainfall intensity. However, as rainfall intensifies, GI ceases to reduce peak flow. For integrated gray‐green infrastructure planning, our results suggest that physical constraints of the stormwater networks and the range of rainfall intensities must be considered when implementing GI.

. Second, although peak flow reduction and base flow increase are also observed in catchment studies (Bhaskar et al., 2016;Jarden et al., 2016;Page et al., 2015;Pennino et al., 2016;Shuster & Rhea, 2013), the cumulative effects of GI are nonlinear and only become observable when there is sufficient density of GI within the study catchment (Jefferson et al., 2017).Third, the watershed-scale research needed for studying the effects of GI in the field is costly to fund and difficult to coordinate, often requiring both financial and political support from property owners, city residents, and regulatory agencies (Torres et al., 2021).
Finally, upscaling the effects of GI requires careful consideration for the environmental and spatial contexts of GI.For example, at the catchment level, GI's marginal runoff reduction depends on external factors (e.g., the rainfall regimes, site conditions, GI types), and the benefits of GI may plateau (i.e., diminishing return) at higher GI coverage (Ahiablame et al., 2012;Kumar et al., 2022;Moore et al., 2016).In more intense storms, higher coverage may be needed before the effects of peak flow reduction saturate (Tansar et al., 2022;Zellner et al., 2016).Additionally, the spatial placement of GI within the catchment will also play a role during upscaling.For example, when dispersed across the watershed, GI's storage capacity may better capture runoff, whereas a clustered placement may under-utilize GI's excess storage capacity (Zellner et al., 2016).When placed in downstream areas closer to the outlet, GI may reduce the peak flow (Chang et al., 2008;Giacomoni & Joseph, 2017;Vittorio & Ahiablame, 2015) at the expense of worsening flood risk (Fiori & Volpi, 2020).In practice, however, municipality's decision of where and what type of new GI to construct is limited by competing stakeholder interests and land availability.Due to these constraints, it is difficult to untangle the relationship between GI and gray infrastructure's spatial arrangement through empirical data alone.
To bypass these challenges associated with upscaling the impacts of GI, we adopt a computational approach that precisely controls for the range of local variability that may confound catchment-level outcomes.By exploring the impact of gray and GI's spatial configurations, we aim to untangle the processes through which GI's placements affect watershed hydrology while considering spatial variations of gray infrastructure networks.In this study, we have developed a simple and flexible gray-green infrastructure computational model with Python-using bioretention cells to represent GI and stormwater pipe networks to represent gray infrastructure, understanding that additional stormwater storage (e.g., BMPs) and a more diverse selection of GI (e.g., rain barrels, green roofs) will contribute greatly to the hydrological outcomes in an urban watershed.The gray infrastructure, represented by urban stormwater drainage networks, is informed by network GIS data from the Minneapolis-St.Paul Twin Cities Metropolitan Area, MN, USA.We use tools from graph theory to classify and reconstruct these networks, which can be stripped down to their most basic components-nodes and edges-that store essential attributes for stormwater modeling.These reconstructed graphs do not rely on any specific network layouts, which are difficult to access from municipalities, yet are flexible enough to capture a wide range of networks with similar spatial configurations represented by flow paths.We stochastically generate stormwater networks based on these generalizable layouts, then model their hydrological responses with different number and placement of GI within the network under a range of rainfall scenarios.
The simulations from this model allow us to explore the spatial interactions between green (i.e., bioretention cells) and gray infrastructure (i.e., stormwater pipe networks) for different rainfall scenarios with increasing intensity.In particular, we focus on three aspects of catchment level gray-green infrastructure planning: (a) the layout of the stormwater networks onto which GI is placed, (b) the density of GI in the network, and (c) the spatial placement of GI within the network.We hypothesize that (a) increasing stormwater network flow path meandering decreases peak flow at the expense of increasing inland flooding, (b) the benefits of GI will decrease with increasing storm intensity, and (c) downstream placement of GI will more efficiently reduce peak flow at the outlet.We anticipate that these results will facilitate stormwater infrastructure planning by explicitly accounting for GI placement within the spatial configuration of stormwater networks.

Methodology
With tools from stochastic processes, we generate a range of graphs to represent synthetic stormwater networks based on real networks derived from shapefiles of Twin Cities Metropolitan Area stormwater infrastructure.We then place GI on these networks.These combined gray-green infrastructure networks vary across three design parameters: (a) the flow path meandering H(s), which measures how much the network's actual stormwater flow paths deviate from the shortest paths, (b) the coverage of GI in the watershed, measured by the proportion of The summary of experimental setup is shown in Figure 1.We first reconstruct existing stormwater networks from municipalities onto square lattice grids to get realistic representation of flow paths H(s) values for our study.The simulation study is then conducted through three experiments (Figure 1a).Experiment 1 investigates the effects of varying network flow path meandering on hydrological outcomes under two different land covers.Experiments 2 and 3 investigate the effects of GI coverage and placement respectively, in low and high H(s) networks.For each experiment, we then simulate each network's hydrological response using the Environmental Protection Agency's Storm Water Management Model (SWMM) across a range of rainfall intensities.We use two intermediate functional variables (Figure 1a) to explain the hydrological outcomes observed from SWMM modeling: (a) network storage capacity, which represents the total stormwater storage volume in both gray and GI; (b) transport efficiency, which relates to the runoff flow path lengths and measures how quickly runoff is removed.Changes in the functional variables will directly affect our desired hydrological outcomes.It is important to note that network storage capacity in GI represents the above-ground storage within treatment practices, and that network storage capacity in gray infrastructure is the storage within the pipe network itself.In particular, increase in any of the functional variables will increase peak flow and decrease flooding volume, resulting in a tradeoff between managing these two hydrological outcomes.For example, pipe conveyance induced storage capacity increase reduces flooding due to reduction in runoff, but it may increase peak flow by allowing more flow to be transported through the network, if no proper rate control devices were installed.Higher transport efficiency will reduce flooding because water can be more quickly conveyed to the outlet, but it will increase the peak flow rate for the same reason.The workflow diagram (Figure 1b) outlines the individual steps for programming the simulation study, with the permanent version of the accompanying Python scripts available on Zenodo (Chen & Feng, 2023).Below, we elaborate on the methods used to generate urban stormwater networks, the design parameters for the experiments, and the specific implementation in SWMM.

Network Flow Path Meandering
We represent stormwater networks mathematically using directed tree graphs on square lattice grids (Text S1 in Supporting Information S1).We use a spanning tree s to represent the realized stormwater network, where its flow paths are unique.In a stormwater network, we are interested in the flow paths to the network outlet v 0 .Here, we denote d G (u, v 0 ) to be the flow path length of the shortest path from node u to the outlet v 0 in a connected square lattice graph G, from which the spanning trees are constructed (Figure 2).In a spanning tree s, if the actual path between nodes u and v 0 , (u, v 0 ), follows the shortest path, then d s (u, v 0 ) − d G (u, v 0 ) = 0. Seo and Schmidt (2012) and Troutman and Karlinger (1992) quantify the network flow path variation by summing all the flow path differences within a network: Although intuitive, this parameter H T (s) does not account for the total size of the lattice grids.As a result, H T (s) is not normalized when comparing flow paths among networks generated from grids of different sizes.We address this shortcoming by normalizing the flow paths by the shortest flow path lengths.Specifically, for a single path from each node in a network to the network outlet, we define flow path variation as the difference between the actual path in s and the shortest path in G, d s (u, v 0 ) − d G (u, v 0 ), normalized by the shortest path, d G (u, v 0 ).We can quantify the network flow path meandering across all nodes in a spanning tree s as H(s), by normalizing the flow path variation from any node on the lattice to the network outlet (Equation 2): The meandering in flow path H(s) is a design parameter that quantifies the spatial configurations of the stormwater network.By using H(s), we can explore the effects of network flow paths on hydrological outcomes, without using any particular municipal networks.We hypothesize that higher H(s) will decrease transport efficiency while increasing network storage capacity.Because the changes in these two functional variables will have the opposite effects on controlling peak flow and flooding, we expect to see a tradeoff between managing these two hydrological outcomes (Hypothesis 1).

Stochastic Generation of Networks
The spanning trees, which are tree graphs that share the same nodes as a n-by-n square lattice graph G, are sampled following the simulated annealing method using a stochastic Markov chain (Aldous, 1987;Seo & Schmidt, 2012;  (Troutman & Karlinger, 1992) and network flow path meandering H(s), and spanning tree t with higher H T (t) and H(t).In the two spanning trees, the nodes represent access points where runoff enters the stormwater network, the larger node at the bottom left corner is the outlet node, the arrows on the edges designate the flow direction, and the thickened edges highlight the flow path from an arbitrary node to the outlet node.The arbitrary node is at the same location on the lattice graph G.The shortest path from such arbitrary node to the outlet is highlighted in graph G. Spanning tree s follows shorter paths and has lower H T (s) and H(s), as the highlighted flow path follows the shortest path.Spanning tree t has higher H T (t) and H(t), as, for example, the highlighted path is much longer than the shortest path.Troutman & Karlinger, 1992).This method results in a stationary probability distribution for the occurrence of a spanning tree that follows a Gibbs distribution (Equation 3).In the Markov chain, the transition probability during each sampling step is defined by the flow path of the current spanning tree H T (s 1 ), the flow path of the proposed spanning tree H T (s 2 ), and a prescribed parameter β (Seo & Schmidt, 2012;Troutman & Karlinger, 1992).The spanning trees with desired H(s) or H T (s) values across multiple prescribed β parameters are retained for further analysis.Network generation is done in Gibbs.py(Chen & Feng, 2023), where the outputs are networks without stormwater attributes (Figure 1b).
The frequency distribution of H(s) and H T (s) of the generated spanning trees s of graph G is verified to follow the Gibbs distribution (Figures S1 and S2 in Supporting Information S1) given by: where H T (s) is the flow path meandering calculated for each graph generated (Equation 1), and C(β) is a normalization coefficient so that the probability density function integrates to 1.More commonly seen in statistical mechanics, Gibbs distribution gives the probability that a system will be in a certain state as a function of that state's energy and the temperature of the system (Landau & Lifshitz, 1980).In the case of spanning trees, the graph G is the system, each spanning tree s generated from G is analogous to a state of the system, and the flow path meandering H(s) of the tree are analogous to the energy of the state (Troutman & Karlinger, 1992).
Here, the single parameter β is used as a proxy for the inverse of thermal energy, where a "warmer" system is associated with a lower β.In our case, β controls the tendency to which the realized flow paths can deviate from the shortest flow paths on G.When β is 0, the sampling distribution of H T (s) is uniform, which means the "cooling power" is low, resulting in a larger range of deviation in the network (or the energy in the network state) with equal likelihood of high and low deviations.When β is large, the resulting samples are more likely to have lower H(s), because high β lowers the overall energy of the system and the degrees of deviation (Troutman & Karlinger, 1992).

Reconstructing Networks From Real Municipal Stormwater Networks
To estimate a reasonable range of H(s) values, catchments from municipal stormwater networks from the Twin Cities Metropolitan Area were delineated and reconstructed on lattice grids.The geospatial information system (GIS) data were provided by local watershed management organizations, including Washington Conservation District, South Washington Watershed District, Capitol Region Watershed District, and Mississippi Water Management Organization, who had previously compiled databases of stormwater networks from constituent municipalities.In particular, we worked with stormwater network shapefiles from the Cities of Stillwater, Woodbury, Lauderdale and Saint Paul, MN.
First, we delineated 20 smaller catchments (i.e., 0.2-1.5 km 2 ) from the stormwater networks, often by city blocks in urban areas or by subdevelopments in suburban neighborhoods.When fitting these networks onto grids, the network nodes are determined by street intersections where applicable, and the distance between nodes determines the edge lengths of networks.When the networks do not obviously follow perpendicular grid structures, they are placed on lattice grids in a way that preserves the physical proximity among upstream branches.For example, if two upstream nodes from separate branches are close to each other on a physical map, the flow paths will meander in such a way that preserves their geographical proximity.The coordinates of the grid-based networks are then used to calculate the flow meandering parameter H(s).

Assigning Stormwater Network Attributes
Attributes for stormwater nodes and edges are populated after the synthetic spanning trees are generated, following the schematic outlined in Figure 1b.These attributes are tabulated in Table 1 and assigned in script hydro_ network.py(Chen & Feng, 2023), and the output is the stormwater network that will be used in SWMM simulation.In particular, nodes in the network represent both access points in the stormwater network and the subcatchments where the runoff drains from, and the corresponding subcatchment attributes are stored at the node.The directed edges represent the stormwater pipes and their flow directions.We assume all subcatchments to have the same design attributes, such as drainage area, impervious area and slope.To keep this exploratory model simple, the nodes do not have any ponding area or other above-ground storage.Therefore, when water overflows a node, the overflow will not return to the network after pipes regain their capacities, and the amount lost to overflow is counted toward flooding loss.Although more conservative, our flooding estimates still inform the relative magnitudes of flood control outcomes across different simulation scenarios.
Some nodes are designated as GI nodes, and all the GI nodes in the network share the same design attributes.The parameters for GI are selected from the self-reported design parameters for bioretention cells from the International BMP Database (Clary et al., 2020).Each GI node only receives direct runoff from its associated subcatchment, and its sizing is not impacted by cumulative upstream stormwater flow in the stormwater network.In the subcatchments that contain GI, the GI practice occupies 5% of the total subcatchment area, sized uniformly at 404 km 2 .
Although all pipes share some physical properties (e.g., pipe slope, roughness, and pipe length), sizing stormwater pipes is a non-trivial consideration as it determines the total capacity of the stormwater networks, and in turn affects how much and how efficiently stormwater flow can be carried downstream.We use the rational method to size pipe diameter, so that the pipes can contain the peak cumulative runoff generated by a design storm.The rational method is given by Q = CiA, which calculates the peak runoff according to a runoff surface coefficient (C), rainfall intensity (i), and the subcatchment area (A).To account for cumulative peak flows upstream to any stormwater pipe, we use  =  ∑  =1  , where there are m pipes upstream of the design pipe.C increases with higher proportion of impervious drainage area.We use a 10-year 24-hr storm based on Atlas 14 precipitation frequency estimates (National Oceanic and Atmospheric Administration, 2013) to size the network, and round up the pipe diameters to the nearest common pipe size.Networks with higher H(s) will have larger overall network storage capacity due to higher flow accumulation throughout the network.

Spatial Configuration of Green Infrastructure on Stormwater Networks
Once the network structure, with its design parameter H(s), and the network attributes are determined, we investigate the effects of GI placements.We expect that adding GI to existing gray infrastructure will modify our functional variables in a number of ways.It will increase useable network storage capacity by reducing the amount of runoff entering the stormwater network, but adding GI reduces the transport efficiency as added storage effectively elongates runoff flow paths.
We also consider the effects of GI density and GI's distance to network outlet.In Experiment 2, we place GI onto the network nodes at different degrees of coverage.A given percentage of nodes are converted to GI nodes and randomly placed across the stormwater network, and this percentage of coverage is increased at fixed increments.In our study, we start with a base case scenario with no GI, and gradually add 10% up to a maximum of 50%.In Experiment 3, we place the GI, which is 10% of the total nodes, in a cluster at different distances from the outlet to test the effects of GI's spatial configuration.We start by placing the GI cluster near the stormwater network outlet, and gradually move them upstream.The GI nodes are clustered together to minimize the confounding effects of spread and network flow path lengths.In our model, we fix the number of GI nodes to 10% of the total nodes.We have hypothesized that downstream clusters of GI may reduce peak flow in the network (Hypothesis 3).

SWMM Implementation and Rainfall Scenarios
This study uses Stormwater Management Model (SWMM), a distributed hydrologic-hydraulic model (United States Environmental Protection Agency, 2022), for all three experiments outlined in Figure 1.By accounting for both the spatial heterogeneity of GI and the hydraulics of the stormwater network, SWMM can accurately represent the transport of urban runoff from inland to streams in a gray-green infrastructure network.We use the bioretention cell setting in SWMM Low Impact Development (LID) module to parameterize the GI nodes.Although the SWMM LID module uses a one-dimensional infiltration model and faces constraints when accounting for transfer between surface water and groundwater, the bioretention model generally generates agreeable results with calibration (Avellaneda et al., 2017;Gülbaz & Kazezyilmaz-Alhan, 2017;Shaneyfelt et al., 2021).
The spatial and hydraulic attributes of the stormwater networks and the environmental conditions of the simulation are then populated to an .inpinput file, which is readable by SWMM engine, and the report summary files are parsed for data analysis.We use the dynamic wave routing method, which can account for the flow rate increase when stormwater flow in pipes goes from open channel flow to pressurized pipe flow.Although SWMM automatically adjusts routing time steps when the dynamic wave routing method is used, the initial routing step for the analysis is 2.5 s.The template for the input file is included in the code repository (Chen & Feng, 2023).
We evaluate the network performance using two hydrological outcomes (Figure 1a): the max flow rate at the stormwater network outfall to represent outlet peak flow, and the percent of flooding loss within the network with regards to total precipitation.In SWMM, flooding refers to all water that overflows a node (i.e., water exceeding storage capacity that would temporarily be stored above the ground surface).Due to the lack of above-ground storage in our exploratory model (previously mentioned in Section 2.4), stormwater overflow does not reenter the network even after the stormwater flow has subsided, resulting in a more conservative (higher) flooding estimate.Using the inland flooding volume and total precipitation volume from the SWMM summary report files, we can calculate the percent flooding loss as a proportion of total precipitation.These outcomes were analyzed using three SWMM output variables from the summary report files: the max flow rate in the outfall loading summary, flooding volumes in the node flooding summary, and total precipitation volume in the runoff quantity continuity summary.To ensure the calculated flood volume realistically represents the flow condition at each node and does not result from a numerical continuity error, the node flooding volume is only included if the flood duration exceeds 0.3 hr.SWMM implementation is done in script make_SWMM_inp.py(Chen & Feng, 2023).
In Experiment 1, we generated 1,000 networks from a 10-by-10 grid, and the H(s) values for each network is rounded down to the nearest 0.5 (i.e., 0, 0.5, 1, and 1.5).These 1,000 networks are sized based on two different land covers.The first case represents an urban business district with higher impervious area in the watershed.We used a runoff coefficient of 0.8 in sizing the pipes and 80% imperviousness in the SWMM simulation.This configuration is typical of a built-out business district in midsize U.S. city (e.g., the City of Minneapolis).The second case is a typical suburban residential neighborhood, which has a rational coefficient of 0.5% and 40% imperviousness (Gulliver et al., 2015;Mays, 2010).For Experiments 2 and 3, we generated 475 networks with low H(s) (H(s) ≈ 0.2) and 349 networks with high H(s) (H(s) ≈ 0.02), all sized to the urban, highly impervious land cover.In Experiment 2, GI nodes are randomly assigned to the network, at an increment of 10 nodes up to 50 nodes.In Experiment 3, clusters of 10 GI nodes are placed from downstream to upstream.
We hypothesize that high-intensity precipitation may reduce GI's effectiveness (Hypothesis 2).To test this hypothesis, we run the model under a few uniform rainfall scenarios for all three experiments.Although green and gray infrastructure is typically designed for a 24-hr storm, we selected the 2-hr storms to reduce the computational time, while preserving the same rainfall total as the two most intense hours of a typical 24-hr design storm, based on either a Soil Conservation Service Type II storm or a MSE 3 storm typical for Central Minnesota (Table S2 in Supporting Information S1).The storm amounts are assigned according to Atlas 14 point precipitation frequency estimates (National Oceanic and Atmospheric Administration, 2013), and they are tabulated in Table 2.The 2-year storm is selected because it is a typical design storm for GI practices (New Jersey Department of Environmental Protection, 2023), whereas the 10-, 25-, and 100-year are typical design storms for stormwater conveyance and flood control structures (Minnesota Pollution Control Agency, 2023).We uniformly distribute the rainfall amount over the 2-hr period to reduce computational instability and simulation time (Figure S4 in Supporting Information S1).
In the input hyetographs, we increase the rainfall intensity (i.e., return periods of 2-, 10-, 25-, 100-year) while maintaining the rainfall duration (i.e., 2 hr), with the precipitation amount uniformly distributed over 120 one-minute simulation intervals.The short hyetograph intervals closely match the routing steps to buffer the inputs and minimize the instability from any sudden large rainfall input to the stormwater network.Over a simulation period of 24 hr, the rainfall hyetograph begins 2 hr after the simulations begins, and lasts for 2 hr.The simulation ends after 20 hr of dry time to allow for steady state.

Results
We first select and reconstruct catchments from municipal stormwater networks onto lattice grids to validate the range of network flow path meandering H(s) we used in our experiments.In Experiment 1, we see that a tradeoff exists between balancing peak flow reduction and flood reduction as we increase network flow path meandering (Figure 4).The results from Experiments 2 and 3 show that at low rainfall intensities, GI and gray infrastructure work together to regulate the tradeoff between peak flow and flood volume.Green infrastructure's ability to mitigate the peak flow response saturates at high rainfall intensities, but GI remains effective for flood reduction (Figures 5 and 6).In the following sections, we elaborate on how rainfall scenarios interact with the spatial configuration of gray-green infrastructure to influence these hydrological outcomes (e.g., peak flow, flooding volume), through directly exploring the rainfall-dependent effectiveness of which our design parameters influence the two function variables (e.g., network storage capacity, transport efficiency).

Reconstructing H(s) From Real Stormwater Networks
Figure 3 shows the map of a subset of catchments selected for the analysis, and examples of reconstructing two of these catchments-one in downtown Stillwater, MN, one in a residential area of Woodbury, MN-onto lattice grids.Table 3 shows the selected catchments' geospatial properties.The size of these catchments ranges from 0.2 to 1.6 km 2 .Three catchments are from urban catchments that have streets with grid layouts, and eight catchments are from suburban neighborhoods that have more cul-de-sacs.The calculated parameters of all 20 selected catchments are included in Table S1 in Supporting Information S1.All the networks except for one are found to  Following these realistic ranges of watershed areas and H(s) values, the areas of our synthetic watersheds are all 0.81 km 2 , and our subsequent analysis on GI placements (Experiments 2 and 3) will compare networks with low and high H(s).We decide to use H(s) ≈ 0.02 (5-10th percentile, as in Figure S3 in Supporting Information S1) and H(s) ≈ 0.2 (around 55th percentile), because they realistically represent stormwater networks' flow path meandering, and they are sufficiently distinct (i.e., a magnitude of 10 apart).

Increasing Network Flow Path Meandering Decreases Peak Flow
When considering the impact of stormwater network flow paths, tradeoff exists between managing outlet peak flow rate and flooding loss.As shown in Figure 4, increasing network meandering H(s) leads to lower peak flow at the expense of higher flooding loss.We use two functional variables-transport efficiency and network storage capacity-to illustrate the effect of network structure, wherein high transport efficiency or high network storage capacity reduces flooding but simultaneously increases peak flow.However, increasing the design parameter H(s) will simultaneously increase network storage capacity while reducing its transport efficiency, and the relative impact of H(s) on the two functional variables will determine its overall impact on flooding and peak flow.Because we know from Experiment 1 that increasing H(s) reduces peak flow while increasing flooding (Figure 4), we surmise that across all rainfall intensities, H(s)'s impact on transport efficiency dominates over its impact on network storage capacity.
The land cover characteristics of the watersheds largely determine gray infrastructure's resilience to intensifying storms, which are represented by longer storm return periods (Figure 4 columns).Overall, both peak flow and flooding in suburban networks are smaller in magnitude compared to those in the urban networks.While all networks are designed for the peak flow of the same 10-year 24-hr design storm, the flooding in suburban networks is also  less sensitive to higher rainfall intensity.Because watersheds with larger impervious surfaces are less resilient to the increase in rainfall intensity, they may benefit more from the installation of GI.Thus, we focus our subsequent analysis on networks in urban watersheds with 80% impervious cover.

Green Infrastructure Coverage: Marginal Peak Flow Reduction Decreases With Storm Intensity
As rainfall intensity increases, higher GI coverage continues to decrease flooding but fails to control peak flow (Figure 5).This suggests that rainfall intensity moderates the marginal hydrological benefits achieved by increasing GI coverage.As explained in Section 2.5, adding GI reduces the transport efficiency (due to more storage and longer flow paths) but increases the network storage capacity (due to more storage).However, rainfall intensities affect how the functional variables influence the hydrological outcomes.
In a 2-year storm, the increase in network storage capacity and the reduction in transport efficiency work together to achieve both reductions in peak flow and flooding.The increase in network storage capacity will sufficiently The tradeoff between peak flow reduction and flood reduction can be understood as H(s)'s effect on transport efficiency reduction dominating over its effect on increasing network storage capacity.
contain the additional runoff generated from a low-intensity storm and reduce flooding.Under this condition, the reduced transport efficiency is able to simultaneously increase the lag time to the outlet, thus reducing peak flow.However, as rainfall intensifies, the increase in storage capacity from additional GI can no longer completely eliminate flooding.In such an overwhelmed network, the reduction in transport efficiency cannot sufficiently reduce peak flow.At this point, the networks with longer flow paths lead to lower peak flow rates, as previously seen in Figure 4. Therefore, at high rainfall intensity, the network structure, represented by network flow path meandering H(s), will control the peak flow (Figure 5) more than the number of GI nodes.

Green Infrastructure Placement: Rainfall Intensity Moderates the Effect of Downstream Green Infrastructure Placement
The effectiveness of GI placement strategy also depends on the rainfall intensity.Figure 6 shows that the placement of GI has negligible flood control benefits.Similarly, the placement also has minimal impact on peak flow reduction.However, it is worth noting that near-outlet placement results in opposite peak flow outcomes in different rainfall intensities: while placing GI downstream reduces peak flow in a 2-year storm, the same placement leads to a wider range of peak flow response, with a higher median value, as rainfall intensity increases.
During low-intensity rainfall, near-outlet placement of GI increases the runoff travel time in the downstream pipe segments, thus reducing the transport efficiency.Meanwhile, upstream placement effectively increases useable network storage capacity by reducing upstream runoff, with negligible effect on transport efficiency.Therefore, placing GI near the outlet during low intensity rainfall results in lower peak flow due to the dominant effects of reduced transport efficiency.
However, during high intensity rainfall, the main mechanism to control peak flow shifts from GI decreasing transport efficiency to gray infrastructure's flow path meandering as the networks become overwhelmed.Placing GI upstream in an overwhelmed network reduces the total upstream runoff, but the reduction is not enough to have an effect on the flow transported downstream.Thus, although downstream placement of GI lowers peak flow during low intensity rainfall, doing so during high intensity rainfall does not have an effect on peak flow control.

Discussion
Scaling the site-level effects of GI to the watershed is a key knowledge gap in stormwater planning.To investigate the impacts of gray and GI on our desired hydrological outcomes at the watershed-scale, our study sets up three experiments that investigate the effects of gray infrastructure network structure and of the interaction between green and gray infrastructure, narrowly represented by bioretention cells and stormwater pipe networks.
In Experiment 1, a baseline scenario without the presence of GI, we find that the degree of imperviousness in the associated drainage watershed land cover predominantly controls the hydrological outcomes, and tradeoff exists between managing peak flow and flooding when varying stormwater network flow path meandering.In Experiments 2 and 3, we find that although GI retains its site-level hydrological benefits when scaled up to the watershed level under light rainfall (e.g., a 2-year storm), its peak flow and flood control benefits often do not increase monotonically or linearly with increasing rainfall intensity.These findings are consistent with other field and modeling studies demonstrating the dependence of GI's effectiveness on rainfall intensity (Alves et al., 2019;Hathaway et al., 2014;Liu et al., 2014;McPhillips et al., 2021;Versini et al., 2018;Webber et al., 2020).
To better understand the complex interactions at play, we used two intermediate functional variables-transport efficiency and network storage capacity-to represent the mechanisms through which green and gray infrastructure may influence hydrological outcomes.The relationship between the design parameters and the functional variables are summarized in Figure 7. Increase in any of the functional variables will increase peak flow and reduce flooding, with an exception that the increase in above-ground network storage capacity due to GI will not increase peak flow.Rainfall intensity shifts GI's relative impact on the intermediate functional variables, resulting in rainfall-dependent hydrological outcomes.
Although the literature agrees that rainfall characteristics influence GI performance, with GI becoming less effective at high-volume, high-intensity rainfall, the specific characteristics of interest and other design parameters have led to a broad range of conclusions.In our study, even with a low-storage GI type such as a bioretention cell, we are still able to find that 40% of GI coverage, as in 40 of all 100 nodes in our study networks are assigned to be bioretention cell, can eliminate flooding in a 2-year storm (Figure 5).But the exact thresholds vary across studies.Moore et al. (2016)'s study concludes that biofiltration practices implemented at 10% of the watershed can moderate but cannot completely eliminate flooding for a 10-year storm.When considering the effect of GI placement, our study finds that placing a cluster of GI downstream can lower peak flow only in low-intensity rainfall (Figure 6), which is also consistent with results from other studies (Chang et al., 2008;Fiori & Volpi, 2020;Giacomoni & Joseph, 2017;Vittorio & Ahiablame, 2015).
Other investigations into GI placement strategies based on different assumptions nevertheless reached different conclusions.Whereas we define our parameter of GI coverage to be across a network, Zellner et al. ( 2016) designated 5%, 10% and 15% of total infrastructure coverage to GI within a subcatchment before the runoff drains to stormwater network.In that setup, they found clustered placement of GI downstream within a subcatchment even at 5% of infrastructure coverage is more effective than upstream cluster, but not as effective as random dispersed placement.Tansar et al. (2022) found that the effects of GI placement on peak flow reduction only become significant at high coverage (>50%) and only in a three-to-five year storm.But even so, rainfall intensity mitigates the effectiveness of different coverage and placement strategies.For example, an 80% uniform placement performs better than clustered placements in a three-to-five-year storm, but underperforms in storms with return periods greater than 50 years.Our results along with other studies all conclude that GI's performance and placement strategy changes based on rainfall intensity.

Integrated Planning of Gray and Green Infrastructure
At a watershed-scale, rainfall intensity controls the effectiveness of GI's coverage and placement by shifting GI's relative control on transport efficiency.At low rainfall intensity, GI can effectively lower peak flow rates by reducing transport efficiency, and it works together with gray infrastructure to lower peak flow rate at the outlet.
Although increasing GI density results in a wider range of peak flow response as rainfall increases (Figure 5, 100-year storm), it still contributes as additional network storage capacity for flood reduction.
However, at high rainfall intensity, the meandering of gray infrastructure flow paths and the land cover of associated drainage area substantially determine the peak flow and flood control outcomes.This finding is consistent with other studies that consider the effect of underlying watershed characteristics on GI performance, and how rainfall facilitates the transition from green-gray infrastructure working together to gray infrastructure dominating the hydrological outcomes.One study finds in large storms, the basin-scale imperviousness largely attenuates GI's (i.e., green roof) impact on peak flow (Versini et al., 2015).Others have also concluded that the effectiveness of GI depends on the degree to which the diverted overland runoff can actually be infiltrated and stored in the pervious areas (Lim & Welty, 2017;Miles & Band, 2015).Coined "watershed capacitance," the input-storage-release mechanism in a watershed depends on the land cover characteristics of the watershed (e.g., land cover types, antecedent soil moisture) and the rainfall characteristics (e.g., rainfall intensity and duration).
The management of hydrological outcome transitions from green-gray together to gray infrastructure only as rainfall intensifies.This rainfall-dependent transition suggests that the common practice to separately plan for gray and GI may not realize the full range of benefits from GI.In reality, stormwater infrastructure exists on a continuum of gray to GI (Bell et al., 2018), and it is important to consider the integration of all stormwater infrastructure for watershed level stormwater management.Many existing studies have already a gray-green infrastructure integration in face of extreme events (Möller & Hansson, 2008;Y. Kim et al., 2017;Denjean et al., 2017), but without explicitly considering the effects of underlying stormwater network structure.Our study provides a scenario-exploration tool that allows researchers and planners to scale up the effectiveness of GI, while considering varying rainfall scenarios, flow path structures of the associated stormwater networks, and land cover types of the associated watersheds.By quantifying the variability in the stormwater network structure (H(s)), this exploratory model unveils the underlying mechanisms through which green and gray infrastructure may influence hydrological outcomes, and it allows us to upscale the gray-green infrastructure's relationship to a watershed-level without being constrained by any specific watershed.Although the spatial optimization and the type of GI is limited in practice, we hope our model results shed light on future integrated gray-green infrastructure planning.

Further Scaling Challenges and Future Directions
Our study scales up the effect of GI by offering insights to the rainfall-dependent spatial relationship between GI and its associated gray infrastructure network.However, our model only investigates one type of GI with fixed model parameters.A typical municipal stormwater control project mixes and matches BMPs to achieve more diverse remediation goals, but it will result in a wider range of hydrological outcomes and make the watershed-level scaling even more challenging (Asleson et al., 2009;Clary et al., 2020).To account for the intrinsic heterogeneity of BMPs into our models, we can incorporate sensitivity analysis of the contributing hydrological parameters (Avellaneda et al., 2017;Worthen et al., 2021;Zhang & Chui, 2020) into our spatially-integrated gray-green infrastructure model, and investigate the placement strategy of different types of BMPs.Based on the management tradeoff we have already observed with stormwater network structure, we hypothesize that the additional storage from BMPs may help us address the tradeoff and simultaneously lower peak flow rates and flooding in the gray infrastructure network.
Our model may be extended from a 10-by-10 small lattice graph-based networks to larger stormwater networks.
Although increasing the grid size of the stochastically generated networks will significantly increase the computation time, tools from machine learning (e.g., Deep Convolutional Generative Adversarial Networks) can be applied to generate large stormwater networks more efficiently (S.E. Kim et al., 2021).We hypothesize that in larger networks, we will continue to observe rainfall intensity moderating gray-green infrastructure's effects on hydrological outcomes, which may have a wider range than in our 100-node networks.

Figure 1 .
Figure 1.(a) A summary of experiments with three design parameters (in circles), two functional variables (in transparent squares) and hydrological outcomes of interest (in shaded squares).Experiment 1 focuses on the effects of flow path meandering H(s), which serves as a proxy for the flow path structure of the stormwater networks.Experiment 2 focuses on the combined effects of flow path meandering and green infrastructure (GI) coverage.Experiment 3 focuses on the combined effects of flow path meandering and green infrastructure (GI) placement.For each experiment, we generate green-gray infrastructure networks based on the design parameters, and run these networks in the Storm Water Management Model (EPA-SWMM) under different rainfall scenarios to get the hydrological outcomes.For each experiment, we examine how the design parameters influence the functional variables, and in turn affect the hydrological outcomes.(b) A workflow diagram of the simulation scripts (Chen & Feng, 2023) used for study simulations.The script Gibbs.py is used to generate networks with different degrees of flow path meandering.The script hydro_network.py is for assigning networks hydraulic attributes and GI attributes.The script make_SWMM_input.py is for running SWMM and recording results for further analysis.

Figure 2 .
Figure 2. Left to right: connected lattice graph G from which spanning trees are generated, spanning tree s with lower network flow path variation H T (s)(Troutman & Karlinger, 1992) and network flow path meandering H(s), and spanning tree t with higher H T (t) and H(t).In the two spanning trees, the nodes represent access points where runoff enters the stormwater network, the larger node at the bottom left corner is the outlet node, the arrows on the edges designate the flow direction, and the thickened edges highlight the flow path from an arbitrary node to the outlet node.The arbitrary node is at the same location on the lattice graph G.The shortest path from such arbitrary node to the outlet is highlighted in graph G. Spanning tree s follows shorter paths and has lower H T (s) and H(s), as the highlighted flow path follows the shortest path.Spanning tree t has higher H T (t) and H(t), as, for example, the highlighted path is much longer than the shortest path.
s) under 0.4, suggesting the actual flow paths of the networks are close to the shortest paths and are rather efficient paths.The network with the highest meandering (H(s) = 1.648) is in an area with many cul-de-sacs.

Figure 3 .
Figure 3. Selected catchments from the Twin Cities Metropolitan Area stormwater municipal networks were reconfigured onto grid lattice grids.H(s) values of these catchments were calculated.Two of the selected catchments (3 and 7) are highlighted to show the process of lattice grid reconstruction from the original shapefiles (left) to the simplified grid-based graphs (right), where the red nodes are the network outlets.

Figure 4 .
Figure 4. Effect of network structure: peak flow and inland flooding in networks with varying degrees of network path meandering H(s), under 2, 10, 25, and 100-year 2-hr storms.Increasing H(s) lowers peak flow rate, but it does so at the expense of increasing inland flooding.The increase in inland flooding due to network structure is more obvious in high-intensity storm and in watersheds that have higher proportion of impervious drainage area.Higher flow path meandering H(s), the design parameter, increases (solid line) network storage capacity while decreasing (dashed line) transport efficiency.The tradeoff between peak flow reduction and flood reduction can be understood as H(s)'s effect on transport efficiency reduction dominating over its effect on increasing network storage capacity.

Figure 5 .
Figure 5. Interaction between green infrastructure (GI) and network structure: peak flow and inland flooding in networks with two fixed networks-H(s) ≈ 0.02 and H(s) ≈ 0.2-and varying percentages of GI nodes, under 2, 10, 25, and 100-year 2-hr storms.Adding GI increases (solid line) the network storage capacity.In low-intensity rainfall, it reduces (dashed line) the transport efficiency (due to more storage and longer flow paths).

Figure 6 .
Figure 6.Effect of green infrastructure (GI) placement shown as the distance to the network outlet: peak flow and inland flooding in multiple networks with H(s) ≈ 0.02 and H(s) ≈ 0.2 and varying placement of GI nodes, under 2, 10, 25, and 100-year 2-hr storms.In low-intensity rainfall, placing GI downstream reduces (dashed line) the transport efficiency.

Figure 7 .
Figure 7. Conceptual diagram that demonstrates the experiments conducted in our study, and the relationship among gray infrastructure and green infrastructure's design parameters, functional variables and hydrological outcomes.

Table 1
Network Design Attributes(SWMM Input)

Table 2
Rainfall Intensities of a 2-hr Storm CHEN ET AL.