Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1, pp. 46–57 (2012) SEMI-EXPLICIT MODELLING OF WATERSHEDS WITH URBAN DRAINAGE SYSTEMS Benjamin J. Dewals *, Pierre Archambeau, Bruno Khuat Duy, Sébastien Erpicum and Michel Pirotton Department ArGEnCo, University of Liege, Chemin des Chevreuils 1, B52/3+1, B-4000 Liege, Belgium * E-Mail: B.Dewals@ulg.ac.be (Corresponding Author)  ABSTRACT: In rainfall-runoff modelling of urbanized and semi-urbanized watersheds, the urban drainage systems considerably influence runoff propagation time. In small scale watersheds, the drainage network may be modelled explicitly. In contrast, for larger watersheds, most hydrological models are based on a rough representation of the effects of drainage systems, thus failing to represent the rapidly-varying real flow dynamics. Therefore, a trade-off methodology has been developed to account for impervious surfaces and drainage effects accurately, without the need for modelling the entire drainage network in detail. Undrained impervious areas have been distinguished from drained ones. Rain falling on the former has been discharged as overland flow, whereas flow on the later has been routed separately using “virtual pipes”, which enable a simplified process-oriented modelling of the drainage network. The methodology has been applied to a 130 km² Belgian catchment, resulting in the simulation of fast flow peaks, which do not appear when the effect of the drainage network is neglected. Keywords: hydrology and water resource, sewers & drains, floods and floodworks 1. INTRODUCTION To account for impervious areas, different methodologies have been developed, depending Flooding in urbanized areas may result from a mainly on the type of hydrological model used variety of causes, including intense storms, and on the extent of the watershed being studied: prolonged rainfall over the urban area, incapacity of drainage systems, tidal surges, rise in  Statistical hydrological models are based on groundwater level, as well as failure or discharge measurements, and as such they overtopping of flood defences along rivers. incorporate the effect of impervious surfaces In the latest case, a thorough analysis often (Fleming, 1975), but they do not take into requires a complete modelling chain (e.g., Khuat account future landuse changes. In particular, a Duy et al., 2010; McMillan and Brasington, 2008), number of recent developments in hydrological combining rainfall-runoff modelling, simulation modelling have been based on soft computing of flood wave propagation in the river network, as techniques (e.g., Chau et al., 2005; Cheng et al., well as detailed analysis of inundation flows 2008; Cheng et al., 2002; Lin et al., 2006; using depth-averaged 2D simulations, possibly Wang et al., 2009; Wu et al., 2009) involving sediment transport or pressurized flows.  Other approaches, such as the rational method In such complex studies, rainfall-runoff modelling (Chow et al., 1988), enable a simple in urbanized areas often plays an important part. representation of urban areas using runoff The present contribution focuses on a semi- coefficients representative of the landuse. explicit approach for modelling watersheds with These coefficients represent mean values based urban drainage systems. on an estimation of the urban density in the From a hydrological perspective, urbanization has considered area (Chow et al., 1988; Thorndahl two main consequences: an increase in et al., 2006). impervious surfaces and the set up of drainage and sewer systems in the catchment. Both need to  One of the most widely used methods be considered in hydrological modelling. worldwide, namely the SCS Curve Number Impervious areas (Beighley et al., 2009; Endreny method, uses runoff coefficients depending on and Thomas, 2009) can be classified as drained or soil type, landuse and rainfall volumes (Garen undrained, depending on whether they are and Moore, 2005; Hjelmfelt et al., 2004). For connected or not to a sewage network (Alley and urban areas, specific coefficients are applied, Veenhuis, 1983; Brabec, 2009). which result from averaging the curve number values for the impervious surfaces and their Received: 12 Feb. 2011; Revised: 5 Aug. 2011; Accepted: 15 Aug. 2011 46 Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1 (2012) surrounding areas (often considered as gardens) falling on undrained areas has been discharged as according to the local urban density overland flow. Rain falling on drained areas has (NRCS/ARS Curve Number Work Group and been routed separately using a simplified Moody, 2004).  modelling of the drainage network, called here  In the model Topurban (Valeo and Moin, 2000), virtual pipes. This separate modelling of the the runoff is routed separately depending on drained impervious areas constitutes our main whether it is produced by pervious or original contribution. impervious areas: rain falling on non saturated pervious areas is considered as throughflow, 2. RAINFALL-RUNOFF AND while runoff from impervious areas is HYDRAULIC MODELS considered as surface flow (with an optional The rainfall-runoff model used here is process-routine accounting for storage ponds). oriented and spatially distributed (Khuat Duy et However, none of these approaches explicitly al., 2010). Using a multi-layer approach based on considers drainage systems, which may layer-averaged equations (Figure 1), it computes significantly modify the flow paths and shorten the main hydrological processes, including propagation times. Modelling the sewage systems overland flow, interflow, infiltration and deep is thus necessary if the rapidly-varying flows percolation. from drained impervious areas need to be represented. Detailed modelling of flows in drainage networks, such as sewers, is generally implemented only as far as urban hydrology is concerned. This refers to small-scale studies consisting mostly of impervious areas (Aronica and Cannarozzo, 2000; Rodriguez et al., 2008; Schmitt et al., 2004). In (a)contrast, when broad-scale analyses of large watersheds are conducted, the impact of sewage Rainfall Evapotranspirationsystems is usually neglected, since most of the Overland flow: computed from Eq. (1) and (2)water travel time is spent in the river network. A Inflow River flow: watershed is considered here as “large” with Infiltration: computed from Green‐Ampt formula to the  computed fromrespect to the impervious areas of interest, if the river Eq. (6) and (7)propagation time in the river network Interflow: computed from Eq. (3) to (5)significantly exceeds the propagation time as Deeppercolationoverland flow in these impervious areas and to the river. (b) Aquifer  In mid-size watersheds, the effect of the drainage Fig. 1 (a) Flow layers and main processes computed network on discharge may not be negligible, in the hydrological model; (b) Flow chart of while it may also not be feasible to take the the main computation steps in the rainfall-drainage network into account with the same level runoff model. of detail as in urban hydrology (Lhomme et al., 2004). In addition, the extent of drained areas The overland flow is modelled with the diffusive remains hard to estimate (Walsh et al., 2009) and wave approximation. This model is derived from the opportunity of conducting direct the fully dynamic shallow-water equations (SWE) measurements is restricted to limited cases mainly by assuming that inertia terms may be neglected in urban hydrology (Alley and Veenhuis, 1983). compared to gravity, friction and pressure terms. In this paper, we present an original methodology This assumption has been widely recognized as to compute impervious surfaces accurately and to generally accurate for modelling overland flow take drainage effects into account without (Jain and Singh, 2005). The hyperbolic fully modelling the entire drainage network in detail. dynamic shallow-water model may then be An existing process-oriented and spatially replaced by the following parabolic equation: distributed hydrological model has been adapted to accommodate the new methodology. The h uh  vh   S              with impervious surfaces have been estimated on the t x ybasis of landuse maps in vector format. Next, they hhave been classified as drained or undrained. Rain S fi  sini  cosi i  x, y , (1) i47 Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1 (2012) where h [m] is the water depth, u and v [m/s] the the infiltration coefficients, as suggested by velocity components along the x and y axis, Nearing et al. (1996). S [m/s] the source terms (rainfall and infiltration), The rainfall runoff model computes lateral Sfi (i=x,y) [-] the friction slopes, θx [-] and θy [-] inflows to the rivers, which are next routed the projected bottom slopes. Using Manning- through the river network by means of the Strickler formula, the velocity components can be hydraulic model Wolf 1D. A suitable shock related to the friction slopes as follows: capturing scheme is used to solve the S S conservative form of the 1D Saint-Venant u 1 h2/3 fx  ; v 1 h2/3 fy equations, enabling the simulation of flow regime n S 2  S 2 1/4 n 2 1/4fx fy S fx  S 2fy  changes and hydraulic jumps (Kerger et al., 2011).  (2) The following set of equations is solved: where n [s /m1/3] is a spatially distributed A  q    Manning coefficient.   t q   x  q2  The infiltration is calculated using Green-Ampt     g cos p A A formula (Chow et al., 1988), while the interflow  q is modelled with the depth-averaged Darcy L  0equations. It is therefore computed using a   hgAsin  gAJ   g coshl b  g cos p  b x 0diffusive wave equation similar to the overland  x flow equation:  (6) h u h   v h  where A [m²] is the river cross section,  q [m³/s] p sub  sub sub  sub sub  Ssub  (3) the discharge, h [m] the water height, qL [m³/s] the t x y lateral inflow, J [-] the energy slope (accounting S for bottom roughness and internal losses), θ [-] u  K fsub,xsub s 1/4 ; the mean bottom slope and lb [m] the channel S 2 2fsub,x  S fsub, y  bottom width. The pressure terms pA [m³] and px [m²] are defined as: S  (4) vsub  Kfsub,x hs S 2 2 1/4 pA (h)  h  l(x, )d ; fsub,x  S fsub, y  0h h  l(x, ) (7) S sfsub,i  sinsub,i  cossub,i i  x, y  (5) px (h)   h   di 0  xwhere p [-] is the soil effective porosity, h  [m] where l [m] is the channel width and a bound subthe subsurface water depth, usub and vsub [m/s] the variable. The energy slope is evaluated using interflow velocity components, Ssub [m/s] the Manning formula as follows: source terms (infiltration, deep percolation), Ks q2[m/s] the lateral hydraulic conductivity (spatially J  2 2 4/3  (8) distributed), Sfsub,i (i=x,y) [-] the subsurface K A RHfriction slopes, and θsub,i (i=x,y) [-] the projected with RH [m] the hydraulic radius and K [m1/3/s] slopes of the soil layer. the Strickler coefficient of the river. Pre-processing of raw data has been carried out The spatial discretization of equations (1), (3) and using the GIS interface of Wolf. The Digital (6) is based on the finite volume method Elevation Model (DEM) has been processed to combined with a self-developed flux-vector remove depressions using an algorithm developed splitting technique. Upwind evaluation of the by Martz and Garbrecht (1999). The DEM-based fluxes is performed, according to the sign of the flow paths have been made consistent with the flow velocity (Erpicum et al., 2010b). The real ones obtained from on-site surveys (Callow, stability of the scheme has been demonstrated 2007; Saunders, 1999). The roughness coefficient using a linear stability analysis (e.g., Kerger et al., n has been deduced from landuse maps, while the 2011). Efficiency, simplicity and low soil type obtained from pedologic maps has been computational cost are the main advantages of used to evaluate the soil conductivity based on this scheme. The resulting ordinary differential pedotransfer functions (Rawls and Brakensiek, equations are integrated in time using an explicit 1989). The influence of landuse on infiltration has Runge-Kutta scheme or an implicit algorithm been taken into account using effective values for based on the Generalized Minimal RESidual 48 Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1 (2012) method (GMRES). An original treatment of the junctions by means of Lagrange multipliers enables modelling of large river networks within a single simulation. The hydrologic inflows are treated as lateral inflows (source term qL) for the 1D simulations and have to be computed separately beforehand. The hydrologic and the river flow equations are therefore decoupled. The herein described models are parts of the modelling system Wolf, developed at the University of Liege. Wolf combines a set of complementary and interconnected modules for simulating free surface flows: process-oriented hydrology, 1D, 2D depth-averaged, 2D-vertical and 3D flow models. Their validity and efficiency (a)have already been verified in numerous  applications (Dewals et al., 2006; Dewals et al., 2008; Dufresne et al., 2011; Erpicum et al., 2009; Machiels et al., 2011), including inundation modelling (Dewals et al., 2011; Ernst et al., 2010; Erpicum et al., 2010a; Roger et al., 2009). 3. METHODOLOGY 3.1 Estimation of impervious surfaces Existing methods to estimate the impervious surfaces can be classified as indirect or direct (Moglen, 2009). The former approach is the most standard one and consists in estimating the surface imperviousness based on the urban density, deduced from the (b)  corresponding landuse class, and using lookup Fig. 2 Computation of impervious areas: (a) tables (Beighley et al., 2009). This technique may Impervious landuse; (b) Impervious part of however prove inaccurate, since the real amount cells. of impervious areas for a given landuse class shows a high variability. Moreover, Thorndahl et In our spatially distributed model, the impervious al. (2006) highlighted that hydrological reduction fraction of the surface has been computed in each factors calculated from measured rainfall and cell based on landuse maps in vector format runoff have been found inconsistent with (Figure 2a). The impervious landuse types (such corresponding values for residential areas as roads, houses, car parks or concrete structures) reported in literature. The indirect method is have been selected and their corresponding nonetheless often used, due to the wide data surfaces have been summed for each cell, using a coverage available from remote sensing landuse specific algorithm dedicated to the computation of recognition techniques. overlay surface between a polygon and a regular In contrast, direct methods identify the mesh. As a result, a raster containing the impervious areas by means of field surveys and impervious fraction of each cell has been obtained aerial imagery, leading to a more accurate (Figure 2b). estimation of impervious surfaces (Beighley et al., 2009; Chabaeva et al., 2009; Jones and Jarnagin, 3.2 Drainage network 2009). The main drawbacks of direct methods are After the computation of the impervious surface the limited availability and the cost of data. in each cell, drained and undrained cells have Endreny and Thomas (2009) introduced a been distinguished. If data on the pipe network combined approach, in which landuse maps were are available, our methodology assumes that all enhanced based on data on impervious areas cells located within a given distance from the pipe extracted from road networks. network are drained. A cell is thus considered as 49 Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1 (2012) drained if the distance from its centre to the nearest pipe is below a threshold value. The choice of this value is further analysed in the case study (section 4.2). In contrast, if the drainage network characteristics are unknown, the percentage of drained areas can be estimated following the approach suggested by  Boyd (1993) for urban catchments. The (a) (b) estimation is based on selected low intensity rainfall events, during which the flow in the rivers Fig. 3 Transformation of the sewage network. Small is essentially generated from the drained circles represent examples of drained cells impervious areas. The corresponding runoff outlets. (a) Original tree-shaped network; (b) coefficients are evaluated and their lower bound Equivalent pipe. provides the part of drained areas in the 2 2catchment. This methodology applies provided a i n u 4/3  (9) sufficiently high number of rainfall events are Rhconsidered. Otherwise, the results may be affected by errors arising from measurement inaccuracies where i [-] is the pipe slope, u [m/s] the flow and from the selection of rainfall events during velocity in the virtual pipe and Rh [m] the which not only drained impervious areas hydraulic radius. The circular pipes are assumed contribute to runoff. This method can also be used to be designed for a constant rainfall of 50 mm/h in combination with the previous one, either to over the drained area, with an occupancy rate of validate the results or to calibrate the threshold 80%. These are proposed default values, which distance from the drained areas to the pipe may be adjusted on a case by case basis. The network. Manning coefficient n of the virtual pipes has -1/3Real sewage systems are tree-shaped networks, been evaluated at 0.015 m s, which is a standard with a high number of small pipes. The complete value for most concrete channels and pipes. This modelling of such systems for a whole watershed results in the following relation: requires large amounts of data, usually partly S 3/8 imp unavailable. Moreover, the density of the network Deq  0.005 1/2   (10) implies using very fine cells (down to a few  i meters) to represent the numerous individual where D  [m] is the virtual pipe diameter, and pipes, reducing consequently the time step and eqS  [m²] the impervious surface drained by the dramatically increasing the computational cost. imppipe. Our methodology consists in merging parts of this tree-shaped network into equivalent pipes, called 3.3 River network hereafter virtual pipes, in order to reduce the model complexity while still representing the key In a process-oriented and spatially distributed processes leading to fast flow propagation to the hydrological model, rivers serve as outlets for the river. As sketched in Figure 3, the runoff from runoff and baseflow, which are next propagated each drained cell is discharged into a through the network of rivers. To this end, our corresponding virtual pipe, at a location enabling methodology enables to determine the river to preserve the same length of flow path as in the network characteristics by combining three real network. complementary data sources. The slope of the virtual pipe is computed by a  First, the DEM enables to generate a river weighted average of the slopes of the network based on flow paths and a convergence corresponding real pipes. Its cross section varies treshold (Renssen and Knoop, 2000; Tarboton as a function of the drained impervious surface. et al., 1991). This network covers the whole Since this cross-section depends on the detailed watershed, but contains limited information on structure of the real network and on the the cross-sectional geometry of the rivers. dimensions of the real pipes, usually unknown, we have developed a relation to estimate default  Second, data from field surveys are exploited to values for the equivalent pipe diameter. Uniform incorporate a more detailed description of the flow conditions are assumed in the pipe, as river geometry and cross sections. However, described by Manning-Strickler formula: such detailed data are generally available only on a limited part of the basin. 50 Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1 (2012)  A third source of data is obtained from the downstreamsewage network, generated using the method described in section 3.2. The overall procedure for generating the river network involves the following six steps:  The DEM is first modified using a dedicated algorithm to force the topography-based river Pipes (sewage system) paths to follow the real river courses where Network from DEM such data are available. Engineering structures Network from on‐site surveys affecting the flowpaths (such as road or railway Splitting points   dikes) are taken into account at this stage of the procedure. Fig. 4 Splitting of networks.  Next, a first river network is generated based downstreamon this modified DEM.  In parallel, a second network is created using the field survey data, including every relevant hydraulic structure such as culverts, footbridges and pipes.  Finally, these two networks and the additional Removed splitting pipes from the sewage system are merged. This points is performed in three sub-steps: Confluences  1. The river branches are first split into Fig. 5 Merging of network appropriate parts. multiple reaches between characteristic points, such as junctions (Figure 4) Rainfall stations 2. Next, the DEM-based river reaches are Discharge gauging sites replaced by field survey-based reaches where these are available. Possible inconsistencies between both networks, such as bed level discontinuities at the junctions, are handled at this stage of the process. 3. Eventually, all river reaches are merged again to form the complete river network (Figure 5). This method automatically generates a 1D network, which (i) covers the whole watershed, (ii)  includes detailed data from field surveys (including cross sections and hydraulic structures, Fig. 6 DEM of watershed (elevation in metres). such as weirs and culverts) and (iii) incorporates 4.1 Watershed description the virtual pipes representing the drainage network. The watershed is mainly covered by meadows (70%) and includes several low or medium-4. CASE STUDY density towns. The DEM of the watershed is characterized by a mean slope of 7.2% (Figure 6). The methodology has been applied to the Hourly discharge measurements are provided by Berwinne watershed located in Belgium. River two gauging stations located on the river, while Berwinne is a tributary of river Meuse and has a rainfall is measured daily at eight stations and catchment area of 132 km². hourly at one station. These stations are located either inside or nearby the catchment. A disagregation process has been applied to the daily series to provide hourly rainfall intensities (Koutsoyiannis, 1994; Koutsoyiannis, 2003). 51 Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1 (2012)  4.2 Application of the methodology According to the procedure described above, the impervious surfaces have been computed based on landuse maps (section 3.1). These landuse maps were available in vector format on 92% of the watershed (source: Belgian National Geographic Institute). Application of this method has resulted in a mean value of 6.3% of impervious surfaces throughout the watershed. Modelling of the drained impervious areas has been performed based on the real sewer system, provided in the “Plan d'Assainissement par Sous-bassin Hydrographique (PASH)” (SPGE 2007). (a)These maps include the coordinates of the existing pipes in vector format (Figure 7a). Using the method described in section 3.2, the sewage network trees have been merged into virtual pipes. Virtual pipes with nearby outlets have been combined in order to further reduce the total number of pipes directly connected to the river. To obtain a reasonable number of virtual pipes in the final network, a minimal distance of 1000 m between two adjacent outlets was imposed. With this value, the propagation time lag between two adjacent outlets remains far below the time step at which rainfall data is available. Applying this methodology has resulted in a set of 30 virtual (b)pipes representing the entire drainage network. Fig. 7 Modelling of sewage system: (a) Existing The virtual pipes can be represented in the sewage network; (b) Locations of outlets of network by arbitrary lines of adequate length and virtual pipes. starting from their outlet point situated along the river (Figure 7b). The cells are considered as drained by the sewage system if their centre lies at a distance below 150 m from the pipes. This distance represents the runoff path on impervious surfaces, in the gutters and in small pipes which do not appear in the drainage network. A sensitivity analysis has been carried out to assess the influence of this parameter on the drained surface. To this end, the distance from each cell to the nearest pipe was first computed (Figure 8). Using this information in the part of the map covered by vector data, the total drained impervious surface (sum of impervious areas from drained cells) has next been plotted as a  function of the cell-pipe distance threshold Fig. 8 Distance (in metres) from each cell to nearest (Figure 9). pipe (zoom on a town located in watershed). For small values of threshold distance from cell to pipe (below 50 m), the surface of drained near the sewage network. Above this distance, the impervious areas is found to increase significantly rate of increase of the drained area with the with the threshold value. This results from the threshold distance becomes significantly lower tree-shape of the network. As can be seen in (Figure 9). This residual increase rate corresponds Figure 8, the chosen distance (150 m) enables to to the impervious areas located outside urbanized consider as drained most urbanized areas located  areas, such as roads or isolated houses and structures. 52 Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1 (2012)  20 018 216 414 612 810 108 126 144 162 180 200 20 40 60 80 100 120 140Time [h]Measured discharge Discharge with the drainage network Discharge without drainage Rainfall Fig. 9 Drained impervious area and rate of change of drained impervious area (with moving Fig. 11 August 1996 flood: Rainfall intensity and average) as a function of the cell-pipe comparison between measured and computed distance. discharge at downstream gauging station. 4.3 Validation To assess the benefits of modelling the drainage network, simulations have been conducted for a real rainfall event (August 2006). Results obtained both with and without modelling of the drainage network have been compared with measured discharges at the most downstream gauging station. In both cases, the considered total impervious area was the same.  Hourly rainfall values, disaggregated from daily Fig. 10 Analysis of rainfall-runoff values for low values at the weather stations, were distributed intensity events. over the catchment using the Thiessen polygons method. The program MudRain developed by In the Berwinne watershed, the computed drained Koutsoyiannis (2003) was used to perform impervious area was found to cover 2.5% of the multivariate disaggregation of rainfall data. Since catchment surface. A complementary analysis of the considered rainfall event occurred during a discharge measurements from low intensity low-flow period, the baseflow has been shown to rainfall events has been carried out to estimate the be negligible compared to the subsequent flood drained impervious areas, as detailed in section discharge. 3.1. Figure 10 shows a rainfall-runoff plot, in For the simulation without consideration of the which the baseflow component was subtracted drainage network, the spatially-distributed friction from the measured discharge. All rainfall events coefficient and soil conductivity were calibrated considered in this analysis took place during low- manually. The default parameter values, which flow periods, for which the baseflow was found were deduced from landuse maps and pedologic almost constant and below 0.5 m³/s. maps, have been multiplied by a different The line shows a runoff coefficient of 2.1%. Most coefficient depending on the landuse category and measured coefficients are bounded by this soil type. minimum value. This is a lower bound, since In Figure 11, the rainfall intensities are plotted as rainfall events during which undrained areas spatially weighted means, while the dotted line contributed to the runoff production have a higher represents the measured hydrograph at the runoff to rainfall ratio. This suggests that the gauging station. Measurement uncertainty is drained impervious areas produce about 2.1% of limited since the gauging curve has been the total runoff. Considering a mean runoff extrapolated for high discharges using local 2D coefficient of 0.85 for these areas, to account for hydraulic modelling, instead of a mathematical losses (evapotranspiration, surface storage), the extrapolation. impervious surface of the catchment is therefore Modelling the drainage network adds a fast flow estimated at 2.1% / 0.85 = 2.5% of its total component, which significantly modifies the surface. This value is fully consistent with the computed hydrograph. This leads to a better value estimated previously on the basis of the representation of the observed flow rates, since landuse map and the drainage network. the fast flow component is indeed similar to the 53 Discharge [m³/s]Rainfall [mm/h]Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1 (2012) measurements. In contrast, the model predictions, but may not be considered as the only disregarding drainage failed to represent the first phenomenon determining the rapidly varying peaks, since the runoff production from flow components of river hydrographs. impervious areas does not reach the river fast The proposed methodology constitutes an original enough to produce the first peak. For this event, approach for the modelling of drained areas. The the Nash-Sutcliffe efficiency coefficient has been proposed simplified modelling of the drainage found to improve from 0.81 to 0.92 (Nash and network can be usefully incorporated into existing Sutcliffe, 1970). rainfall-runoff models and can be run Nevertheless, some differences remain between automatically, requiring limited extra effort from the two curves. Although runoff on drained areas the modeller. Further verification using additional generates a rapidly varying flow component, rainfall events and different catchments is still other processes also influence the discharge at an needed to assess more generally the performance hourly time steps. For instance, saturated areas of the methodology and the sensitivity of the located close to the river may similarly induce results. rapidly varying flows (Cosandey, 1996). The In future developments of the concept of virtual disaggregation of daily rainfall to hourly series is pipes, the main original contribution of the another source of uncertainty. present research, the focus should be set on the sizing of the pipes, in order to compute more 5. CONCLUSIONS accurately the flow propagation through the real network. Additionally, it would be of key Impervious surfaces can produce a significant part practical interest to identify objective criteria to of the surface runoff. In small and mid-size assess a priori the need for explicitly modelling watersheds, the flow dynamics of the water the drainage network, based on the watershed size falling on such impervious areas differs and its mean degree of urbanisation. significantly depending on whether these areas are connected or not to a drainage network, such as a sewage system. A methodology has been NOMENCLATURE developed to accurately compute the contributions of impervious surfaces and to take A river cross section [m²] into account the main effects of the drainage Deq equivalent pipe diameter [m] network without the need for its detailed h water depth [m] modelling. The impervious part of each cell hsub subsurface water dept [m] surface is computed as the sum of the impervious i slope of virtual pipe areas delimited in a detailed landuse map in J energy slope (accounting for bottom vector format. The real drainage network is roughness and internal losses) modelled by means of a simplified network, called virtual pipes, which routes the rain falling K Strickler coefficient of the river 1/3on drained cells to the river. The separation of [m /s] drained and undrained areas is based on a Ks lateral hydraulic conductivity [m/s] maximal distance criterion between the cells and l channel width [m] the nearest pipe from the drainage network. lb channel bottom width [m] In the case study of river Berwinne, the model not n Manning coefficient [s /m1/3] taking into account the drainage network failed to p p  pressure terms defined in (7) [m³] simulate the fast evolution of the river discharges, ω, xeven after a calibration process. In contrast, the [m²] improved capacity of the developed model to q discharge [m³/s] represent the river dynamics at an hourly time qL lateral exchanges (lateral inflow and step has been demonstrated in this case study. exchanges between mainstream and Nevertheless, the computed river discharges at floodplains) [m³/s] such small time scales may also be affected by R  the hydraulic radius [m] other hydrological processes (e.g. rainfall falling hon saturated areas near the rivers) and by S source terms (rainfall and infiltration) uncertainties due to data resolution (e.g. daily [m/s] rainfall time series disaggregated into hourly Sfi friction slopes (i=x,y) series). Therefore, accounting for the drained Sfsub,i subsurface friction slopes (i=x,y) impervious areas is necessary to improve flow Simp impervious surface drained by a 54 Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1 (2012) virtual pipe [m²] chaos. Water Resources Management Ssub source terms (infiltration, deep 22(7):895-909. percolation) [m/s] 11. Chow VT, Maidment D, Mays LW (1988). u, v velocity components along x and y Applied Hydrology. McGraw-Hill. 12. Cosandey C (1996). Surfaces saturées, axis [m/s] surfaces contributives : Localisation et usub, vsub interflow velocity components [m/s] extension dans l'espace du bassin versant. θ mean bottom slope Hydrological Sciences Journal 41(5):751-761. θsub,i projected slopes of the soil layer 13. Dewals BJ, Erpicum S, Archambeau P, (i=x,y) Detrembleur S, Pirotton M (2006). Depth-θ θ  projected bottom slopes integrated flow modelling taking into account x, y bottom curvature. J. Hydraul. Res. 44(6):787-795. REFERENCES 14. Dewals BJ, Kantoush SA, Erpicum S, Pirotton M, Schleiss AJ (2008). Experimental 1. Alley WM, Veenhuis JE (1983). Effective and numerical analysis of flow instabilities in impervious area in urban runoff modeling. rectangular shallow basins. Environ. Fluid Journal of Hydraulic Engineering Mech. 8:31-54. 109(2):313-319. 15. Dewals BJ, Erpicum S, Detrembleur S, 2. Aronica G, Cannarozzo M (2000). Studying Archambeau P, Pirotton M (2011). Failure of the hydrological response of urban dams arranged in series or in complex. catchments using a semi-distributed linear Natural Hazards 56(3):917-939. non-linear model. Journal of Hydrology 16. Dufresne M, Dewals BJ, Erpicum S, 238(1-2):35-43. Archambeau P, Pirotton M (2011). Numerical 3. Beighley E, Kargar M, He Y (2009). Effects investigation of flow patterns in rectangular of impervious area estimation methods on shallow reservoirs. Engineering Applications simulated peak discharges. Journal of of Computational Fluid Mechanics 5(2):247-Hydrologic Engineering 14(4):388-398. 258. 4. Boyd MJ (1993). Pervious and impervious 17. Endreny TA, Thomas KE (2009). Improving runoff in urban catchments. Hydrological estimates of simulated runoff quality and Sciences Journal 38(6):463-478. quantity using road-enhanced land cover data. 5. Brabec EA (2009). Imperviousness and land- Journal of Hydrologic Engineering use policy: Toward an effective approach to 14(4):346-351. watershed planning. Journal of Hydrologic 18. Ernst J, Dewals BJ, Detrembleur S, Engineering 14(4):425-433. Archambeau P, Erpicum S, Pirotton M (2010). 6. Callow JN (2007). How does modyfying a Micro-scale flood risk analysis based on DEM to reflect known hydrology affect detailed 2D hydraulic modelling and high subsequent terrain analysis? Journal of resolution land use data. Nat. Hazards Hydrology 332:30-39. 55(2):181-209. 7. Chabaeva A, Civco DL, Hurd JD (2009). 19. Erpicum S, Meile T, Dewals BJ, Pirotton M, Assessment of impervious surface estimation Schleiss AJ (2009). 2D numerical flow techniques. Journal of Hydrologic modeling in a macro-rough channel. Int. J. Engineering 14(4):377-387. Numer. Methods Fluids 61(11):1227-1246. 8. Chau KW, Wu CL, Li YS (2005). 20. Erpicum S, Dewals BJ, Archambeau P, Comparison of several flood forecasting Detrembleur S, Pirotton M (2010a). Detailed models in Yangtze River. Journal of inundation modelling using high resolution Hydrologic Engineering 10:485. DEMs. Engineering Applications of 9. Cheng CT, Ou CP, Chau KW (2002). Computational Fluid Mechanics 4(2):196-208. Combining a fuzzy optimal model with a 21. Erpicum S, Dewals BJ, Archambeau P, genetic algorithm to solve multi-objective Pirotton M (2010b). Dam-break flow rainfall-runoff model calibration. Journal of computation based on an efficient flux-vector Hydrology 268(1-4):72-86. splitting. Journal of Computational and 10. Cheng CT, Wang WC, Xu DM, Chau KW Applied Mathematics 234(7):2143-2151. (2008). Optimizing hydropower reservoir 22. Fleming G (1975). Computer Simulation operation using hybrid genetic algorithm and Techniques in Hydrology. American Elsevier Publishing Co., Inc. 55 Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1 (2012) 23. Garen DC, Moore DS (2005). Curve number 35. Mcmillan HK, Brasington J (2008). End-to-hydrology in water quality modeling: uses, end risk assessment: A coupled model abuses, and future directions. Journal of the cascade with uncertainty estimation. Water American Water Resources Association Resources Research 44(W03419):14. 41(2):377-388. 36. Moglen GE (2009). Hydrology and 24. Hjelmfelt AT, NRCS/ARS Curve Number impervious areas. Journal of Hydrologic Work Group, Moody HF (2004). Estimation Engineering 14(4):303-304. of direct runoff from storm rainfall. In: Part 37. Nash JE, Sutcliffe JV (1970). River flow 630 Hydrology National Engineering forecasting through conceptual models part I - Handbook. Ed. USDA. A discussion of principles. Journal of 25. Jain MK, Singh VP (2005). DEM-based Hydrology 10(3):282-290. modelling of surface runoff using diffusion 38. Nearing MA, Liu BY, Risse LM, Zhang X wave equation. Journal of Hydrology 302(1- (1996). Curve numbers and Green-Ampt 4):107-126. effective hydraulic conductivities. Water 26. Jones JW, Jarnagin T (2009). Evaluation of a Resources Bulletin 32(1). moderate resolution, satellite-based 39. NRCS/ARS Curve Number Work Group, impervious surface map using an independent, Moody HF (2004). Hydrologic soil-cover high-resolution validation data set. Journal of complexes. In: Part 630 Hydrology National Hydrologic Engineering 14(4):369-376. Engineering Handbook. Ed. USDA. 27. Kerger F, Archambeau P, Erpicum S, Dewals 40. Rawls WJ, Brakensiek DL (1989). Estimation BJ, Pirotton M (2011). A fast universal solver of soil water retention and hydraulic for 1D continuous and discontinuous steady properties. In: Unsaturated Flow in flows in rivers and pipes. Int. J. Numer. Hydrologic Modeling. Ed. Morel-Seytoux HJ. Methods Fluids 66:38-48. Kluwer Academic Publishers. 28. Khuat Duy B, Archambeau P, Dewals BJ, 41. Renssen H, Knoop JM (2000). Global river Erpicum S, Pirotton M (2010). River routing network for use in hydrological modelling and flood mitigation in a Belgian modeling. Journal of Hydrology 230(3-catchment. Proc. Inst. Civil. Eng.-Water 4):230-243. Manag. 163(8):417-423. 42. Rodriguez F, Andrieu H, Morena F (2008). A 29. Koutsoyiannis D (1994). A stochastic distributed hydrological model for urbanized disaggregation method for design storm and areas - Model development and application to flood synthesis. Journal of Hydrology 156(1- case studies. Journal of Hydrology 351(3-4):193-225. 4):268-287. 30. Koutsoyiannis D (2003). Rainfall 43. Roger S, Dewals BJ, Erpicum S, disaggregation methods: Theory and Schwanenberg D, Schüttrumpf H, Köngeter J, applications. In: Workshop on Statistical and Pirotton M (2009). Experimental und Mathematical Methods for Hydrological numerical investigations of dike-break Analysis, Rome. induced flows. J. Hydraul. Res. 47(3):349-31. Lhomme J, Bouvier C, Perrin JL (2004). 359. Applying a GIS-based geomorphological 44. Saunders W (1999). Preparation of DEMs for routing model in urban catchments. Journal use in environmental modelling analysis. of Hydrology 299(3-4):203-216. 1999 ESRI User Conference, San Diego, 32. Lin JY, Cheng CT, Chau KW (2006). Using California. support vector machines for long-term 45. Schmitt TG, Thomas M, Ettrich N (2004). discharge prediction. Hydrological Sciences Analysis and modeling of flooding in urban Journal 51(4):599-612. drainage systems. Journal of Hydrology 33. Machiels O, Erpicum S, Archambeau P, 299(3-4):300-311. Dewals BJ, Pirotton M (2011). Theoretical 46. SPGE (2007). Plans d'Assainissement par and numerical analysis of the influence of the Sous-bassins Hydrographiques - document de bottom friction formulation in free surface travail - : mise à jour juin 2007. flow modeling. Water SA 37(2):221-228. 47. Tarboton DG, Bras RL, Rodriguez-Iturbe I 34. Martz LW, Garbrecht J (1999). An outlet (1991). On the extraction of channel networks breaching algorithm for the treatment of from digital elevation data. Hydrological closed depressions in a raster DEM. Processes 5(1):81-100. Computers & Geosciences 25:835-844. 48. Thorndahl S, Johansen C, Schaarup-Jensen K (2006). Assessment of runoff contributing 56 Engineering Applications of Computational Fluid Mechanics Vol. 6, No. 1 (2012) catchment areas in rainfall runoff modelling. Water Science & Technology 54(6-7):59-56. 49. Valeo C, Moin SMA (2000). Variable source area modelling in urbanizing watersheds. Journal of Hydrology 228(1-2):68-81. 50. Walsh CJ, Fletcher TD, Ladson AR (2009). Retention capacity: A metric to link stream ecology and storm-water management. Journal of Hydrologic Engineering 14(4):399-406. 51. Wang WC, Chau KW, Cheng CT, Qiu L (2009). A comparison of performance of several artificial intelligence methods for forecasting monthly discharge time series. Journal of Hydrology 374(3-4):294-306. 52. Wu CL, Chau KW, Li YS (2009). Predicting monthly streamflow using data-driven models coupled with data-preprocessing techniques. Water Resources Research 45(8).  57