Space variability of hydrological responses of Nature-Based Solutions and the resulting uncertainty

During the last decades, the urban hydrological cycle has been strongly modified by the built environment, resulting 10 in fast runoff and increasing the risk of waterlogging. Nature-Based Solutions (NBS), which apply green infrastructures, have been more and more widely considered as a sustainable approach for urban stormwater management. However, the assessment of NBS performance still requires further modelling development because of their hydrological responses sensitively depends on the representation of multiscale space variability of both the rainfall and the NBS distribution. Indeed, we initially argue this issue with the help of the multifractal intersection theorem. To illustrate the importance of this question, the spatial 15 heterogeneous distributions of two series of NBS scenarios (porous pavement, rain garden, green roof, and combined) are quantified with the help of their fractal dimension. We point out consequences of their estimates. Then, a fully-distributed and physically-based hydrological model (Multi-Hydro) was applied to consider the studied catchment and these NBS scenarios with a spatial resolution of 10 m under two different types of rainfall: distributed and uniform, and for three rainfall events. These simulations show that the impact of spatial variability of rainfall on the uncertainty of peak flow of NBS scenarios 20 ranges from about 8 % to 17 %, which is more pronounced than those of the total runoff volume. In addition, the spatial variability of the rainfall intensity at the largest rainfall peak responds almost linearly to the uncertainty of the peak flow of NBS scenarios. However, the hydrological responses of NBS scenarios are less affected by the spatial distribution of NBS. Finally, the intersection effects of the spatial variability of rainfall and the spatial arrangement of NBS seem more pronounced for the peak flow of green roof scenarios and the total runoff volume of combined scenarios. 25


Introduction
The increasing of extreme flooding risks seem strongly related to two essential drivers: rapid urbanization and climate change (Lovejoy and Schertzer, 2013). Adapting to climate change and mitigating urban flooding constitute now significant societal https://doi.org/10.5194/hess-2020-468 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License.
The rainfall module (MHRC) can treat different kinds of rainfall data (from radar or rain gauge). In order to adapt the rainfall inputs for Multi-Hydro, the intersection between the pixels of the model (with a 10 m spatial resolution in this study) and the pixels of the X-band radar data (with a 250 m spatial resolution) were performed by the QGIS interface using the following 160 equation (Paz et al., 2018): where 10 , 10 is the rainfall rate computed for each pixel of the model; 250 , 250 is the rainfall rate of X-band radar data in each pixel; 10 , 10 , 250 , 250 is the area of each pixel of the model and the X-band radar data, respectively.
The surface module (MHSC) of Multi-Hydro uses the code of the Two-Dimensional Runoff Erosion and Export (TREX) 165 model that computes the surface runoff and infiltration at each pixel, depending on the land use classification (Velleux et al., 2008). The diffusive wave approximation of the Saint-Venant equations is used for calculating the overland flow, following the conservation of mass and momentum equations. The groundwater module (MHGC) is based on the Variably Saturated and 2-Dimensional Transport (VS2DT) model, which simulates the process of infiltration in the unsaturated subsurface zone. For a simplification of the NBS scenarios, this module is not applied in this study. The drainage module (MHDC) in Multi-Hydro 170 uses the code of 1D SWMM model (James et al., 2010) to simulate the sewer network. This model represents the flow computed by 1D Saint-Venant equations in conduits and nodes.
The high spatial resolution of Multi-Hydro allows an easy implementation of small-scale controlled measures, like the rain garden, green roof, bio-retention swale, porous pavement, and rainwater tank, by locally modifying the land use parameters to link the size and shape of the corresponding NBS infrastructures, with their infiltration and storage capacities. 175

Multi-Hydro implementation in the Guyancourt catchment
Based on the fully distributed character of Multi-Hydro, users can choose a specific spatial resolution. In this study, Multi-Hydro was implemented with a 10 m spatial resolution (the grid system creates square grids with a cell size of 10 m), and a temporal loop of 3 min. The 10 m resolution was performed because it sufficiently represents the heterogeneity of the 180 catchment, and also for saving the computation time.
The implementation of Multi-Hydro in a new catchment starts with the conversion of the original GIS data (e.g., land use, topography) into the standard rasterised format with the desired resolution by using the MH-AssimTool (Richard et al., 2013), a supplementary GIS-based module for generating the input data for Multi-Hydro). During this process, a unique land use class was assigned to each pixel, specifying its hydrological and physical properties. In order to attribute a unique land use class to 185 each pixel, the following priority order was used in this study: gully, road, parking, house, forest, grass, and water surface.
Because the gully is the only land use class able to connect the surface module and the drainage module, it has the highest priority (i.e., if a raster pixel contains gully and the other land use classes, the whole pixel will be considered as gully).
Generally, this order considers the impervious land use classes have higher priority than the permeable land use classes, which result in an overestimation of impervious land uses (see Ichiba et al., 2017, for an alternative approach). After the rasterization 190 process, the impervious land uses occupied 54 % of the Guyancourt catchment (Fig. 4). In this study, all the standard model parameters related to the land use classification were selected from the Multi-Hydro manual (Giangola-Murzyn et al., 2014), as they are shown in Table 2. Besides the land use, the elevation is also assigned to each pixel, which was obtained by performing the interpolation on the raw DEM information, freely available at the 25m-resolution.
Green roof is a special NBS measure that can be simulated by a specific module in Multi-Hydro (Versini et al., 2016). 195 Accordingly, five physically-based parameters are defined for the green roof. They are based on the experimental site of Cerema (Ile-de-France) where several green roof configurations were monitored (see Versini et al., 2016). In detail, the chosen configuration is the following: substrate thinness of 0.03 m and characterized by a porosity of 39.5%, an initial moisture condition of 10 %, a field capacity of 0.3, and a hydraulic conductivity of 1.2 m h -1 .

Simulation scenarios 200
For achieving the purpose of the study, a series of NBS scenarios were created and simulated under both different types of rainfall (described in Sect. 2.2). The baseline scenario is considered as the current configuration of the Guyancourt catchment, without implementing any NBS (Fig. 2 left). The baseline scenario will be used later on for the model validation. All details concerning the scenarios of the NBS implementations, including a detailed description of each NBS and the percentage of the space required for its implementation, are presented in Table 3, while the maps of the resulting land use are illustrated on Fig  205 5.
The first set of NBS scenarios includes porous pavement (PP1), rain garden (RG1), green roof (GR1), and their combined scenario (Combined1). For each scenario, the corresponding NBS are implemented heterogeneously over the catchment, while respecting the local catchment conditions and stormwater management requirements. They are applied to assess the impact of spatial variability of rainfall on the hydrological responses of NBS scenarios. 210 The second set of NBS scenarios (PP2, RG2, GR2, and Combined2) was proposed with a different arrangement to assess the potential effects of a heterogeneous implementation of NBS at the urban catchment scale. Regarding these two sets of NBS scenarios, the ones related to the same type of NBS (e.g., PP1 and PP2) require the same percentage of the space for their implementation over the whole catchment. But, both scenarios significantly differ in terms of spatial distributions of the considered asset. This distribution is characterized by the across-scale indicator, called fractal dimension presented in details 215 in the following paragraph.

Fractal dimension of NBS scenarios
To quantify the multi-scale space heterogeneity of NBS in each NBS scenario, we applied the concept of fractal dimension (DF), which was initially introduced to describe the scale invariance of some irregular geometric objects (Mandelbrot, 1983). https://doi.org/10.5194/hess-2020-468 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License.
Namely, a similar structure can be observed in any scale. DF has been often used in catchment hydrology (e.g., Schertzer and 220 Lovejoy, 1984;Schertzer and Lovejoy, 1987;Schertzer and Lovejoy, 1991;Lavallé e et al., 1993;Gires et al., 2013;Gires et al., 2016;Ichiba et al., 2017;Paz et al., 2020;Versini et al., 2020). In this study, a standard box-counting technique was applied to estimate the DF of each NBS scenario (Hentschel and Procaccia, 1983;Lovejoy et al., 1987). The DF of a geometrical set A (here represented by the non-overlapping pixels of NBS embedded in a 2-D space) is obtained with the following power-law: where , is the number of non-empty (containing NBS) pixels to cover the set A, at the resolution , which is defined as the ratio between the outer scale L and the observation scale (λ = ). The symbol ≈ means an asymptotic relation (i.e. for large resolution and possibly up to a proportionality prefactor).
Based on Eq.
(2), we count the number of pixels containing at least one NBS, by starting with the smallest pixel size ( = 10 m in this study), then continuously increasing the pixel size by simply merging the 4 adjoined pixels. This procedure is repeated 230 until reaching the largest pixel size (L). Thus, , is counted at different resolutions, and the results are plotted in the log-log plot (see Fig. 6). Corresponding to Eq. (2), the fractal dimension DF of each NBS scenario is defined as follows: Here, for each scenario, a square area of 128 x 128 pixels was extracted from the catchment to make the fractal analysis (see the example of the PP1 scenario in Fig. 5). In order to avoid the "no data" areas, which would bias the fractal dimension 235 estimate, the selected square area is the greatest possible size characterized by a multiple of two in the studied catchment.
As shown in Fig. 6, all the NBS scenarios are presented with two scaling behaviour regimes, with a scale break roughly at 80 m. For each regime, the scaling is robust, with linear regression coefficients (R 2 ) around 0.99. For the first regime corresponding to the small-scale range (10 m -80 m) that related to the assets implementation level, the dimension DF is around 1 for most of NBS scenarios. It is in contrast with the second regime, the large-scale range (from 80 m to 1280 m) that 240 exhibits a scaling behaviour with a DF ranging from about 1.75 to 1.98. We also applied fractal analysis on the impervious surface of the baseline scenario in the same selected area, and we also found the same scale break at 80 m (the DF of the baseline scenario in each regime are presented in Fig. 6). Therefore, it rather confirms that the spatial distribution of NBS is strongly constrained by the urbanisation level of the catchment.
The DF of each NBS scenario is summarized in Table 3. The dimension DF measures the implementation level of NBS across 245 scales. For instance, the DF (large-scale) of the two combined scenarios (Combined1 and Combined2) is close to 2, showing that NBS are rather homogeneously distributed. However, it is important to notice that, in spite of initially identical percentage at a given scale of the NBS implementation over the catchment, the resulting DF could be quite different. It is simply because the percentage of the space is a scale dependent quantity, while DF quantifies the propagation of the spatial heterogeneity for each of NBS scenarios, from the smallest scale to the outer scale of the catchment. This propagation remains scenario 250 dependent and hence a subject to its optimisation. https://doi.org/10.5194/hess-2020-468 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License.

Multifractal intersection theorem
We would like now to illustrate and emphasise why it is so much indispensable to take into account the multiscale space variability of both the rainfall and the NBS distribution. For instance, both "hot spots" (extremes) of the rainfall and NBS are scarce and therefore could rarely coincide, i.e., rainfall spikes may fall more often elsewhere than on NBS. Similar questions 255 can occur for less extreme events. The effective NBS performance could be therefore biased with respect to their potential performance due to this problem of intersection between rainfall intensity and NBS. It reminds us of the so-called multifractal intersection theorem applied to the intersection of a rainfall with extreme space variability and a rain gauge network that provides quantitative estimates of this intersection (Tchiguirinskaia et al., 2004). Figure 7, adapted from this paper, schematically represents the intersection at a given time of a (multifractal) rainfall, displaying quite variable pixel intensities 260 ranging from light blue to dark brown (e.g., from 1 to 100 mm h -1 ), with a heterogeneous rain gauge network (light brown pixels). The resulting measured rainfall field M is simply the product of the rainfall intensities R by the gauge characteristic function N (=1 if there is a gauge in this pixel, 0 otherwise). The intersection theorem states that for fractal objects, like for the usual (Euclidean) geometric ones, the codimensioni.e. the complement = − of the dimension to the embedding space dimensionof the measured field above a given intensity threshold is the sum of the codimensions of the network 265 ( = − ) and of the "real" field ( = − ) above the same intensity threshold: For instance, the intersection in a plane ( = 2) of two straight lines ( = = 1; = = 1) corresponds to a point ( = 2, = 0). Of particular interest is the case where the intersection is so scarce that its codimension cM is larger than the embedding dimension d, i.e. has a negative dimension DM (Schertzer and Lovejoy, 1987). Due to Eq. (4), the codimension 270 of the network is thus the critical dimension of the (real) field under which the rainfall intensity is rarely measured by the network: More precisely, the smaller DR is with respect to cN , the rarer the real field R is measured. Let us mention that Paz et al., (2020) used this intersection theorem to determine when the adjustment of radar data by a rain gauge network becomes misleading 275 instead of improving the data.
The assessment of the performance of an NBS network cannot be reduced to the binary question of presence or not of an NBS, like done for a rain gauge of a network. However, we can immediately state they will be more and more ineffective for rainfall intensity whose fractal dimension is more and more below the codimension of the network. This is already an important information that can be used to design NBS and their networks. This also explains why we estimated in the previous subsection 280 the fractal (co-) dimension of the NBS network, as well as to compare in section 4.3 simulations resulting from spatially uniform rainfalls ( = , = 0) and spatially heterogeneous rainfalls ( < , > 0). https://doi.org/10.5194/hess-2020-468 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License.

Modelling experiments
The overall target of the study is to investigate whether the spatial variability of rainfall and the spatial arrangement of NBS have an impact on the hydrological responses of NBS scenarios at the urban catchment scale. For this purpose, three sets of 285 modelling experiments were prepared, and two indexes (PEQp, percentage error on peak flow; PEV, percentage error on total runoff volume) were used for quantifying the uncertainty associated to rainfall and NBS spatial distribution in the hydrological response of the catchment. Figure 8 presents the flow chart of the three sets of modelling experiments. In addition, the corresponding descriptions are presented as follows: The first set is used to investigate the impact of spatial variability of rainfall on the hydrological responses of NBS scenarios. 290 In this first set, we employed the following scenarios: baseline, PP1, RG1, GR1, and Combined1. These five scenarios were simulated under the distributed and uniform rainfall events. Then, we compute the ratio on peak flow (Eq. (6)), and the PEQp and PEV (Eq. (7) and Eq. (8)) indexes for each scenario under two different kinds of rainfall.
The second set is used to analyse the impact of the spatial distribution of NBS on the hydrological responses of NBS scenarios.
In this experiment, we compared the two groups of NBS scenarios mentioned in the Section 3.2 (GR1 vs GR2 for instance). 295 The eight scenarios were only simulated with the uniform rainfall in order to avoid the impact of spatial variability of rainfall and to focus on the uncertainty associated with the spatial arrangement of NBS.
The third set is used to analyse the intersection impact of spatial variability of rainfall and the spatial distribution of NBS on the hydrological responses of the catchment. In this experiment, the eight mentioned NBS scenarios were simulated under the distributed and uniform rainfall, respectively. Then, the PEQp and PEV of each NBS scenario were computed by comparing the 300 results obtained for the two different kinds of rainfall events (distributed and uniform). Finally, we compare the difference of PEQp and the difference of PEV between the NBS scenarios characterized by the same solutions/measures. The peak flow ratio and the two indexes are especially calculated for the sum of four highlighted conduits connected to the catchment outlet (the right side of Fig. 2) with Eq. (6), (7) and (8): where 1 and 2 refer to the peak flow of scenarios under distributed rainfall and uniform rainfall respectively for the first and third modelling experiment. For the second experiment, they represent the peak flow of the first set of NBS scenarios and the second set of NBS scenarios, respectively. Correspondingly, for the first and third modelling experiments, 1 and 2 refer 310 to the total runoff volume of scenarios under the distributed and uniform rainfall respectively. For the second modelling experiment, they represent the total runoff volume of the first set of NBS scenarios and the second set of NBS scenarios, respectively.

Validation
Before the simulation of NBS scenarios, Multi-Hydro was validated with the water levels of the storage basin by applying the 315 baseline scenario under the three distributed rainfall events. The simulations were then repeated with the three uniform rainfall events, respectively. The model performance was evaluated through two indicators: Nash-Sutcliffe Efficiency (NSE) and percentage error (PE). The Nash-Sutcliffe Efficiency (NSE ≤1) is an indicator generally used to verify the quality of the hydrological model simulation results, described as follows: where refers to simulated values, refers to observed values, and represents the average of the observed values. The NSE closer to 1 indicates that the model is more reliable, whereas NSE closer to 0 indicates that the simulation does not better than that of the average observed value , which means the simulation performance is rather poor. If NSE is far less than 0, it means the simulation is even less performing than .
The percentage error (PE) represents the difference between observed values and simulation values, which reflects the 325 reliability of the simulation values.
With respect to the observed and simulated water levels in the baseline scenario, the NSE coefficients and PE values are summarized in Table 4. For the three distributed rainfall events (Fig. 9), the NSE are larger than 0.9, and PE are lower than 5 %. For the uniform rainfall event of EV2, the model represents the water levels with NSE equal to 0.95, and PE equal to 330 1.96 %: only a slight overestimation of the observed water levels is observed between hours 4 to 7. For the uniform rainfall of EV1 and EV3, the temporal evolutions of simulated water levels slightly underestimate the observed ones, with NSE around 0.8, as well as PE around 7 %. Regarding the temporal evolutions of simulated water levels under the distributed rainfall of EV1 and EV3, they are more consistent with the observed ones. The reason is that the rainfall intensities of the distributed rainfall are generally higher than those of the uniform rainfall at the storage basin location. Namely, in uniform rainfall events, 335 the accumulated water levels in the storage basin are less than that of in distributed rainfall events. Overall, the distributed rainfall gives slightly better results, and the simulated water levels using uniform rainfall also match sufficiently well the observed ones to validate the Multi-Hydro implementation in the Guyancourt catchment.
Hence, the Multi-Hydro is suitable and sufficiently reliable to investigate the impacts of spatial variability, either of rainfall and/or of NBS arrangements, on the hydrological responses under various NBS scenarios.

Impacts of spatial variability of rainfall
The impact of spatial variability of rainfall on the hydrological responses of each NBS scenario over the whole catchment was evaluated integrally in terms of the sum of flow in four conduits (highlighted in the right side of Fig. 2). These four conduits are chosen because they collect the runoff from the whole catchment and finally merge into the storage unit representing the 345 outlet of the drainage system. To be more specific, the PEQp and the PEV computed for the first set of modelling experiments (described in Sect. 3.3) are presented in the following section.

Baseline scenario
Before going on, it is important to evaluate the 'baseline' scenario under both distributed and uniform rainfalls, by using the simulations already performed to validate the Multi-Hydro implementation in the Guyancourt catchment. As shown in the 350 hydrographs ( Fig. 10a, b, c), the higher peak flow was generated by the distributed rainfall in EV1 and EV2. Hence, the peak flow ratio computed by comparing distributed rainfall and uniform rainfall is larger than 1 (see the first column of Fig. 12c), but this ratio is around 0.9 in EV3. The reason is that during the largest rainfall peak of EV1 and EV2, the rainfall intensity of all radar pixels in distributed rainfall is higher than those of uniform rainfall. While in EV3, the rainfall intensity of around 30 % radar pixels in uniform rainfall is about 28 mm h -1 higher than that of the distributed rainfall.
According to the SD of the rainfall intensity at the largest rainfall peak of each event (Table 1), the spatial variability of the rainfall intensity of EV2 is more pronounced than that of EV1 and EV3. Accordingly, the PEQp of baseline scenario in EV2 is the highest. Regarding the total runoff volume (Fig. 12b), the PEV of the baseline scenario for the three rainfall events range from 1 % to 3.7 %. Contrary to the PEQp, the PEV of the baseline scenario is not correlated to the SD of the total rainfall depth. 360 For the baseline scenario, it is noticed that the PEQp is more pronounced than PEV for all rainfall events. These results can be explained by the fact that the spatial variability of rainfall intensity at the largest rainfall peak is strong in all three rainfall events, while the spatial variability of the total rainfall depth is relatively uniform. Figure 11 presents the simulated flow of the first set of NBS scenarios under the three distributed rainfall and three uniform 365 rainfall events. The results are generally consistent with the results of the baseline scenario. Indeed, as shown in Fig. 12c, the peak flow ratios between distributed rainfall and uniform rainfall simulations for the four NBS scenarios are larger than 1 for EV1 and EV2, and around 0.8 for EV3 for the reason mentioned in the previous section.

NBS scenarios
As shown in Fig. 12a, the results of PEQp for PP1, RG1, and Combined1 scenarios are generally in agreement with the baseline scenario: PEQp is the lowest for EV1, and the highest for EV2. For these three NBS scenarios, PEQp range from about 8 % to 370 PEQp of each NBS scenario (Fig. 13a) show that PEQp (the uncertainty related to the peak flow) computed for PP1, RG1, and Combined1 scenarios increase simultaneously with the increase of the SD of the rainfall intensity. The results computed for GR1 scenario do not depict the same tendency: PEQp computed for EV3 is higher than those computed for the two other events.
The reason is related to various factors. Namely, it may be affected by the intersection effects of the spatial variability of 375 rainfall and the spatial arrangement of green roofs in the catchment. As already mentioned, in EV3, the rainfall intensity of 30 % radar pixels in uniform rainfall is about 28 mm h -1 higher than that of the distributed rainfall. Coincidentally, in the GR1 scenario, the green roofs are mainly implemented on the locations with low distributed rainfall intensities. As demonstrated by many previous studies (Qin et al., 2013;Palla and Gnecco, 2015;Ercolani et al., 2018), GR are usually more effective for low rainfall intensity. In the case of the GR1 scenario under the distributed rainfall of EV3, GR measures stored more runoff 380 than in the uniform rainfall during the main rainfall peak. This enlarges the variability of the hydrological response in terms of peak flow.
Regarding the percentage errors on total runoff volume, it is noticed that the computed PEV are lower than 6 % for all NBS scenarios under the three rainfall events, especially in EV3, where they are lower than 2 %. This demonstrates that the resulting uncertainty on the total runoff volume is little influenced by the spatial variability of the rainfall. The reason is that the spatial 385 variability of total rainfall depth is less pronounced with respect to the spatial variability of the rainfall intensity, and also there is no highly localized storm cell in studied events. As illustrated in Fig. 13b, the relationship between the SD of total rainfall depth and the PEV of NBS scenarios is nonlinear. This can be explained by the fact that the three rainfall events are relatively long, and the hydrological performances of NBS are gradually changed during the event (e.g. they can efficiently infiltrate or store water at the beginning, and be saturated after a long rainfall period). Comparing the PEV of each NBS scenario for all 390 three rainfall events (Fig. 12b), those computed for GR1 and Combined1 appear to be the highest for EV2. It could be also related to the intersection effects of spatial location of GR measures and the spatial variability of rainfall. Indeed, these GR measures (considered in the GR1 and Combined1 scenarios) are mainly located in the north side of the catchment. In this area, the first distributed precipitation of EV2 (1-3.5 h), is relatively weak and variable (i.e., there is no rainfall or the rainfall with very low intensity in some localization pixels). Furthermore, as the initial moisture condition of GR measures are considered 395 as unsaturated in both distributed and uniform rainfall, the GR measures are more efficient at the beginning of the distributed rainfall than in the uniform rainfall, and finally enlarge the uncertainty associated with precipitation variability (i.e., the corresponding PEV). More discussion about the intersection effects is presented in Section 4.3.

Impacts of the spatial distribution of NBS
In order to analyze the impacts of the spatial distribution of NBS on the hydrological responses of NBS scenarios, the results 400 of the second set of modelling experiment (described in Section 3.3) are presented as follows. As shown in Fig. 14a, the PEQp of all NBS scenarios are lower than 5 %, and the PEV of all NBS scenarios are lower than 7 %, which indicates that the hydrological responses of NBS scenarios are little affected by the spatial distribution of NBS in the catchment. This result is generally consistent with the observation of Versini et al., (2016), who pointed out that the impact of the spatial distribution https://doi.org/10.5194/hess-2020-468 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License. of green roofs on the catchment response is minimal. However, comparing the PEQp of each NBS scenario, those computes 405 for PP and GR scenarios range from about 2 % to 5 %, which are slightly higher than those related to other scenarios, especially for EV1 and EV3. The reason can be explained by two factors: (i) the infiltration or detention capacity of PP and GR measures are less effective for rainfall characterized by strong intensity and long duration (Qin et al., 2013;Palla and Gnecco, 2015), whereas the RG measures are artificial depressed green areas (simulated with a 0.3 m substrate) with higher retention capacity (Dussaillant et al., 2004); (ii) the differences of DF (large scale; i.e., the second regime) between PP1 and PP2 scenarios as 410 well as between GR1 and GR2 scenarios are larger than that of the other NBS scenarios (Table 3). Figure 15a shows the difference of DF between the same types of NBS scenarios is proportional to the corresponding PEQp. It is found that the larger the difference of DF, the higher the PEQp is. Regarding the PEV of NBS scenarios for the three uniform rainfall events (Fig.   14b), those comparing PP1 and PP2 scenarios (which ranges from about 4 % to 7 % for the three rainfall events, especially higher for the two strong and long events) are slightly higher than those related to the other scenarios. Because porous 415 pavements are infiltration-based measures (they do not retain water for a limited period), their performances are more related to the heterogeneity of their performed location. Namely, some PP measures implemented in drained areas may suffer more from surface runoff, are therefore more easily saturated (see Fig. 5 for a comparison of the spatial arrangement of PP measures for two PP scenarios). As shown in Fig. 15b, the difference of DF between the same types of NBS scenarios has a moderate positive correlation (r =0.61) with the corresponding PEV. Our study hypothesizes that this rather weak correlation is related 420 to the complexity of rainfall with several peaks and dry periods, the retention/infiltration capacity of NBS changes with the rainfall intermittency.

Intersection effects of spatial variability of rainfall and spatial arrangement of NBS
In the following, we present the results of the third modelling experiment set described in Sect. 3.3. The aim is to analyse the potential intersection effects of spatial variability of rainfall and spatial distribution of NBS on the hydrological responses of 425 NBS scenarios.
The resulting uncertainty on the peak flow and total runoff volume (PEQp and PEV) of the third set of modelling experiments are shown in Fig. 16. Firstly, we found that the spatial variability of rainfall has a certain extent impact on the peak flow of each scenario, with the PEQp ranging from about 8 % to 17 %. With the exception of GR1, all the NBS scenarios have a similar tendency: the PEQp are the lowest for the first event, and the highest for the second one. Namely, for most of NBS scenarios, 430 the PEQp (uncertainty on peak flow) increases with the increase of the spatial variability of rainfall intensity. As shown in Fig.   16a, comparing the PEQp between scenarios of PP1 and PP2, RG1 and RG2, as well as Combined1 and Combined2 for the three rainfall events, the maximum difference is less than 3 %. However, comparing the PEQp between GR1 and GR2, the difference is larger, especially in EV3 (> 5 %). For the GR1 scenario, PEQp range from about 8 % to 17 % in all three rainfall events, and those of GR2 range from about 10 % to 15 %. Furthermore, for GR1, the largest PEQp is in EV3, but for GR2, the 435 largest PEQp is computed for EV2. The difference of PEQp between GR1 and GR2 scenarios demonstrated that the spatial variability of rainfall and the spatial arrangement of GR measures have some intersection effects on the peak flow of GR https://doi.org/10.5194/hess-2020-468 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License.
3. The fractal dimension DF appears as a useful tool to quantify the spatial heterogeneity of NBS across a range of scales.
The DF of each NBS scenario is associated with the urbanization level of the catchment, which confirms that the level of implementation of NBS is reasonable to match the catchment conditions. The fractal dimension combined with the fullydistributed model is an innovative approach that is easily transportable to other catchments.
4. The spatial distribution of rainfall and the spatial arrangement of NBS have intersection effects on the hydrological 475 responses of NBS scenarios, especially significant for the peak flow of GR scenarios (with a maximum difference between the scenario of GR1 and GR2 reaching about 5 % on peak flow). The intersection effects on the total runoff volume of each NBS scenario is quite variable because the chosen NBS present some limitations in terms of infiltration or detention capacity during a long rainfall event with high intermittency. However, the RG scenarios appear to be less affected by the intersection effects, with a difference lower than 3 % on peak flow and lower than 1 % on total runoff volume. 480 5. The fully-distributed model takes into account the spatial rainfall variability and the spatial distribution of NBS at very small scale, which presents the uncertainty of hydrological performances of NBS scenarios (e.g. peak flow and total runoff volume) is related to these two factors. Whereas the lumped and semi-distributed models are not able to consider these heterogeneities, and the level of spatial resolution incorporate in this type of models may result in some critical heterogeneous information be missing that can influence model predictions. 485 In our specific case, the GR scenarios are more sensitive to the spatial variability of rainfall and the spatial arrangement of GR measures, while the performances of RG scenarios and combined scenarios are more stable under any condition. Apparently, these findings already give some incites to decision-makers on Why they need to prioritize given NBS within the urban planning process.
However, larger rainfall samples, including extreme rainfalls, will be needed to provide a better knowledge towards the 490 somewhat universal solutions and to give an answer on How to prioritize these NBS. Under this perspective, the obtained results already demonstrated that new scale-independent indictors, like the fractal dimension DF applied in this study, will be essential for more profound quantitative evaluation of the diversity of combined impacts, including for other heterogeneous catchments.

495
Data availability. Datasets for this study are under preparation for public release by the Chair Hydrology for Resilience Cities.
They can already be freely obtained via contact.hmco@enpc.fr.           Percentage error on total runoff volume between the same types of NBS scenarios under the three uniform rainfall events.