EHSMu: a New Ecohydrological Streamflow Model to Estimate Runoff in Urban Areas

A conceptual lumped ecohydrological streamflow model (EHSMu) is presented as a promising tool to simulate runoff in urban catchments. The model, based on the interaction between a soil bucket and two linear reservoirs, enables also evapotranspiration and aquifer recharge to be estimated. Notwithstanding its minimalism, EHSMu describes interactions among soil moisture dynamics, hydrological fluxes and ecological processes. The model was calibrated and validated within two densely urbanized sub-basins in Charlotte (US). A Monte Carlo procedure is used to investigate the efficiency of random sets of 8 model parameters. Results show the high model performance (NSE = 0.72). The influence of land use change is evaluated, by varying the imperviousness and crop coefficients. Synthetic experiments show that increasing urbanization triggers a linear decrease in evapotranspiration and aquifer recharge, while it increases the fast runoff. An opposite response is achieved by installing vegetation with higher potential evapotranspiration, which would contribute to the actual evapotranspiration making up 50–55% of the total water balance.

Urban areas, however, are complex systems and, due to the high spatial variability and heterogeneity, they are difficult to model (Fletcher et al. 2013;Salvadore et al. 2015). Several models have been proposed to estimate the urban runoff, trying to account for spatial and temporal variability of rainfall and catchment (Cristiano et al. 2017). Most of these models focus on the identification and representation of the stormwater drainage network (Zoppou 2001), and they combine it with surface dynamics, in order to estimate the amount of rainfall that reaches the sewer system (Bermúdez et al. 2018;Rubinato et al. 2013;Versini et al. 2016). Semi-and fully-distributed physically-based models have been developed for many urban areas, showing good performance in the runoff estimation (Aronica and Cannarozzo 2000;Pina et al. 2016). However, in most of the cases, these models require a long computational time and also spatially and temporally distributed input data (Cristiano et al. 2017). The choice of the most suitable model is constrained by the data availability and final goals. Berne et al. (2004) suggested a criterion based on the relation between available rainfall resolution and catchment area: if the rainfall resolution is larger than the catchment area, a conceptual modelling approach is suggested, otherwise, a distributed approach could be more efficient.
In the case of a lack of distributed data or a necessity of fast simulations, conceptual models would be preferable. Sarma et al. (1973) compared different conceptual rainfall-runoff model performances in urban areas and showed that, for small basins, the single linear reservoir leads to good runoff estimation. These models, however, do not account for other components, such as evapotranspiration, subsurface flow and aquifer recharge.
The combination of hydrological and ecological elements within ecohydrological models has been explored in depth for natural environments and the importance of taking into account both components has been largely underlined (Botter et al. 2007;Laio et al. 2001;Rodriguez-Iturbe et al. 1999). In urban areas, this aspect has recently become more popular: to guarantee a sustainable management of soil and water resources, a better understanding of the whole urban ecosystem is required (Marchionni et al. 2019). Investigating hydrological and ecological processes and their interactions provides a strong support for urban design and water management (Li 2012). Recently, Meili et al. (2020) proposed an urban ecohydrological model, Urban Tethys-Chloris (UT&C), which couples energy and water balance and combines principles of ecosystem modelling with an urban canopy scheme. UT&C explicitly accounts for the biophysical and ecophysiological characteristics of roof vegetation, ground vegetation and urban trees in simulating water and energy fluxes.
Many ecohydrological urban models focus on small scale, investigating the interactions among impervious surfaces, vegetation, and atmosphere in urban areas (Grimmond et al. 2010;Marchionni et al. 2019;Vico et al. 2014). At the moment, however, only a few ecohydrological models have been developed at catchment scale (Liu et al. 2013;Shields and Tague 2015) and the integration of ecohydrological processes to define the runoff generation at urban scale needs to be further investigated.
In this context, we introduce a new parsimonious lumped conceptual ecohydrological streamflow model for urban areas (EHSMu) to estimate the runoff generation in an urban environment, taking into account water-soil dynamics, vegetation types, evapotranspiration fluxes and aquifer recharge. The model, described in Section 2, combines a soil bucket and two linear reservoirs, to simulate the runoff at the catchment outlet and to estimate evapotranspiration and aquifer recharge. Two catchments in the Charlotte metropolitan areas (US) were chosen as case studies to calibrate and validate the model, using rainfall and potential evapotranspiration at hourly scale as input for the model. Rainfall and streamflow data have been used in previous studies to investigate the urban hydrological response in relation to storm characteristics and flood frequencies (Cristiano et al. 2019;ten Veldhuis et al. 2018;Wright et al. 2014), and they have shown good reliability, with high rainfall-runoff correlation. Results are discussed in Section 3, where the impact of land use and land cover changes is investigated through a sensitivity analysis.

Model Description
The EHSMu model is based on a simple schematization of physical elements within an urban catchment and on several hypotheses regarding water fluxes. Urban areas are described as a lumped ensemble of impervious (roads, roofs) and pervious (gardens, parks) surfaces. When rain falls over impermeable surfaces, it is supposed to be directly transferred towards the drainage system, in a time depending on the catchment area and drainage network geometry and materials (Al-Janabi et al. 2019). The rain falling over permeable areas drives soil moisture dynamics, which frame evapotranspiration rates, excess runoff generation and leakage occurrences. The latter are supposed to partially feed deep aquifers, thus introducing a new element to the water balance.
The hydrological scheme of the urban basin, implemented in EHSMu and presented in Fig. 1, is defined as the combination of three interconnected elements: a soil bucket below the permeable area and two linear reservoirs, named fast and slow reservoir respectively. The fast reservoir collects both the runoff from impermeable areas (Q imp ) and runoff generated by the saturation excess in the permeable areas (Q exc ). The slow reservoir receives as input a fraction c 1 of the percolation leakage L originated from the soil bucket below the permeable area, while the component (1 − c 1 )L infiltrates the aquifer. The total outflow runoff Q out is defined as the combination of the fluxes deriving from the fast and slow reservoir.

Soil Water Balance
Hourly rainfall depth R (mm/h) is considered as the main climatic forcing and it is used as input for the model as a time series spatially averaged over the catchment. As mentioned in the Model description, rainfall is partitioned between an impermeable fraction c 0 and a permeable fraction (1 − c 0 ), where water fluxes between atmosphere and soil bucket may occur through rainfall infiltration and evapotranspiration processes, affecting the balance of soil moisture s. Interception by vegetation is assumed to be negligible, since the contribution over a long time series in urban areas is often quite small compared to the other fluxes (Zabret and Šraj 2019). The permeable fraction of soil is characterized by the soil porosity n and by the root depth Z r , whose product nZ r (mm), often referred to as "active soil depth", is a parameter of the model.
The assessment of relative soil moisture variations Δs in the permeable part of the soil bucket during an hourly temporal step is achieved with the numerical solution of the following water balance equation carried out through a forward finite differences method: where soil moisture s is defined as the dimensionless ratio between the soil water volume and the soil control volume, Q exc (mm) is the runoff generated by saturation excess and ET (mm) and L (mm) represent the soil water losses due to evapotranspiration and leakage, which will be described in the next section. In particular, Q exc is produced when the rainfall exceeds the maximum storage capacity of the permeable soil.

Loss Components: Evapotranspiration and Percolation Leakage
In ecohydrology, evapotranspiration is one of the most critical processes to represent: its interdependency with the soil moisture plays a significant role in the modelling (Asbjornsen et al. 2011) and it has often been modelled as a stepwise function (Laio et al. 2001;Pumo et al. 2008;Viola et al. 2008). Maximum evapotranspiration E max generally occurs when the soil water content is above the soil field capacity s l : in these conditions the plant living processes are guaranteed. The complex relations between evapotranspiration and soil water content are often modelled through a simplified approach by assuming that, as soil moisture decreases to the hygroscopic point s 0 , evapotranspiration linearly decreases from E max to zero (Hale and Orcutt 1987). Under this simplified framework, in the EHSMu model, evapotranspiration is assumed to be governed by the following Eq. 2 as a function of soil water content s.
Following the FAO-56 Penman-Monteith method described in Allen et al. (1998) and widely used in literature, the maximum evapotranspiration rate E max was derived as a function of the Fig. 1 Schematic representation of the hydrological processes implemented in EHSMu reference evapotranspiration ET0. This value depends on several factors, such as solar radiation, air temperature, humidity and wind speed. The maximum evapotranspiration rate can be obtained as the product of the reference evapotranspiration ET0 and the crop coefficient K c , which is a characteristic value of the vegetation type and can be easily derived from agronomy literature (Allen et al. 1998). This approach has been integrated in EHSMu, where the crop coefficient K c is a parameter to be calibrated, due to the high variability in urban areas.
In addition to evapotranspiration, losses derived from percolation leakage are mainly controlled by the hydraulic conductivity, which varies depending on the soil moisture. If the soil water content is below the soil field capacity, the hydraulic conductivity is null and, from that point, we assume that it linearly increases to a maximum in correspondence with the soil saturation (Clapp and Hornberger 1978). This schematization enables the presence of macro porosity and fractures in the soil to be represented, which may lead water directly to deep aquifers. According to this simplified framework, the percolation leakage can be estimated with the following equation: It is worth mentioning that, although the leakage process is here represented as a pulse function, with instantaneous transfer from the soil bucket to the slow reservoir at each hourly time step, the exponential recession response of the reservoir itself mitigates errors that could arise from this approximation (Viola et al. 2014).

Runoff
In EHSMu, the runoff generated at the catchment outlet Q out derives from the combination of the contribution of fast and slow linear reservoirs. This choice was supported by several studies, which largely discussed the number of linear reservoirs required to model properly the catchment outflow (Chow 1964;Jakeman et al. 1990). The runoff from impervious surface, Q imp , and the one from the excess saturation, Q exc , originating from the soil bucket surface, are routed in the fast reservoir, characterized by a residence time 1/K f . The outflow from this reservoir, Q f , is assumed to be a linear function of the water storage h f , and it is defined as: This reservoir represents the quick hydrological response of the basin, where the residence time 1/K f is expected to be close to the basin concentration time, generally less than a day. The second linear reservoir is characterized by a longer residence time 1/K s , and it is fed by a fraction c 1 of the leakage from the soil bucket. The outflow from the slow reservoir Q s can be expressed, again, as a linear function of the water stored h s : Considering that the second reservoir is introduced to describe the slow hydrological response of the catchment, the residence time 1/K s is expected to be in the order of a few days. Hence, the total runoff at the catchment outlet can be written as:

Study Case
The Charlotte metropolitan area (North Carolina, US) was chosen as a case study. In particular, Little Sugar Creek (#507; 111 km 2 ) and Little Hope (#470; 7 km 2 ) were selected to investigate the applicability of the EHSMu on two basins with different extensions (Fig. 2). Land use and cover data have been derived from high-resolution (30 m) gridded datasets, available from the American National Land Cover Dataset NLCD. The investigated basins are highly urbanized, with almost 30% of the area covered with medium or high intensity urbanization (Fig. 2b). Rainfall data were obtained exploiting a network of 72 tipping-bucket type rain gauges, maintained by the United States Geological Survey (USGS), which covers the Charlotte metropolitan area and provides data with a 15-min temporal resolution. Eight stations from the mentioned network (Fig. 2a, red dots) were selected, choosing the ones that have a higher spatial influence and can provide a continuous data set. Streamflow data, measured at the hydrological outlet of each sub-catchment (Fig. 2a, blue square), can be accessed through the USGS website and can be downloaded with a temporal resolution of 15 min. A 15-year time series, from 1st January 2001 to 1st October 2015, was selected for both rainfall and streamflow data and aggregated to an hourly time step to feed the model.
Potential evapotranspiration was estimated with the Thornthwaite equation (Thornthwaite 1948), using monthly average temperature T i , downloaded from the Global Historical Climate Network (GHCN) (Menne et al. 2012) at daily scale. The weather station located at the Charlotte airport was chosen to obtain the temperature data needed to estimate the potential evapotranspiration rate.

Parameter Calibration
The input time series were divided into two periods: the first one used for calibration and the second one for validation. For the smallest basin, #470, the first 40% of the time series was used as the validation period and the remaining 60% as the calibration time series. This peculiar fractioning of hydrologic dataset into calibration and validation periods was motivated by the presence of a period of missing data in runoff records for the basin #470. Since these gaps are approximately after 40% of the considered period, the longer period was used for calibration purposes and the shorter one for validation.
EHSMu requires a set of 8 parameters to represent the hydrological behaviour of the catchment: the active soil depth nZr, the soil moisture triggering the leakage s l and the hygroscopic point s 0 , the crop coefficient K c , the slow and fast reservoir constants K s and K f , the percentage of impervious surface c 0 and the fraction of leakage that contributes to the total runoff c 1 . For each parameter a range of possible values was defined, based on literature and on a priori knowledge. The active soil depth nZr is quite difficult to identify without detailed information from pedological profiles and for this reason a large range is assumed for this parameter, which can vary between 1 mm and 1000 mm. The term s l strongly depends on the soil characteristics and it is assumed to be between 0.3 and 0.95, while the hygroscopic point s 0 is limited between 0 and s l . The crop coefficient K c represents the potential capacity of the plants to sustain the evapotranspiration processes, under optimal weather conditions, ensuring no water stress for the vegetation. This coefficient is assumed to be equal to 1 for standard grass (Allen et al. 1998) and it is generally lower for vegetation with a Crassulacean acid metabolism, such as sedum, and higher for plants that require more water for survival (Farg et al. 2012;Guerra et al. 2016). The crop coefficient is assumed to vary between 0 and 2.5. The reservoir constants K s and K f represent the inverse of the residence time, which in the case of fast reservoir can be considered a proxy of the basin concentration time. The residence time in the slow and fast reservoir is assumed to range from 1 h to 240 h and 10 h, respectively. The fraction of imperviousness c 0 varies between 0, if there is no urbanization, and 1, a condition in which all the surfaces are impermeable. The coefficient c 1 changes from 0, if the all the percolated water contributes to the aquifer recharge, to 1, if all the soil bucket leakages supply the slow reservoir. To calibrate the model, different parameter sets θ[nZr, s l , s 0 , K c , K s , K f , c 0 , c 1 ] were randomly generated through a Monte Carlo simulation, selecting the parameters from the previously described ranges. Among the different sets, only the first 20,000 with a Nash-Sutcliffe efficiency NSE (Nash and Sutcliffe 1970) higher than zero in the calibration period have been retained and considered behavioural (Beven and Binley 1992). NSE is defined as: where Q outi and Q obsi are simulated and observed runoff at i-th time step, while Q obs is the average of the observed runoff and N is the total number of time steps. NSE can vary between −inf and 1, where 1 indicates the perfect model representation of the observations. Considering only the simulations with NSE > 0 ensures that only sets of parameters that lead to a high model performance are analysed.

Model Calibration and Reference Set of Parameters
Through the Monte Carlo behavioural calibration, it was possible to obtain a set of parameters that ensure the highest performance in terms of NSE. For the investigated basins the maximum NSE during the calibration period (NSE = 0.72 in both cases) was achieved for the following parameter sets: However, as often reported in hydrologic literature, equifinality conditions are present (Beven and Freer 2001;Foulon and Rousseau 2018), implying that similar model performance (e.g. in terms of NSE) can be achieved by using different sets of parameters. To deal with this source of parametric uncertainty, out of the 20,000 behavioural Monte Carlo simulations, the best 5%, corresponding to the 1000 simulations with the highest NSE, were retained and used with two aims. The first was to identify a reference set, taking into account the parameters' a posteriori distribution. At the same time, we used this restricted parameter ensemble to provide model output with parametric uncertainties. Each model parameter cumulative probability distribution and the correspondent cumulative distribution weighted on the model performance NSE have been plotted in Fig. 3. The terms nZr, s l and s 0 were combined in one factor, nZr(s l − s 0 ), which represents the amount of water that can be stored in the soil bucket. The analysis of Fig. 3 reveals that parameter cumulative distributions are quite different for the two considered basins. Once the influence of NSE is considered, the cumulative distribution weighted on the model performance (orange) departs from the unweighted one (red) if the model has relatively higher performances for some parameter ranges. Otherwise, the two cumulative distributions are close to each other, as is the case for the fast reservoir constant K f . The model parameter probability distribution was derived from the cumulative distribution weighted on NSE: reference parameter sets were chosen considering their modal values. The following reference sets were selected for the investigated basins: These data sets lead to an NSE equal to 0.7 for #470 and to 0.69 for#507 during the calibration period. In the validation period, the model performance is, as expected, a bit lower, corresponding to 0.65 and 0.33, respectively. The lower performance during the validation period in the catchment #507 could be due to the high climatological difference between the calibration and validation period.

Model Results
Using the model outputs associated with the 5%-best behavioural simulations, uncertainty bounds were evaluated and used to define an envelope of the 5%-best model outputs. Fig. 4 illustrates, for both basins during calibration and validation periods, an example rainfall event (blue) and corresponding observed streamflow (red), simulated streamflow (black), obtained from the reference parameter set, and the 5%-best simulation model output envelope (green shadow). The plots highlight that most of the observed streamflow falls within the envelope of the 5%-best simulation model outputs: in particular 94% (#470) and 96% (#507) of the observed streamflow used for calibration is inside this envelope. For the validation period, 91% and 96% of the observations are contained in the 5%-best simulation interval, respectively. Different envelopes, built from the 2.5%-and 7.5%-best simulations, have been investigated, confirming the increase in the percentage of points falling within the envelopes as the percentage of simulations rises. Simulated outflow results have also been compared with the one obtainable from the ecohydrological streamflow model EHSM proposed by Viola et al. (2014), where aquifer recharge processes are not considered and the entire leakage derived from the soil bucket feeds the slow reservoir. The percentage of points that fall inside the envelopes is almost the same for both models, highlighting the good reliability of EHSMu. Moreover, EHSMu shows a higher model performance in terms of NSE than EHSM: during the calibration period for EHSM, the maximum NSE is equal to 0.61 and 0.64 for #470 and #507 respectively, while using EHSMu the best simulation led to 0.72. The higher performance obtained with EHSMu underlines the importance of evaluating the aquifer recharge as a fundamental component in the water balance for urban catchments.

Long Term Water Balance
EHSMu gives as outputs the runoff at the basin outlet, the actual evapotranspiration and the aquifer recharge at hourly time scale. We aggregated the water balance components over long time intervals to compare flux magnitude among calibration and validation Fig. 4 Example of rainfall events and corresponding observed (red) and simulated streamflow obtained with the reference set (black) and 5%-best simulation confidential interval (green areas) for the investigated basins, during calibration and validation period Fig. 5 Water Balance plotted for both basins, during the calibration and validation period using the reference parameter sets periods and between the basins. Fig. 5 shows the water balance over the calibration and validation period for the investigated basins, highlighting evapotranspiration, slow and fast streamflow and aquifer recharge as a percentage of the cumulated rainfall. Evapotranspiration corresponds to about 45% in all considered scenarios, 15% is composed of aquifer recharge and the runoff at the basin outlet is 40%, divided in 10% slow streamflow and 30% fast streamflow.
Water balance during the calibration period for basin #507 is slightly different, since it presents a lower (40% instead of 45%) total evapotranspiration and consequently a higher aquifer recharge. This is probably due to the statistical properties of rainfall: for the basin #507, in the calibration period there are less but more intense rainfall events than in the validation period. This implies higher leakages from the soil bucket, with a consequent increase of the aquifer recharge and slow streamflow components and at the same time less evapotranspiration.  100, 0.91, 0.3, 1.7, 0.55, 0.39, 0.29, 0.38] and θ 507 [nZr, s l , s 0 , K c , K s , K f , c 0 , c 1 ] = [100, 0.78, 0.35, 1.8, 0.25, 0.18, 0.3, 0.35]. Potential water balance variability is plotted for both basins, during the calibration and validation period

Land Use and Land Cover Changes
In this section, we analyse and discuss the potential impact of land use and land cover changes on the water balance of the Charlotte metropolitan area. The reference parameter sets have been used as zero scenario and variations of the imperviousness and crop coefficients are investigated.
The presence of houses, roads and car parks within an urban catchment is represented by the EHSMu with the imperviousness coefficient c 0 , which was chosen as equal to 0.29 for #470 and to 0.30 for #507 in the reference parameter sets θ 470 and θ 507 . These values are a good representative of the actual land use scenario presented in Fig. 2b, where about 30% of the surface presents a medium-or high-intensity development. To estimate the impact of land use changes on the urban basin water balance, a sensitivity analysis of c 0 was performed. Figure 6a illustrates the percentages of each component of the water balance, as discussed in Section 3.3, estimated using the following reference parameter sets where c 0 varies from 0 to 1: In Fig. 6a, the blue dashed line indicates the zero scenario. Values higher than 0.3 represent a potential future scenario, where the urbanization and the presence of high-density buildings are increased. For these scenarios, a higher percentage of fast runoff is expected, with a corresponding decrease of evapotranspiration and aquifer recharge. Results highlight how aquifer recharge, slow streamflow and evapotranspiration processes linearly decrease with higher values of the imperviousness fraction. Aquifer recharge is slightly higher in the larger basin (#507), where it could constitute 20-25% of the water balance, in the case of a no impervious surface: this can be due to the lower c 1 coefficient, which defines the percentage of leakage that feeds the slow linear reservoir.
The potential impact of vegetation coverage changes is evaluated by investigating different values of the crop coefficient K c . Fig. 6b investigates the potential impact on the water balance of possible future variations of the vegetation type, represented through a perturbation of K c : Lower values of K c represent vegetation with a low evapotranspiration rate, such as cacti and succulent plants. In this scenario, evapotranspiration decreases with a corresponding increase in aquifer recharge and slow streamflow. Fast streamflow, instead, is not affected by crop coefficient and variations in the total outflow are caused only by increase in the slow runoff component. This is clearly observable in all the cases investigated in Fig. 6a: slow streamflow generally varies from being more than 20% of the total water balance, when K c = 0, to 7-8%, for K c = 2.5. Aquifer recharge is also strongly influenced by the crop coefficient, varying between 35% to 10% for catchment #470, and from 40% to 15% for the larger basin.

Conclusions and Future Works
This work presented a new lumped conceptual ecohydrological streamflow model for urban areas that combines hydrological and ecological processes to estimate the runoff generation at the basin outlet at hourly scale. EHSMu was calibrated for two densely-urbanized catchments in the Charlotte metropolitan area. The 8 parameters that characterize the model were calibrated through Monte Carlo simulations, considering a subset of 20,000 runs with a NSE higher than zero. The model shows good performance, reaching a NSE equal to 0.72 in both basins: this underlines the importance of including the aquifer recharge component in the hydrological analysis and of evaluating in depth soil moisture dynamics in relation to evapotranspiration.
The water balance, analysed considering a reference set of parameters, highlighted that 45% of rainfall is lost due to evapotranspiration processes and 15% is directed to the aquifer. The runoff at the basin outlet, corresponding to 40% of the total rainfall, is composed of 25% of slow runoff derived from the soil bucket leakages and 75% of fast runoff from impermeable and saturated permeable surfaces.
The sensitivity analysis on the impacts of potential land use changes highlighted the drastic implication of an increase of imperviousness degree. Evapotranspiration, aquifer recharge and slow runoff linearly decrease for high imperviousness values, while the fast runoff largely increases. On the other hand, installing plants with a higher crop coefficient and faster metabolism will lead to an increment in evapotranspiration losses, with a consequent reduction in aquifer recharge and slow runoff.
This work aimed to present a conceptual model to estimate the runoff generation in urban areas. The model performance, however, could be further investigated using multiple indicators. Moreover, several different catchments should be evaluated in order to analyse the potential applicability of the EHSMu model to basins with different land use and land cover and located in different climatological regions. appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.