HydroSCAPE : a multi-scale framework for streamflow routing in large-scale hydrological models

Introduction Conclusions References


Introduction
The emerging of socio-hydrology as a research area of hydrological sciences, addressing multiple scale feedbacks between human activities and hydrological processes, fostered new developments in large-scale hydrological models as a component of more Figures

Back Close
Full comprehensive, and complex, Earth System Models (ESMs).Manabe (1969) was the first to add a land component in a climate model.His work prompted the development of a first generation of Global Circulation Models and a parallel development of Land Surface Models (LSMs) for hydrological applications (see Haddeland et al. (2011) and Prentice et al. (2015) for a comprehensive review).Land Surface Models are developed with the intent of providing a realistic and detailed representation of vertical water, energy and CO 2 fluxes, with the perspective to facilitate coupling with atmospheric models.On the other hand, Large-scale Hydrological Models (LHMs) have been developed with the perspective of a realistic representation of water resources and horizontal water transfer (Haddeland et al., 2011).VIC (Liang et al., 1994), LaD (Milly and Shmakin, 2002), H08 (Hanasaki et al., 2008a, b), Noah-MP (Niu et al., 2011), WEHY-HCM (Kavvas et al., 2013), and CLM (Oleson et al., 2013), are examples of LSMs, while MacPDM (Arnell, 1999), WBM (Vörösmarty et al., 1998), WGHM (Döll et al., 2003), WASMOD-M (Widén-Nilsson et al., 2007), PCR-GLOBWB (Van Beek et al., 2011), LISFLOOD (van der Knijff et al., 2010), and mHM (Samaniego et al., 2010) can be classified as belonging to the category of hydrological models.This classification notwithstanding, the boundary between these two categories is blurred since in their recent developments most of these models are converging to a comprehensive representation of the terrestrial processes, in an attempt to increase realism, though this is often achieved at the expense of reliability and robustness (Prentice et al., 2015).However, at the current stage of the developments, both category of models suffer from a discretization often too coarse to represent routing in the river network with enough detail to capture geomorphological dispersion (Rinaldo et al., 1991).
The hydrological component of both categories of models, which for simplicity we indicate here as LHMs, rely on simplified conceptualizations and empirical upscaling procedures (Nazemi and Wheater, 2015), when dealing with heterogeneities that characterize hydrological fluxes across a hierarchy of scales, ranging from the hillslope to the catchment and the continent.In addition, most of the available models inherit the grid approach from LSMs, which works fine for the vertical fluxes, but renders grid de-Figures Back Close Full pendent the surface routing.In most cases routing is performed by solving either the kinematic wave or the de Saint-Venant equation (or one of its simplifications) by using the same discretization adopted for resolving the vertical fluxes, thereby leading to scale-dependent inaccuracies in the representation of horizontal fluxes, unless a very fine discretization is used, which is untenable at large scales also for the currently available high performance computational resources.The introduction of improved routing schemes to adequately resolve horizontal fluxes with an acceptable computational effort is therefore indicated as one of the priorities in ESMs, and in LHMs as well, (Clark et al., 2015).In fact, grid-based models are typically applied with a spatial resolution ranging from 0.1 to 0.5 • .This resolution has been proven to be insufficient to capture geomorphological dispersion and travel time distribution at the level of accuracy needed to model horizontal fluxes, in particular at intermediate scales, i.e. scales of the order of thousand or tens of thousand km 2 (Gong et al., 2009;Verzano et al., 2012), which are relevant in modeling flood events.
Hyper-resolution LHMs relying on global digital drainage networks at fine scales, such as HydroSHEDS at 90 m resolution (Lehner et al., 2008), represent a possible strategy to overcome the above limitations and obtain reliable estimates of horizontal fluxes (Wood et al., 2011).However, applying a LHM at a such fine discretization is untenable for the large computational cost associated to it, which becomes unbearable when inversion is applied to infer model parameters from observational data.This burden is currently too high for LHMs adopting explicit hydrodynamic routing through the numerical solution of the mass and momentum conservation equations (i.e. the de Saint-Venant equations), but also for models adopting cell-by-cell routing algorithms based on mass conservation and relationships between river-channel storage and streamflows (Yamazaki et al., 2011;De Paiva et al., 2013).If modeling high flow events is the objective of the analysis a hourly, or smaller, time scale should be adopted, thereby increasing the computational burden, for the same spatial discretization, with respect to the daily or monthly time scale typically adopted in large-scale simulations of water resources.Introduction

Conclusions References
Tables Figures

Back Close
Full Mixed schemes in which routing is separated from runoff have also been employed (see e.g.Gong et al., 2009Gong et al., , 2011;;Wen et al., 2012;Lehner and Grill, 2013).Hydro-ROUT is a vector-based routing scheme fully integrated with ArcGIS (ESRI, 2011) developed by Lehner and Grill (2013), in which streamflow is obtained by routing to the catchment outlet the surface runoff (as provided by an external runoff simulator to be coupled with the routing scheme) accumulated at the nodes of the network.Routing is performed by using a "plug-flow" routing scheme similar to that implemented by Whiteaker et al. (2006) and applied to the HydroSHEDS river network (Lehner et al., 2006).Wen et al. (2012) proposed a multi-scale routing framework employing a pdf distribution for the overland flow path lengths lumping the effect of unmodeled sub-grid variability in the parameters describing the pdf.The kinematic wave routing method is then employed for both overland and channel flow simulations, thereby resulting in a highly computational demand, as already discussed above.In addition, the routing scheme generates streamflow only at the outlet of the catchment without the possibility to simulate streamflows at internal nodes, such as lakes, reservoirs or other infrastructures, while the use of response functions aggregated at the daily time scale limits the applicability to flood events occurring in large river basins, with the approach proposed by Gong et al. (2009Gong et al. ( , 2011)).To overcome the above limitations without resorting to hyper resolution hydrological models (when not needed to better reproduce spatial variability of soil water storage and transmission), we propose a multi-scale approach for streamflow routing based on the travel time approach.More specifically, we propose a scheme based on the Width-Function Instantaneous Unit Hydrograph (WFIUH) theory (see e.g.Rodríguez-Iturbe and Rinaldo, 1997), applied to a hybrid raster-vector data structure, according to the definition proposed by Lehner and Grill (2013).To emphasize that the model is designed to reproduce the effect of landscape on horizontal water transfer, we coined the name "HydroSCAPE", where "SCAPE" also recalls that the model is inherently SCAlable and highly ParallelizablE.In fact, by design, routing is independent from the Introduction

Conclusions References
Tables Figures

Back Close
Full grid size adopted to simulate water storage and runoff generation (we call this part the land component of the model), which as in all LHMs depends on the resolution of the climate or weather forecasting model used to simulate the meteorological forcing.In particular, geomorphological dispersion (Rinaldo et al., 1991) is invariant with respect to the grid size of the land component and the network response is perfectly upscaled to the computational grid.Scale-invariance is an important feature characterizing our modeling approach not enjoyed by grid-based LHMs, which rely on empirical upscaling procedures in order to represent unmodeled geomorphological dispersion.
An additional crucial aspect is the computational time efficiency, which stems from the fact that the most demanding step of the procedure is the computation of the geomorphological width functions, which is performed only once as an offline pre-processing procedure.This characteristic, coupled with easiness to parallelization and parsimony, makes HydroSCAPE an appealing tool for uncertainty assessment of the predictions, and in general for simulations conducted in a Monte Carlo framework.This paper is organized as follows: Sect. 2 describes the multi-scale hydrological conceptual model, details of the pre-processing module and derivations of node specific width functions are provided in Sect.3, with reference to the Upper Tiber case study.An example of application for two flood events is discussed in Sect.3.3, and finally concluding remarks are drawn in Sect. 4.

Model development
As stated in the Introduction, our aim is to develop a simple, parsimonious and computationally efficient method for modeling streamflow (and particularly floods) in large basins.To this aim, we adopt different modeling strategies for the river network and the associated hillslopes: a lumped non-linear formulation for the latter and a linear, geomorphologically-based approach for the former.The proposed modeling framework reflects the current understanding of hydrological processes at the hillslope and watershed scales (see e.g.Sivapalan, 2003).Introduction

Conclusions References
Tables Figures

Back Close
Full The sketch of Fig. 1 displays the conceptual model adopted here.The modeled domain A, which can be of any size and may include any number of disconnected river networks (for example, the sketch of Fig. 1 contains two river networks), is subdivided into N macrocells, each one characterized by a contributing area A (i ) , such that N i =1 A (i ) = A. We emphasize that the generic A (i ) does not need to coincide with the macrocell area, for instance this happens when the macrocell contains parts of neighbor watersheds (see the macrocells at the boundary of Fig. 1).There are no constraints on A, except that it must cover all the watersheds of interest.For easiness of representation, in Fig. 1 macrocells are represented as squares, but this is not by any means a constraint of the model and irregular macrocells can be used if convenient in the application of the model.However, size and geometry of the macrocells can be set to coincide with the gridding of the weather or climate model used to provide the meteorological forcing.
Drainage characteristics of the basins are obtained from the Digital Elevation Model (DEM) of the area of interest.The spatial resolution of the DEM should be fine enough to adequately capture the spatial structure of the drainage basins and the embedded river networks.Following procedures widely adopted for the identification of drainage direction and hillslope-channel separation (Tarboton et al., 1991;Montgomery and Foufoula-Georgiou, 1993) (see Sect. 3.2), the river networks are extracted and the hillslopes identified.
The area of the hillslopes changes through the domain, unless identification of the river networks is performed by using a constant threshold area.The link between the hillslope and the channel is denoted as hillslope-channel transition site.As an example, a synthetic DEM grid is shown in Fig. 1 (25 DEM cells per macrocell) together with the identified river networks.The figure shows also a few hillslopes (colored areas), each of them associated to the corresponding hillslope-channel transition site (colored squared symbols).Notice that in the example the brown hillslope is divided between two macrocells (i = 1, 2), a common situation in our approach since the discretization Introduction

Conclusions References
Tables Figures

Back Close
Full of the domain into macrocells may be arbitrary (i.e. it does not necessarily use topographic information).The next step in the construction of the model is the identification of the N k nodes where streamflow is desired.In the example of Fig. 1, N k = 5 nodes, identified by red bullets, are distributed within two networks.No limitations are imposed to the number and position of these nodes, which represent locations were streamflow is computed either to compare it with observational data, as part of the inversion procedure, or for other purposes, such as to verify flood protection structures (or the need to build them).Each macrocell i may feed one or more nodes; for instance, the macrocell i = 3 in Fig. 1 includes contributing areas to nodes k = 1 and k = 2, with different routes.It may also occur that a macrocell contributes to the same node k with different routes, as for the case of macrocell i = 2 which contributes to node k = 1, through both hillslopes highlighted in green and brown.
We denote with The streamflow generation processes occurring at the hillslope scale can be modeled by using schemes of different level of complexity, from the simple SCS-CN method (US Soil Conservation Service, 1964) to methods based on the solution of the Richards equation, or one of its simplification (Clark et al., 2015), depending on the objectives of the analysis.Whatever the hillslope model, hereafter we indicate with η ] the water discharge per unit area produced by the hillslope of area a , which belongs to the macrocell i and contributes to the streamflow at the closest downstream node k along the river network (for instance, the hillslope highlighted in green in Fig. 1  the network system through the hillslope-channel transition site, and is subsequently routed through it.From this kinematic scheme, it follows that the streamflow contribution of the hillslope , belonging to the macrocell i , to node k can be written as follows where a is the travel time from the hillslope-channel transition site of the hillslope to node k, and t [T] is the current time.
Under the hypothesis that the stream velocity V c [L T −1 ] is constant through the network, the travel time assumes the following expression: τ k is the distance, measured along the network, from the hillslope-channel transition site of the hillslope to the first downstream node k.The assumption of constant V c is crucial for the linearity of the processes at the watershed scale, as stated at the beginning of this section.Equation ( 1) is consistent with the general conceptual framework used to derive the Width Function Instantaneous Unitary Hydrograph by rescaling the geomorphological width function through a suitable constant velocity (Gupta et al., 1986;Mesa and Mifflin, 1986;Gupta and Mesa, 1988;Rodríguez-Iturbe and Rinaldo, 1997).
In agreement with the WFIUH theory, stream hydrodynamic dispersion is neglected, owing to its small to negligible effect on the hydrological response, which has been demonstrated to be dominated by geomorphological dispersion already embedded into the rescaled width function (Rinaldo et al., 1991(Rinaldo et al., , 1995;;Rodríguez-Iturbe and Rinaldo, 1997).
The streamflow q ] generated by macrocell i and contributing to the node k is then obtained by summing up all the contributions stemming from the hillslopes of the macrocell draining to the node k: Introduction

Conclusions References
Tables Figures

Back Close
Full (2) Under the hypothesis that the hillslope water discharge per unit area is constant through the cell i , i.e. by assuming that η (2) simplifies to: where k weighted by the relative hillslope area a (i ) k .In Eq. ( 3) the asterisk denotes convolution.
Finally, water discharge Q k (t) [L 3 T −1 ] at the node k is given by the sum of the direct contribution of each macrocell i to the node plus the contribution of the nodes upstream of k: where τ j k = D j k /V c is the travel time from the node j , located upstream to the node k, to the node k, D j k is the distance between the two nodes, and N up k is the number of nodes upstream of k.In the first right hand term of Eq. ( 4) summations are extended over all the macrocells, with the convention that q for all macrocells not contributing, i.e. not connected through the network, to the node k (the same applies for index j ).Notice that, in the second right hand term of Eq. ( 4) the streamflow Introduction

Conclusions References
Tables Figures

Back Close
Full computed at the node j is rigidly translated to the node k with the delay τ j k depending on the distance between the two nodes, thereby neglecting again the effect of stream hydrodynamic dispersion, as typically done in the WFIUH approach (see e.g.Gupta et al., 1986;Van Der Tak and Bras, 1990;Botter and Rinaldo, 2003;Giannoni et al., 2005).
The above method is simple and computationally effective.The underlying principle is similar to the block-effective dispersion employed in groundwater hydrology (see e.g.Rubin et al., 1999;De Barros and Rubin, 2011), with the travel time pdf used to represent the effect on streamflow of morphological heterogeneity at scales smaller than the macrocell, with a lower cutoff given by the DEM's scale.The linearity of the transfer processes at the scale of the network makes the algorithm easy to parallelize.In principle, streamflow generated by a macrocell can be elaborated by a single processor (with the convolution of Eq. ( 3) being the most demanding step), for the whole duration of the simulations, independently from the other processors.Then, the streamflow Q k (t) at the nodes of interest is further processed with Eq. ( 4) when all processors terminated the elaboration of their macrocell.For the same reasons, this model is particularly suited to multiple Monte Carlo runs (e.g. for parameter estimation, uncertainty analysis or multi-model or multi-scenario analysis).The main advantage of the method is that, by design, it preserves the global geomorphological dispersion of the basin, as calculated by the fine grid DEM, no matter the size of the macrocells.Consequently, the upscaling of river network dispersion is perfectly resolved, without resorting to hyper-resolution numerical grids.This point shall be illustrated in the ensuing Section.In addition, with the proposed approach we address one of the main limitation of the WFIUH formulation: the inadequacy of the original method (generally based on a single WFIUH for the whole river basin) to properly account for spatial variability of precipitation.It is in fact recognized that at intermediate spatial scale (typically beyond a few thousands of km 2 ) the spatial distribution of rainfall plays a fundamental role in shaping the hydrological response of river basins (see e.g.Nicótina et al., 2008;Volpi et al., 2012;Sapriza-Azuri et al., 2015).With our approach this limitation is overcome by assembling in the right Introduction

Conclusions References
Tables Figures

Back Close
Full hand term of Eq. ( 4) the convolution of Eq. ( 3) between the macrocell-specific rescaled width functions f k and the discharges per unit area η (i ) , with the latter embedding the spatial variability of rainfall patterns according to the macrocell resolution (a similar approach, but based on a partition of the catchment into sub-basins, can be found in Rinaldo et al., 2006;Rigon et al., 2015;Bellin et al., 2015).
Routing requires the definition of only a parameter, the channel velocity V c , which is a very parsimonious, yet effective, parametrization with respect to grid-based routing schemes.If the DEM resolution is high and the total domain A where the model is applied is large, the preprocessing step can be time consuming; the effort is however compensated in the application of the model, particularly if the modeling activity is performed in a multiple run framework.

Study area
In this section we describe an application of HydroSCAPE to the Upper Tiber river basin, providing a practical example of model characteristics and performances.The study area covers the upper portion of the Tiber river basin, located along the Apennine ridge (central Italy) between 42 • 36 and 43 • 51 latitude North and 11 • 48 and 12 • 55 12 longitude East (see Fig. 2a).The basin drains an area of approximately Intense precipitation events are typically associated with humid frontal advection from the Mediterranean sea and condensation due to the orographic uplift.Because of strong topographic gradients, headwaters experience intense rainfall events, mostly occurring from autumn to spring, associated with frequent flood events.Relevant flood events have been also observed in the floodplain of the river (southern part) where most of population and economical activities are clustered (Manfreda et al., 2014).Topography is represented through a 20 m resolution DEM provided by the Istituto Geografico Militare (IGM, available online at http://www.igmi.org/).Digital maps of land use and lithological characterization were supplied by the European Environmental Agency (Corine Land Cover project) and by Italian Agency for Environmental Protection and Research (ISPRA), respectively (maps not shown).Furthermore, land use classes from Corine classification and infiltration capacity estimates were used to associate at each DEM cell a value of the Curve Number parameter (CNII, see Fig. 2c), which shall be used in the SCS-CN runoff model as described in Sect.3.3.

Macrocell discretization, width functions derivation and perfect upscaling
The control sections adopted for multi-site model calibration (see Sect. 3.3) are located at 5 stream-gauge stations: Santa Lucia (SL), Ponte Felcino (PF), Ponte Bettona (PB), Ponte Rosciano (PR) and Ponte Nuovo (PN), with the latter being the outlet of the river basin (see Fig. 3).Drainage area (ranging between 1000 and 4000 km 2 ), longest flow path, and other geomorphic characteristics of the sub-catchments identified by the 5 control nodes are reported in Table 1.The control nodes are located along the main course of the Tiber river and its two major tributaries (Chiascio and Topino rivers), and they are equipped with gauges registering water levels at 30 min time step.Stage measurements in the period 1990-2000 together with validated stagedischarge relationships have been provided by the Hydrographic Service of Umbria Region (http://www.idrografico.regione.umbria.it).The meteorological forcing is described with half-hourly precipitation at 32 meteorological stations managed by the same institution.Figure 3 shows a map with the locations of the meteorological and stream 9067 Introduction

Conclusions References
Tables Figures

Back Close
Full gauging stations together with the subdivision of the watershed into the 5 inter-basin areas.
In the following the effect of spatial discretization on the hydrologic response is analyzed with reference to macrocells of different dimensions.In particular, the study area was overlaid with macrocells of increasing size, from 1 to 150 km (the latter including the whole Upper Tiber river basin within a single macrocell), and corresponding to about 37 and 1 • 36 , respectively.Domain discretization with macrocells of 5, 10 and 50 km are shown as an example in Fig. 4.
Identification of the drainage network and associated geomorphic metrics was performed by adopting standard DEM pre-processing techniques.In particular, the identification of the flow path lengths involved the following steps: (i) pit and flat area removal following the procedure of Tarboton et al. (1991); (ii) determination of the drainage directions by using the standard single direction D8 algorithm (O'Callaghan and Mark, 1984); (iii) identification of the space filling tree network, which connects each site to the outlet; (iv) definition of channel initiation, by adopting a combination of the thresholdslope area and the threshold-support area criteria (Di Lazzaro, 2009).
For a given resolution of the macrocell grid, it is thus possible to derive the frequency distribution f (i ) k of the flow path lengths L (i ) k pertaining to macrocell i and connecting the hillslope-channel transition site of the hillslope to the control node k (i.e. the macrocell-node specific width functions introduced in Sect.2), which is the best possible approximation of the width function, given the scale of the DEM.
Figure 4 shows as an example width functions constructed at the Santa Lucia node (upstream node on the main river-course, see Fig. 3), for macrocell size of 5 km (Fig. 3a), 10 km (Fig. 3b) and 50 km (Fig. 3c), respectively.We observe that the distribution of path lengths is wider for larger macrocells, because more hillslopes contributing to the node are included into the same macrocell, such that a larger portion of available path lengths is sampled.When a single macrocell is used, containing the entire catchment, all the hillslopes are included and the macrocell width function coincides with that of the catchment.On the other hand, when the macrocell coincides with the cell of Introduction

Conclusions References
Tables Figures

Back Close
Full the DEM the distribution of the lengths reduces to a Dirac delta.This is consistent with the underlying rationale of the model, which is intended to keep the geomorphologic component of dispersion as it is derived from the finest scale description of topography at hand (i.e. the DEM scale), even when runoff generation processes are represented at a larger scale (e.g. to comply with the output of climate models).When all the width functions of the macrocells contributing to the SL node are combined (i.e.57, 19 and 2 macrocells for 5, 10 and 50 km, respectively), the global width functions are exactly the same, thereby confirming that routing is insensitive to the size of the macrocell as desired.Hence, upscaling of geomorphological dispersion is by construction free of errors or scale effects, in our approach.This is not the case for models in which the geomorphological description of the drainage network is performed at the same scale of the macrocells grid (see the discussion in Sect.1).To show this, in Fig. 5 we compare the width function calculated at Ponte Nuovo, the outlet of the catchment, derived from the original 20 m DEM (grey line), which in our formulation is perfectly preserved for any choice of the macrocell size, with the analogous width functions derived after DEM aggregation at 5 km (blue), 10 km (orange), and 50 km (red), respectively.The figure shows, in particular, how aggregating digital topographic information over macrocells of progressively increasing size results in a deterioration of the width function, and, as a consequence, of the geomorphologic response of the watershed.This is not the case in our formulation in which the width function maintains all the information derived from the spatial resolution of the DEM.

Application example
In this section we present an example of application of HydroSCAPE for flood prediction in the Upper Tiber basin, with the purpose to illustrate its major computational and functional features.To focus on routing, the rainfall-runoff model was kept as simple as possible.In particular, we combined the widely used SCS-CN method (US Soil Conservation Service, 1964) for runoff simulation with a linear reservoir model describing the travel time distribution within the hillslope (Rodríguez-Iturbe and Rinaldo, 1997, Introduction

Conclusions References
Tables Figures

Back Close
Full 3).This is consistent with the notion that travel times in hillslopes are important in shaping the hydrologic response of a watershed and cannot be neglected even for large river basins, where the channeled lengths are usually much larger than the mean hillslope size (Botter and Rinaldo, 2003;D'Odorico and Rigon, 2003;Di Lazzaro and Volpi, 2011).Subsurface contribution to streamflow was not explicitly considered for this specific model configuration, which is focused on floods.As an alternative to the linear reservoir model, a hillslope scale rescaled width function can be used to represent the travel time distribution at the hillslope scale.In this case, rescaling may be obtained by using a hillslope specific velocity V V c .At the hillslope scale runoff is computed by using the classic SCS-CN scheme: where P i (t) [L] and R i (t) [L] are the cumulative rainfall and the cumulative runoff, respectively, at time t, both assumed uniform within the macrocell i .In addition, S i [L] is the soil potential maximum infiltration (identified on the basis of the map of CNII shown in Fig. 2c), I a = α c S S [L] is the initial abstraction, with α < 1 [-] introduced to represent the initial abstraction as a fraction of the maximum infiltration, and c S [-] is a multiplicative factor accounting for uncertainty in the identification of S. Therefore, the effective rainfall intensity p i [L T −1 ] at time t can be computed as follows: and Eq. ( 5) is applied at discrete times, according to the time step ∆t.
The specific water flux produced by the hillslopes of the macrocell i can be obtained by applying mass conservation at the hillslope scale by considering the effective precipitation as inflow and runoff η as the only outflow (evapotranspiration can be neglected Introduction

Conclusions References
Tables Figures

Back Close
Full since the model is applied at the flood event temporal scale): where λ [T] is the mean residence time of the linear reservoir and the left hand term is the first order approximation of the time derivative of runoff η (i ) .Parameters α, c S and λ were assumed uniform through the river basin, i.e. all the macrocells share the same coefficients.On the basis of preliminary analysis α was found not to be a sensitive parameter and was set to 0.08 (which is in agreement with the values found by D'Asaro and Grillone, 2012), while c S and λ are calibration parameters, together with the channel velocity V c .We emphasize that this simplified version of HydroSCAPE is event-based since it does not include a continuous soil-moisture budget; however, this is enough for the purpose of this example application, mainly focused on flood events, which aim is to show how HydroSCAPE implements routing.As explained in Sect.2, HydroSCAPE is not limited to this simplified implementation, yet effective for the purpose of flood forecasting, and can work with more sophisticated runoff generation schemes.
In order to illustrate model performance we selected two major rainfall events within the decade 1990-2000, that generated significant, yet not extreme, floods.The streamflow triggered by these rainfall events was compared with observational data at the 5 nodes described in Sect.3.1.The two events occurred in December 1992 and February 1999, respectively.In both cases precipitation was caused by humid frontal advection from the Mediterranean Sea followed by condensation due to orographic uplift (Calenda et al., 2005).The spatial distribution of the precipitation was obtained by interpolation of the measurements at the available rain gauges (18 and 32 for the events of December 1992 and February 1999, respectively) by means of Kriging with External Drift (see e.g.Goovaerts, 1997).The precipitation was interpolated separately at each time step by using the same exponential semivariogram which has been obtained by analyzing offline the available data.In particular, precipitation was first calculated over 9071 Introduction

Conclusions References
Tables Figures

Back Close
Full a 1 km resolution grid and successively aggregated at the macrocell scale, according to the resolution adopted in the simulations.
In order to test the computational efficiency of HydroSCAPE, model calibration was performed generating a large number (i.e. 100 000) of model parameter sets using the Latin Hypercube Sampling technique (McKay et al., 1979) with the following boundaries: V c ∈ [0.5, 4] m s −1 , c S ∈ [0.3, 3], and λ ∈ [0.01, 1] d.The optimization procedure was based on the maximization of the Nash-Sutcliffe Efficiency (NSE) index (Nash and Sutcliffe, 1970) for streamflow evaluated at the outlet of the basin (i.e.Ponte Nuovo, see Fig. 3).The model was run with four spatial resolutions (i.e. 5, 10, 50 and 150 km, see Sect.3.2) with a computational time step of ∆t = 0.5 h, and calibrated on the event occurred in the period 6-12 February 1999 at Ponte Nuovo (PN) station.Results in terms of optimized parameter sets, NSE index, and computational time for the entire set of 100 000 runs are summarized in Table 2.
In all cases, the NSE index at the calibration section (PN) assumes high values, close to one, indicating a very good model fit to the observed streamflow data.Optimal parameter sets assume similar values at all the scales, suggesting that the model is able to preserve geomorphological dispersion when the domain is discretized with macrocells of increasing dimension.This is verified also when a single macrocell of 150 km resolution is used, though in this case the impossibility to reproduce the spatial variability of the rainfall (given that only a single macrocell is used the precipitation is considered uniform over the entire basin) resulted in an inaccurate description of inter-basin propagation of fluxes, as emphasized by the negative values of the NSE index averaged over all nodes.Conversely, for all the other spatial resolutions, overall NSE values between 0.56 and 0.69 was obtained.Notice that all cases with the average NSE > 0.5 are with a macrocell dimension equal or smaller than the integral scale of the precipitation, which is about 36 km (E.Volpi, Modello di struttura spaziale del campo di precipitazione, unpublished technical report).It is therefore clear that the inaccuracies encountered with large macrocells are due to the inaccurate spatial description of the precipitation.Finally, the computational cost for 100 000 runs and for a Introduction

Conclusions References
Tables Figures

Back Close
Full single processor (Intel(R) Xeon(R) W5580 at 3.20 GHz core), the code being written in Fortran 90, is shown to increase from a few seconds in the case of one single macrocell to about 10 min for the finer resolution (1 km, corresponding to 476 macrocells).We emphasize that the computational efforts can be reduced considerably by implementing parallel computing techniques, to which HydroSCAPE is particularly suited thanks to its inherently parallel formulation (see also Sect.2).
Model validation was carried out by means of a combined multi-site, multi-event approach, and was coupled with a Monte Carlo based uncertainty analysis performed on a subset of parameter combinations sampled during calibration at PN with the 1999 flood event as observational data.This subset was identified according to a model efficiency rejection criterion that classifies as behavioral all sets of parameters with a NSE index greater than zero, resulting in 14 958 parameters sets and model realizations.Successively, the 95 % uncertainty bands associated with the retained simulations were evaluated using the standard likelihood weighted procedure proposed by Freer et al. (1996).Results obtained for the 10 km spatial resolution configuration are presented in Fig. 6, which shows 95 % prediction uncertainty bands and observed streamflow at the 5 nodes shown in Fig. 3. Simulated hydrographs obtained adopting the optimal parameter set reported in Table 2 are also shown with a continuous black line.Figure 6a-e show uncertainty bands for the February 1999 event (the calibration event) considering all the gauging stations including PN, which was the only used in calibration.Other indicators of goodness are P and R factors (see e.g.Abbaspour et al., 2009), which are defined as the percentage of data bracketed by the confidence band, and the ratio between the average width of the band and the standard deviation of observations, respectively.Computed water discharge at all nodes provided high P factor values (80 % for SL and 100 % for all the others), and moderate R factor values (1.99, 2.56, 1.62, 2.87, and 3.21 for PN, PF, SL, PR, and PB, respectively).The somewhat suboptimal performance with respect to the R factor is in part due to the decision of considering behavioral all the models with NSE > 0, instead of the typical choice of setting the threshold at NSE = 0.5, with the consequent reduction of the uncertainty band Introduction

Conclusions References
Tables Figures

Back Close
Full thus of the R factor.Visual inspection of Fig. 6 and the above performance factors indicate that the model is able to encompass most of the observed discharges, while retaining reasonable uncertainty band amplitudes.The same analysis was performed also for the event occurred between 4 and 7 December 1992 (multi-site, multi-event validation).Results are presented in Fig. 6f-i, which suggest reasonable model prediction capability (P factor equal to 100, 39, 17, 100 %, and R factor equal to 1.63, 1.47, 1.06, and 2.28, for PN, PF, SL, and PR, respectively; we note that no water discharge data were available at PB during this event), although a general tendency to underestimate the flood volume is evident.This is likely due to inherent differences between precipitation conditions (e.g.intensity, spatial distribution) during the two events and in the preceding days, which reflect into different initial soil moisture conditions that cannot be fully captured with the simple event-based SCS-CN model used here.

Conclusions
This work presents an innovative, multi-scale streamflow routing method based on the travel time approach.The principal aim is to develop a simple, parsimonious and computationally efficient method for modeling streamflow (and particularly floods) in large basins.The model, coined as HydroSCAPE, aims at correctly reproduce the relevant hydrological fluxes across the scales of interest, from a single catchment to the whole continent.The method is based on the WFIUH theory applied to a hybrid rastervector data structure, that allows to derive localized information on travel times and flow characteristics without the need of narrowing the resolution of the computational grid adopted for the study area.The relevant features of the model are illustrated through the modeling of two flood events in the Upper Tiber river basin (Central Italy), with 4 different domain discretizations, i.e. different dimensions of macrocells.
The main results of the present work can be summarized as follows: -HydroSCAPE employs a strategy for modeling cell-scale runoff dispersivity such that the catchment response is independent from the grid size, which in turn is 9074 Introduction

Conclusions References
Tables Figures

Back Close
Full function of the resolution of the atmospheric model or the integral scale of observed precipitation, in case ground-based rainfall measurements are used.In particular, geomorphological dispersion is invariant with respect to the grid size, and the upscaling of the river network response at the cell scale is automatically taken into account.This "perfect upscaling" characteristic of HydroSCAPE is particularly important in all cases when the catchment response needs to be accurately represented, e.g. when dealing with extreme events like floods and inundations.
-The above "perfect upscaling" characteristic allows adopting large cells, making the model suitable to large-scale models, up to the continental scale.The overall response function of the river networks will anyway be preserved, no matter the discretization.
-Computational efficiency is another relevant feature of the proposed approach.Efficiency stems from the fact that the demanding calculation of the width functions is a pre-processing, one-time effort.Furthermore, the model is prone to parallelization, stemming from the linearity of routing and independency of the streamflow generation modules at the cell scale.These features make HydroSCAPE an appealing tool for uncertainty assessment of the predictions, and for simulations conducted in a Monte Carlo framework.
-The routing component of the model (including hillslope routing) depends on two parameters, while the model inherits the parameters introduced in the conceptual model of runoff generation at the hillslope scale.While no limitations are posed to the latter conceptualization we are in favor of a pragmatic "downward" approach, which limits the total number of parameters, to reduce uncertainty and overparametrization.Parsimony is important for a meaningful and reliable parameter estimation procedure and uncertainty analysis, especially when dealing with large-scale and complex basins.We believe that all of the above characteristics make of HydroSCAPE an appealing routing tool to be implemented in LHMs, particularly suitable for climate change impact studies where the accuracy of the streamflow routing may be significantly affected by the spatial resolution adopted.Introduction

Conclusions References
Tables Figures

Back Close
Full  Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | contributes to node k = 1, which is the first node encountered moving downstream).The resulting water flow depends on the partitioning of hydrological fluxes, triggered by rainfall or snowmelt, at the hillslope scale and may include the contribution of groundwater.According to the above conceptual scheme, water flow produced by the hillslope enters Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | q Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | about 25 % of the entire Tiber basin at the river mouth in the Tyrrhenian sea.The basin is predominantly mountainous, with elevation ranging from 145 to 1560 m a.s.l.(see Fig. 2b) and it is aligned to the North-West South-East direction, with the Apennine ridge-line representing an important physical boundary at East that causes variability in precipitation and temperature.From a geological point of view most of the catchment is underlined by low-permeability formations, chiefly flysch, sandstone clay, and limestone clay.However, high-permeability formations (calcareous lithology) are found in the upper part of the basin and on the eastern divide.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ch. 7.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1 .
Figure 1.Sketch of basin conceptualization: subdivision of the study area into macrocells and nodes (red dots).River network is subdivided into hillslope-channel transition sites (colored squared symbols) each associated to a pertaining hillslope area a areas).For this simple case we consider N = 9 macrocells and k = 5 nodes.An example of two pathways characterized by the same length L (i ) k and travel time τ

Figure 3 .Figure 4 .
Figure 3. Map showing the subdivision of the watershed into 5 inter-basins, each one identified by a node where water discharge is computed (black triangles).The location of the meteorological stations are also shown as colored dots.

Figure 6 .
Figure 6.Comparison between 95 % confidence band, water discharge simulated with the optimal parameter set and observed at all the gauging stations and events: (a)-(e) 6-12 February 1999, and (e)-(i) 4-7 December 1992.Shaded areas identify the period considered in the evaluation of P and R factors.Model parameters are estimated by using the discharge at Ponte Nuovo (PN) measured during the 1999 event (a).

Table 1 .
Main geomorphic characteristics of the inter-basin drainage areas within the Upper Tiber river basin (CV: coefficient of variation).