Coupling a Distributed Time Variant Gain Model into a Storm Water Management Model to Simulate Runoffs in a Sponge City

Abstract

The storm water management model (SWMM) has been used extensively to plan, implement, control, and evaluate low impact development facilities and other drainage systems to solve storm-related problems in sponge cities. However, the calibration of SWMM involves a variety of sensitive parameters and may bring significant uncertainties. Here we incorporated the distributed time variant gain model (DTVGM), a model with a simple structure and few parameters, into the SWMM (called DTVGM-SWMM) to reduce the complexity but keep the mechanistic representation of the hydrological process. The DTVGM runoff module parameters were calibrated and validated using the Nash–Sutcliffe efficiency (NSE), based on measured data and the results of SWMM. It was then coupled with the SWMM routing module to estimate catchment runoffs and outflows. Finally, the performance was evaluated using NSE (0.57~0.94), relative errors of the flow depth (−7.59~19.79%), and peak flow rate (−33.68~54.37%) under different storm events. These implied that the DTVGM-SWMM simulations were generally consistent with those of the control group, but underperformed in simulating peak flows. Overall, the proposed framework could reasonably estimate the runoff, especially the outflow process in the urban catchment. This study provides a simple and reliable method for urban stormwater simulation.

1. Introduction

Rapid urbanization has caused various problems in constructed areas [1]. Increased surface runoff, frequent flooding, and water deterioration have significantly impacted the environment and the economy [2,3]. In the last decade, the Chinese government has launched sponge city construction [4], which often uses low impact development facilities (LIDs) as source control measurements [5]. A sponge city is a city with low impact development infrastructure, which aims to make the catchment hydrological response approach the pre-development status with minimum cost. In sponge cities, the runoff response to rainfall has significant nonlinear characteristics, which can be depicted by the storm water management model (SWMM), widely used in urban hydrology [6,7,8]. However, SWMM requires many parameters and much measured data; moreover, it operates in a complicated way.

On the other hand, the distributed time variant gain model (DTVGM) has few parameters and a simple structure for representing the rainfall-runoff response [9,10]. It can simulate complex nonlinear hydrological processes, but is generally used in large natural watersheds with low-intensity human activities under prolonged rainfalls. Thus far, DTVGM has rarely been used to simulate the rainfall runoff process in sponge cities as urban areas always have small areas, complicated land use and land cover, and, in most cases, short and intense storms. It would be questionable to apply it to a study subject with a small spatiotemporal scale [11]. Moreover, a direct reason for its rare usage is that DTVGM lacks a routing module for the drainage networks, which the SWMM can provide (written in C).

To this end, we aim to incorporate the runoff module of DTVGM into the routing module of SWMM for developing the DTVGM-SWMM to predict the runoff and outflow in a sponge city. The proposed model will contribute to the modeling and management of the rainfall-runoff process in sponge cities. This article is organized as follows: Section 2 presents the study area and data; Section 3 introduces the SWMM, DTVGM runoff module, coupling method, and performance indicators; Section 4 presents the results and a discussion from four perspectives; Section 5 summarizes the conclusions.

2. Study Area and Data

Fengxi New City is located between Xi’an and Xianyang (Figure 1b), Shaanxi Province, China (Figure 1a), with a total area of 143 km². Since 2013, the study area (22.5 km², Figure 1d) has been developed with a LID-based stormwater management system (Figure 1c). Fengxi was designated one of China’s pilot sponge cities in 2015 and a UNESCO global eco-hydrology demonstration site in 2019 (http://ecohydrology-ihp.org/demosites/view/1220 (accessed on 16 January 2023)), for featuring a temperate continental monsoon climate and nonself-weight collapsible loess land cover. The average annual temperature is about 13.6 °C, and the average annual rainfall depth is about 520.0 mm. Precipitation is mainly concentrated from July to September (50–60% of the annual amount), resulting in severe summer floods. The land cover includes built areas, roads, grasses, trees, bare grounds, water, etc. The study area has an average slope of 0.5%, with an elevation of 378.0–392.0 m, and outflows from outfalls are drained to the Weihe River and Fenghe River.

Most of the study area data, including measured rainfall and runoff data and the input file for SWMM, was provided by the Fengxi New City Management Committee; other data were collected or measured by our group.

Table 1. Statistics of storm events used to calibrate and validate the SWMM and DTVGM-SWMM.

Storm Event Name	Return Period (Year)	Duration (min)	Time-to-Peak Coefficient	Depth (mm)	Mean Intensity (mm/min)	Usage
Chicago-1 1	1	120	0.35	20.67	0.17	Calibrate DTVGM-SWMM
Chicago-2	2	120	0.35	27.92	0.23	Validate DTVGM-SWMM
Event-20170820 2	n/a	96	0.57	13.40	0.14	Calibrate SWMM and Validate DTVGM-SWMM
Event-20170909	n/a	968	0.47	16.00	0.02	Calibrate SWMM and Validate DTVGM-SWMM
Event-20170916	n/a	771	0.36	11.40	0.01	Calibrate SWMM

¹ For simplicity, Chicago denotes a rainfall time series calculated using Equation (1) of Yang et al. [12]. ² Event denotes a rainfall time series and the corresponding outflows recorded at that time.

lists the statistics of the storm events used to calibrate and validate the models.

3. Methodology

We ran the DTVGM and SWMM in the study area to account for rainfall-runoff heterogeneity between subcatchments (i.e., grid cells) within the catchment. A framework for coupling the DTVGM and SWMM was proposed to simulate runoffs and outflows in sponge cities (Figure 2). The flowchart consists of four parts: (1) the preprocessing and splitting of the observed storm events, including design storms and corresponding simulated flows, into a calibration group and a validation group; (2) calibration and validation of the SWMM model to simulate the rainfall runoff in subcatchments under storm events; (3) calibration of the DTVGM runoff module using a genetic algorithm [13]; (4) coupling the calibrated DTVGM runoff module with the SWMM routing module to develop the DTVGM-SWMM; (5) evaluating the coupling model using performance matrices.

3.1. Storm Water Management Model

SWMM is a widely used rainfall-runoff model developed by the US Environmental Protection Agency (https://www.epa.gov/water-research/storm-water-management-model-swmm (accessed on 16 January 2023)). The software consists of hydrological, hydrodynamic, and water quality modules and helps predict runoff quantity and quality from LIDs and drainage systems. It is free, open-source, and allows for the easy achievement of secondary development [14].

Based on the topography and underlying conditions, the study area was divided into 268 subcatchments (Figure 1c). In this paper, the SWMM model has 603 nodes, 603 sections of pipeline, and 25 outfalls (Figure 1d). Table 2 lists the values of the options and main hydrological parameters of SWMM [12,15]. Note that the SWMM used here was calibrated and validated with measured data, i.e., Event-20170820, Event-20170909, and Event-20170916.

Table 2. Values of options and main calibrated parameters for the storm water management model.

Options or Parameters	Value 1
Infiltration model	Horton [16,17,18]
Routing model	Dynwave [6,19,20,21,22,23,24,25,26,27]
Reporting time step (minute)	1
Routing time step (second)	10
Catchment slope (%)	0.5
Imperviousness (%)	65.0
Percent of the impervious area with no depression storage (%)	35.0
Depression storage in impervious areas (mm)	2.1
Depression storage in pervious areas (mm)	3.6
Conduit roughness (s/m1/3)	0.013
Surface roughness for overland flow in impervious area (s/m1/3)	0.013
Surface roughness for overland flow in pervious area (s/m1/3)	0.150
Minimum infiltration rate on Horton curve (mm/h)	3.56
Maximum infiltration rate on Horton curve (mm/h)	25.40
Decay rate constant of Horton curve (1/h)	7
Drying time (day)	7

¹ Data from Yang et al. [12,28]. For more details on SWMM, please refer to Rossman [29,30].

Table 2. Values of options and main calibrated parameters for the storm water management model.

Options or Parameters	Value 1
Infiltration model	Horton [16,17,18]
Routing model	Dynwave [6,19,20,21,22,23,24,25,26,27]
Reporting time step (minute)	1
Routing time step (second)	10
Catchment slope (%)	0.5
Imperviousness (%)	65.0
Percent of the impervious area with no depression storage (%)	35.0
Depression storage in impervious areas (mm)	2.1
Depression storage in pervious areas (mm)	3.6
Conduit roughness (s/m1/3)	0.013
Surface roughness for overland flow in impervious area (s/m1/3)	0.013
Surface roughness for overland flow in pervious area (s/m1/3)	0.150
Minimum infiltration rate on Horton curve (mm/h)	3.56
Maximum infiltration rate on Horton curve (mm/h)	25.40
Decay rate constant of Horton curve (1/h)	7
Drying time (day)	7

¹ Data from Yang et al. [12,28]. For more details on SWMM, please refer to Rossman [29,30].

3.2. Runoff Module of Distributed Time Variant Gain Model

The DTVGM divides the catchment into cell grids (i.e., subcatchments) based on GIS/DEM and applies the time variant gain model, a conceptual hydrological model, to each cell grid for runoff and routing calculations. Previous studies have shown that the DTVGM can perform well in large watersheds under prolonged rainfalls. In this section, we only explain the runoff module of DTVGM because the routing module of DTVGM is not used.

The number of grids in the DTVGM equals the number of subcatchments in the SWMM. For the ith subcatchment at time t, the runoff (R_i(t), mm) is the sum of the surface runoff ($R_i^S\left(t\right)$, mm) and subsurface runoff ($R_i^{SS}\left(t\right)$, mm). The surface runoff yielded in the subcatchment is calculated thus:

$$R_i^S\left(t\right)=G_i\left(t\right)P_i\left(t\right)$$

(1)

where P_i(t) is the rainfall depth in the ith subcatchment at time t, mm. G_i(t) is the time variant gain factor (determined by soil moisture) in the ith subcatchment at time t:

$$G_i\left(t\right)=g_i^1+g_i^{2}AP{I}_i\left(t\right)$$

(2)

where $g_i^1$ and $g_i^2$ are dimensionless time variant gain parameters of surface runoff for the ith subcatchment. API_i(t) is the antecedent precipitation index of the ith subcatchment at time t, which represents the recent soil moisture state of the subcatchment, mm [31,32,33]:

$$AP{I}_i\left(t\right)={\displaystyle \sum }_{j=1}^t\left(1-e^{-\frac{1}{K_i}}\right)e^{-\frac{j}{K_i}}{P}_i\left(t-j+1\right)$$

(3)

where j is the time and t is the number of times; K_i, the dimensionless decay factor and is related to the soil property and evaporation within the ith catchment.

Lastly, for the ith subcatchment at time t, the subsurface runoff is calculated as:

$$R_i^{SS}\left(t\right)=g_i^{3}AP{I}_i\left(t\right)$$

(4)

where $g_i^3$ is the dimensionless time variant gain parameter of subsurface runoff for the ith subcatchment.

3.3. Coupling Model

Though a monitoring system has been developed at a few sites in the catchment, there were no observed rainfall, runoff, and outflow data at a subcatchment scale. Thus, we calibrated the DTVGM runoff module using the genetic algorithm (https://ww2.mathworks.cn/help/gads/ga.html (accessed on 16 January 2023)) based on the subcatchments’ runoff time series calculated by SWMM. For the ith subcatchment, the decision variables of the optimization problem are $g_i^1$, $g_i^2$, $g_i^3$, and $K_i$. Thus, we set the number of iterations to 8000, namely, 2000 times the number of decision variables.

The Nash-Sutcliffe efficiency (NSE) [28,34,35,36,37,38,39] was used to quantify the overall consistency between the hydrographs by DTVGM and observed values (including some hydrographs by SWMM). Note that NSE ranges from $-\infty$ to 1, being equal to 1 if a perfect agreement exists. It is calculated as follows:

$$NS{E}_i=1-\frac{{{\displaystyle \sum }}_{t=1}^N{\left[R_i^{SWMM}\left(t\right)-R_i^{DTVGM}\left(t\right)\right]}^2}{{{\displaystyle \sum }}_{t=1}^N{\left[R_i^{SWMM}\left(t\right)-\overline{R_i^{SWMM}}\right]}^2}$$

(5)

where N refers to the number of time steps of the runoff time series, min; $Q_i^{SWMM}\left(t\right)$ and $Q_i^{DTVGM}\left(t\right)$ are the runoff time series simulated by SWMM and DTVGM for the ith subcatchment at time step t, respectively, m³/s; $\overline{Q_i^{SWMM}}$ is the average value of runoff time series simulated by SWMM for the ith subcatchment, m³/s.

The objective function is to maximize the NSE value for the ith subcatchment. Here, the calibration was conducted 30 times, thus obtaining 30 sets of optimal parameters (decision variables) with corresponding NSE values. Thus, the optimization model for the calibration can be expressed as follows:

$$\begin{array}{c}{f}_i\left(g_i^1,g_i^2,g_i^3,K_i\right)=min\left(-NS{E}_i\right)\end{array}$$

(6)

subject to:

$$\begin{array}{c}-1\le g_i^1\le 1\\ 0\le g_i^2\le 2\\ 0\le g_i^3\le 1\\ 0\le K_i\le 100\end{array}$$

(7)

Generally, developing a DTVGM-SWMM consists of three steps: (1) call the MATLAB engine (https://www.mathworks.com/help/matlab/calling-matlab-engine-from-c-programs-1.html?s_tid=CRUX_lftnav (accessed on 16 January 2023)) in Visual Studio (https://visualstudio.microsoft.com/ (accessed on 16 January 2023)) to run the runoff module of DTVGM that was developed in MATLAB by our group; (2) use the runoff results of DTVGM to replace the values of vOutflow in SWMM, which is a variable that refers to the rate of flow leaving from the catchment; (3) run the SWMM routing module to obtain catchment outflow time series.

3.4. Performance Criteria

We used three indicators, namely, the NSE (Equation (5)), the relative error of runoff depth (${\delta }_R$, Equation (8)), and that of peak runoff (${\delta }_P$, Equation (9)) to evaluate the coupling model. A higher NSE and a lower δ_R and δ_P indicate better performance. As defined in Equations (8) and (9), the values of δ_R and δ_P are positive when the DTVGM overestimates and are negative when it underestimates.

$${\delta }_R\left(i\right)=\left[\frac{R_i^{DTVGM}-R_i^{SWMM}}{R_i^{SWMM}}\right]\times 100\%$$

(8)

$${\delta }_P\left(i\right)=\left[\frac{\underset{1\le t\le N}{max}{R}_i^{DTVGM}\left(t\right)-\underset{1\le t\le N}{max}{R}_i^{SWMM}\left(t\right)}{\underset{1\le t\le N}{max}{R}_i^{SWMM}\left(t\right)}\right]\times 100\%$$

(9)

4. Results and Discussion

4.1. Simulation of Storm Water Management Model

The relevant parameters of SWMM were adjusted; for details on the calibration and validation, please refer to Yang et al. [12]. In addition, we used a runoff coefficient-based method to verify the SWMM parameters. Theoretically, the integrated runoff coefficient for the study catchment is 0.45~0.60 because the site is densely built-up [40]. In fact, under Event-20170820, the subcatchment runoff coefficients calculated by the SWMM simulation were 0.52–0.57 for the no LIDs scenario and 0.18–0.28 for the LIDs scenario. Thus, the calibrated SWMM is satisfactory for calibrating and validating the DTVGM runoff module and the coupling model.

4.2. Performance of Distributed Time Variant Gain Model

Figure 3a shows the performance of the DTVGM-SWMM based on the NSE values of each subcatchment in the calibration phase (i.e., Chicago-1) and validation phase (i.e., Chicago-2 and Event-20170820). We can see that most subcatchments under the three storms have an NSE > 0.6, with only five exceptions: subcatchments 13, 31, 171, and 238 under Event-20170820. Figure 3b further illustrates these cases, displaying a positive correlation trend between the calibration and validation NSE values. For most subcatchments, especially the underperformed five subcathments, NSE values under Chicago-2 are greater than those under Event-20170820 because the former pattern is more similar to Chicago-1, which was used for calibration. Figure 3c presents the violin plot to analyze the statistical characteristics of the NSE values. We found that the average (median) NSE values for the runoff simulation under the three storm events were 0.92, 0.71, and 0.84 (0.91, 0.64, and 0.85), respectively, indicating an acceptable model performance. Overall, the DTVGM could simulate the runoff over 263 subcatchments (98.1% of all).

Figure 4 shows the calibrated parameters (g₁, g₂, g₃, and K) of the DTVGM runoff module. The values of g₁ and g³ for different subcatchments displayed more variance than the values of g₂ (Figure 1a,b); the values of K had a considerable variation range (Figure 1c,d). In a subcatchment, the value of the four parameters reflected the complicated runoff generation characteristics in urban areas because g₁ and g₂ are gain parameters of surface runoff, and g₃ is that of subsurface runoff. For example, in most subcatchments, a greater g₁ value denotes a higher conversion from rainfall to runoff; for g₂, a greater value denotes a greater impact of soil moisture on the surface runoff generation. However, to some extent, the four parameters were not entirely based on physical mechanisms, but were obtained using a statistical method. Further work should be conducted on the parameter explanation based on a detailed investigation of the hydrological conditions in each subcatchment.

Table 3. Statistics of calibrated parameters of distributed time variant gain model for all subcatchments.

Variable	Range of Variation	Calibrated Value
		Average	Median	Minimum	Maximum
g1	[−1, 1]	−0.241	−0.281	−1.000	−0.004
g2	[0, 2]	1.326	1.575	0.000	2.000
g3	[0, 1]	0.169	0.138	0.078	0.583
K	[0, 100]	62.469	61.906	2.619	100.000

lists the statistics of the calibrated parameters of DTVGM for all subcatchments. We found that all the parameters’ values were within the range of variation, i.e., subjected to the constraints in Equation (7). Moreover, these calibrated parameters were consistent with other studies on DTVGM [10]. Overall, these results suggest that the calibration of DTVGM was reasonable.

4.3. Catchment Runoff and Outflow Simulations

The catchment runoff time series was calculated as the sum of the subcatchments’ runoff time series. For a subcatchment under a given storm event, two runoff time series were obtained by the SWMM simulation and the DTVGM-SWMM. Then, they were used to calculate the NSE with Equation (5), the relative error of runoff depth with Equation (8), and the relative error peak runoff with Equation (9). The calibrated DTVGM runoff module, determined based on the NSE values in Section 4.2, was used to generate the subcatchment runoff time series under different rainfalls. These series were finally used to drive the SWMM routing module to obtain the catchment outflow series.

We analyzed the applicability of DTVGM-SWMM by comparing the simulation results with those of SWMM in Figure 5, which illustrates the time series of rainfall, and the corresponding catchment runoff and outflow (simulated by DTVGM and SWMM) in calibration and validation phases. We found that the outflow time series were flattened after routing, and the times to peak were delayed. In addition, we found that the simulated hydrographs by DTVGM or DTVGM-SWMM demonstrated a good fit to those by SWMM under all storm events, especially in the rising and receding limbs.

However, there were gaps in the magnitudes of the peak runoffs of two simulations. Specifically, the DTVGM runoff module underestimated the peak runoff in the calibration phase of Chicago-1 (Figure 5a) and the validation phase of Event-20170909 (Figure 5d), but overestimated them in the validation phase of Chicago-2 and Event-20170820 (Figure 5b,c); similar phenomena were also found for the catchment outflow estimation by DTVGM-SWMM. No pattern of underestimation or overestimation was observed from these simulations. That is, the calibrated DTVGM demonstrated considerable uncertainty in predicting the peak rate of runoff. This can be explained by the fact that only the NSE, which reflects the accuracy and precision of the simulation for the entire time, was chosen as the objective function for calibrating the DTVGM runoff module. To conclude, the peak outflows of DTVGM-SWMM under three validation storms were very different from those of SWMM, one main reason being the poor simulation of peak runoff by DTVGM.

4.4. Performance of Coupling Model

Table 4. Depths and peak flow rates of runoff and outflow simulated by storm water management model and the coupling model under four storm events.

Storm	Chicago-1	Chicago-2	Event-20170820	Event-20170909
Precipitation (mm)	20.67	27.92	13.40	16.00
SWMM model
Runoff depth (mm)	2.41	4.07	1.06	2.75
Outflow depth (mm)	2.45	4.12	1.08	2.65
Peak runoff rate (m3/s)	10.90	23.54	1.99	8.72
Peak outflow rate (m3/s)	3.62	7.76	0.98	3.72
DTVGM-SWMM model
Runoff depth (mm)	2.53	4.67	1.27	2.54
Outflow depth (mm)	2.47	4.53	1.22	2.45
Peak runoff rate (m3/s)	7.23	24.65	3.07	5.92
Peak outflow rate (m3/s)	3.50	10.63	1.12	2.52

lists the statistics of the runoff and outflow time series simulated by SWMM and DTVGM-SWMM under different storms (illustrated in Figure 5). The depths and peak flow rates of runoffs and outflows were as given for the subsequent evaluation. We observed that the outflow depths might be greater than the runoff depths under a given storm for the SWMM model, but not for the DTVGM-SWMM model because SWMM considers the rainfall-derived inflow and infiltration into sewage infrastructure [41,42], while the coupling model neglects this.

Figure 6 summarizes the performance criteria of DTVGM-SWMM for estimating the runoff and outflow time series at a catchment scale under storm events (described in Table 1). We observed that:

First, the NSE values of the runoff and outflow simulations during the validation phase (Chicago-2, Event-20170820, and Event-20170909) were smaller than those for the calibration phase (Chicago-1), which is to be expected. The NSE minimum was 0.57 (runoff simulation under Event-20170909), suggesting a satisfactory accuracy for runoff and outflow simulations. In any case, for all storm events, the NSE value of the outflow simulation was higher than that of the runoff simulation. During the validation phase, the runoff simulation under Event-20170820 outperformed that under Chicago-2 and Event-20170909. The outflow simulation was similar. The NSE values of outflow simulation showed a strong correlation with the values of runoff simulation. In short, the NSE values indicate an acceptable performance for overall catchment runoff and outflow simulation.

Second, during validation, the absolute maxima of the relative error of flow depth (δ_R) were 19.79% for the runoff simulation and 13.37% for the outflow simulation; they were both under Event-20170820, which had the largest NSE values. We found that when the calibrated DTVGM-SWMM simulated runoff well in terms of overall fit, the runoff depth would be poorly estimated. A similar pattern could be found for the outflow simulation.

Third, δ_R values under Chicago-1, Chicago-2, and Event-20170820 were positive, indicating that estimations by DTVGM-SWMM were higher than those of SWMM; under Event-20170909, it was the opposite. For any storm event, the outflow depth simulation was better than the runoff depth simulation. Overall, the runoff and outflow depth simulations were reasonable and acceptable.

Fourth, the absolute maxima of the relative error of peak flow (δ_P) were 54.37% for runoff (under Event-20170820) and 36.98% for outflow (under Chicago-2). In most cases (7 out of 8), ${\delta }_P$ was greater than ${\delta }_R$ except in the runoff simulation under Chicago-2 (i.e., δ_P = 4.70, δ_R = 14.53).

Overall, the DTVGM-SWMM had a sizable relative error in runoff depth but a small relative error in the outflow depth after the routing process. The reason is that only the NSE was considered in the objective function for the model calibration without considering the accuracy of simulating the peak runoff where the duration is short. The DTVGM runoff module was well-calibrated for estimating the runoff process at the catchment scale, and the DTVGM-SWMM showed high accuracy in outflow simulation.

5. Conclusions

We proposed and validated the DTVGM-SWMM by combining the DTVGM runoff module with the SWMM routing module to develop a simple and efficient model for rainfall-runoff simulation in a sponge city. We conclude that: (1) the NSE and relative error of flow depth are good when using the calibrated coupling model to simulate the catchment runoff and outflow processes, especially for estimating outflows; (2) the relative errors of peak flow rate were significant for runoff and outflow simulation, indicating that the calibrated DTVGM-SWMM is bad in predicting peak flow rates for urban catchment. Overall, the proposed model can characterize the nonlinear characteristics of rainfall-runoff with few parameters; thus, it can be used for effective and efficient prediction, like machine learning models.

This research demonstrates that the DTVGM-SWMM excels in simulating catchment runoff and outflows in a sponge city. Future studies can be conducted on: (1) packaging the DTVGM-SWMM into software, like Mat-SWMM [14], and developing a feature that models the runoff for different underlying surfaces and the LIDs constructed in sponge cities; (2) improving the calibration by taking the peak flow rate as another objective function, as well as using more performance indicators; (3) investigating and quantifying the uncertainty (input data, structure, and parameters of the DTVGM and SWMM used), sensitivity, efficiency, and ease of development, and comparing them with other distributed models; and (4) developing and integrating a multi-objective decision-making tool for LIDs implementation [43,44].