Robust Vegetation Parameterization for Green Roofs in the EPA Stormwater Management Model (SWMM)

In increasingly expanding cities, roofs are still largely unused areas to counteract the negative impacts of urbanization on the water balance and to reduce flooding. To estimate the effect of green roofs as a sustainable low impact development (LID) technique on the building scale, different approaches to predict the runoff are carried out. In hydrological modelling, representing vegetation feedback on evapotranspiration (ET) is still considered challenging. In this research article, the focus is on improving the representation of the coupled soil–vegetation system of green roofs. Relevant data to calibrate and validate model representations were obtained from an existing field campaign comprising several green roof test plots with different characteristics. A coupled model, utilizing both the Penman–Monteith equation to estimate ET and the software EPA stormwater management model (SWMM) to calculate the runoff, was set up. Through the application of an automatic calibration procedure, we demonstrate that this coupled modelling approach (Kling–Gupta efficiency KGE = 0.88) outperforms the standard ET representation in EPA SWMM (KGE = −0.35), whilst providing a consistent and robust parameter set across all green roof configurations. Moreover, through a global sensitivity analysis, the impact of changes in model parameters was quantified in order to aid modelers in simplifying their parameterization of EPA SWMM. Finally, an improved model using the Penman–Monteith equation and various recommendations are presented.


Introduction
Increasingly expanding cities lead to larger areas of sealed surfaces, offering less space for infiltration and evapotranspiration. This results in increased flooding with stronger peaks, which may be intensified during high rainfall events. Additionally, climate change may contribute to the increase in flood risk [1,2]. In an effort to counteract the negative impacts of urbanization on the water balance and reduce flooding through retention, roofs are still underutilized areas. Various concepts focusing on more sustainable approaches to improve urban stormwater management are currently under development, such as low impact development (LID) technologies, which are ecologically based approaches that combine urban and natural systems [3].
Green roofs, as a sustainable LID technique on the building scale, are becoming an alternative to impermeable roofs and aim to integrate the built infrastructure into the landscape [4]. They basically consist of a substrate layer of growing medium covered by vegetation, a drainage system and a waterproof membrane with root penetration resistance [1,5]. Green roofs are often classified into extensive or intensive, depending on the thickness of the substrate layer. Extensive green roofs have a maximum substrate layer thickness of approx. 150 mm and are usually covered by drought resistant plants (i.e., Sedum), whereas intensive green roofs can present a greater thickness and are often covered by perennial plants (i.e., shrubs and trees) [1,6].
Green roofs provide important benefits, as they reduce peak stormwater runoff and runoff volume, improve water quality of runoff and lower the risk of frequent urban flooding. Moreover, they provide important ecosystem services, including regulation of building temperatures and reduction in urban-heat island effects [7][8][9]. The benefits of green roofs are related to their evapotranspiration (ET) performance and retention capacity, which, apart from the climatic conditions, can be linked to the vegetation and the growing medium's characteristics (i.e., the available amount of water on the substrate layer, the plant's stomatal resistance, etc.) [10].
So far, different approaches have been adopted to predict the runoff from green roofs. Reviews on published findings on green roof hydrology reveal that antecedent conditions, growing medium and plant species are among the major factors that affect a green roof's performance [11]. Because plants contribute to the ET processes (i.e., through the surface resistance) and higher ET rates may suggest a reduction in surface runoff, an accurate estimation of the parameters describing the vegetation should be provided.
In this context, the U.S. Environmental protection agency stormwater management model (EPA SWMM) is well capable of simulating the water balance of green roofs (i.e., [12]), especially when the ET processes are considered within the model [13]. Compared to other simpler models, EPA SWMM seems to be better suited for modelling detailed green roofs, due to its accuracy, algorithms and scales that are behind the different processes involved in the model [14]. Furthermore, EPA SWMM is able to carry out continuous simulations, which offer major advantages over event-based simulations when estimating design floods [15,16].
However, modelling of green roofs is still subjected to uncertainties. For instance, Burszta-Adamiak and Mrowiec [17] highlight a number of limitations of the EPA SWMM model, mainly due to its weakness to assign an initial depth in the substrate layer and the low number of parameters used to characterize the surface layer on green roofs. The release of newer versions of the EPA SWMM model helped to overcome some of the limitations faced in the past, for example, by introducing a specific green roof module [18]. Still, the performance of the model could be improved by including improved ET computations in long-term simulations and through calibration, especially considering the parameters that describe the retention capacity of water in green roofs [19]. Alfredo et al. [4] also express the need to use available experimental data to calibrate and validate the model, due to the unknown nature of some model parameters.
Niazi et al. [20] conclude that, overall, EPA SWMM can be used as a flexible tool for simulating urban hydrology, but special attention should be paid to the current gaps within the model, for instance, the lack of guidance on parameter estimation, sensitivity analysis and calibration, especially while considering LID structures for water-quality simulations. Moreover, storm event-based simulations on green roofs show that water availability in the soil becomes the most important aspect for ET processes and, thus, a good estimate on the antecedent dry period should be calculated [21].
One challenge is then the representation of the special reaction to dryness of the plant genus Sedum, which is usually used for roof greening. Because Sedum species switch quickly from C3-(basic type of photosynthesis) to Crassulacean Acid Metabolism (CAM) in response to drought stress, they are known as C3-CAM intermediates or facultative CAM plants [22]. The facultative CAM plants maximize their performance during C3 photosynthesis and ideal conditions, while switching temporally to CAM under water stress. While C3 plants lose water as by-product during photosynthesis and CO2 uptake at daytime (transpiration), CAM plants absorb CO2 in the night, while evaporation stress is lower [23].
Therefore, in the following research, the improved representation of the model-based vegetation-water balance feedback of Sedum is investigated. Similar to what was done by Feng and Burian [24], the Penman-Monteith scheme is introduced in the EPA SWMM model to determine ET. This approach allows heterogeneous and sub-daily potential ET inputs, which simple methods cannot represent (i.e., the Hargreaves equation integrated in EPA SWMM). This study aims to improve the simulated runoff from extensive green roofs, with special attention on the ET processes, based on an existing field campaign. The added value of this study is the robust estimation of vegetation parameters, especially the surface resistance in the Penman-Monteith equation, which was not subject to calibration in previous studies on modelling green roofs in EPA SWMM (i.e., [24]).

Field Campaign
The field campaign served to comprehensively examine some roof greening systems, their production, material use and optimization as well as the effects of the different material arrangements (construction methods) on perennial plant development. The product properties were tested in the laboratory. In the following, the most important results of the entire study are summarized, because they have not been published in this form so far. This makes it possible for later simulations to be carried out with other models that can process further parameters. Details can be found in the supplement. Furthermore, the daily runoff measurements are available online (see data availability) just as the meteorological data (www.muk.uni-hannover.de/258.html). For the simulations in this paper, only some parts of the original studies were used (e.g., water storage, daily runoff volumes). From 12 construction methods number #10 was not considered in the simulations because it was tested as a special construction method and was not designed with conventional roof greening substrate according to the FLL (Forschungsgesellschaft Landschaftsentwicklung Landschaftsbau e.V.) Roof Greening Guidelines [25].
The field campaign was carried out from July 2002 to the end of August 2006, to determine the vegetation characteristics and the average annual runoff coefficients of thinlayer green roof structures. The campaign was led by Professor Dr. Hans Joachim Liesecke ( †) in Hanover-Herrenhausen, at the Institute for Landscape Architecture, at the Department of technical-constructive basics of open space planning (Professor Dipl. Ing. Gilbert Lösken). The campaign examined the properties of the materials used, load assumptions and water storage of the construction methods, the development of the vegetation, the annual water retention and the determination of annual runoff figures. Results compiled up until 2004 have been published by Liesecke [26,27].
The average temperature in Hannover-Herrenhausen is 10.3 °C and the average amounts to 636 mm a −1 . The climate corresponds to the cfb climates according to the Köppen-Geiger climate classification [28], i.e., warm-temperate humid climates with cool seasons with mean temperatures of the coldest month between +18° and −3 °C and the mean temperature of the warmest summer month is below 22 °C (summer-warm climates).
The field campaign involved a total of 12 superstructures in triple repetition ( Figure  1 and Supplementary Figure S1). The individual test plots were 2 m × 2 m with a slope of 2% and a drainage opening in the middle of the lowest point of the slope. The outflowing water was collected in non-weighable lysimeters (rain barrels) that were read and emptied at 8 A.M. every day. On weekends, readings were taken the following workday. Figure 1. Setup of the field campaign in Hanover-Herrenhausen with 36 green roof test plots in total (12 superstructures used). The picture shows rows 2-6. Runoff from green roofs is collected in rain barrels.
The rain barrels have a readable, transparent plastic tube on the outside that indicates the water level (principle "communicating tube") as shown in Figure 2. They were previously calibrated in the laboratory. The marks corresponding to 1 liter filling level, approx. 0.25 mm per 4 m 2 test field. The maximum capacity is 17 mm, i.e., 68 liters. Further runoff leads to overflowing of the rain barrels. For maintenance, a build-up-saturating irrigation and a partial moistening in an initial dry period (2002)(2003) as well as an annual fertilization with 30 g m −2 of a complex long-term fertilizer were administered. Therefore, the undisturbed period from 2004 to 2005 was chosen for the modelling experiments in this study. Figure 2. Transparent plastic tube with water level. Scale was calibrated in the laboratory on the 4 m 2 test field (4 L = 1 mm; accuracy (± 2 mL).

Laboratory Tests
For the substrates used, the grain size distribution, volume weight, maximum water capacity, total pore volume, air volume at maximum water capacity and at pF 1.8, pH, organic material, carbonate content and salt content were determined (Supplementary Table S1 to Table S7 and Figure S2). The investigations of the vegetation mats, the protection, storage and drainage mats were carried out on three sections of 40 cm × 40 cm each and the results were extrapolated to 1 m 2 .
Supplementary Tables S8 and S9 summarize the materials used, the respective thickness, the load assumptions dry and at maximum water capacity, as well as the water storage capacity of the construction methods. Due to the thickness and fine-grained composition of the substrate, the vegetation mat has a relatively high water storage capacity. With fleeces, water storage increases significantly while increasing thickness and dry weight. In relation to their thickness, both storage fleeces have a considerably higher water storage capacity than the vegetation mat. The water storage of the drainage and filter mat is very low and insignificant due to its function and structure. In the present order, the thickness of the structure increases essentially from about 3 cm to about 10 cm, which, with some deviations, also results in higher load assumptions and higher water storage.  Figure  S3) stabilized in the years 2005 and 2006. Except for variants 1 and 2, which have gaps and a total coverage of 90%, it is complete for all variants. The total cover is determined exclusively by the succulents (except variant #1). In summary, the transformation of the herb aspect present at the beginning of the experiment into a sedum-determined aspect was completed by 2006.

Perennial Plant Development
Although Sedum album has become the dominant species, the visual aspect is significantly influenced by the variety of other species with different structure and growth height. Supplementary Tables S13 and S14 compare the results obtained for the five Sedum types during the classification days in July and for the experimental years 2002 to 2006. This shows that the intended transformation of the Sedum stand with herbs into a Sedum stand was already achieved after one year and the further development took place in the Sedum group (see Supplementary Tables S13 and S14). The goal to set the foundation for a permanent Sedum population by cultivating the vegetation mat was achieved. The intensity of rooting, the stock density and the shear strength were enhanced by the addition of herbs. While the herbs decreased or completely disappeared after placing the vegetation mat, further development towards a fully developed Sedum population started. This completed development, suggesting a constant vegetation layer, is important for the modelling experiment of the current study, in order to avoid transient vegetation characteristics throughout the simulation period.

Runoff
The runoff was recorded in the period from August 16, 2002 to August 31, 2006. In winter, when temperatures reached freezing point and below, measurements were interrupted, in order to avoid damaging the measuring equipment by ice formation. The downtimes were excluded from further analyses. Based on the precipitation measurements, determined by the nearby meteorological measuring station Hannover-Herrenhausen, operated by the Institute of Meteorology and Climatology (approx. 200 m distance from the study site) and the recorded runoff quantities, the percentage of water retention and the annual runoff coefficient were determined (see Table A1 in the Appendix).
To illustrate the influence of weather conditions, the percentage of water retention for different seasons can be summarized as follows: • Summer period (from 1st of May to 30th of Sep) between 58% and 72%, subject to a decrease by 1/3 in cooler weather in • Spring/autumn (from 16th of Mar to 30th of Apr and 01st of Oct to 15th of Nov) between 41% and 49%.

•
Winter (from 16th of Nov to 15th of Mar) between 5% and 12% On average, over the year, a water retention of between 42% and 51% was determined.

Water Balance Simulation
In order to predict the water balance of green roofs, the Storm Water Management Model EPA SWMM (version 5.1.014) [29] software was used. This model, which is opensource and freely available, is a rainfall-runoff model, widely used to simulate and estimate the hydrological response of urban areas. The LID editor implemented in EPA SWMM makes it possible to define different horizontal layers for green infrastructure and nature-based solutions (NBS) in a sub-catchment. The aim of implementing these tools is to infiltrate, store and evaporate rainwater runoff close to the source, at the expense of runoff. The module "green roof" makes it possible to describe infiltration and evaporation from green roofs in the model. In terms of their area, green roofs are represented as a fraction of their sub-catchments. The empirical Hargreaves approach is the standard method in EPA SWMM to compute ET, based on temperature recordings. It is expressed by the following formula and gives ET0 in mm d− 1 [30]: where is the theoretical maximum extraterrestrial radiation for the considered day of year (mm d −1 ), TC is the mean daily temperature (7-day running average) [°C], and are the maximum and minimum daily temperatures (°C), respectively. This representation of the Hargreaves model (Equation (1)) has fixed parameters, which cannot be modified by the user.
Since this approach does not provide any vegetation-specific parameterization, this study employs a different approach.

Physically-Based ET Simulation
The Penman-Monteith approach describes the diffusion of energy by sensible and latent heat in a network of resistances. It is derived from first principles considering the radiation balance and energy balance of the surface, and therefore it is a physically-based ET model. The Penman-Monteith equation can be written in terms of the latent heat flux LE (W m −2 ) [31]: Where ∆ is the slope of the saturation vapor-pressure curve (hPa °C −1 ), is the net ra- , is the mean air density at con- is the bulk surface resistance (s m −1 ) and the number of seconds within a time step (daily in this case) [31].
Because it is a physically-based approach, the climatological boundary conditions (i.e., time series of temperature, humidity, wind speed and downwelling shortwave radiation) and vegetation parameters have to be known and understood. Evaporation, using the Penman-Monteith approach, is calculated externally using the open source Python tool "PyETo" [32]. This includes a standardized formula from Allen et al. [33] for ET, from short grass as reference (i.e., the crop reference ET). In this study, the simplified formula for the crop reference ET is adopted, in order to allow adjustments to individual parameters that are fixed in case of the crop reference ET. This, in turn, enables a tailored parameterization of different types of vegetation. In this context, the surface resistance is an important plant parameter that governs ET. It is the resistance of vapor flow from within plant leaves and from beneath the soil surface. It includes both the canopy (leaf area index) and the stomatal resistance, which describes the average resistance of an individual leaf. This resistance is crop specific and differs among crop varieties and crop management.
Combined with the weather data from the test site, suitable ET time series are calculated for arbitrary values of surface resistance [33].

Coupled Modelling Approach
For the simulation, hourly recordings of precipitation are used to force EPA SWMM (Section 2.2.1). As EPA SWMM requires daily values for ET, the calculated ET, using the Penman-Monteith approach (Section 2.2.2) was introduced as an external daily time series. The results obtained from the simulations are aggregated to achieve daily totals, in order to make them comparable with the daily observations of runoff. As a final step, a comparison between simulated and observed runoff is performed to estimate the predictive skill of the calibration schemes. The adjustable parameters that were used to perform the initial model runs are listed in Table A2 in the Appendix (under column "default"). The properties of the coupled model set up are described in Table 1:

Property Set up
Sub-catchment A sub-catchment represents a roof of 4 m 2 (0.0004 ha) that is 100% impermeable (covered with PVC sheets; N-Imperv PVC: 0.01 [34]), has a slope of 2% and drains entirely to an outfall. It is covered by different low impact development (LIDs) (green roofs).

LID
The sub-catchment is covered by the LID type "Green Roof". Different LID controls represent the different structures of the test setups. As long as the parameters are unknown, they were estimated and calibrated.

Precipitation
Hourly precipitation data are available. The data for the calibration period ( The Penman-Monteith approach is applied to calculate daily evapotranspiration rates. The data were introduced as external time series to EPA SWMM. The daily evaporation data for the calibration period (

Calibration and Validation
Since various parameters in the coupled modelling approach are unknown, they are subject to the calibration approach described below. To determine the accuracy of the model outside the calibration period, an independent validation is also considered, following a split sample test [35]:  Table 2 both periods differ in terms of weather: the first season (calibration) has a lower average temperature and a higher rainfall total compared to the second season (validation). Thus, ET is higher in the validation period.
The current research employs the open source Python tool spotpy [36] for model calibration. It is a statistical parameter optimization tool that can be used for calibration, uncertainty and sensitivity analysis. In order to customize spotpy to any arbitrary hydrological model, it is necessary to develop an interface to the externally used model. The interface for SWMM was adopted from [37] in order to read different files of ET, based on the sampled resistance values. In addition, the hourly simulated runoff is aggregated to daily data (daily totals from 8 A.M. to 8 A.M. of the subsequent day) after each iteration of spotpy, as observation data are provided in daily resolution.
The Shuffled Complex Evolution (SCE-UA) global optimization algorithm [38] is used in this study. It was found to be a robust algorithm for hydrological model calibration and is widely used in hydrological applications [39]. As an objective function, the Kling-Gupta model efficiency (KGE) is employed [40]. A wise use of the objective function is important because of its high sensitivity. For this case, the KGE seems to be a robust function, as it considers different types of model errors (error in the mean, the variability, and the dynamics) [41]. In addition to KGE, the Nash-Sutcliffe Efficiency NSE [42] and the root mean square error (RMSE) are also computed for each run. The simulated results are compared to the mean of observed values. For all performance measures, days without observations (e.g., weekends) were dropped and those with a runoff higher than 17 mm (This applies to 7th and 10th of May 2004) were set to 17 mm as this was the measurement limit (i.e., the capacity of the rain barrels that collect the runoff). This modification for computing the objective function (i.e., KGE) is necessary to give less emphasis on uncertain peak values in the calibration scheme. However, all display items in this paper (see, Section 3) show the original time series of the model without this modification. Table A2 in the Appendix lists all parameters that have been adjusted in the framework of the calibration procedure. Special emphasis is put on the surface resistance, since we intend to find a robust parametrization valid for all test sites, which show the same vegetation characteristics. This is motivated by the fact that, despite the different structures, the different green roof variants underwent a similar evolvement of the vegetation coverage as evidenced by the uniform vegetation development in the field campaign (Section 2.1). Hence, it is assumed that vegetation characteristics in our modelling experiments are uniform in space across all test sites. In turn, calibration involves the following steps: • 1st iteration: Run spotpy, applied to the coupled modelling approach with 5 repetitions for each test plot (to see how the stochastic nature of sampling affects results), in order to find a unique representative surface resistance value. Various evapotranspiration time series are calculated (with a surface resistance range from 60 to 300 s m −1 ). These files, or more specifically the file path in the EPA SWMM's input file, are available as adjustable parameters to the calibration scheme. For this purpose, an existing EPA SWMM model driver in spotpy has been adjusted (https://github.com/mmmatthew/swmm_calibration (last access 24 Nov 2020)). This way, the surface resistance is estimated directly as a yearly average, suggesting that a seasonal variation is neglected.

Global Sensitivity Analysis
To consider parameters with a high impact on the model results, an additional global sensitivity analysis, also using the Python tool spotpy [36], complements the previous calibration procedure. The algorithm "Extended FAST" based on [43] is selected for this purpose. Besides "Sobol", this is the only method that takes parameter interaction into account [44]. The method is based on the Fourier amplitude sensitivity test (FAST) and allows quantifying the total contribution of each input factor to the variance of the output to be calculated. These methods are complex and to keep the explanations concise, we will only provide a brief overview of how Extended FAST works in principle. At first, the algorithm creates matrices for parameter sampling [45]. Then a bootstrap technique is applied to the matrices which relies on replacement of elements. Based on that, real model parameters are derived from sampled values and the corresponding model output is then used to calculate the parameter sensitivity [44] based on the estimator proposed by Jansen [46].
The main advantages of "Extended FAST" are its robustness, especially with small sample sizes, and its computational efficiency [43]. Important for this algorithm is the small number of iterations in order to obtain reliable information about the parameters. It can be calculated according to [44]. Finally, the algorithm provides a sensitivity index for each parameter. The sensitivity index is a reliable measure to rank parameters according to their relative importance [44]. Hence, comparing the sensitivity index among all parameters allows for ranking the importance of parameters.
First, the literature (default) values for common value ranges of the EPA SWMM parameters are defined as the initial parameter set (see "default" column in Table A2 in the Appendix). This parameter set includes those parameters that are adjusted in each iteration of the algorithm and generally unknown (e.g., the conductivity slope). Indeed, known characteristics of the field campaign, e.g., the dimension of each green roof variant and the layer structure remain unchanged during parameter sampling. As a simplificationexcept for the surface resistance-the Penman-Monteith parameters that describe vegetation characteristics remain fixed (Albedo: 0.2, crop height: 0.12 m, height wind exponent a: 0.5, rs: 78 s m− 1 ). This is assumed to be feasible, since the surface resistance is the most sensitive parameter for plants subject to limited water availability [47], which applies in particular to green roofs. To avoid possible negative impacts arising from the selected value ranges for each parameter, each range is extended by ±20 and ±50%. Since "Extended FAST" also involves KGE as an objective function, it serves as another independent approach for model calibration.

Representative Surface Resistance
In the first step of the calibration procedure, the identification of a unique, representative value of surface resistance is addressed. The SCE-UA algorithm in spotpy is applied five times to each test site (green roof variant). Figure 3 shows a histogram including all surface resistance values and box plots, visualizing the best estimate of surface resistance for each green roof variant. The results indicate that SCE-UA yields consistent estimates of surface resistance across all test sites. Since we seek to understand surface resistance as a vegetation property, which is assumed to be homogeneous in all experiments, an average value is computed (first iteration) and applied in each following run of SCE-UA algorithm in spotpy (second iteration).

Model Runs
The surface resistance found to be representative for all test plots is then used as a fixed value in the 2nd iteration of the calibration procedure (see Section 2.

Model Run with Default Parameters
For the uncalibrated model, the default parameter values from the literature review were applied [21,29]. Simulation and observation show similar patterns in terms of timing and variability of the time series, but the simulation overestimates the runoff. Results for green roof variant #6 are visualized in Figure 4. The value for the objective function, more specifically, the Kling-Gupta efficiency (KGE), is 0.57.

Calibrated Model
For the calibrated model, the calibrated LID parameters for the surface, soil and drainage mat (second iteration) are considered, together with the surface resistance value from the first iteration (Table A2 in the Appendix). Modelled and observed time series coincide well, even though the simulation sometimes slightly underestimates observed runoff (except for days when observed runoff exceed the capacity of the rain barrel of 17 mm). Results for green roof variant #6 are presented in Figure 5. The Kling-Gupta efficiency (KGE) reaches 0.88, which underlines the added value of calibration, compared to the uncalibrated model (KGE = 0.57).   The good model performance is also reflected by the statistical summary in Table 3 (first column). Both observed and modeled average runoff compare well. However, the longest dry period (i.e., consecutive days with zero runoff) is overestimated by the model, which is also obvious from comparing the time series in Figure 5. The runoff event in late June 2004 is not reproduced by the model. However, the majority of events is captured well in the simulations.

Model with the Integrated Hargreaves Method for ET
For reasons of comparison, the model was also calibrated and run using the original Hargreaves method, which is already implemented in EPA SWMM. Similar to the previous runs, runoff time series for green roof variant #6, obtained from the EPA SWMM setup with the Hargreaves method, are presented in Figure 7. The simulation does not only overestimate total runoff but also runoff variability. The Kling-Gupta efficiency (KGE) is −0.35, suggesting that the Penman-Monteith approach clearly outperforms the Hargreaves method in case of our test site.

Validation
The validation of the calibrated model (Section 3.2.2.) shows similar patterns between observation and simulation, but the simulation sometimes underestimates the runoff. Results for green roof variant #6 are presented in Figure 8. The Kling-Gupta efficiency (KGE) is in this case equal to 0.84, which is a value close to the corresponding value achieved for model calibration (KGE = 0.88). Table 3 also confirms that runoff characteristics are presented well by the model, even though the longest dry period is also overestimated in the validation period (the event in late July is not captured by the model). However, it is assumed that the model shows predictive skill outside the calibration period, since the performance measures of the validation period reflect a similar accuracy as those values achieved in the calibration period (Table A2 in the Appendix)

Global Sensitivity Analysis
The results of the global sensitivity analysis are presented in terms of the total sensitivity index (with a maximum of one, see Section 2.4). According to the sensitivity analysis, the model is most sensitive towards the following parameters: surface resistance and berm height. The computed sensitivity index is smaller for soil thickness, soil porosity, drainage mat thickness, and drainage mat void fraction, indicating that these parameters are less important in our experimental setup (Figure 9). This allows for re-running the calibration scheme with a smaller number of adjustable parameters. As the soil thickness was given by the field campaign, this value was not calibrated and another simulation is carried out calibrating only the remainder of the most sensitive parameters (i.e., berm height, soil porosity and drainage mat void fraction). The calibrated surface resistance of 78 s m− 1 is also used here. By calibrating only these sensitive parameters, the simulation yields a Kling-Gupta Efficiency (KGE) of 0.79 (run not shown here), which is consistent with the calibration utilizing SCE-UA (KGE = 0.88).

Synthesis of Model Runs
To summarize the different model runs described previously, the total runoff and the Kling-Gupta efficiencies (KGE) obtained in each model are compiled in Figure 10. Especially when using the Hargreaves method implemented in EPA SWMM, the total runoff is clearly overestimated while ET is underestimated. The calibrated model and the model calibrated with only sensitive parameters (referred to as "sensitivity" run) show the best efficiencies (KGE 0.88 and 0.79), which is confirmed by similar model skill values achieved for the validation period (see Table A2 in the Appendix). In contrast, the efficiencies using the model with the Hargreaves method are −0.35 and −0.52, for the calibration and validation periods, respectively. This mismatch, expressed by the negative KGE values, is mostly determined by the too low values of ET computed by the Hargreaves method, which in turn exaggerates the hydrological response in terms of runoff, which is more than twice that of the observed value. While all previous results pertain to green roof variant #6, Figure 11 compiles corresponding metrics for the entire set of eleven green roof variants (see Section 2.1 for reference to the green roof variants). Besides total runoff in the calibration period, the Kling-Gupta efficiencies (KGE) of calibration and validation periods are compiled likewise. The performance measures for all green roof variants are listed in Table A3 in the Appendix. Except for green roof variants #1, #7 and #12, the total runoff is always slightly underestimated (percent bias (PBIAS) = 6% on average (Please note that positive PBIAS values indicate underestimation of the model, while negative values suggest a model overestimating observed values [48].)). On the one hand, this underestimation is acceptable for streamflow simulations (very good performance if PBIAS ≤ ±10% [48]). On the other hand, the model has been calibrated against observed runoff while maximizing KGE, which considers not only the model bias but also the correct timing and magnitude of peaks, thus representing a compromise of several aspects of calibration. In green roof variant #12, however, runoff is overestimated by 20%, in green roof variant #7 by 26%. For all green roof variants, the Kling-Gupta efficiencies (KGE) vary from 0.77 to 0.91 for the calibration period and from 0.7 to 0.89 for the validation period, respectively. Only the KGE for green roof variant #12 is significantly lower with 0.21. It is also noted that the Kling-Gupta efficiencies of calibration and validation periods are similar (with exception of green roof variant #12), which justifies using the calibrated parameters outside the calibration period.

Discussion
To investigate an improved representation of the hydrological response of green roofs, a combined model using the Penman-Monteith approach and EPA SWMM was set up. Detailed data from a field campaign served as a basis and the measurements were used to calibrate and validate the model.
A multiple parameter calibration of LID parameters (surface, soil and drainage) shows distinct improvement to the default model using literature values (increase in KGE from 0.57 to 0.88). The calibration of the surface resistance seems to be the best method to take the vegetation feedback into account, since it also helps to improve the representation of the water balance closure of the model, through more representative ET estimates. ET is computed with daily temporal resolution, which is consistent with the expected format of time series and can be used to force EPA SWMM with external ET simulations. The model, however, was integrated using hourly time steps, which is considered a reasonable approach, since EPA SWMM is mostly applied to urban drainage studies that involve subdaily time steps due to short times of concentration in urban areas [49].
The added value of the coupled model, involving both EPA SWMM and the Penman-Monteith approach, is obvious, given that the uncalibrated model with KGE = 0.57 even outperforms a calibrated EPA SWMM setup with Hargreaves method to compute ET. The latter model, shows a poor model performance due to ET underestimates, yielding negative KGE values. As long as the coefficients of the Hargreaves equations in EPA SWMM cannot be adjusted, it is not possible to calibrate ET simulations in the original approach of EPA SWMM, employing the Hargreaves equation [50]. An external calculation with calibrated coefficients might yield better results, but this goes beyond the scope of the research question, since this study is not about comparing different evapotranspiration methods, but rather about showing an improvement in the modelling of runoff from green roofs with EPA SWMM.
Extreme events and rainfall events after long dry periods are represented with lower accuracy in our modelling experiments. Since the interest of the study is on the vegetation feedback rather than on extreme events, this inaccuracy is accepted here. In addition, observed values are subject to a maximum value of 17 mm per day due to the limited capacity of the rain barrels used to collect runoff from the test plots, rendering the accurate representation of extreme events impossible.
A global sensitivity analysis shows that in case of our modelling experiments, only a few soil parameters are sensitive in EPA SWMM. One of them is the berm height that describes how much water can pond, therefore suggesting an additional storage. The construction thickness and porosity also describe the storage volume. Accordingly, the model reacts in a particularly sensitive way, as a consequence of changing water storage. Other descriptive parameters are less influential and might be neglected in the calibration according to the investigations. Here, we describe the possibility to consider a simplification of the parameterization of EPA SWMM, suggesting a lower number of parameters that need to be adjusted in the calibration procedure. The coupled model is particularly sensitive to the surface resistance. This selected vegetation parameter influences ET significantly and thus governs the simulation of the water balance of green roofs. Due to its high sensitivity, it is recommended to focus on this parameter, especially in cases of different vegetation.
Regarding the parameters analyzed in the framework of the automatic calibration ( Table A2 in the Appendix and Supplementary Table S17), similarities are found among the green roof configurations, which can be related to their characteristics (Supplementary  Tables S8 and S9). For instance, the porosity is higher for green roof variants #7, #8, and #9 when compared to green roof variants #1 to #6. This sampling of parameters is feasible, given that the green roof variants #7, #8, and #9 consist of a single-layer substrate with higher porosity than the geotextiles. Likewise, for green roof variants #3 and #5, the drainage mat void fraction is higher than for #1, #2, #4, and #6, which is supported by the construction details (e.g., green roof variants #5 have two drainage layers instead of one). The drainage mat of green roof variants #6, however, shows a much higher density of the material, suggesting that the void fraction is much lower (Supplementary Table S17). All these parameters were found through optimization, suggesting that they are also effective parameters rather than physical constants [51]. For instance, the berm height is a calibrated parameter loosely linked to the physical characteristics and thus representing model uncertainty and scaling issues as well. However, we can identify some groups of green roof variants that show similar parameters (i.e., porosity, field capacity), e.g., green roof variants #1 to #6 with geotextiles, green roof variants #7 to #9 which represent single layer substrates and green roof variants #11 and #12 that are multi-layer substrates, whereby #12 has the highest vertical extent.
The validation for several green roof variants confirms the accuracy of the model, suggesting that it can be applied to other time periods as well as to other green roof variants. Moreover, the results confirm that a unique, representative estimate of the surface resistance parameter highlights that our model setup reflects a robust model parameterization for vegetation, even though we neglect important mechanisms in our model setup, e.g., the CAM metabolism of Sedum and seasonality of surface resistance [22]. Still, an open question refers to the lower model performance found for green roof variant #12 (higher storage). Given that the setup clearly differs from the remaining green roofs, several reasons could in principle hold responsible for this mismatch. To gain insight into this question, another objective function or a different optimization might help to see if better model results are possible when applying EPA SWMM to green roof variant #12. However, this goes beyond the scope of this study.

Conclusions
In this research article we present a coupled modelling approach comprising EPA SWMM and PyETo. In this method, the hydrological model EPA SWMM is complemented by a more comprehensive representation of ET, which is computed using the Penman-Monteith approach. Based on numerous green roof test plots, which are represented in the coupled model, we show that: (i) the model confirms that shallow layered extensive green roofs show good retention capabilities with a minimum roof load, as witnessed by the lower average value and longer dry spell period observed for runoff (compared to rainfall, see Table 1 and  Table 3); (ii) the coupled modelling approach outperforms the standard EPA SWMM model with Hargreaves ET computation, even without any further calibration of the ET component; (iii) a robust parametrization of the vegetation (or more specifically surface resistance) is possible for different green roof configurations with similar vegetation cover; (iv) only a subset of LID parameters in EPA SWMM are sensitive for continuous simulations, namely surface resistance, berm height, soil porosity and the drainage mat void fraction.
Since these findings mostly apply to 11 independent green roof tests sites (three repetitions for each), for which we identified a unique representative surface resistance parameter governing ET, the parameterization of the coupled modelling approach is viewed as being robust in terms of transferability in space (different test plots) and time (independent validation period).
The coupled modeling approach allows for the determination of the water balance components on green roofs by applying a more accurate continuous simulation, which shows major advantages for design flood estimation when compared to event-based approaches [52], for example, by the implicit incorporation of antecedent soil moisture conditions [16]. However, long and continuous rainfall series are not always available and present a low resolution. Thus, stochastic rainfall models could be used in combination with the coupled model to overcome this problem [15,53] and perform continuous simulations in urban hydrology.
However, it is worth mentioning that several uncertainties have so far not been considered. Among other aspects, these include: (i) the design of the existing field campaign was not tailored to provide a basis for modelling studies. Hence, hydrological quantities that might complement constraining the model for dry and wet periods, such as soil moisture measurements, have not been addressed; (ii) The non-stationarity of vegetation characteristics was not considered in our study.
For instance, Sedum tends to recover within one week after drought stress [54]. Moreover, surface resistance might undergo seasonal changes and might be subject to changes as a response to drought stress, as idealized irrigation experiments after drought might suggest [55]; (iii) Equifinality (i.e., different parameter sets yield similar results or accuracy in terms of objective function) is still a source of uncertainty in studies aimed at parameter identification [56]. Some parameters are effective parameters, suggesting that their original physical meaning is not always given due to model uncertainty and scaling issues [51]. However, since the field campaign comprises three independent samples for each green roof configuration-the differences between which are small (working with an average for each configuration is feasible)-our results are considered robust in terms of parameterization. Even though some uncertainties arise from the limit of 17 mm per day, reflecting the maximum runoff observed on a specific day, there is only limited evidence on the predictive skill for extreme events, especially since rainfall observations from the weather station are also subject to uncertainties.
Apart from these limitations, our coupled model setup is a substantial step forward in modelling green roofs in urban drainage studies, especially since the components involved here (EPA SWMM, PyETo, and spotpy) are available free of charge and coupling the models is straightforward-without the need to modify the EPA SWMM code-as it is possible to force the model with external ET files. The full set of meteorological quantities needed to utilize Penman-Monteith (instead of temperature only) is also available from re-analysis data, e.g., ERA5 (European Centre for Medium-Range Weather Forecasts. European Reanalysis 5), making it possible to drive the coupled modelling approach for arbitrary sites. Even with limited information on the detailed characteristics of the green roof, we demonstrate that the coupled modelling, even without any further calibration, is capable of outperforming the original EPA SWMM ET simulation, based on the Hargreaves equations. Moreover, some practical recommendations on constraining the number of adjustable parameters are given.
As an outlook for future research, it is recommended to better address the non-stationarity of vegetation characteristics of green roofs. In the present study, we demonstrate that the assumption of a constant surface resistance value across different green roof variants with similar vegetation coverage is valid. This is an important finding, since the surface resistance is hence transferable in space. However, a seasonally varying surface resistance value might help to better predict season onsets and drought stress recovery (e.g., by involving weather data to adjust parameters depending on the season [57]). Thus, a time-varying parameterization might complement the spatially uniform parameterization of vegetation across the test sites. In addition to a constant surface resistance applied to all sites (as studied here), a uniform method that describes the temporal variation of the surface resistance parameter should be investigated in the future.
For the design of field campaigns for green roofs, we also recommend sampling hydrological quantities other than runoff, to better constrain the parameter optimization of the hydrological model, through employing a multi-objective calibration. Moreover, field campaigns with similar characteristics should be carried out in different climate zones, to investigate the effect of climate on parameter optimization.

Supplementary Materials:
The following are available online at www.mdpi.com/2306-5338/8/1/12/s1, Figure S1: Arrangement of the 12 green roof variants used in the field campaign, Figure S2: Grain-size distribution curves of the materials used as substrates, Figure S3 Table S1: Grain size distribution of the materials used as substrates, Table S2: Volume weight (dry, at delivery and max. water) and water content at delivery after compression of the different materials used as substrates, Table S3: Maximum water capacity and infiltration rate of the materials used as substrates, Table S4: Total pore volume, air volume and volume reduction after compression of the materials used as substrates,   Acknowledgments: We would like to cordially thank Professor Hans Joachim Liesecke ( †) who conducted the field work in the years 2002-2006. We gratefully acknowledge that he allowed other researchers at the Institute of Landscape Architecture to proceed with further research based on the data of the field campaign. Without the well archived data and this permission, this work would not have been possible. Last but not least, we also thank the reviewers for their helpful comments and Larissa van der Laan for proof reading.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Construction technique and material description (including thickness of each substrate, total thickness and water retention) for the different green roof variants and runoff coefficients obtained.