The study area is in the Aricanduva River basin, a left-bank tributary of the Tietê River, which is highly urbanized. It is situated in the eastern zone of the city of São Paulo and is part of a cluster of 39 municipalities that make up the urban conurbation of the São Paulo metropolitan region (Figure 1A). The basin stretches for approximately 20 km in the SSE-NNW direction and covers an area of 100.62 km². The elevation ranges from 994 m above sea level (asl) upstream to 718 m downstream (Figure 1B). The climate in the study area is characterized by a humid summer from October to March and a dry winter from May to August. Precipitation levels vary seasonally, with a climatological average of 211 mm in February during the high wet season, and 38 mm in August during the dry season (EM, 2020).

The classification of land use and land cover, based on physical aspects, shape, and surface texture, obtained using high spatial resolution satellite images (SPOT 2.5 m, Rapid Eye 5 m) and orthophotos (1 m) (Rossini-Penteado et al., 2007; Rossini-Penteado and Gilberti, 2008), resulting in eight distinct classes (Table 1).

The Aricanduva River basin is occupied by approximately 62% of small residential, commercial, and service buildings, as well as streets and sidewalks. It is occupied by 16% of large commercial buildings such as malls, commercial warehouses, and others. Vegetation covers 19% of the area, including urban forests, shrubs, grasslands, and trees in public spaces. Approximately 3% of the area consists of exposed soil without buildings or vegetation, transitioning into urban development (Table 1; Figure 1C).

The boundary conditions and meteorological forcing for the hydrological modelling were defined using a database containing topographic, geological available at: http://geosampa.prefeitura.sp.gov.br/PaginasPublicas/_SBC.aspx, meteorological and hydrological information available at: https://www.saisp.br/estaticos/sitenovo/home.html.

These following conditions are essential to set up the model simulations in the study area:

i) Topography and soils: Digital elevation model with 0.1 m resolution using LiDAR 3D data in the urban area (Figure 2), and soil classes used to estimate the Curve Number (CN) of the basin for the HEC-HMS model, obtained from GEOSAMPA/Secretaria Municipal de Urbanismo e Licenciamento (SMUL) of the city of São Paulo.

The flow and water level data observed at 10-min intervals during flood events underwent quality control procedures, including the application of filters to identify faults, sensor malfunctions, spurious records, and extreme values. Consistency checks were performed in three stages, following Feng et al. (2004) and Tomaz (2016): 1) Removal of values outside prescribed thresholds for maximum and minimum values; 2) Analysis of data from the previous and subsequent days for each event in the station series to detect inconsistent sequences; 3) Exclusion of inconsistent data by comparison with nearby stations and 4) Gap filling of missing flow and level data using Simple Linear Regression, constrained by strong correlation with one or more nearby stations (Bertoni and Tucci, 2020).

The spatial distribution of precipitation in the basin was determined using data from the four rain gauges (Figure 3; Table 3) calculated with the Thiessen Polygons method using the Analysis Tools—Proximity tool in ArcGIS, which provides the percentage of precipitation for the centroids of the sub-basins (Pokojski and Pokojska, 2018). For the hydrological modelling case studies, seven intense precipitation events were selected from a time series of 236 events between years 2010 and 2020. The events were chosen accordingly to classification of emergency or overflow conditions (Table 3), that were predominantly in the range of 45–115 mm (Figure 4). The measured rainfall was spatially distributed within the contributing area at a 10 min resolution.

Our study involved three main steps in the calculations: 1) a comprehensive dataset of boundary conditions and hydrometeorological forcings was prescribed to configure the hydrological contributing area with high-resolution and quality control. This dataset was then used to run the computational codes of the HEC-HMS and HEC-RAS 2D models; 2) the model was calibrated and validated using multiple case studies of flood events, with significant statistical performance and sensitivity analysis, and 3) hypothetical scenarios were built to evaluate the impact of introducing Green Roofs (GR), or permeable pavement (PP), on runoff and peak flow during a flood event. These scenarios involved replacing impermeable surfaces with permeable alternatives of GR or PP. The analysis focused on evaluating the changes in runoff and peak flow resulting from the implementation of those infrastructure.

The HEC-HMS (Hydrologic Engineering Center—Hydrologic Modelling System) version 4.10 was chosen to simulate the rainfall-runoff processes in our urban basin to quantify the flow during rapidly changing water level events. The computational code is integrated with database information and geoprocessing tools, resulting in various outputs of soil moisture and flow (Rouf, 2015; USACE, 2016a). This integration occurs through a database management platform known as HEC-DSS, designed to store and retrieve time series and paired data generated in hydraulic projects within the HEC family. HEC-DSS was built to standardize data transfer between the programs. All communication between the HEC-HMS and HEC-RAS models is conducted through these DSS-format files. Technical information on database generation can be accessed online at https://www.hec.usace.army.mil/confluence/dssdocs/dssvueum/introduction/overview-of-the-hec-data-storage-system.

The basin was defined by seven hydrological elements: 1) Subbasins: the physical subdivisions of individual micro basins within the larger basin. 2) Reaches: segments of the river where discharge occurs from upstream to downstream. 3) Junctions: nodes representing the confluence of rivers and their tributaries. 4) Sources: a flow introduced at a specific point in the system through the drainage area. 5) Sinks: water withdrawal from a specific micro basin. 6) Reservoirs: elements of water storage and loss within the system. 7) Derivations: the transfer of flow between different elements of the system (USACE, 2016a-b). With the definition of these elements, meteorological forcings were spatially distributed within the basin, that enabled the simulation for a temporal interval. The storage and discharge functions of the attenuation reservoirs were considered in the simulations (Decina, 2012). The simulations employed the simple methods of initial and constant loss and the unit hydrograph method, to estimate the runoff resulting from net precipitation in the subbasins (USACE, 2016a), see Table 4. The methods assumed a single hypothetical soil layer, considering changes in moisture content without accounting for detailed subsurface characteristics. Soil moisture uptake was neglected, supposed to be suitable for short-duration and intense rainfall events.

The estimation of urban flooding was performed using the hydrodynamic model HEC-RAS 2D version 6.3.1 (USACE, 2016b; Rangari et al., 2019). The 2D computational mesh was generated in the pre-processing step using a high-resolution digital terrain model (DTM) with 10 cm resolution. This DTM raster file was input into RAS-MAPPER to generate the computational mesh, allowing for the estimation of water level elevation and flow velocity and the temporal and spatial evolution of the flood. The floodplain was modeled using a 2D flow mesh consisting of 2,853 cells, each measuring 40 m × 40 m. Manning’s roughness values (n) were spatially prescribed in the model based on the land use classification (Table 1; Figure 1C) of the floodplain. The values were estimated as the average for the left and right banks and the main channel of the Aricanduva River. In addition, flow hydrographs and water surface profiles (hydrograph) were used as boundary conditions to simulate the spatial and temporal evolution of flow in the basin. The simulations were conducted using a time step of 1 min, enabling a detailed analysis of the flood dynamics.

To ensure the consistency and robustness of the modelled scenarios, we previously built the calibration and validation of the HEC-HMS and HEC-RAS 2D models (Cunderlik and Simonovic, 2004; USACE, 2016a). For the HEC-HMS model, the initial parameter values for Initial Loss (I_L) and Constant Rate (CR) were predefined based on USACE (2016a). The automatic calibration was performed using the Univariate method with a maximum of 100 iterations and a tolerance of 0.01 in the optimization tests. Four precipitation events (Table 3) were used for calibration. Manual calibration (Golmohammadi et al., 2014; Ouédraogo et al., 2018) was only required in cases where convergence was not achieved during the automatic calibration. The values of the Impervious (I_M) parameters were determined based on the Land Use classification (Table 1; Figure 1C) (USDA, 2017), and the Lag time (LT) parameters were determined in the pre-processing step using ArcHydro and HEC-GeoHMS extensions (USACE, 2016a). The optimized parameter values are shown in Supplementary Table S1, which were found to be comparable to estimates of permeable pavements (Collins et al., 2008; Monrose and Tota-Maharaj, 2018; Zheng et al., 2021).

In the HEC-RAS model, flow, and water level observations at P2 and the observed water level at P3 were used as boundary conditions in the internal 2D grid to evaluate the lateral contributions from streams. The observed water level at P4 (Figure 2) was also used for validation. Initial values of the Manning’s roughness coefficient were estimated in the 2D flow mesh area (Chow, 1959) (Figure 5). The maximum elevation of the event was recommended to adjust the Manning’s roughness coefficient (n) in the mesh, as described in the application guide provided by USACE (2016b). The first step of the calibration estimated the Manning’s roughness coefficient (n) considering different land use types on the right and left banks of the river, and hydraulic characteristics of the channel, with values of 0.034, 0.035, and 0.036, respectively. To ensure the best condition of stability and numerical accuracy, a spatial resolution of 40 m × 40 m was used for the 2D flow mesh. Subsequently, numerical stability tests were performed for different time steps, and instability issues were observed for computational intervals above 5 min. Therefore, a time step of 1 min was chosen for the final simulations.

The performance was assessed based on the following statistical indices: PEV, PEPF, PEE, R², RSR, NSE, and PBIAS, as listed in Supplementary Table S2.

To enhance the understanding of the model’s functioning as a predictive tool for this study, we assessed the sensitivity of key parameters. Sensitivity analyses assess how variations in the value of a parameter can result in significant disparities between observed and simulated volumes, capturing the isolated impact of each variation (Cunderlik and Simonovic, 2004; USACE, 2016a). The following parameters were selected: 1) Initial loss (in mm); 2) Constant rate (in mm/h); 3) Imperviousness (in %); and 4) Lag time (in minutes). The optimized values obtained after the calibration of each parameter were perturbed from −25% to 25%, with increments of 5%, while keeping the remaining parameters constant.

In our simulations, we assumed that all activities following the pre-development period in the upstream basin of observation point P2 influenced the hydrological properties of the soils, leading to a shift towards soil types D and D/D with Curve Numbers (CN) ranging from 77 to 98, as suggested in TR-55 (USDA, 1986) and estimated during preprocessing using ArcHydro and HEC-GeoHMS extensions (USACE, 2016a) (Figure 5). We also assumed that all impervious surfaces are hydraulically connected to the micro and macro drainage network upstream of observation point P2.

The two GIs considered in this study, namely, Extensive Green Roofs (GR) and Permeable Pavement (PP), form the basis of our simulation scenarios. Technical information on the hydraulic behavior of extensive green roofs and permeable pavement used in this study was obtained from literature sources, drawing upon various microscale experiments conducted worldwide (Bateni et al., 2021; Zheng et al., 2021). The retention of precipitation by green roofs correlates inversely with the volume of rainfall, with retention efficiency varying considerably based on rainfall event volume. On average, green roofs can retain approximately 62.2% of rainfall and 69.3% of peak flow (Berndtsson, 2010). Carter and Rasmussen (2006) found that the percentage of rainfall retained for small storms (<25.4 mm) reached 88%, retention for medium storms (25.4–76.2 mm) was around 54%, and for large storms (>76.2 mm), it could be as high as 48% of precipitation. For rainfall events with precipitation less than 10 mm, green roofs can retain nearly all the rainfall, as reported by Zheng et al. (2021). For a 12 mm precipitation event, green roofs showed retention ranging from 26% to 88%, depending primarily on precipitation intensity, residual moisture in the substrate, vegetation species, geometric properties (such as coverage area and slope), substrate characteristics (type, depth, porosity, and density), and drainage layer (type and depth), as observed by Simmons et al. (2008). Larger precipitation events of 28 mm and 49 mm showed retention ranging from 43% to 8% and 44% to 13%, respectively (Berndtsson, 2010; Liu et al., 2020a; Zheng et al., 2021).

In contrast to traditional pavements that generate runoff from almost all rainfall, permeable pavements encourage rainfall infiltration, increasing moisture within the soil profile (Bateni et al., 2021) and has multiple applications in commerce, and residences, reducing runoff volume by over 40% and decreasing peak flow by 50% to 90%. Commercially, there are five categories: permeable asphalt (PA), permeable concrete (PC), permeable interlocking concrete pavement (PICP), and concrete grid pavement (CGP) or plastic grid pavement (PGP), filled with grass or sand and gravel (Collins et al., 2008; Bateni et al., 2021). Gomez-Ullate et al. (2011) demonstrated that the hydrological performance of all permeable pavements (PPs) in terms of peak flow reduction is very similar, with runoff reduction potentially exceeding 98%. Weiss et al. (2015) found that variability in peak flow and infiltration rates among different types of PPs was primarily influenced by design, material quality, aggregate mix, periodic maintenance, and other factors (Lin et al., 2014; Barszcz, 2015; Bateni et al., 2021).

The evaluation of the effectiveness of GIs involves the use of values from the Soil Conservation Service—Curve Number (SCS-CN) method to simulate runoff based on local precipitation, land use, and soil data (Ahiablame et al., 2012; USDA, 2017). The CN method is an empirical two-parameter procedure (CN and initial abstraction, S) widely used in stormwater management as well as in complex watershed models to determine how much of a given rainfall event becomes direct runoff (USDA, 1986; Sample et al., 2001). The initial abstraction, which describes all precipitation losses before runoff begins (interception, infiltration, surface storage, and evaporation), is a function of CN and is calculated as (USDA, 1986): S = 25400 C N − 254 (1)

Considering that the precipitation, P _h (mm), is greater than 0.2S, the direct runoff, Q _h (mm), is estimated as: Q h = P h − 0.2 S 2 P h + 0.8 S (2)

Q _h = 0 when P _h ≤ 0.2S Q v = Q h x A (3) where Qv is the water volume, and A is the area of interest.

In this context, simulations of GRs and PPs with the HEC-HMS model involve the use of SCS-CN values to estimate runoff. CN is a key parameter common to both GIs used in this study. The two green infrastructures (GRs and PPs) are represented with initial CN values suggested in USDA (1986), Sample et al. (2001), MDE (2009), and Ahiablame et al. (2012). In the simulations, we will assume that all activities following the pre-development period in the watershed have influenced the soil’s hydrological properties, leading to a shift to hydrological group D, with CN ranging from 77 to 98, as suggested in TR-55 (USDA, 1986) and estimated in the preprocessing using ArcHydro and HEC-GeoHMS extensions (USACE, 2016b). In the study, we also accept that all impervious surfaces are hydraulically connected to the micro and macro drainage network in the watershed.

The CN values used for green roofs are 85 for all four Hydrologic Soil Groups A, B, C, and D, corresponding to low, moderately low, moderately high, and high runoff potential, respectively (Sample et al., 2001). For permeable pavement, the CN values used are 70, 80, 85, and 87 for the 4 Hydrologic Soil Groups A, B, C, and D (Sample et al., 2001). These CN values were recommended to adjust the model values to characterize the effects of GIs on runoff, allowing for comparison between hydrological conditions before and after GI implementation. Therefore, for the construction of simulation scenarios with green infrastructures (GIs) in the Aricanduva River basin, we considered extensive green roofs (substrate thickness <15.24 cm) and a reduction in the Curve Number to 77, along with permeable pavements (PPs) with a Curve Number of 79. The hydraulic conductivity for both green roofs and permeable pavements was set at 25.4 mm/h, following the guidelines of MDE (2009) and USDA (2017). Simulations incorporating green roofs and permeable pavements adhered to the following criteria: i) the grey roofs of buildings were replaced with green roofs; ii) areas with established vegetation in the basin were not replaced by green infrastructure; iii) paved areas (streets, sidewalks, parking lots, and open areas on private properties) were replaced by permeable pavements.

The hypothesis of this investigation aims to evaluate the hydrological response of the urban watershed under interference conditions with green infrastructure, such as Green Roofs (GR), or gray infrastructure, such as Permeable Pavement (PP). We assume significant attenuation of the hydrological pulse in both conditions but with different performances. The interference of GR or PP is compared to the Control simulation, which serves as the model calibration.

Interference boundary conditions were applied to the 13 sub-basins within the domain (Table 5). However, hydrological results were effectively derived only from the upstream contributing area of gauge P2, which encompassed 9 sub-basins. These sub-basins consisted of built-up areas, green spaces, and other small spaces, totaling an area of 67.01 km² (Supplementary Table S3). We used the precipitation event of 16/02/2019 when the river overflowed at gauges P2 and P3.

The Control simulation was configured with a uniformly distributed conventional (gray) roof across the watershed. The average built-up area on impermeable gray roof surfaces within the watershed was 36%, varying from 19% (sub-basin W880), with extensive vegetated areas, to 48% (sub-basin W550), with high building density.

For streets and sidewalks in public areas, the average built-up area was 22% across the watershed, ranging from 19% (sub-basins W810 and W750) to 26% (sub-basin W530). In private areas, the ground plot area was defined as the open space available for the insertion of PP or internal open space (plot area minus roof area), accounting for 41% of the total area across the watershed. The variation in ground plot area ranged from 26% (sub-basin W530) to 60% (sub-basin W880). Detailed occupancy information can be found in Supplementary Table S3.

Thus, we constructed 7 simulation scenarios for GR or PP (Table 5). Starting from scenario 1 (Control), scenarios 2, 3, and 4 introduced GR roofs, replacing gray roofs in increasing order of individual plot area (A), i.e., conditioned on A > 500 m² (scenario 2), A > 200 m² (scenario 3), and A > 50 m² (scenario 4), respectively, thereby increasing the surface coverage with GR. In the remaining scenarios, PP was introduced in 50% of the ground plot area for scenario 5, in streets and sidewalks of public areas for scenario 6, and in scenario 7 as the combination of scenarios 5 + 6.

The calibration and validation of the model generally showed high performance in simulating runoff volume and peak flow (Supplementary Table S4), within the acceptable range of indices: R2 > 0.50, NSE ≥0.65, PBIAS ≤ ±10%, and RSR ≤0.60 (Moriasi et al., 2007; Ouédraogo et al., 2018). Figure 6 shows the good agreement between observed and simulated runoff in the model validation events. Overall, we noticed that the fit was slightly better in the rising limb of the hydrograph compared to the falling limb. There was a small delay in the model’s peak flow timing in some cases (Figures 6B, C, F), no more than 15 minutes, while in other cases, the simulation of the peak timing was very close (Figures 6A, D, E).

In the event of 16/02/2019, the water level at gauge P2 (Figure 7A) caused flooding when it reached the maximum channel elevation (737.04 m asl) at approximately 17:50 h. It then rose an additional 0.76 m by 18:40 h, reaching a peak discharge of 306 m³/s. Finally, it returned to the overflow level 2 h later at 19:50 h, leaving flooded areas downstream of the gauge. At gauge P3, the hydrological pulse was slightly ahead of P2 but followed a similar pattern. It started overflowing at 17:30 h, rose 1.29 m, reached its maximum elevation at 18:40 h, and then gradually receded back to the channel level after approximately 2 h (Figure 7B). At gauge P4, the river reached a maximum elevation of 727.89 m at 19:30 h, below the overflow level (Figure 7C). The calibration of water levels at the three gauges showed excellent performance in this event and equally in the validation simulations (not shown), as indicated by the obtained statistical indices (Supplementary Table S5).

The 5 attenuation reservoirs played a role in delaying the inflow of water into the main channel and mitigating both the peak flow and total runoff volume in the channel, as evident from the flow hydrograph pattern (Figure 6F). However, they were unable to prevent the occurrence of overflow. The reservoirs had a maximum retention capacity of 2.02 million m³ and operated upstream of gauges P2 and P3. In the event of 16/02/2019, all reservoirs operated below their maximum capacities. When calculating the ratio between the peak storage during the event and the maximum storage capacity, we observed values ranging from approximately 27% in reservoir 2 to 46% in reservoir 4. The total volume attenuated in the five reservoirs was approximately 2.35 million m³, which retained 49.7% of the observed volume (4.67 million m³) in P2, with peak discharges in 4 reservoirs occurring after 20:00 h (Table 6).

For a more detailed verification of the role of the reservoirs, we conducted additional simulations of the event by sequentially disconnecting the reservoirs from the drainage pattern, one by one cumulatively, so as not to assimilate the surface runoff water. In summary, Table 7 shows that the rate of change in runoff volume increases by less than 2% compared to the peak flow, which increases at more significant rates of up to 26.4% because of retention.

The model simulated the flooding over a spatial extent around the overflow of the Aricanduva River canal, as observed at gauges P2 and P3 (Figure 8A) and noted in records from the Flood Alert System of São Paulo (FCTH, 2022). For a pictorial illustration, we sketched the elevation level reached on the side of the river channel, with the water level above the level of vehicles on the avenue (Figure 8B).

We conducted a simple sensitivity analysis to assess the scale dimension of the hydrological response, as it depends on the variation of key model parameters that characterize the degree of surface impermeability. With a 25% increase in impermeability, we observed a 3.7% increase in runoff volume, and with a 25% reduction, there was a corresponding 3.7% reduction in volume (Figure 9A). According to the Percent Error in Peak Flow (PEPF) criterion, the peak flow showed higher sensitivity to the Lag time parameter, exhibiting an inversely proportional relationship. A 25% increase in Lag time resulted in a 9.8% reduction in peak flow, while a 25% reduction led to a 9.8% increase (Figure 9B).

The other parameters showed lower sensitivity in the model. For example, the Initial loss parameter resulted in a variation of approximately ±0.4% in runoff volume and ±2.0% in peak flow. Similarly, the Constant rate parameter, associated with saturated hydraulic conductivity, exhibited low sensitivity, with the runoff volume varying approximately between ±1.4% and the peak flow between ±1.1% (Figures 9A, B) (Supplementary Table S6).

The Manning’s roughness coefficient (n), which varies spatially, showed a more pronounced inverse relationship with the depth of inundation (in percentage terms). The variation in n by approximately ±25% resulted in a depth change of approximately ±5.3% (Figure 9C). On the other hand, perturbing the cell area of the computational grid had a very small effect on the depth of inundation (in percentage terms), with a variation of approximately ±0.9%.

We identified the Impervious parameter as the most sensitive based on the Percent Error in Volume (PEV) criterion.

We conducted simulations for scenarios 2, 3, and 4, gradually expanding the coverage of GR in the contributing area of the watershed. In scenarios 2 and 3, 1,164 and 6,444 conventional roofs were replaced upstream of P2, totaling 1.6 and 2.9 km² of area, resulting in a percentage of permeable area of 3.5% and 6.3%, respectively. Scenario 4 prescribed a significant increase in permeable surface, with 125,843 conventional roofs covering 13.9 km² with GR, resulting in a permeable area percentage of 30.4%, that showcased approximately 9 times more than scenario 2 and 5 times more than scenario 3.

The introduction of GR confirmed the hypothesis of reduced runoff impact. In scenarios 2 and 3, the peak flows were significantly reduced compared to the Control scenario and were very close to each other throughout the runoff event (red and green lines in Figure 10). Despite reducing peak flow by around 15% and runoff volume by approximately 23%, these scenarios fell short in preventing canal overflow, as the runoff remained above the critical level for about 1 h. In contrast, scenario 4 showed a significantly different outcome. It achieved a reduction of approximately 37% in peak flow and approximately 31% in runoff volume, and importantly, it successfully prevented canal overflow (orange line in Figure 10).

The simulations of scenarios 5, 6, and 7 described, in sequence, an increasing area of the watershed covered by PP instead of traditional gray pavement. In scenario 5, only public areas such as streets and sidewalks were replaced, while in scenario 6, only private areas (50% of the opened ground) were replaced. Scenarios 5 and 6 both entailed a replaced area of 10.1 and 9.8 km² respectively, with a permeable area percentage of 22.1% and 21.4% upstream of P2, indicating comparable extents of permeable areas in both. Scenario 7 represented the combination of scenarios 5 and 6, resulting in a total replaced area of 19.9 km² and a permeable area percentage of 43.5%.

The results of scenarios 5, 6, and 7 confirmed the hypothesis of attenuating the hydrological pulse with the introduction of PP. In scenarios 5 and 6, the reduction in total volume was approximately 43% and 44% respectively, and the reduction in peak flow of 31% and 32% respectively (Figure 11). In both cases, calculated the peak flow was approximately equal to the critical overflow level. Scenario 7 excelled, reducing runoff volume by 52% and peak flow by 44%, successfully averting canal overflow unlike the previous scenarios.