A look at the links between drainage density and flood statistics

We investigate the links between the drainage density of a river basin and selected flood statistics, namely, mean, standard deviation, coefficient of variation and coefficient of skewness of annual maximum series of peak flows. The investigation is carried out through a three-stage analysis. First, a numerical simulation is performed by using a spatially distributed hydrological model in order to highlight how flood statistics change with varying drainage density. Second, a conceptual hydrological model is used in order to analytically derive the dependence of flood statistics on drainage density. Third, real world data from 44 watersheds located in northern Italy were analysed. The three-level analysis seems to suggest that a critical value of the drainage density exists for which a minimum is attained in both the coefficient of variation and the absolute value of the skewness coefficient. Such minima in the flood statistics correspond to a minimum of the flood quantile for a given exceedance probability (i.e., recurrence interval). Therefore, the results of this study may provide useful indications for flood risk assessment in ungauged basins.


Introduction
Drainage density (D d ) was defined by Horton (1945) as the ratio of the total length of streams in a watershed over its contributing area. It describes the degree of drainage network development and was recognised by many authors to be significantly effective on the formation of flood flows (see, for instance, Gardiner and Gregory, 1982). D d is higher in arid areas with sparse vegetation cover and increases with increasing probability of heavy rainstorms (Gregory, 1976; Correspondence to: A. Montanari (alberto.montanari@unibo.it) Woodyer and Brookfield, 1966). D d is also higher in highly branched drainage basins with a relatively rapid hydrologic response (Melton, 1957).
The drainage density exertes on flood peaks significant controls which can be broadly divided between direct and indirect effects (Merz and Blöschl, 2008;. Among the most significant direct effects there is the control associated with the length of the stream network and hillslope paths. Because flow velocity is higher in the river network, D d significantly affects the concentration time and therefore the peak flow magnitude. It follows that an increasing drainage density implies increasing flood peaks. Moreover, a long concentration time implies more opportunities for water to infiltrate. Therefore a decreasing D d generally implies decreasing flood volumes. Among the indirect effects there are those that can be ascribed to the role of D d as an index of geology. A reduced D d can be attributed to the presence of impervious rocky hillslopes, and therefore reduced storage volumes and high flood peaks. But a low D d may also be due to the presence of karstic areas, highly weathered bedrock, and/or highly permeable fluvial deposits in the valley floors, which all may imply large storage volumes and response times and hence small flood peaks and volumes. Another indirect control may be given by the interaction of landform evolution, soil formation, erosion and floods (driven by climate and modulated through geology). Over time scales of centuries, large floods may shape catchments through a positive feedback loop to increase topographic gradients, increase D d and decrease storage volumes which in turn may increase flood peaks and volumes (Merz and Blöschl, 2008). Finally, a significant indirect control, that can be very effective in semiarid areas, is that exerted by D d through its connection with the vegetation cover. Bare soils are much prone to soil erosion and therefore characterised by high drainage density and high runoff production, therefore implying large flood peaks and volume. This effect is clearly visible in watersheds of central Italy, where arid badlands are in fact characterised by significantly larger drainage density (Grauso et al., 2008).
The scientific literature proposed numerous contributions dealing with the links between D d , geomorphological behaviours of the basin, climate and river runoff regime (see, for instance, Morgan, 1976;Gurnell, 1978;Rodriguez-Iturbe and Escobar, 1982). In our view, significant contributions within this respect were provided by Murphey et al. (1977) who analysed the relationship among flood hydrograph characteristics, like rise time, base time, mean peak discharge and flood volume, and basin parameters such as, among the others, contributing area and drainage density. Plaut Berger and Entekhabi (2001) and Humbert (1990) proved that drainage density is significantly related to the basin runoff coefficient. Gresillon (1997) studied 100 African catchments and concluded that D d and watershed area are two very important variables to explain the shape of the recessing limb of the hydrograph, while in tropical regions D d is effective on the runoff coefficient during flood events. Yildiz (2004) carried out a numerical simulation showing that streamflows simulated by a rainfall-runoff model steadily increase with increasing D d . Recently, Grauso et al. (2008) showed that D d is a fundamental explanatory variable in regression models for estimating the annual amount of suspended sediment yield. The explaining capability of D d was found to be extremely relevant, therefore suggesting that a significant link indeed exists between D d and the river flow regime. A recent and interesting review is provided by Wharton (1994) who studied the usefulness of D d in rainfall-runoff modelling and runoff prediction.
Nowadays, there is an increasing interest in the potential role of geomorphic parameters to effectively describe and represent some peculiar behaviours of the rainfall-runoff transformation. Such interest is motivated by the increasing attention that the hydrologic community is dedicating to the problem of prediction in ungauged basins (PUB). Many researchers are involved in the activity of the PUB initiative, which was launched in 2003 by the International Association of Hydrological Sciences (Sivapalan et al., 2003). Furthermore, the continuous development of remote sensing techniques makes the estimation of a number of geomorphic parameters more accessible (Jena and Tiwari, 2006).
Within the context of PUB, it is extremely important to investigate the potential utility of geomorphic parameters in the identification and calibration of hydrological models. Such parameters can be useful "hydrological signatures", that is, large scale markers of intrinsic local scale hydrological processes. While the connection between geomorphic parameters and river runoff regime has long been recognised (see the brief review presented above), expressing such link in a quantitative form is still an open problem.
The analysis presented here aims at shedding more light on the links between drainage density and streamflow regime. In fact, we believe the topology of the river network is poten-tially able to convey much more information on the formation of the river flows than what is currently known. In particular, the present analysis focuses on the potential interrelationships between D d , the mean (µ), coefficient of variation (CV ) and coefficient of skewness (k) of the annual maximum peak flow. Given that µ, CV and k summarise the frequency regime of flood flows, the possibility to derive indications about their value in ungaunged basins from catchment and river network descriptors (such as D d ) is extremely important in the context of PUB.
A major problem that limits the possibility to assess the dependence of flood statistics on D d is the scarcity of historical data that makes the estimation of the statistics themselves highly uncertain. Moreover, a lot of additional forcings, other than D d , are effective on the rainfall-runoff transformation, like the climate, the geology of the contributing area, the presence of preferential flow and many others. Therefore, the main objective of our study, i.e. assessing the explanatory capability of D d , is a difficult task which is complicated by the presence of a significant noise. Ideally, we would like to assess under the same climatic conditions how the flood frequency regime of a given basin changes when changing the D d . This is clearly not possible when dealing with real world case studies, for which D d , as well as the flood frequency regime, is a direct expression of all forcings (e.g., climatic, geological, etc.) that characterise a given catchment.
To overcome the above problem, we first focus on an ideal case study of a river basin for which we generate arbitrarily long time series of synthetic river flow data. We perform many different simulations for different values of D d , while keeping all other forcings unchanged. This first numerical study enables us to obtain a first indication on the effects of D d on flood statistics and points out the possible existence of a optimal value of D d for which the minimum CV and k are attained. To elaborate a sound physical interpretation of this result, we provide an analytical derivation of the relationship between D d and CV . The analytical derivation is based upon a number of simplifying assumptions, but proves indeed that (1) the pattern of CV that emerges from the numerical study is feasible and (2) an optimal value of D d exists under certain hypotheses. Finally, we assess and discuss the results obtained through the simulation study against empirical evidences for a set of river basins located in northern Italy.
We believe that recognising the pattern of the D d −CV and D d −k relationships is an important issue. In fact, if such patterns were known, and the critical value for D d was confirmed and quantified, one could estimate to what extent a river basin is prone to extreme floods on the basis of its drainage density. Our study is not conclusive in defining the patterns above in a quantitative way, but provides a first indication about their possible shape.

Links between the drainage density and flood statistics -a numerical simulation through a spatially distributed rainfall-runoff model
The links between the drainage density and the statistics of the annual maximum peak flows are first inspected through a numerical simulation. With reference to the Riarbero River basin, which is a right tributary to the Secchia River (Italy) with a drainage area of 17 km 2 , we simulate many 100-year long time series of hourly river flows for different D d values, setting all of the other external forcings and model parameters constant. River flow simulation is carried out by using a spatially distributed rainfall-runoff model, while the meteorological input is obtained by generating synthetic rainfall and temperature series at hourly time scale.
The rainfall-runoff model applied here is AFFDEF (Moretti and Montanari, 2007). It determines the catchment hydrologic response by composing the two processes of hillslope runoff and channel propagation along the river network. A description of the structure of the model is provided in the Appendix. AFFDEF is based on the application of a conceptual model for the continuous time simulation of the infiltration and runoff formation processes at local scale, as well as energy and mass balance concepts for simulating runoff concentration and propagation. Therefore it is capable of fully simulating the direct controls of D d on flood statistics. It can also account for the indirect controls exerted by D d , by using information at local scale about watershed morphology, soil type and vegetation cover to compute runoff formation and propagation (Moretti and Montanari, 2007).
AFFDEF has been calibrated for the Riarbero River basin during a previous study (see Moretti and Montanari, 2008). For performing the hydrological simulation, a 100-year hourly rainfall time series for the Riarbero River basin was first generated by using the Neyman-Scott rectangular pulses model. The model represents the total rainfall intensity at time t as the sum of the intensities given by a random sequence of rain cells active at time t. Extensive details about the model and the simulation we performed can be found in Cowpertwait (1996), Rodriguez-Iturbe et al. (1987), Burlando (1989) and Moretti and Montanari (2008).
Then, a 100-year record of hourly temperature data was generated by using a stochastic model, namely, a fractionally differenced ARIMA model (FARIMA). This kind of model has been shown in many applications to be able to well fit the autocorrelation structure of temperature series which, for increasing lag, is very often affected by a slow decay that may suggest the presence of long term persistence. FARIMA models, which are characterised by a high flexibility in their autocorrelation structure, are capable of fitting long term persistence by means of the fractional differencing operator. More details on FARIMA models and the simulation procedure herein applied for the temperature data can be found in Montanari (2003) and Montanari et al. (1997).
The synthetic rainfall and temperature data have been subsequently routed through AFFDEF, therefore obtaining a 100-year long sequence of river flows relative to Riarbero River outlet. The simulation was performed for varying values of the critical support area A 0 that is used to distinguish rill flow from channel flow (see Appendix A), therefore obtaining basin configurations that correspond to different D d values, while keeping any other forcing unchanged. In particular, the D d values were computed as the ratio between the total length of the river network and the drainage area. The total length of the river networks depends on A 0 and was automatically retrieved from a Digital Elevation Model (DEM) of the catchment, as described in the Appendix A. Then, for each simulation (i.e., each D d value) we extracted the annual maximum river flows and computed the sample values of the following statistics: µ, σ , CV and k. Figure 1 shows the patterns of µ, σ , CV and k depending on D d . For the sake of generality, Fig. 1 reports standardised values of D d , CV and k. Values of CV and k were standardised by their empirical minima while the D d values were standardised by the D d of the Riarbero river basin (0.54 km −1 ), which was computed as the ratio between the total length of the stream network, as described by the blue lines on 1:25 000 topographic maps, and the drainage area of the catchment. The results show that µ and σ increase for increasing D d . This outcome was expected, because D d significantly affects the concentration time. As a matter of fact, the surface flow is much slower on the hillslopes than in the river network and therefore higher D d values imply lower concentration times, which in turn imply amplified peak river flows because of the fast response of the catchment.
A more interesting result is obtained by looking at the patterns of CV and k. In fact, the simulation provided by AFFDEF supports the existence of an optimal value of D d for which a minimum is attained for the above statistics. This outcome has a significant implication: the optimal value of D d would correspond to a situation for which the river basin is less prone to heavy floods, while significant departures from the optimum would be signatures of a higher flood risk.
Of course the numerical simulation is affected by uncertainty and therefore its outcome may not be fully representative of what actually occurs in nature. Therefore one may wonder whether or not the patterns found for CV and k are realistic. In order to better investigate the possible physical explanations for the critical value of D d we carried out an analytical derivation of the D d −CV relationship, that is presented in the following Section.

Analytical derivation of the links between the drainage density and flood statisticsa conceptual approach
The statistics of the flood flows can be analytically computed by means of a derived distribution approach. This analysis is meaningful in order to confirm and better understand the patterns highlighted by the numerical simulation. In order to make the analytical computation possible, simplifying assumptions need to be introduced in the schematisation of the rainfall-runoff transformation.

Gross rainfall
We assume that the extreme mean areal gross rainfall intensity r l is constant during the storm event and is described by a scale-invariant intensity-duration-frequency (IDF) curve of the type r l =a(T )d n−1 , where n is fixed, a(T ) is a coefficient that depends on the return period T and d is the storm duration. For more details about scale invariant IDF curves see Burlando and Rosso (1996). It is also assumed that a(T ) is distributed according to a Gumbel distribution (Kottegoda and Rosso, 1998) and therefore can be estimated through the relationship where µ * and σ * are mean and standard deviation of the annual maximum rainfall for storm duration of one hour (Burlando and Rosso, 1996).

Net rainfall
It is assumed that the mean areal net rainfall intensity can be estimated as where r p is the mean areal infiltration rate during the storm. We adopted Horton's equation for infiltration (Horton, 1933) to describe r p , which provides a reasonable schematisation for the Riarbero River basin, where steep slopes with a reduced permeability dominate. If one assumes that r l >r p,0 , where r p,0 is the initial infiltration rate, the average rate of rainfall loss during the storm event can be computed by with r p,c and denoting the asymptotical infiltration rate and a parameter, respectively.

Rainfall-runoff
We assume that the rainfall-runoff transformation can be schematised by using a unit hydrograph rainfall-runoff model. Accordingly, under the assumption of constant rainfall introduced above, the flood hydrograph can be estimated through the linear relationship for t≤t 0 +d, where t 0 is the initial time of the hydrograph, A is the drainage area and h(t−τ ) is the travel time distribution of the catchment, which is supposed to be a stationary function of time. Let us denote with t c the concentration time of the catchment. We assume that the travel time distribution can be written in the form where c is a constant parameter.

Computation of the annual maximum flood statistics
Finally, we assume that the flood flow with recurrence interval T , Q(T ), is induced by a rainstorm event with return period T , duration d=t c and net rainfall intensity given by Eq.
(2). Under the assumptions introduced in Sect. 3.1.1, 3.1.2 and 3.1.3, the T -year flood Q(T ) can be written in the form while µ, σ and CV of annual maximum floods can therefore be written as:

Link between the annual maximum flood statistics and the drainage density
We suppose that the main impact of D d on the frequency regime of annual maximum flood, summarised by the statistics computed in Sect. 3.1.4, is due to its monotonically decreasing relationship with the concentration time t c of the catchment. We assume this relationship to be expressed by where M is constant. The simple conceptual model introduced above can explicitly account for the direct controls of D d on the flood statistics. In fact, Eq. (8) assumes the existence of a direct link between the drainage density and the concentration time, and therefore the duration of the critical storm (direct controls). Soil use, basin morphology and vegetation cover (indirect controls) are considered in a lumped form and only implicitly, through the parameter . It is worth remarking that we assumed in this study that the return period of the annual maximum flood coincides with the return period of the rainfall event generating the flood. This assumption is rarely providing a satisfactory schematisation in real world cases (Viglione and Blöschl, 2009). Nevertheless, it was deemed appropriate for the scope of this study as it simplifies the casual relationship between the rainfall forcing and the induced flood.  The results show that the patterns obtained through the hydrological simulation are confirmed. In particular, the analytical derivation confirms the presence of an optimal value for D d with regard to CV , which is due to the interplay between the depth-duration-frequency curve for rainfall and Horton's infiltration equation.

Results of the analytical computation
Of course, the outcome of the conceptual approach is strongly affected by the underlying assumptions. In particular, it is worth noting that a different relationship between D d and t c (see Eq. 8) would affect the patterns in Fig. 2. Nevertheless, the monotonically decreasing relationship between D d and t c assures the existence of an optimal D d regardless of the mathematical expression assumed to describe this relationship. Therefore, the hypothesis of describing the relationship between D d and t c through (Eq. 8) could be relaxed and still an optimal D d could be found under the remaining assumptions.
Furthermore, it is interesting to note that the simple conceptual model that was used to perform the analytical derivation, with the parameters given above, would have provided a constant value of CV =0.333 over an arbitrary range of drainage densities if the rainfall losses were neglected. The adoption of Horton's equation for explaining the losses along with the other simplifying assumptions mentioned above results in the presence of the optimal value of D d .
Finally, it is important to observe that the analytical model described above, in view of its linear structure, leads to obtaining k=k * for any value of D d , where k * is the skewness of the annual maximum rainfall corresponding to storm duration of one hour. Therefore, the analytical model cannot explain the pattern simulated by AFFDEF for the third order moment. There are many possible reasons for the presence of a critical value of D d with respect to k, that nevertheless imply relaxing some of the assumptions underlying the conceptual framework introduced above, with the consequence that an analytical derivation could be no longer possible. For instance, if one keeps the assumption of linear rainfall-runoff model, the variations of k depending on drainage density could be explained by an analogous variation of k * over dif- For what concerns the AFFDEF simulation the increase of k for increasing drainage density is dominated by a corresponding increase of the rainfall skewness, while the causes for the increase of k for very low drainage densities are less clear. We believe the complexity of the spatially distributed infiltration pattern dominates in this case as a consequence of the increased residence time of water in the hillslopes.

Links between the drainage density and flood statistics -case studies.
The patterns identified for the D d −CV and D d −k relationships are displayed by the outcome of mathematical models that are based on assumptions that were described in detail above. Although the models are based on reasonable conceptualisations, their output is affected by a relevant uncertainty, which is induced by many causes including input uncertainty, model structural uncertainty (also induced by the underlying assumptions) and parameter uncertainty. Therefore, there is no guarantee that the identified patterns actually correspond to what can be observed in the real world. The analysis so far described has been performed through mathematical models in view of the aforementioned impossibility to study the dependence of the hydrological response on the drainage density for a real world watershed. However, an analysis of real world data can be carried out by considering different watersheds characterised by different drainage densities. We believe such type of analysis is necessary in order to provide a consistent support to what has been previously found.
However, the real world analysis is affected by a relevant uncertainty as well, because different catchments have not only different drainage densities, but also different climate, geomorphology and so forth. Therefore the patterns induced by the drainage density could be masked by patterns induced by the additional forcings mentioned above and focusing on different regions may lead to different results (Merz and Blöschl, 2009). Thus, we expect that the patterns obtained through the real world analysis are affected by a significant scatter of the empirical data, which may make the interpretation of the results difficult and to some extent subjective. We do not expect a clear confirmation of the previous outcomes, but at least a result allowing us not to reject their plausibility. In order to limit the ambiguity induced by the additional forcings, it is important to focus on catchment with similar behaviours.

The study area and data set
We focus on 44 subcatchments of the Po River basin with drainage area ranging from 20 to 10 000 km 2 . The Po river flows in northern Italy, with a total contributing area of about 70 000 km 2 . The main stream is about 652 km long. The mean discharge at its mouth is around 1560 m 3 /s while the mean annual inflow into the Adriatic sea is around 47 billions m 3 . The river regime in the alpine part of the watershed is sensitive to temperature, with an important contribution by snowmelt to runoff during summer. For the tributaries from the Apennines, rainfall is the most significant contribution to river runoff. Mean annual rainfall over the Po river watershed is around 1200 mm/y, with minimum and maximum values around 600-800 mm/y and 1600-1800 mm/y, respectively. The rainfall regime is predominantly continental over the Alpine portion of the basin, with a maximum during summer and a minimum during winter. The remaining portion of the basin is mainly characterized by two rainfall maxima in Spring and Autumn (Autorità di Bacino del Fiume Po, 1986).
Annual maximum series (AMS) of river flows are available for the considered 44 catchments for the period between 1918 and 1970. These data were collected by the Italian Hydrometric Service. The series are affected by missing values, so that the sample size of the AMS ranges between 8 and 44, with an average size of 23. For each series the sample mean, standard deviation, CV and k were computed.

Computation of the drainage density
For the 44 considered catchments, D d is calculated as the ratio between the total length of the stream network, as described by the blue lines on 1:25 000 topographic maps, and the drainage area of each catchment. Figure 3 shows the pattern of the computed D d over the Po River watershed. Figure 4 shows the progress of reduced µ and σ depending on D d for the considered subwatersheds of the Po River Basin. Reduced µ and σ are considered instead of µ and σ to remove the dependence on drainage area and to enable one to compare Fig. 4 with the corresponding diagrams of Figs. 1 and 2. The reduced values of the moments are obtained by dividing the empirical moments by A β , where β is a regional scaling exponent. Empirical values of β vary from 0.2 to 0.8 for various regions of the world (see e.g., Castellarin et al., 2005;Castellarin, 2007;Gaume et al., 2009) and a value of 0.2 appears to best suit the considered study area. Figure 4 also reports the progress of CV and k as a function of D d . In order to aid the interpretation of the plots, a moving average (window length=9) of the statistics is also shown, along with a linear interpolation for µ and σ .

Results
As expected, Fig. 4 shows that empirical data present a significant scatter, therefore making the interpretation of the results highly uncertain and the identification of patterns subjective. However, the moving average curves reported in the diagrams do not contradict what was inferred through the previous numerical and analytical investigations. Although the linear regression models reported in panels a and b of Fig. 4 are not statistically significant, we decided to include them in the figure to better illustrate the general increasing tendency of reduced µ and σ with D d . Also it seems useful to remark here that the progression of µ and σ in Figs. 1 and 2 can be effectively represented through a power law (i.e., linear model of the log-transformed variables) and that a linear regression of the log-transformed case study data is significant at the 10% level for both reduced µ and σ .
The considerations on the analysis of real world basins and their flood statistics enable us not to reject the hypothesis of the existence of a critical value of the drainage density with respect to CV and k. Although this conclusion could certainly be questioned in view of the above mentioned scatter, the practical effect of the detected critical value is evident if one looked at the progression of the dimensionless flood quantile (DFQ) as a function of D d . The DFQ is the ratio between the flood quantile and the corresponding mean value of the annual maximum flood. If one assumes that CV and k are given by the moving averages shown in Fig. 4, the DFQ can be easily computed by fitting a Generalised Extreme Value (GEV) distribution (Kottegoda and Rosso, 1998) to the available dimensionless flood data. The GEV distribution has been widely shown to satisfactorily reproduce the sample frequency distribution of hydrological extremes (in particular, precipitation depths and flood flows) observed in different geographical contexts around the world (Castellarin, 2007). Figure 5 shows the progression of the 100-year DFQ depending on the corresponding value of D d . Indeed, it can be seen that the variation of DFQ depending on D d is significant and certainly affects the magnitude of the peak flows.
It is worth pointing out that one should not be surprised by the differences in terms of drainage density in Figs. 1, 2 and 4. In fact, D d was estimated by using different information in each application and therefore the obtained values are not directly comparable.

Conclusions
Drainage density is a classical descriptor of catchment morphology which is known to control the formation of river flows. As such it may influence significantly the frequency regime of flood flows. In this study, we addressed the analysis of the relationships between drainage density and flood frequency regime by performing a numerical simulation, an analytical study based on a conceptual hydrological model and an investigation of real world data. We found that we cannot reject the hypothesis of the existence of a critical value of the drainage density for which a minimum is attained in the coefficient of variation and in the absolute value of the skewness coefficient for the annual maximum sequences of flood flows. This finding is relevant from a practical viewpoint. A minimum for the flood statistics mentioned above may imply a more pronounced minimum in the estimated design-flood (the flood flow associated with a given exceedance probability, or recurrence interval). Therefore drainage density could provide interesting indications about the flood risk for an ungauged river basin.
Previous studies (Blöschl and Sivapalan, 1997) found that the coefficient of variation of annual peak flows seems to reach a maximum at a certain threshold area of the upstream river basin, while we found here a minimum at a certain threshold drainage density. These findings are apparently in disagreement if one assumed that the concentration time is the main control in both cases. However, it is important to note that Blöschl and Sivapalan (1997) found the maximum of the CV for a catchment area of about 100 km 2 . In this study we showed that a the progression of CV with D d may exhibit a minimum independently of drainage area (see Sects. 2 and 3). Concerning the case study (see Sect. 4), we found that the critical value of D d , i.e. ∼0.26 km −1 (see Fig. 4), characterizes a group of basins with highly variable drainage areas, which span approximately over the entire range of areas of the case study. Therefore the above results are not in disagreement, they are rather complementary to each other and likely to postulate that the progress of CV with the basin concentration time might be fluctuating. We  postulate that such progress is governed by the interplay between the depth-duration-frequency relation for rainfall and the infiltration curve and we are currently investigating this issue in more detail.
Our study is an initial effort to understand how descriptive the drainage density is in the representation of the flood frequency regime for a given basin. For this reason we did not provide in this study any quantitative indication for identifying the critical value of the drainage density for a single ungauged basin. At this stage of our knowledge, we believe that an extended set of comparative evaluations are needed in order to gain additional insights. These further investigations should be performed for climatically and geologically homogeneous areas, at least, in order to keep unchanged the most significant forcings on the flood statistics other than the drainage density.
Finally, we would like to remark once again that the flood statistics are related to many other forcings other than drainage density like, for instance, climate, geomorphology and so forth. The probability distribution of the peak flows is actually the result of many influencing processes and catchment behaviours and the drainage density is not the most effective in general. Therefore the patterns detected here might be not there in other contexts.

Appendix A AFFDEF: a spatially distributed rainfall-runoff model
AFFDEF discretises the basin in square cells coinciding with the pixels of the Digital Elevation Model (DEM). The river network is automatically extracted from the DEM itself by applying the D-8 method (Band, 1986;Tarboton, 1997), which allows one to estimate the flow paths and the contributing area to each cell. In detail, the network determination is carried out by first assigning to each DEM cell a maximum slope pointer and then processing each cell in order to organise the river network. Digital pits are filled in a preprocessing step, before extracting the channel network from the DEM of the catchment. Each cell receives water from its upslope neighbours and discharges to its downslope neighbour. For cells of flow convergence, the upstream inflow hydrograph is taken as the sum of the outflow hydrographs of the neighbouring upslope cells.
Distinction between hillslope rill and network channel is based on the concept of constant critical support area (Montgomery and Foufoula-Georgiou, 1993). Accordingly, rill flow is assumed to occur in each cell where the upstream drainage area does not exceed the value of the critical support area A 0 , while channel flow occurs otherwise. The value of A 0 conditions the drainage density of the river basin. Decreasing values of A 0 induce an increase of the linear extension of the drainage network and therefore a corresponding increase of D d .
The interaction between soil, vegetation and atmosphere is modelled by applying a conceptual approach. The model firstly computes the local gross rainfall P l [t, (i, j )], for each DEM cell of coordinates (i, j ), by interpolating the observations referred to each raingauge through an inverse distance approach. Then, for each cell a first rate of rainfall depth is accumulated in a local reservoir (interception reservoir) which simulates the interception operated by the vegetation. The capacity of such interception reservoir is equal to C int S(i, j ), being C int a parameter, constant in space and time, and S(i, j ) the local storativity. The latter is computed depending on soil type and use according to the Curve Number method (CN method, Soil Conservation Service, 1987;Chow et al., 1988).
Once the interception reservoir is full of water, the exceeding rainfall reaches the ground. Then, surface and subsurface flows are computed according to a modified CN approach that is able to simulate the redistribution of the soil water content during interstorm periods. In detail, it is assumed that a linear reservoir (infiltration reservoir), which collects the infiltrated water, is located in correspondence of each DEM cell at the soil level. The local surface runoff and the infiltration are computed according to the relationship where P [t, (i, j )] is the rainfall intensity that reaches the ground at time t, P n [t, (i, j )] is the intensity of surface runoff, F [t, (i, j )] is the water content at time t of the infiltration reservoir located in correspondence of the cell (i, j ), and H ·S(i, j ) is the capacity of the infiltration reservoir itself, computed by multiplying the calibration parameter H by the soil storativity previously introduced. The quantity represents the intensity of the infiltrated water. The outflow W [t, (i, j )] from the infiltration reservoir to the sub surface river network, which is assumed to coincide with the surface one, is given by the linear relationship where H S is a calibration parameter. H and H S are assumed to be constant with respect to both space and time.
The hourly intensity of potential evapotranspiration E P [t, (i, j )] is computed at local scale by applying the radiation method (Doorenbos et al., 1984). When some water is stored in the interception reservoir, the effective evapotranspiration E [t, (i, j )] is assumed to be equal to E P [t, (i, j )] and is subtracted from the water content of the interception reservoir itself. When this latter is empty, or is emptied while subtracting the evapotranspiration rate, the remaining part of E P [t, (i, j )] is subtracted from the water content of the infiltration reservoir. In this case, it is assumed that E [t, (i, j )] is varying linearly from 0 when F [t, (i, j )]=0, to E P [t, (i, j )] when F [t, (i, j )] =H ·S(i, j ). Evapotranspiration is the only source of water losses in the model, which primarily depends on the capacity of the interception reservoir and hence on the parameter C int . Therefore a first estimation of C int can be obtained by comparing observed and simulated runoff coefficients. Figure A1 shows the scheme of the interaction among soil, vegetation and atmosphere operated by the model.
The continuity equation applied to the infiltration reservoir can be written as: By combining Eqs. (A1), (A2) and (A3) and taking the effective evapotranspiration into account, the mass balance equa- Surface and sub surface flows are propagated towards the basin outlet by applying the variable parameters Muskingum-Cunge model. Extensive details can be found in Cunge (1969) and Orlandini et al. (1999) for surface and sub surface propagation, respectively. For the surface flow, the kinematic celerity is computed by considering rectangular river cross section with fixed width/height ratio. The latter parameter and the channel roughness can assume different values along the river network and on the hillslopes. In particular, the channel roughness in the river network is allowed to vary from a minimum to a maximum value depending on the contributing area. For the subsurface flows, the kinematic celerity is instead computed as a function of the saturated hydraulic conductivity of the soil. It is interesting to note that the model describes in a simplified manner the dynamics of the sub surface flows. In particular, it does not distinguish between near surface and deep water flow, and assumes that the calibration parameters H and H S are constant with respect to both space and time. This simplified description has been used in order to reduce the number of model parameters and, consequently, the amount of historical data required for model calibration. On the other hand, one may expect a significant approximation in the simulation of the low river discharges, especially when referring to highly permeable basins.
Moreover, the formation of the surface runoff is modelled according to a scheme that is very similar to the one adopted by the CN method, which is considered by many authors as an infiltration excess approach (Beven, 2000). Therefore one may expect that the proposed model is better suited for basins characterised by low permeability and prevalently impervious hillslope, where the surface runoff is more likely to be given by excess of infiltration instead of excess of saturation.