Design and Simulation of Stormwater Control Measures Using Automated Modeling

Stormwater control measures (SCMs) are decentralized technical elements, which can prevent the negative effects of uncontrolled stormwater flow while providing co-benefits. Optimal SCMs have to be selected and designed to achieve the desired hydrological response of an urban catchment. In this study, automated modeling and domain-specific knowledge in the fields of modeling rainfall-runoff (RR) and SCMs are applied to automate the process of optimal SCM design. A new knowledge library for modeling RR and SCMs, compliant with the equation discovery tool ProBMoT (Process-Based Modeling Tool), was developed. The proposed approach was used to (a) find the optimal RR model that best fits the available pipe flow measurements, and (b) to find the optimal SCMs design that best fits the target catchment outflow. The approach was applied to an urban catchment in the city of Ljubljana, Slovenia. First, nine RR models were created that generally had »very good« performance according to the Nash–Sutcliffe efficiency criteria. Second, six SCM scenarios (i.e., detention pond, storage tank, bio-retention cell, infiltration trench, rain garden, and green roof) were automatically designed and simulated, enabling the assessment of their ability to achieve the target outflow. The proposed approach enables the effective automation of two complex calibration tasks in the field of urban drainage.


Introduction
Stormwater control measures (SCMs), also known as Sustainable urban drainage systems (SUDS), Low impact development (LID), Best management practices (BMP), Water sensitive urban design (WSUD), and Nature-based solutions (NBS) [1,2], are technical elements that are designed to prevent and mitigate negative effects of uncontrolled stormwater flow (e.g., urban floods or excessive combined sewer overflows) [3,4]. Common to all of them is a decentralized (or on-site) approach to rainwater management by using processes typical for the natural water cycle (e.g., infiltration, retention . . . ) while providing ecosystem services and additional benefits for the city. Many studies have investigated the influence of individual SCMs on the quantity and quality of stormwater in urban sub-catchments, i.e., on a micro-scale [5][6][7][8]. Moreover, models are being used to design SCMs and predict SCM performance for different design and weather scenarios [9][10][11]. To investigate the impact of SCMs at meso-scale and macro-scale these models need to be integrated into or coupled with catchments rainfall-runoff (RR) models. Models incorporate the previously acquired knowledge on water-related physical phenomena (e.g., water percolation [12,13], surface runoff [14] . . . ) in a form of mathematical formulations.
In this context, the open-source Stormwater Management Model (SWMM) [15] is one of the most widely used models [16][17][18][19] enabling RR and hydraulic modeling of 1.
To develop a library of components for modeling RR and SCMs, compliant with the equation discovery tool ProBMoT (Process-Based Modeling Tool) [32]; 2.
To apply the proposed automated model discovery approach to find the optimal RR model for an urban sub-catchment in the city of Ljubljana, Slovenia, which best fits the available pipe flow measurements; 3.
To find an optimal design of SCMs that would best fit the target catchment outflow.

Case Study Area
The case study area is located in the western part of the city of Ljubljana, Slovenia, and covers about 30 ha. The predominant land use in this area is family houses with gardens ( Figure 1). The area has a temperate continental climate, with a mean long-term  annual rainfall of about 1380 mm [33]. The largest part of the area is served by a mixed sewer system with a length of approx. 5.4 km.

Data
Precipitation data were provided by the Slovenian Forestry Institute, for a rain gauge station located on the north side of the case study area (location: 46 • 03 06.82 N, 14 • 28 47.58 E, 306 m above mean sea level (MAMSL) [34]. The measurements were performed with a Davis ® (0.2 mm) Rain Gauge Smart Sensor [35]. Data on extreme rainfalls were obtained from the IDF curves for the Ljubljana-Bežigrad weather station [36]. These served for the simulation of design rain events with a duration of 1 h, and return periods of 5, 10, 25, and 50 years, with total precipitation depth of 38 mm, 44 mm, 53 mm, and 59 mm, respectively.
The local public utility company (JP VODOVOD KANALIZACIJA SNAGA d.o.o.) provided the information on the combined sewer network and the flow measurement data. The flow measurements were performed between 1 March 2019 and 30 September 2020 in a pipe (of the combined sewer system) collecting all contributing water from the case study area (location: 46 • 02 44.56 N, 14 • 29 01.82 E, 295 MAMSL) ( Figure 1). The flow rate was calculated based on the combination of non-contact radar velocity measurements and ultrasonic water level measurements [37]. As the focus was only on rainfall-runoff (RR) (i.e., stormwater) modeling and the measurements were conducted in a combined sewer system, the dry weather flow, which is independent of the RR and cannot be explained with the RR model, was deduced from the total measured flow. For the preparation of precipitation and flow data, a 5-min time step was used.
Within the entire flow measurement period, rainfall data were collected corresponding to the measured flow. To test the usefulness of the proposed approach for single events and continuous simulations, as defined by the EPA SWMM [15], two types of rainfall periods were selected: shorter (approximately 1 day) and longer (approximately 6 days). Four rainfall periods were selected for model calibration (C1-C4) and another four for model validation (V1-V4) ( Table 1). Hence, two short (C1 and C2) and two long (C3 and C4) rainfall periods were selected for model calibration. The longer calibration periods include multiple rainfall events of different intensities and durations so that the ProBMoT tool could learn from the different surface runoff responses of the catchment. The validation periods were selected based on the following criteria: (a) duration, with shorter (V1 and V2) and longer (V3 and V4) periods; (b) seasonal distribution; (c) similarity in peak flows; (d) stability of measured flows; and (e) antecedent weather conditions [38,39].

Rainfall-Runoff Model
The rainfall-runoff (RR) model is based on the principles and equations used by the EPA Storm Water Management Model (SWMM) [15]. Variables and constants included in Equations (1)-(8) are described in Tables 2 and 3. The catchment is represented as a nonlinear reservoir, governed by surface storage mass balance, i.e., conservation of mass: Evaporation was not included in the model, due to its limited influence on the water mass balance within the short rainfall events used for SCM design. The runoff flow rate per unit of the surface area is based on Manning's equation [15]: Afterward, the total runoff flow from the catchment (L/s) is calculated by multiplying the catchment surface area (m 2 ) and q (m/s).
There are several well-known alternative methods for modeling the infiltration process. Three alternatives were considered in this study: the Soil Conservation Service Curve Number (SCS CN) method, the Variable UK runoff equation, and the UK Water Industry Research runoff equation. Thus, the complete RR model can have different structures, depending on the infiltration method selected. In addition, different infiltration methods can be applied in different sub-catchments, resulting in many plausible model structures for a given catchment.  This method assumes that the total infiltration capacity of a soil is related to the soil's tabulated Curve Number (CN). The CN value is determined by the hydrologic soil group, land use, and hydrologic condition. CN values range from 30 to 98. The latter value is assigned to paved roadways, roofs, and other impervious surfaces. At higher CN values (i.e., closer to 98), precipitation is mainly translated into surface runoff. On the other hand, at lower CN values (i.e., closer to 30), rainfall is mainly infiltrated and is thus not translated into a runoff [14]. A modified version of the SCS CN equation was used, as described in the SWMM Reference Manual-Hydrology [15]. In the modified version, the initial abstraction (I a ) is not included, as it is already included in the depression storage (d s ; see Equation (2). The following three equations (Equations (3)-(5) are used to calculate the infiltration rate:

Variable UK Runoff Equation
The Variable UK runoff equation (VARUK) has three components: runoff from impervious areas, runoff from pervious areas, and initial losses [40]. It is based on data from 11 UK catchments and 112 rain events. The VARUK equation is as follows [41]: Rainfall can be converted into infiltration by using the following equation:

UK Water Industry Research Runoff Equation
The UK Water Industry Research runoff equation (UKWIR) [42] was developed to overcome some of the limitations of VARUK, which are described in detail in the report Development of the UKWIR Runoff Model (2014). As VARUK, it has a fixed runoff component for paved surfaces (IF n × PIMP n ). It was upgraded with a variable runoff component for paved surfaces (1 − IF n ) × PIMP n . Additionally, a component for pervious surfaces was added (1 − PIMP TOTAL ), which enables differentiation between winter and summer runoff (i.e., negative antecedent precipitation index (NAPI)). The UKWIR equation is then as follows:

Stormwater Control Measures
In this study, SCMs are based on the principles and equations used by SWMM for modeling low impact development (LID) controls [15]. These are bio-retention cell (BRC), rain garden (RG), green roof (GR), and infiltration trench (IF). Based on the principles and equations typical for these measures, two additional SCMs were defined, namely the detention pond (DP) and the storage tank (ST). The processes that characterize individual SCMs are described in equations 9-30. Variables and constants that are included in these equations are described in Tables 4 and 5.
To make the computation of SCMs hydrologic performance less complex, the following simplifications are assumed by SWMM [15]: (1) the cross-sectional area of the unit remains constant throughout its depth; (2) flow through the unit is 1D in the vertical direction; (3) inflow to the unit is distributed uniformly over the top surface; (4) moisture content is uniformly distributed throughout the soil layer; (5) matric forces within the storage layer are negligible so that it acts as a simple reservoir that stores water from the bottom up.

Bio-Retention Cell
Bio-retention cells (BRCs) are depressions that contain vegetation grown in an engineered soil mixture placed above a gravel storage bed. They provide storage, infiltration, and evapotranspiration of both direct rainfall and runoff captured from surrounding areas. As such, they can be used as a generic SCM model, which can be then customized to describe the behavior of other SCM types. Due to the short length of the design rain events (i.e., 1 h); evapotranspiration was not included in the model, as it would not have a significant influence on the results. A conceptual model of a typical BRC can be represented with three overlaying horizontal layers ( Figure 2). The inflow to the surface layer includes both direct precipitation (i) and captured runoff from contributing areas (q 0 ). On the other hand, the outflow can occur through infiltration into the soil layer below (f 1 ) or by the overflow of ponded surface water (q 1 ). The latter may occur if the level of ponded surface water exceeds the freeboard height of the surface layer (D 1 ). The soil layer receives infiltration (f 1 ) from the surface layer and loses water through percolation (f 2 ) into the storage layer below. The storage layer receives percolation (f 2 ) from the soil layer above and loses water by exfiltration (f 3 ) into the underlying natural soil. Materials with high porosity (e.g., gravel) are used to provide water storage capacity. Variables and constants included in the SCM equations are described in Tables 4 and 5. Along with their original (EPA-SWMM) [15] names, the newly defined names used in ProBMoT are provided. Furthermore, Table 5 provides the minimum and maximum values of the constants.
The current amount of water within every layer can be calculated by solving mass balance equations that take into account all the inflows and outflows (Equation (9)-(11)). The infiltration rate (f 1 ) (Equation (12)) is based on the Green-Ampt equation, whereas for the percolation rate (f 2 ) (Equation (13)), Darcy's law is used [15]. The overflow rate (q 1 ) (Equation (14)) is simply calculated as any ponded water that exceeds the maximum freeboard (or depression storage) height within each time step.
Water fluxes between layers are limited by the current conditions in the connected layers. Two types of conditions need to be fulfilled for infiltration to occur: (1) a sufficient amount of water is needed in the upper layer (Equations (16) and (18)), (2) a sufficient amount of empty pore space is needed in the lower layer (Equations (15) and (17)).

Rain Garden
SWMM defines a rain garden (RG) as a bio-retention cell without a storage layer [15]. Its governing equations are, therefore, Equations (9) and (10). The nominal soil percolation rate f 2 (Equation (13)) is limited to the smaller of the following two values: the amount of drainable water available in the soil layer (Equation (16)) and the saturated hydraulic conductivity of the native soil (K 3S ).

Green Roof
SWMM's green roof (GR) is also similar to a bio-retention cell, except it uses a drainage mat instead of gravel aggregate in its storage layer and has an impermeable bottom (K 3S = 0). Therefore, the bottom exfiltration (f 3 ) is replaced by the drainage mat flow rate (q 3 ) (Equation (23)). Furthermore, there is no term for captured runoff from impervious areas (q 0 ) (Equation (19)). Both the runoff rate from the soil layer surface (Equation (22)) and the drainage mat flow rate (q 3 ) (Equation (23)) are computed by using the Manning equation for uniform overland flow.

Infiltration Trench
SWMM's infiltration trench (IT) is also similar to a bio-retention cell, except that it has only a surface and a storage layer, where f 1 assumes that any ponded surface water (d 1 ) drains into the storage layer over the time step. Furthermore, the surface void fraction (φ 1 ) does not appear in the surface layer equations (Equation (24)) as a gravel-filled trench would have no vegetation growing above it.
Both the infiltration and exfiltration rate are limited by the amount of water currently in the storage layer (Equations (27) and (28)):

Detention Pond and Storage Tank
Based on the principles and equations presented with the previously defined elements, two additional elements were defined, namely a detention pond (DP) (Equation (29)) and a storage tank (ST) (Equation (30)). They both include only a surface layer, which is deeper than with the previously defined elements. The inflow for DP includes direct rainfall and runoff from contributing areas, whereas ST receives only runoff from contributing impervious areas. For both elements, the outflow (i.e., overflow) is calculated by using q 1 (Equation (14)).

Equation Discovery and Process-Based Modeling
The proposed automated modeling approach is based on the Process-Based Modeling Tool (ProBMoT), developed byČerepnalkoski et al. [43]. ProBMoT allows for the integration of domain knowledge, formalized as template components for the construction of process-based models, into the procedure of equation discovery from measured data. It automatically identifies both the structure and parameter values of an appropriate process-based model, given: (a) a knowledge library (i.e., a mathematical formulation of the selected domain) in the form of model components, or, more specifically, template entities and processes, (b) a conceptual model of the observed system, and (c) measurements ( Figure 3). Candidate model structures are generated from the knowledge library and a user-specified conceptual model of the observed system. The candidate models are transformed into equations, calibrated against provided data, and ranked according to their errors. The latter is calculated as the root-mean-squared-error (RMSE), i.e., the discrepancy between the model simulation and measured data. To use ProBMoT within the presented case study of RR modeling and SCM design, the following steps were taken: (A.1) The RR models and SCM processes were encoded in a modeling library; (A.2) Conceptual models of the case study and SCM scenarios were elaborated; (A.3) ProBMoT was used to discover the optimal model structure and parameters among viable RR models, following the conceptual model of the case study and flow measurements; (B.1) The RR model with the best performance was used to simulate catchment outflow for rainfall events with a duration of 1 h and return periods of 5, 10, 25, and 50 years; (B.2) Three design events (i.e., data sets) that represent target catchment outflows were defined. The target outflow of precipitation (with a return period of 10, 25, or 50 years) was set to be the outflow typical for a 5-year return period; (C.1) ProBMoT was used to determine the optimal SCM parameters (i.e., SCM design), following the conceptual model for each SCM scenario and target catchment outflow (i.e., outflow reduction); (C.2) To determine the best SCM design, the conceptual models for SCMs were iteratively changed, based on the preliminary results of SCM design.

Library of Components for Modeling Rainfall-Runoff and Stormwater Control Measures
The library consists of entity templates, process templates, and compartment templates. Each template captures general knowledge that can be applied to different cases and can be reused when dealing with a specific task. The dynamic system to be modeled, i.e., the catchment, can be structured by using compartments. Compartments are organized in a nested, tree-like structure. Each compartment contains entities and processes and can also contain other sub-compartments (e.g., sub-catchments) [29]. In the urban hydrology domain, processes calculate the change of water fluxes (e.g., surface storage) within a time step and entities aggregate these changes over the simulated time.
The equations presented in Section 2.3 were encoded in the knowledge library (see Supplementary Material S1: Knowledge library) as template processes named general processes, as they are common to all the models. These include D (surface storage), Q (surface runoff), OF1 and OF3 (outflow), and F (infiltration), which can be modeled with three alternative methods: F_SCS, F_UKWIR, and F_VARUK. The SCMs equations presented in Section 2.4 were encoded in the knowledge library as template processes named after SCMs (e.g., detention pond processes). Names of processes and entities are derived from variable names and SCM type abbreviations (e.g., D1_DP). Uppercase letters are used for variable names related to processes (e.g., D1_DP) and lowercase letters for variable names related to entities (d1_DP). Outflow processes that sum the discharges from different sub-compartments (i.e., sub1, sub2, etc.) were introduced. Namely, the OF3 (i.e., the sum of discharges from an impervious and pervious area with no implemented SCMs (i.e., OF1)), the OF4_SCM (i.e., the sum of the overflow from SCMs (i.e., OF2_SCM)), and the OF5_SCM (i.e., the sum of OF3 and OF5_SCM). The time step in the presented equations is 1 s, so ProBMoT is calculating the infiltration rate in m/s and discharge in L/s. However, the actual time step of the input data (e.g., precipitation, measured flow) is 5 min, which is also the reporting time step.
Furthermore, the parameters that appear in the RR model (Table 3) and SCMs (Table 5) were listed as constants within the template entities »Surface« and »SCM« in the knowledge library, respectively, together with their expected ranges and units.

Conceptual Models of the Case Study Area and Stormwater Control Measures
To apply ProBMoT to a specific catchment, a conceptual model of the observed system must be provided (see Supplementary Material S2: Conceptual model). The conceptual model consists of entity, process, and compartment instances. Each instance is created using a template from the urban hydrology modeling library. The specification of a process instance includes a list of arguments (i.e., names of the entity instances that are involved in the process) and a list of constants with their exact values (e.g., entity: Surface, constant: slope = 0.002). For unknown values of the constant parameters or the parameter that will be calibrated (e.g., SCM design), the parameters are assigned a special value »null«, along with its expected range of values (e.g., entity: SCM, constant: D1 {value: null; fit_range: <0, 1.5>}). The values and ranges of the constants that appear within the RR model (i.e., entity Surface) were adjusted based on the (im)perviousness of the sub-catchment (Table 3) following the values proposed in the literature [14,15,[40][41][42]. The values and ranges of the constants that appear within SCMs (i.e., entity SCM) were defined based on the values proposed by EPA [15] ( Table 5).
The conceptual model of the case study area (Figure 4a) was structured as a single compartment (i.e., catchment), divided into two sub-compartments (i.e., sub1 and sub2), where sub1 represents the impervious part of the catchment (60% of the total surface area) and sub2 the pervious part of the catchment (40% of the total surface area). A new sub-compartment (i.e., sub3) (Figure 4c) was introduced to account for SCMs, including SCM area and contributing areas. SCMs can be placed either within a part of an existing pervious area (DP, IT, RG, BRC) or an impervious area (ST, GR). However, in the case of GR, the measure itself also acts as a contributing area, by receiving direct rainfall. To ensure transparency, a separate conceptual model was defined for each SCM.

Rainfall-Runoff Models
Given the conceptual model and the modeling knowledge library, ProBMoT explored nine alternative structures of RR models (M1-M9) (Appendix A Table A1). Each model chooses a different combination of modeling templates for infiltration for each of the two sub-catchments. Each of the nine models was calibrated and validated against the measured data. All models were calibrated by simultaneously using data (C1-C4) from two short time periods of measured flow, with a duration of approximately one day (C1 and C2), and two long periods of measured flow, with a duration of approximately one week (C3 and C4). The calibrated models were validated on four different rainfall periods (V1-V4) that vary in rainfall duration and intensity. The values of the Nash-Sutcliffe efficiency (NSE) coefficient [44] for the calibration and validation data sets/events are presented in Table 6. The rainfall periods are presented in consecutive order, based on the modeling phase (first calibration-C, then validation-V) and the duration of the modeled rainfall period (e.g., 1D-one day).
Based on the criteria for assessing the goodness of fit for models, proposed by Moriasi et al. [45], the models on average had »very good« performance (i.e., 0.75 < NSE ≤ 1.00) for all calibration periods. In the validation process, all models had »very good« performance for all validation periods, confirming the usefulness and efficiency of the proposed approach for calibration of model parameters. When comparing the performances of the models, differences can be noticed for individual rainfall periods. Nevertheless, they all achieved a similar average NSE value for calibration (i.e., average NSE value of 0.794) and validation periods (i.e., average NSE value of 0.859). Therefore, no general conclusion can be made about which model is better, since model performance results are greatly influenced by the characteristics of calibration and validation periods (e.g., rainfall intensity, number of (sub) rainfall events) [46]. However, the model with the highest average NSE values for the calibration and validation periods (i.e., M3) was used for further analysis. Next, the 1-h rainfall events with return periods of 5, 10, 25, and 50 years were simulated using the calibrated model M3, resulting in a total outflow of 6429 m 3 , 7562 m 3 , 9426 m 3 , and 10,403 m 3 , respectively. Based on these results, three target catchment outflow scenarios and datasets were defined, namely the DE_Q5_P10, DE_Q5_P25, and DE_Q5_P50. These design events represent a combination of an outflow typical for 5-year precipitation (i.e., Q5) and precipitation with longer return periods (i.e., P10, P25, P50). In this way, SCMs are designed to reduce the peak and total volume of the catchment outflow.

Stormwater Control Measures
SCMs were designed to retain the difference between outflows from the catchment provided by precipitation with 10, 25, and 50 years return period (i.e., P10, P25, P50) and 5 years return period (Q5, or the target outflow). More precisely, they were designed to maximize the extent to which each SCM type can achieve the target outflow, using the proposed parameter ranges. The match between the simulated and the target outflow was evaluated using the Nash-Sutcliffe efficiency coefficient (NSE), peak flow ratio (PFR), and total volume ratio (TVR) (Appendix A Tables A2-A7).
In the first round of SCM design, the conceptual model with three sub-compartments ( Figure 4) was used. Furthermore, the surface area of all SCMs was set to 5000 m 2 , except for GR where this value was set to 50% of all the impervious area (i.e., 90,000 m 2 ). The results showed that, in many cases, an SCM with uniform characteristics could not sufficiently follow the event dynamics. Namely, DPs and STs provide detention for a limited amount of time after which the outflow from SCM appears (Figure 5a). Therefore, the initial SCM sub-compartment (i.e., sub3), was further divided into three equally sized subcompartments with SCMs (i.e., sub 3-5) (Figure 5b), enabling the design of three elements of the same SCM type with different characteristics (e.g., dimensions). Moreover, the initial simulation provided information on which SCMs could take up less space (i.e., shallow elements) (e.g., ST, DP, IF) or would need more space for its implementation (i.e., thick elements or significant overflow) (e.g., RG, BRC). Taking into account the above-mentioned observations, the following results were obtained.

Detention Pond and Storage Tank
Due to their similarity, detention ponds (DPs) and storage tanks (STs) provided comparable results. Namely, a near-perfect match between the target and simulated outflow was achieved for both measures (i.e., an NSE value of 0.99) (Appendix A Tables A2 and A3). In the second round of SCM design, the areas of ST and DP were reduced by 80% (i.e., to total an SCM area of approximately 1000 m 2 ). Thus, these are the two most space-efficient measures. As expected, larger and deeper elements were designed for less frequent design events and vice versa. In general, the measures can be divided into those that detain all the inflow from contributing areas, and those that provide outflow. For the latter, at least one of the dimensions was set to (or very close to) the minimum value (i.e., area to 250 m 2 or depth to 0.5 m).

Infiltration Trench
Infiltration trenches (ITs) also achieved a good fit between the target and simulated outflow (i.e., NSE value of 0.99) (Appendix A Table A4). However, when compared to ST and DP, these elements on average take 2.5 times more space (i.e., between 1997 and 3000 m 2 ). This is caused by a smaller maximum depth (i.e., 3.5 m), smaller maximum voids fraction (i.e., 40%) and limited infiltration rate of the existing soil (i.e., 2.6 × 10 −5 m/s). On the other hand, IT can be implemented as a set of elements with tailor-made shapes, based on site conditions (e.g., along a plot border). As expected, larger and deeper elements were designed for less frequent design events and vice versa. For the design event DE_Q5_P25 (Figure 6), two ITs were designed to capture the whole run-off from contributing areas with the maximum area of 1.000 m 2 and storage layer depths of 3.40 m and 2.62 m. The third IT (i.e., sub4) was designed using the min. dimensions (i.e., area 500 m 2 , 0 m of freeboard (i.e., D1) and a storage layer that is 0.9 m deep). Although the freeboard is set to 0 m, surface ponding can still occur. In this case, the maximum ponding depth (i.e., d1) is 0.25 m, highlighting the fact that IT with an even smaller area could be designed. However, this could increase the ponding height unreasonably.

Rain Garden and Bio-Retention Cell
RGs and BRCs yielded similar results (Appendix A Tables A5 and A6), as they are both limited by the soil infiltration rate (f1), which is calculated based on the soil's saturated hydraulic conductivity (K2S). Therefore, both measures performed well for DE_Q5_P10 (Figure 7a) but proved insufficient for DE_Q5_P50 (Figure 7b). Namely, although maximum parameter values, that would increase detention and infiltration, are selected for DE_Q5_P50, significant overflow occurs between the 45 and 75 the minute of simulation. This can be observed also in Figure 8a, where the overflow occurs when d1 exceeds the maximum freeboard (i.e., D1 is 0.30 m). As mentioned before, the infiltration rate is too low to transfer all the ponded water into the soil layer fast enough. Looking at Figure 8b, we can see that the soil layers' moisture content riches maximum values (i.e., 0.60) only after 2 h. Among the SCMs that are placed within existing pervious areas, these two SCM types take up the most space, between 3000 and 6000 m 2 . However, in contrast to other SCMs that do not include a soil layer, they can offer additional ecosystem services (i.e., water filtration, biodiversity, etc.). It can be concluded that taking into account the proposed parameter values, including the SCM area, RGs and BRCs are more suitable measures for more frequent rain events (i.e., DE_Q5_P10).

Green Roof
In comparison to other SCMs, the green roof (GR) acts as a contributing area and an SCM area at the same time. Namely, it does not receive surface runoff (q0) from adjacent areas. GRs are a replacement for conventional roofs that act as impervious areas (i.e., sub1). Therefore, changes in the GR area have a direct influence on the size of the impervious area within sub1. Namely, if the GR area is increased, the impervious area has to be decreased by the same amount. Simulations showed that GRs could sufficiently detain all direct rainfall. Thus, the optimal GR area depends on the outflow from impervious (i.e., sub1) and pervious areas (i.e., sub2), which is in direct correlation with their size. Figure 9a presents results for DE_Q5_P50, where GRs replace 50% of existing impervious areas (i.e., 90,000 m 2 ), resulting in an outflow from impervious and pervious areas (i.e., of3) that is lower than the target outflow (i.e., Q5). Consequently, GRs are designed in a way that they provide additional outflow (i.e., from sub4). As this design is suboptimal, the area of GRs was further reduced to 40% of existing impervious areas (i.e., 72,000 m 2 ) (Appendix A  Table A7). The results for this scenario show that only a small outflow from GRs occurs and that of3 almost matches the target outflow (i.e., Q5) (Figure 9b). For both design events of higher frequencies, the of3 was lower than the target outflow (i.e., Q5). Thus, the GRs were designed in a way to provide additional outflow.

Stormwater Control Measures Performance
In general, all of the SCM scenarios performed well, except for the scenarios 3RG and 3BRC, which exceeded the target peak flow, for the design events DE_Q5_P25 and DE_Q5_P50 (Table 7). As mentioned before, RGs and BRCs were limited by their design characteristics (e.g., soil's saturated hydraulic conductivity) and the proposed area for implementation (i.e., 5% of the pervious area). If the additional area was available, these two scenarios could also perform well. Scenarios 3ST, 3DP, and 3IT slightly exceed the target outflow values for the design event DE_Q5_P50; however, they remain within tolerable limits. In this way, an optimal SCM design is provided, by setting some of the characteristics of the elements to the maximum possible value (e.g., depth of the surface layer, area . . . ). The 3GR scenario, assuming that 40% of the current impervious areas will be replaced with green roofs, achieved the best performance. Table 7. The performance of stormwater control measure (SCM) scenarios for the analyzed design events, evaluated with the total volume ratio (TVR) and peak flow ratio (PFR).

Rainfall-Runoff Models
The development of a validated RR model is a prerequisite for the assessment of hydrologic effectiveness of different SCM scenarios [25]. In this study, an automated equation (model) discovery approach that enabled the automatic generation and calibration of RR models was applied. The derived RR models can be classified as lumped (regarding spatial resolution) and continuous (regarding temporal resolution) [47]. To build fully distributed models, hydraulic processes (i.e., pipe flow) should be included. Due to some limitations of the ProBMoT tool, only hydrological processes were included. Hydraulic (pipe flow) modeling is crucial when modeling complex and large sewer system. In our case, the case study area and the adjacent sewer system were relatively small, thus the hydraulic (pipe flow) modeling would not significantly affect the outflow dynamics.
The developed knowledge library, compliant with the ProBMoT tool, enabled a systematic comparison of three alternative infiltration methods (i.e., SCS CN, VARUK, and UKWIR) and their combinations. This is a unique and novel approach, as normally only one infiltration method is used within one RR model [20].
In general, all RR models were judged to have »very good« performance, with only minor differences across models. Thus, no general conclusion can be made on which model is better, since model performance results are greatly influenced by the characteristics of calibration and validation periods (e.g., rainfall intensity, number of (sub) rainfall events) [46]. However, two of the models that combine infiltration methods performed best (i.e., M3, M4), highlighting the true value and contribution of this research, which is the discovery of new knowledge in the urban runoff-modeling domain by automated modeling.
In contrast to traditional machine learning methods for knowledge and model discovery that are purely data-driven, ProBMoT allows for integrating expert background knowledge from the domain of use in the process of equation/model discovery. The knowledge library, presented in this paper, provides a blueprint for building mechanistic, transparent models of hydrological events as opposed to black-box models that would be obtained by traditional data-driven methods. The blueprint results in a small number of candidate model structures, which significantly reduces the danger of overfitting the model to calibration (i.e., training) data. This robustness to overfitting explains why the resulting models often perform better on validation data and confirms the model generality, i.e., its ability to properly explain unobserved hydrological events and accurately predict their outcomes and implications.
We should distinguish between two major types of users of the proposed approach. The first type of users are domain experts and experienced users, who have a thorough understanding of the underlying physics. These users are allowed to create and modify libraries of domain knowledge, which include templates for modeling different processes that are consistent with the physics of the domain. The second type of users are so-called end-users. They can use libraries of domain knowledge, conceptual models, and measured data to solve specific urban catchment modeling and SCM design problems. Assuming that they lack in-depth knowledge of the underlying physics, they would not be allowed to modify the libraries of domain knowledge and the template formulas for modeling different processes.

Stormwater Control Measures
In this study, SCMs were designed based on the target catchment outflow (i.e., hydrograph). This is typical for auto-calibration of SCMs, where an objective function has to be defined [25]. Otherwise, the performance of SCMs is usually evaluated only with hydrological parameters (e.g., peak flow reduction, total volume) after their predetermined design and model simulation [4,16,19]. The hydrograph-based SCM design can be extremely helpful in situations where not only the peak flow is a limitation, but also the flow distribution (i.e., downstream convergence of outflows from catchments with different dynamics). In this way, this approach can be used to design SCMs and investigate to what extent they can reduce flood hazards. Moreover, in this study, all relevant SCMs parameters were included in the calibration (i.e., design) process, including material characteristics, as opposed to a more typical approach, where only a search for the optimal SCMs dimensions (e.g., surface) is conducted. A comparative analysis with other studies of the obtained SCM performance results cannot be conducted, due to different ratios between the contributing and the SCM area, the difference in implementation surface among SCM types, and specific local conditions.

Conclusions
In this study, an automated equation (model) discovery approach was applied to the field of RR modeling and SCM design. First, a new library of model components, compliant with the equation discovery tool ProBMoT was developed, formalizing the knowledge on RR modeling and SCM design. Next, a conceptual model of the experimental urban sub-catchment within the city of Ljubljana, Slovenia, was defined. The proposed automated model discovery approach was used to find the optimal structure among the viable RR models, based on pipe flow measurements. The best RR model was used to simulate catchment outflow for rainfall events with a duration of 1 h and return periods of 5, 10, 25, and 50 years. The results were used to define three design events (i.e., data sets) that represent target catchment outflows. Namely, by using SCMs, the precipitation with a return period of 10, 25, or 50 years, should generate outflow typical for a 5-year return period. Next, ProBMoT was used to determine the optimal SCM parameters (i.e., SCM design), following the conceptual model for each SCM scenario and target catchment outflow (i.e., outflow reduction).
The main findings of the study are as follows: 1. The proposed automated model discovery approach for finding the optimal RR models proved to be very efficient. Nine RR models were considered that generally had »very good« performance. The best performance was achieved by two models that used a combination of two different infiltration methods; 2. The proposed automated model discovery approach enabled the design of six SCM scenarios and the assessment of their ability to achieve the target outflow; 3.
Detention ponds (DPs) and storage tanks (STs) provided comparable results. Namely, a near-perfect match between the target and simulated outflow was achieved for both measures (i.e., NSE value of 0.99). These are the two most space-efficient measures, with an average area of approx. 1000 m 2 ); 4.
Infiltration trenches (ITs) achieved a good fit between the target and simulated outflow (i.e., NSE value of 0.99). When compared to ST and DP, these elements on average take 2.5 times more space, due to smaller maximum depth, smaller maximum voids fraction and limited infiltration rate of the existing soil; 5.
RGs and BRCs provided similar results, as they are both limited by the soil's saturated hydraulic conductivity (i.e., parameter K2S). Therefore, both measures performed well only for the design event with the lowest intensity. Among the SCMs that are placed within existing pervious areas, these two SCM types take up the most space; 6.
The changes in the GR area have a direct influence on the size of the impervious area. Thus, the optimal GR area depends on the outflow from impervious and pervious areas. The outflow from impervious and pervious areas for the design event with a return period of 50 years almost matches the target outflow (i.e., Q5), if GRs replaced 40% of existing impervious areas.
Finally, this approach is transferable to any catchment and enables the discovery of viable RR models and target-based SCM design. To do so, the user has to specify a conceptual model of the selected case study area, compliant with the knowledge library. Furthermore, the user has to provide rainfall and flow measurement data and catchment characteristics (i.e., land use, topography, soil properties), based on which the values of the model constants can be determined.
The currently defined knowledge library encompasses knowledge defined in SWMM, but it is not limited to it. Namely, it can be further developed and upgraded with new SCM elements, equations, and parameters. Thus, it can serve as an open and flexible vehicle for testing new functionalities. In future research, a search among scenarios that integrate multiple SCMs can be conducted. Furthermore, the automated design of SCMs can be based on additional criteria (i.e., investment costs, biodiversity) that can be included within new objective functions.