Comparison of Nonlinear Reservoir and UH Algorithms for the Hydrological Modeling of a Real Urban Catchment with EPASWMM

Abstract

This paper presents a comparative analysis of two hydrological models in the Storm Water Management Model (SWMM) software, namely, the non-linear reservoir (N-LR) and the unit hydrograph (UH), on the urban catchment of Cascina Scala, Pavia in Italy. The two models were applied for the simulation of the rainfall-runoff transformation in the 42 sub-catchments in Cascina Scala, while flow routing in the underground channels was simulated by means of the De Saint-Venant equations. A dataset of rainfall and runoff for 14 events from 2000 to 2003 was adopted for the calibration and validation of the models. The calibration was performed on 7 out of the 14 events by maximizing the fit of modeled-to-measured hydrographs in the final channel of the system. Prediction performance was assessed through different indices. Results from both models fit measured data well in terms of the total hydrograph. Whereas the time to peak was reliably predicted by both models, the N-LR was found to slightly outperform the UH in terms of total volume and peak flow prediction, though it requires a more detailed knowledge of the system for its implementation. Accordingly, the UH must be preferred in the case of a scarcity of data.

1. Introduction

The high rate of urbanization and the transformation of vegetated areas into impervious ones have significantly impacted water hydrological processes [1]. An increase in urbanization potentially results in increased runoff volumes, flash floods, and intensive water quality degradation [2,3]. Moreover, climate change consequences are putting more pressure on water network systems through extreme precipitation events [4,5]. Therefore, the analysis and design of urban catchments has become a crucial subject to achieve sustainable water resource management. In this regard, an accurate estimation of rainfall-runoff processes is of crucial importance for proper drainage system analysis and planning, water storage design, and flood mitigation measures. However, rainfall-runoff transformation is a complex hydrological task involving several variables related to rainfall patterns, catchment characteristics, soil type, and physiography [6]. Indeed, several physical processes are present in urban hydrology, such as precipitation and interception, infiltration, water movement through saturated and unsaturated soil and groundwater, evaporation, transpiration, and runoff [7]. The degree of complexity varies with the catchment’s spatial and temporal scales. Rainfall-runoff models can be categorized depending on their spatial representation (lumped or distributed) and/or temporal representation (event-based or continuous time) [8,9]. Lumped models consider space-average values of model parameters to represent the entire catchment [10]. Distributed models consider, instead, space-variable parameters in the rainfall-runoff transformation [11]. Although distributed models are helpful in comprehending the physics of hydrological processes, they need accurate calibration of several parameters, which represents a very challenging task [12]. The rainfall-runoff models can be calibrated using event-based or continuous storm event simulation. Event-based models consider a single rainfall event with a total duration of some hours. Continuous simulation models instead use an extended period made up of rain events and inter-event dry periods, with a total duration ranging from months to years. In general, for the event-based approach, losses include only the infiltration process during the simulation, whereas in continuous modeling, the evapotranspiration process is also taken into account [13]. Several research works have used continuous calibration models [14,15,16,17,18]. Continuous calibration delivers an adequate estimation in terms of total runoff volume [19,20]. Event-based calibration instead provides a better estimation of time to peak, peak flow rate, and overall hydrograph shape, which are critical prerequisites for developing an effective pollutant dynamic modeling.

Distributed models typically use conceptual hydrological modeling to simulate rainfall-runoff transformation outside the channel network. The unit hydrograph (UH) [21] and the non-linear reservoir (N-LR) [22] are two methods often adopted in this hydrological modeling. In the N-LR, the surface runoff is generated by modeling the study catchment area as a nonlinear reservoir. In the UH, the surface runoff is reconstructed with the convolution integral, starting from the response of the catchment to an instantaneous input of rainfall with a unitary depth.

The open-source software Storm Water Management Model, created by the United States Environmental Protection Agency (US EPA-SWMM), [23] is one of the most frequently used software programs to simulate the rainfall-runoff transformation and the related hydrological processes in urban catchments [24]. Either the N-LR or the UH can be used in the US EPA-SWMM to model the formation of runoff over external catchment surfaces. Then, the De Saint-Venant equations are used for modeling the flow routing in underground channels, which are fed by external catchment surfaces. EPA-SWMM can be used to perform event-based or continuous simulations and can be applied to combined sewer systems, sanitary sewer systems, and catchments containing storm drains [23,25]. The accuracy and reliability of rainfall-runoff models are highly dependent on how the parameter values are defined [26,27]. Some of these can be physically measured from the catchment characteristics (e.g., area and slope). Others require parameterization based on high-quality rainfall-runoff measured data. This parameterization consists of finding the set of parameters that maximize the fit of model results to experimental observation. Traditionally, this process is carried out manually by using the trial-and-error method. However, manual parameterization is time consuming and subjective, and it depends on expert judgement. Furthermore, it may fail to reach the global optimal solution [28]. Therefore, automatic calibration methods using computer-based tools have been widely developed and used to obtain a more robust and efficient estimation of the best-fit parameter values [29,30,31]. In general, automatic calibration algorithms can be classified into deterministic and stochastic methods. Deterministic optimization algorithms are developed to locate the optimal parameter set close to a starting solution. Therefore, they are generally fast local optimization methods. To avoid getting trapped in local optima, stochastic optimization methods can be used, though they generally feature a higher computational overhead. The presence of many parameters [25] makes stochastic optimization methods preferable in hydrological models [32,33,34]. The genetic algorithm (GA) [35] is one of the efficient optimization approaches known for its capability of detecting global optimal solutions. It is a population-based heuristic algorithm inspired by the concept of natural selection that drives biological evolution. The GA has been applied successfully to several environmental research problems [36,37,38]. It is also widely used to optimize the parameters of rainfall-runoff transformation [39,40].

This research work attempts to develop an automatic calibration method for the event-based rainfall-runoff transformation in a small urban catchment. The main novelty of the present work consists of the application and comparison of the two rainfall-runoff models N-LR and UH in the framework of the EPA-SWMM software. In fact, to the best of our knowledge, no previous work has ever attempted to compare the performance of these models on a real case study.

The following sections describe the case study, methodology and applications, followed by the discussion and conclusions.

2. Case Study

The selected study area is Cascina Scala, which is located in northern Pavia, Italy [41] (Figure 1). It is an urban catchment hosting 1500 residents. The contribution area is 12.7 ha, where 7.9 ha (62%) is impervious, and 4.8 ha (38%) is pervious.

The impervious area is composed of 2.88 ha (22%) of roofs and 5.02 ha (39.6%) of streets and paved surfaces. The area is quite flat, with an average slope of 0.15% inclined from northwest to southeast and connected to the urban drainage system directly. The catchment can be separated into 42 sub-catchments, with a single pipe drain for each sub-catchment.

Table 1. Cascina Scala sub-catchment characteristics. See footnote for description of the symbols.

ID	A (ha)	S (%)	IA		ID	A (ha)	S (%)	IA
			R (%)	S&S (%)				R (%)	S&S (%)
1	0.27	0.5	0	29.2	22	0.50	0.5	20.4	56.3
2	0.20	0.2	21.2	54.8	23	0.13	0.1	0	85.9
3	0.39	0.1	23.1	40.7	24	0.81	0.1	23.5	47.8
4	0.28	0.1	32.8	54.4	25	0.41	0.2	24.2	31.9
5	0.07	0.1	0	100	26	0.13	0.1	0	100
6	0.49	0.1	32.5	40	27	0.39	0.1	22.1	26
7	0.09	0.1	0	100	28	0.21	0.1	37.9	9
8	0.16	0.1	32.4	33.2	29	0.20	0.1	0	77.9
9	0.07	0.1	24.5	43.2	30	0.06	0.1	0	75.4
10	0.57	0.1	5.1	35.6	31	0.10	0.1	0	80.8
11	0.12	0.1	12.4	62.7	32	0.50	0.1	37.6	21.9
12	0.36	0.3	30.8	30.5	33	0.05	0.2	19.5	80.5
13	0.20	0.3	31.3	46.8	34	0.33	0.1	23.1	40.2
14	0.16	0.1	20.3	53	35	0.42	0.1	32	41.5
15	1.35	0.1	26.3	25.6	36	0.27	0.1	16	40.7
16	0.21	0.1	29.3	37.1	37	0.24	0.1	16.3	65.7
17	0.11	0.1	23.9	43.4	38	0.16	0.2	14.4	51.9
18	0.06	0.3	0	100	39	0.05	0.1	0	48
19	0.05	0.1	30.6	54.3	40	0.19	0.6	29.6	43.9
20	0.17	0.1	26.2	42	41	0.12	0.3	25.2	49.6
21	0.47	0.1	26.5	29.3	42	0.22	0.3	26.8	57.3
					42 *	1.55	0.3	25.1	21.1

Footnote: For the 42 sub-catchments, ID: identifier of sub-catchments, A: area, S: slope, IA: impervious area, R: roof, S&S: streets and squares; 42 * indicates the modified catchment 42 after May 2001.

shows the characteristics (area, slope, and impervious area percentage divided into roofs and streets) of each sub-catchment.

After May 2001, a new area (1.33 ha) was connected to the urban drainage system, specifically to sub-catchment number 42. Therefore, the area of sub-catchment 42 was considered to grow from 0.22 to 1.55 ha starting from May 2001 (in Table 1, sub-catchment 42 * is the upgraded version of sub-catchment 42). Two tipping-bucket rain gauges were installed in the catchment area with 0.2 mm accuracy (P1 and P2 in Figure 1) to collect the rainfall data. The two rainfall gauges have a funnel with an area of 1000 cm². The presence of two gauges with a mutual distance of 310 m ensured that the spatial uniformity of the precipitation was controlled. Furthermore, it provided an accurate assessment of rainfall volume. To measure the runoff at the outlet of the channel network, an ISCO (ISCO, Genova, Italy) 4230 bubbler depth meter was placed upstream from a Venturi flume at the final reach of the channel network in the southeastern corner of the catchment (Figure 1). The system of underground channels in the Cascina Scala catchment is made of 42 concrete pipes. The upstream channels are circular conduits with a diameter ranging between 0.4 m and 0.6 m. The downstream channels are egg-shaped conduits with horizontal and vertical sizes ranging between 0.6–0.9 and 0.7–1.05 m, respectively. The conduit slope ranges between 0.15 and 1.01%, with an average value of 0.42 %. The pipe length ranges between 16 m and 75 m, with a total length of 2045 m. Twenty-three rainfall events were monitored for both rainfall-runoff transformation and water quality in the period between June 2000 and October 2003 (more details about the events can be found in [42]). As for the rainfall-runoff transformation, simultaneous measurements of rainfall intensity at the two rain-gauges and of the flow rate in the final catchment channel were performed. A representative subset was selected for rainfall-runoff simulation, meeting the following requirements:

• Regular operation of rain gages and flow meter without any instrumentation failures or pressurized flow conditions;
• Total rainfall depth of at least 5 mm;
• Maximum rainfall intensity equal to or greater than 0.1 mm/min;
• Maximum rainfall depth of at least 2 mm over 15 min.

Fourteen events (N_r,tot) were found to fit the requirements (

Table 2. Main characteristic of the selected rainfall events. The first seven events are for calibration. The following seven events are for validation.

Event	Total Rainfall Depth Vtotal (mm)	Rainfall Duration Ttotal (min)	Peak Flow Rate QMax (m3/s)
3	11.8	197	0.255
5	16.4	108	0.551
7	7	50	0.326
9	10.6	215	0.19
14	15.8	380	0.257
17	23.4	964	0.139
23	39.8	1133	0.149
8	11	64	0.376
11	26.2	478	0.281
12	18.6	443	0.263
13	8.4	111	0.157
19	12.6	248	0.234
20	16.2	231	0.245
21	7	286	0.06

), having the following characteristics: total rainfall depth (V_total) ranging between 7 mm and 39.8 mm, rainfall duration (T_total) ranging between 50 and 1133 min, and peak flow rate (Q_Max) ranging between 0.06 m³/s and 0.55 m³/s.

The events were split into two subsets, i.e., the calibration and validation subsets, which were used for the optimization-based parameterization of the hydrological model and for the testing of the optimal set of parameters, respectively. The splitting was carried out while guaranteeing wide and representative ranges for V_total, T_total, and Q_Max in the two subsets.

3. Materials and Methods

In this study, the distributed modeling approach was used for the rainfall-runoff transformation. For each sub-catchment of the Cascina Scala catchment (see previous Section 2), the N-LR and the UH models were constructed to simulate the runoff formation and routing, upstream from the channel network. The following sections describe the N-LR and UH models and their implementation in EPASWMM, followed by the genetic algorithm adopted for their parameterization.

3.1. Non-Linear Reservoir (N-LR)

This method models the generic sub-catchment as a nonlinear reservoir, for which water storage is a function of inflow (precipitation) and outflow (infiltration, evaporation, and runoff), as shown in Figure 2. During the generic rain event, the reservoir storage starts to fill. Runoff starts only after water storage depth overcomes the depression storage depth associated with surface wetting, interception, and ponding.

In this study, the Horton model was used to model infiltration. The required parameters for Horton’s approach include the maximum infiltration rate, the minimum infiltration rate, the decay coefficient, the maximum infiltration volume, and the regeneration coefficient. Evaporation is ignored in this work, since the temporal scale of event-based modeling is small.

As a result of mass conservation, the net change in water depth d (m) in the sub-catchment per unit of time t (s) is the difference between inflow (rainfall) and outflow (infiltration + runoff) according to the following equation:

$$\frac{\partial d}{\partial t}=i-f-r$$

(1)

where i = rainfall intensity (m/s), f = infiltration rate (m/s), and r = runoff rate per unit of surface area (m/s).

Assuming a uniform flow over the sub-catchment surface, which features a width W (m), a slope S (-), and an area A (m²), r can be computed by Manning’s equation:

$$r=\frac{W{\left(d-d_s\right)}^{5/3}{S}^{1/2}}{n A}$$

(2)

where d_s = depth of depression storage (m), and n = Manning roughness coefficient (s/m^1/3), which depends on surface characteristics, being obtained as a function of the roughness of the various pervious and impervious areas present in the sub-catchment.

In the US EPASWMM, the use of the N-LR requires the following set of parameters to be defined for the external sub-catchments:

• Slope (-);
• Area (m²);
• Percentage of impervious and pervious area (-);
• Manning roughness coefficient for impervious area (roof) (s/m^1/3);
• Manning roughness coefficient for impervious area (street) (s/m^1/3);
• Manning roughness coefficient for the pervious area (s/m^1/3);
• Depth of depression storage on impervious area (mm);
• Depth of depression storage on pervious area (mm);
• Maximum infiltration rate (mm/h);
• Minimum infiltration rate (mm/h);
• Decay coefficient (1/h);
• Maximum infiltration volume (mm);
• Drying time (days);
• Width coefficient (-).

In the applications of the present work, the optimal set of parameters for the N-LR method was derived as follows. Three parameters (1,2,3) were directly assigned based on the geometrical and land use characteristics of the sub-catchments (e.g., see Table 1 and [41]). After preliminary analyses, not reported in this paper, proving their poor impact on the results, the two parameters of the Horton model (12,13) were set to constant values of 0 mm and 7 days, respectively, according to the values suggested in [23]. The remaining nine parameters (4,5,6,7,8,9,10,11,14) were searched for by the genetic algorithm (see Section 3.3), assumed to be equal for all the sub-catchments. For parameter 14 (width coefficient), it was then multiplied by the length of the pipe in which the single sub-catchment drains the runoff, for calculating the equivalent width of each sub-catchment. Indeed, as indicated in [23], the width is a model parameter assumed to be proportional to the length of the drainage channel. For the channel network, the following set of parameters must be defined:

15.

Manning roughness for conduit (s/m^1/3);

16.

Conduit length (m);

17.

Junction elevations (m);

In this case, two parameters (16,17) were directly assigned based on the available information about the sewerage system (see [41]), and the remaining parameter (15) was searched for by the genetic algorithm (see Section 3.3) and assumed to be equal for all the pipes, since they are all made of the same material (concrete).

3.2. Unit Hydrograph (UH)

UH is traditionally used in surface water hydrology to estimate the direct runoff hydrograph resulting from a unit input of rainfall. UH techniques can be adapted to estimate rainfall-dependent infiltration/inflow (RDII) flows into a channel network. One of the widely used and most flexible UH approaches for estimating RDII is the UH-RTK method (the RTK abbreviation comes from the key three parameters that feature the unit hydrographs used by the technique, as shown in Figure 3). Notably:

• R is the fraction of rainfall entering the channel network, also known as runoff volumetric coefficient (-);
• T is the time to peak (h);
• K is the ratio of time to recession to time to peak (-).

The RTK method was developed by [43] and became one of the most commonly used unit hydrograph methods in the storm water management field [44]. This method uses the single triangular unit hydrograph to represent the sub-catchment response to a rainfall event. These unit hydrographs can be employed for any event to generate the corresponding RDII flow rates. Figure 3 shows a single triangular UH assumed to represent the RDII flow induced by one unit of rainfall over a unit of time.

As the integral of the UH is equal to 1, its peak UH_peak at time T is given by the following formula:

$$U{H}_{peak}=2/\left(T+KT\right)$$

(3)

By applying the convolution integral, the runoff r per unit of sub-catchment area at the generic time t can be calculated as:

$$r\left(t\right)={{\displaystyle \int }}_0^{t}R i\left(\tau \right)UH\left(t-\tau \right)d\tau$$

(4)

In the US EPASWMM, the use of the UH requires the following set of parameters to be defined for each external sub-catchment:

• Area (m²);
• Parameter R (-);
• Parameter T (h;
• Parameter K (-);

In the applications of the present work, the optimal set of parameters for the UH method was derived as follows. One parameter (1) was directly assigned based on the geometrical characteristics of the sub-catchments (e.g., see Table 1), and the remaining three parameters (2, 3, 4) were searched by the genetic algorithm (see Section 3.3) for all the sub-catchments, for a total of (42 × 3 + 3 (for the sub-catchment 42 *)) = 129 variables. For the channel network, the following set of parameters must be defined:

5.

Manning roughness for the conduit (s/m^1/3);

6.

Conduit length (m);

7.

Junction elevations (m).

Also in this case, two parameters (6,7) were directly assigned based on the available information about the sewerage system (see [41]), and the remaining parameter (5) was searched by the genetic algorithm (see Section 3.3) and assumed to be equal for all the pipes.

3.3. Genetic Optimization

Inside the GA [35], a population of individuals exists, in which the generic individual is encoded in genes representing the decision variables of the problem. Following an initial random generation, the population evolves thanks to the processes of crossover and mutation, till it converges towards a global optimum in terms of fitness, represented by the objective function (OF) to be optimized.

In the present work, the GA was used to optimize the parameters of the N-LR and UH algorithms for the 42 catchments of the Cascina Scala catchment (see Section 2) to maximize the fit of modeled-to-measured hydrographs in the final channel of the system. The parameters were optimized in a subset of seven rainfall events (calibration set) considering the following OF to be minimized:

$$OF={\displaystyle \sum }_{i=1}^{N_r}{\displaystyle \sum }_{j=1}^{N_{t,i}}{\left(Q_{p,i,j}-Q_{m,i,j}\right)}^2$$

(5)

in which N_r and N_t,i are the number of rain events (seven) considered for the optimization (the first seven of Table 2) and the number of time instants in the generic i-th rain event (discretization of measured water-discharge data 1 min), respectively. Q_m,i,j (m³/s) and Q_p,i,j (m³/s) are the measured and predicted water discharges (in the final channel of the system), respectively, at the j-th time and in the i-th rain event.

Notably, the application of the N-LR requires the calibration of 10 parameters (i.e., nine for the sub-catchments and one for the channels), whereas for the UH, the number of parameters to be calibrated is 130 (i.e., 129 for the sub-catchments and one for the channels).

In this work, the GA toolbox of MATLAB^® was linked to the EPA-SWMM dll in order to perform the optimization for both the N-LR and the UH. For each individual of the GA, the sequence of instructions shown in Figure 4 was executed.

First, the parameters of either the N-LR or the UH were defined based on the decision variables encoded in individual genes. Then, the EPA-SWMM dll was run, i.e., the hydrological module first and then the hydraulic module, to estimate the runoff generated in each sub-catchment and to simulate the flow routing in the underground channels. Finally, the hydrograph at the exit of the Cascina Scala catchment was evaluated and compared with the measured one for each event of the calibration set, leading to the evaluation of the objective function OF.

3.4. Goodness-of-Fit Indices

The performance of the calibration (carried out by minimizing Equation (5)) and validation steps was evaluated using goodness-of-fit indices, as suggested in the literature [45,46,47]. Three efficiency criteria were selected, the coefficient of efficiency (CE) (also known as Nash–Sutcliffe efficiency), the coefficient of determination (R²), and the root mean square error (RMSE), the equations of which are:

$$\mathrm{CE}=1-\frac{{{\displaystyle \sum }}_{\mathrm{k}=1}^N{\left({\mathrm{V}}_{\mathrm{m},\mathrm{k}}-{\mathrm{V}}_{\mathrm{p},\mathrm{k}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{k}=1}^N{\left({\mathrm{V}}_{\mathrm{m},\mathrm{k}}-\overline{{\mathrm{V}}_{\mathrm{m}}}\right)}^2}$$

(6)

$$R^2=\frac{{{\displaystyle \sum }}_{\mathrm{k}=1}^N{\left({\mathrm{V}}_{\mathrm{p},\mathrm{k}}-\overline{{\mathrm{V}}_{\mathrm{m}}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{k}=1}^N{\left({\mathrm{V}}_{\mathrm{m},\mathrm{k}}-\overline{{\mathrm{V}}_{\mathrm{m}}}\right)}^2 }$$

(7)

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum }_{\mathrm{k}=1}^N\left({\mathrm{V}}_{\mathrm{m},\mathrm{k}}-{\mathrm{V}}_{\mathrm{p},\mathrm{k}}\right)}²$$

(8)

where $\overline{{\mathrm{V}}_{\mathrm{m}}}$ is the mean of the measured variable V_m,k, $\overline{{\mathrm{V}}_{\mathrm{p}}}$ is the mean of the predicted variable V_p,k, N is equal to N_t,i or N_r,tot, depending on whether the generic goodness-of-fit index is calculated on the hydrograph or on one of the global variables Total Volume, peak flow Q_Max, and time to peak T_Peak.

4. Results

The results of the calibration step for the N-LR method are listed in

Table 3. The calibrated parameters of the N-LR method.

Parameter	Range of Variation	Calibrated Value
Manning coefficient for roof (s/m1/3)	0.014–0.030	0.014
Manning coefficient for street (s/m1/3)	0.010–0.055	0.011
Manning coefficient for pervious area (s/m1/3)	0.10–0.35	0.321
Depression storage on impervious area (mm)	1.27–2.54	1.27
Depression storage on pervious area (mm)	2.54–5.08	3.08
Maximum infiltration rate (mm/h)	58–170	118
Minimum infiltration rate (mm/h)	1–57	50
Decay coefficient (1/h)	2–7	4.3
Width coefficient (-)	1–2	1.99
Manning coefficient for conduit (s/m1/3)	0.012–0.02	0.012

, in which the range of variation (taken from [23]) and the final calibrated values for the 10 parameters are shown. The results of the calibration step for the UH method are listed in

Table 4. The calibrated parameters of the UH method. (42 *) indicates the modified catchment 42 after May 2001.

N	R (-)		T (h)		K (-)
	Range of Variation	Calibrated Value	Range of Variation	Calibrated Value	Range of Variation	Calibrated Value
1	0.23–0.35	0.24	0.08–0.33	0.33	0.5–2	1.99
2	0.61–0.91	0.63	0.08–0.33	0.32	0.5–2	1.46
3	0.51–0.77	0.51	0.08–0.33	0.08	0.5–2	0.57
4	0.7–1	0.73	0.08–0.33	0.09	0.5–2	0.53
5	0.8–1	0.8	0.08–0.33	0.33	0.5–2	2
6	0.58–0.87	0.58	0.08–0.33	0.33	0.5–2	1.88
7	0.8–1	0.8	0.08–0.33	0.17	0.5–2	0.62
8	0.52–0.79	0.53	0.08–0.33	0.1	0.5–2	0.69
9	0.54–0.81	0.61	0.08–0.33	0.09	0.5–2	0.52
10	0.33–0.49	0.37	0.08–0.33	0.33	0.5–2	2
11	0.6–0.9	0.67	0.08–0.33	0.08	0.5–2	0.57
12	0.49–0.74	0.55	0.08–0.33	0.08	0.5–2	0.56
13	0.62–0.94	0.64	0.08–0.33	0.33	0.5–2	1.97
14	0.59–0.88	0.59	0.08–0.33	0.08	0.5–2	0.56
15	0.42–0.62	0.53	0.08–0.33	0.09	0.5–2	0.55
16	0.53–0.8	0.69	0.08–0.33	0.08	0.5–2	0.52
17	0.54–0.81	0.57	0.08–0.33	0.09	0.5–2	0.56
18	0.8–1	0.91	0.08–0.33	0.08	0.5–2	1.5
19	0.68–1	1	0.08–0.33	0.08	0.5–2	0.5
20	0.55–0.82	0.82	0.08–0.33	0.09	0.5–2	0.5
21	0.45–0.67	0.48	0.08–0.33	0.27	0.5–2	1.91
22	0.61–0.92	0.62	0.08–0.33	0.33	0.5–2	1.68
23	0.69–1	0.72	0.08–0.33	0.09	0.5–2	1.37
24	0.57–0.86	0.62	0.08–0.33	0.08	0.5–2	0.5
25	0.45–0.67	0.58	0.08–0.33	0.08	0.5–2	0.5
26	0.8–1	0.94	0.08–0.33	0.08	0.5–2	0.5
27	0.38–0.58	0.58	0.08–0.33	0.3	0.5–2	1.2
28	0.38–0.56	0.44	0.08–0.33	0.08	0.5–2	0.89
29	0.62–0.93	0.63	0.08–0.33	0.33	0.5–2	1.86
30	0.6–0.9	0.72	0.08–0.33	0.1	0.5–2	0.52
31	0.65–0.97	0.7	0.08–0.33	0.33	0.5–2	1.7
32	0.48–0.71	0.48	0.08–0.33	0.08	0.5–2	1.16
33	0.8–1	0.97	0.08–0.33	0.08	0.5–2	0.52
34	0.51–0.76	0.7	0.08–0.33	0.08	0.5–2	0.51
35	0.59–0.88	0.59	0.08–0.33	0.1	0.5–2	0.79
36	0.45–0.68	0.67	0.08–0.33	0.08	0.5–2	0.58
37	0.66–0.98	0.98	0.08–0.33	0.08	0.5–2	0.59
38	0.53–0.8	0.59	0.08–0.33	0.08	0.5–2	0.59
39	0.38–0.58	0.4	0.08–0.33	0.16	0.5–2	0.5
40	0.59–0.88	0.88	0.08–0.33	0.09	0.5–2	0.6
41	0.6–0.9	0.9	0.08–0.33	0.08	0.5–2	0.83
42	0.37–0.55	0.55	0.08–0.33	0.21	0.5–2	1.5
42 *	0.67–1	0.99	0.08–0.33	0.09	0.5–2	0.6

, in which the range of variation and the final calibrated values for the 129 parameters related to the sub-catchments are shown.

Notably, for the parameter R, since it represents the runoff volumetric coefficient, its range of variation was set equal to ±20% of the ratio between the impervious area and the total area of each sub-catchment. For the parameter T, which represents the time to peak, it was allowed to vary, for all the sub-catchments, between 5 and 20 min, according to the topography and land use of the area. Finally, for the parameter K, according to the characteristics of the analyzed area, the range of variation was allowed to vary from 0.5 to 2 times the time to peak of the hydrograph. For the Manning coefficient of the conduit, the range of variations and the final calibrated value were equal to that of the N-LR model, and specifically n = 0.012 s/m^1/3.

Without loss of generality, the hydrographs (with the corresponding hyetographs) for event 23 (from the calibration set) and event 12 (from the validation set), obtained with the N-LR (grey solid line) and the UH (dark-gray dotted line) models, are plotted in Figure 5a,b, respectively. Measured data are shown with star symbols. Results from both models were found to be in good agreement with the measured data, either in the calibration or the validation step, with a slight underestimation of the total runoff and peak flow, and with a very good prediction of the time to peak. Notably, as can be seen from the values of the R², RMSE, and CE in Figure 5, the N-LR provides, either in calibration and validation, slightly better prediction performance:

-

In calibration, the N-LR method yields R² = 0.99, RMSE = 0.01 m³/s, and CE = 0.93, whereas the UH method yields R² = 0.98, RMSE = 0.01 m³/s, and CE = 0.90;

-

In validation, the N-LR method yields R² = 0.95, RMSE = 0.01 m³/s, and CE = 0.92, whereas the UH method yields R² = 0.89, RMSE = 0.02 m³/s, and CE = 0.87.

Table 5. Simulation results. The yellow color is used to highlight the performance of the better model.

Event		R2 (-)		RSME (m3/s)		CE (-)
		N-LR	UH	N-LR	UH	N-LR	UH
Calibration	3	0.78	0.87	0.03	0.02	0.77	0.85
	5	0.92	0.96	0.03	0.03	0.91	0.94
	7	0.51	0.42	0.07	0.08	0.50	0.30
	9	0.90	0.86	0.02	0.02	0.69	0.58
	14	0.94	0.94	0.02	0.02	0.88	0.86
	17	0.91	0.91	0.01	0.01	0.72	0.74
	23	0.99	0.98	0.01	0.01	0.93	0.90
Validation	8	0.81	0.68	0.04	0.06	0.81	0.57
	11	0.98	0.96	0.01	0.01	0.95	0.95
	12	0.95	0.89	0.01	0.02	0.92	0.87
	13	0.94	0.82	0.02	0.02	0.85	0.78
	19	0.97	0.97	0.02	0.02	0.91	0.92
	20	0.69	0.75	0.03	0.03	0.51	0.58
	21	0.85	0.88	0.01	0.01	0.49	0.44

reports the values of the three goodness-of-fit indices defined in Section 3.4, for all the events (from both the calibration and validation sets) and for both models, calculated considering measured and predicted water-discharge patterns at the catchment outlet. It is evident that the values of the RMSE are generally small, with the highest values for event 7 (in the calibration set) and event 8 (in the validation set), which are also the shortest rainfall events considered (durations of 50 min and 64 min, respectively). Therefore, both models seem to perform satisfactorily. In more detail:

-

In calibration, looking at all the indices, the N-LR has a better performance in three events (7, 9, and 23); the UH prevails in two events (3 and 5); finally, the performance is similar in the remaining two events (14 and 17).

-

In validation, looking at all the indices, the N-LR is better in four events (8, 11, 12, and 13); the UH prevails in two events (20 and 21); finally, the performance is similar in the remaining event (19).

-

Globally, the N-RL seems to perform slightly better, also exhibiting higher performance in validation.

In Figure 6, Figure 7 and Figure 8, the comparison between the prediction performance of the N-LR (light grey) and UH (dark grey) model, in the calibration (values indicated with solid circle) and validation (values indicated with solid triangle) sets, is shown in terms of total volume (Figure 6), peak flow (Figure 7), and time to peak (Figure 8). The three indices described in Section 3.4 were also applied to the total volume, peak flow, and time to peak considering both models on the whole dataset of fourteen rainfall events.

Notably, in Figure 6, the comparison between the measured and the predicted total volume is shown. The scattered dots along the bisector for both models confirm a good prediction capability either in calibration or validation. From the performance indices, it is clear that the N-LR model outperforms the UH in terms of R² (0.92 against 0.91), RMSE (183 m³ against 222 m³), and CE (0.91 against 0.86). Overall, both models overestimate the total volume, and only in one case (for the rainfall event with the largest total volume belonging to the validation set), they provide a remarkable underestimation.

Figure 7 shows the comparison between the measured and the predicted peak flow Q_Max. Also in this case, the scattered dots along the bisector confirm the good agreement of the predicted values with the measured data, either in calibration or validation, with a general slight underestimation of the peak flow, especially for those events featuring higher values of the peak water discharge. From the performance indices, it is clear that the N-LR model generally outperforms the UH in terms of R² (0.92 against 0.84) and CE (0.84 against 0.82), while showing the same average error RMSE = 0.05 m³/s.

Finally, in Figure 8, the comparison between the measured and the predicted time to peak is shown. From the performance indices and the alignment of the dots along the bisector, it is evident that both models, either in calibration or in validation, provide a very accurate estimation of the time in which the peak flow occurs, with an average error of just a few minutes (RMSE = 2.95 min for N-LR and RMSE = 3.9 min for the UH). The R² is excellent for both models (R² = 0.99).

5. Discussion

In the last decades, the Storm Water Management Model (SWMM) software of the EPA has been widely used for modeling the hydrology of urban and rural catchments due to its versatility and flexibility in reproducing the complex hydrological processes in various contexts. Overall, the good performance of this software was fully confirmed in this work on the Cascina Scala catchment. The results obtained considering both the adopted hydrological models, namely, the N-LR and the UH, were found to satisfactorily agree with the measured data in terms of the total hydrograph, total volume, peak flow, and time to peak at the catchment outlet. The N-LR showed a slightly better performance than the UH, at the price of a much more detailed amount of information required on the study area, mainly in terms of geometrical and morphological characteristics of pervious and impervious areas inside each sub-catchment. When this information is unavailable, the UH model can be conveniently adopted, as its parameterization does not require this detailed degree of knowledge. However, the easier applicability of the UH model is paid back by a larger parametrization burden (number of parameters equal to 10 and 130 for the N-LR and UH models, respectively). More specifically, both models provide a slight underestimation of the peak flow, especially for those events that feature the highest values of the peak water discharge and the shortest durations (see Figure 7 and Table 2). This could be because for modeling the catchment, a uniform spatial distribution of rainfall was assumed, and it is plausible that short and intense rainfall events are not uniformly distributed over the catchment. Furthermore, it must be highlighted that the results of the parameterization may be affected by the objective function adopted, as it was pointed out in [48], in which the authors investigated the influence of objective functions on the rainfall-runoff models’ performance. In fact, the analysis of Table 5 points out that the best performance, considering the three goodness-of-fit indices altogether, is obtained for the longest rainfall events, namely, 14, 17, and 23, for which T_total = 380, 964, and 1133 min, respectively. This can be explained by the structure of the objective function used for the optimization (Equation (5)). Indeed, it minimizes the sum of the squared difference between the measured and predicted water discharges for each time instant of each event. Accordingly, this could drive the optimizer to favor parameter sets that fit better to longer rainfall events, which obviously feature a higher number of time instants and a larger number of terms to be summed as a result. By considering the events in the calibration step all together, the optimizer attempts to find a good compromise but implicitly gives a higher weight to longer events. Indeed, thanks to the multi-objective approach, the adoption of other objective functions considering various variables of interest, such as the total volume and peak flow, is expected to yield different sets of optimal parameters in comparison with those obtained in the present work. It is worth highlighting that this aspect is easily implementable in the proposed procedure, and therefore, this issue will be the topic of future investigations.