- Email: [email protected]
A hybrid method for flood simulation in small catchments combining hydrodynamic and hydrological techniques
Journal of Hydrology 540 (2016) 331–339
Contents lists available at ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
A hybrid method for flood simulation in small catchments combining hydrodynamic and hydrological techniques Vasilis Bellos ⇑, George Tsakiris Laboratory of Reclamation Works and Water Resources Management & Centre for the Assessment of Natural Hazards and Proactive Planning, School of Rural and Surveying Engineering, National Technical University of Athens, Athens, Greece
a r t i c l e
i n f o
Article history: Received 6 April 2016 Received in revised form 20 May 2016 Accepted 18 June 2016 Available online 20 June 2016 This manuscript was handled by K. Georgakakos, Editor-in-Chief, with the assistance of Jennifer Guohong Duan, Associate Editor Keywords: 2D flood modelling Small catchments FLOW-R2D model Unit hydrograph Kostiakov equation Green-Ampt model
a b s t r a c t The study presents a new hybrid method for the simulation of flood events in small catchments. It combines a physically-based two-dimensional hydrodynamic model and the hydrological unit hydrograph theory. Unit hydrographs are derived using the FLOW-R2D model which is based on the full form of two-dimensional Shallow Water Equations, solved by a modified McCormack numerical scheme. The method is tested at a small catchment in a suburb of Athens-Greece for a storm event which occurred in February 2013. The catchment is divided into three friction zones and unit hydrographs of 15 and 30 min are produced. The infiltration process is simulated by the empirical Kostiakov equation and the Green-Ampt model. The results from the implementation of the proposed hybrid method are compared with recorded data at the hydrometric station at the outlet of the catchment and the results derived from the fully hydrodynamic model FLOW-R2D. It is concluded that for the case studied, the proposed hybrid method produces results close to those of the fully hydrodynamic simulation at substantially shorter computational time. This finding, if further verified in a variety of case studies, can be useful in devising effective hybrid tools for the two-dimensional flood simulations, which are lead to accurate and considerably faster results than those achieved by the fully hydrodynamic simulations. Ó 2016 Elsevier B.V. All rights reserved.
1. Introduction It is widely accepted that the physically-based models cannot simulate the total complexity of the hydrological cycle using the full set of equations for each phenomenon (e.g. interception, infiltration, overland flow, evaporation, etc.), due to the large number of parameters which should be estimated, the bulk of data needed, and the prohibitive computational cost. Therefore, for the majority of the physically-based models only one or two principal processes of the hydrological cycle are simulated comprehensively, whereas all the other processes are simulated by simplified methods (i.e. SWAT software (Neitsch et al., 2002)). This study focuses on the hydrodynamic modelling for the simulation of the flood in a catchment created by a storm episode and especially on the two-dimensional (2D) modelling in which the full dynamic form of 2D Shallow Water Equations (2D-SWE) are used. It should be noted that in the majority of cases, the hydrodynamic simulation is applied only to the river channel or the floodplain (Abderrezzak et al., 2008; Bellos, 2012; Vacondio et al., 2016). In ⇑ Corresponding author. E-mail addresses: [email protected] (V. Bellos), [email protected] (G. Tsakiris). http://dx.doi.org/10.1016/j.jhydrol.2016.06.040 0022-1694/Ó 2016 Elsevier B.V. All rights reserved.
the following paragraphs, some of the efforts made for the 2D flood modelling at catchment scale, are briefly discussed. The first studies for flood flow simulation employing 2D distributed models, used simplified forms of mass and momentum conservation equations instead of the full dynamic form of 2D-SWE, due to the computational cost of the latter. Indicatively, interesting applications in this category were published by Tayfur et al. (1993), Feng and Molz (1997), Liu et al. (2004), Kazezyilmaz-Alhan and Medina (2007) and Gottardi and Venutelli (2008). However, the dynamics of the very shallow flows, such as the flow created by precipitation, (especially with low intensity rain) are quite complex and are characterised by mixed subcritical and supercritical conditions. According to some researchers, phenomena such as the backwater effect (Liang, 2010; Costabile et al., 2012b) or the river reach confluences, where momentum transfer is significant (Singh et al., 2014), was not able to be simulated satisfactorily by simplified approaches, such as the kinematic or diffusive wave methods. Therefore, the 2D Shallow Water Equations (2D-SWE) were used as an improvement of the simulation. Fiedler and Ramirez (2000) developed their algorithm solving the 2D-SWE, using the Finite Difference Method (FDM) and the McCormack numerical scheme. Infiltration is simulated through
332
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
the Green-Ampt method. In their study, numerical experiments were performed in a micro-scale simulating the rainfall-runoff process. Esteves et al. (2000) have also developed a numerical model solving the 2D-SWE with the same procedures (FDM and the McCormack numerical scheme and Green-Ampt method for the infiltration losses), which was validated using the experimental data derived from an experimental basin at a micro-scale, constructed in the Niger river (West Africa). Liang et al. (2007) combined 2D-SWE for overland flow with the Boussinesq equations for groundwater flow. Costabile et al. (2012a, 2013) have tested their model for various benchmark scenaria. They have also applied their model in a sub-basin of the Reno river (Italy), simulating an extreme rainfall episode. Kim et al. (2012) applied their 2D model, which incorporates precipitation, ice melting, infiltration and groundwater dynamics, in a real world episode which occurred in Peacheater Creek watershed of the Oklahoma state, USA. Singh et al. (2014) simulated a rainfall-runoff episode in Goodwin Creek watershed (Mississippi, USA). For the infiltration simulation the Green-Ampt method was used. Finally, the authors of this paper have proposed an in-house 2D fully hydrodynamic model, called FLOW-R2D, which is based on the solution of the Shallow Water Equations using a modified form of McCormack numerical scheme incorporating artificial viscosity (Tsakiris and Bellos, 2014). The model was tested recently in complex built-up areas (Bellos and Tsakiris, 2015a). A short description of this model is presented in the next section of this paper. Most of the recent applications of hydrodynamic models showed an increased accuracy in 2D flood simulations. However, although the speed of calculations has been considerably increased due to the technological innovations in computing facilities, the computational time for a fully hydrodynamic simulation (particularly in large watersheds) is still high and in many cases still unmanageable. This led researchers to redirect their efforts towards hybrid methods which exploit the merits of the hydrodynamic modelling but at the same time use also some simpler hydrological approaches aiming at reducing the computational time. A recent example in this direction is the hydrologicalhydraulic model proposed by Nguyen et al. (2015), called HiRes Flood-UCI. This model combines the hydrological model HLRDHM as a rainfall-runoff model, but instead of the classical hydrological routing scheme, it uses the 2D hydrodynamic BreZo model. The model has been successfully tested in a catchment of the Illinois river in USA using a decameter resolution. In the present study a novel hybrid method which combines both the hydrodynamic and the hydrological approaches is used, in order to simulate flood events in small catchments. The implementation of the method is performed in three steps: (a) the Unit Hydrographs of the catchment are derived using a 2D hydrodynamic model, (b) the effective rainfall depth of the storm episode is calculated using an appropriate infiltration model, (c) the final flood hydrograph of the event is generated based on the principle of hydrograph superposition. The performance of the proposed method is tested in a flash flood event which occurred in a small catchment of the Halandri stream, in a suburb of Athens-Greece, in February 2013. For the derivation of the Unit Hydrographs, the in-house 2D hydrodynamic model FLOW-R2D is used and for the infiltration process, both the Kostiakov equation and the Green-Ampt model, are used. The results from the proposed hybrid method in the above event are compared directly with the available recorded data and the results derived from the fully hydrodynamic simulation of the entire flood episode. 2. Theoretical background For the hydrodynamic part of the hybrid method, the numerical model FLOW-R2D is used, based on the solution of the twodimensional Shallow Water Equations (2D-SWE) by the Finite
Difference Method, in a cell-centred, non-staggered computational grid (Tsakiris and Bellos, 2014):
@W @F @G þ þ ¼D @t @x @y where
h uh 2 W ¼ uh ; F ¼ u2 h þ 0:5gh ; vh uv h rf D ¼ ghS0;x sb;x =q ghS0;y sb;y =q
ð1Þ vh G¼ uv h ; 2 v 2 h þ 0:5gh
ð2Þ
and h is the water depth, u is the horizontal component of flow velocity at x-direction, v is the horizontal component of flow velocity at y-direction, S0,x, S0,y are the bottom slopes for the x and y-directions, respectively, and q is the fluid density. Shear stresses sb,x and sb,y can be modelled by various methods such as the Manning equation which is used in the present study, as follows:
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
sb;x gn2 uð u2 þ v 2 Þ sb;y gn2 v ð u2 þ v 2 Þ ¼ ; ¼ q q R1=3 R1=3
ð3Þ
where n is the Manning coefficient and R is the hydraulic radius. The two terms, r and f in the Continuity Equation, represent the precipitation rate (source term) and the rate of losses, such as infiltration or drainage (sink term), respectively. It should be noted that infiltration and rainfall are not incorporated into the Momentum Equations of the model because their magnitude is very small in comparison with the other components. A modification of the explicit McCormack numerical scheme (McCormack, 1969) is used to solve the 2D-SWE. It is noted that this scheme has been also used by several researchers (i.e. Kalita, 2016). With this modification, artificial viscosity is added through a diffusion factor. The oscillatory errors are smoothed out while the shock capturing capability of the scheme is retained (the solution has still second order accuracy) and therefore discontinuities of the flow can be simulated. Wet/dry modelling is achieved through a water depth threshold which distinguishes the wet and the dry cells. In the dry cells the water depth and flow velocity are taken as zeros. In the following equations, the discretisation of 2D-SWE using the McCormack numerical scheme is obtained in two steps (predictor strep in Eq. (4) and corrector step in Eq. (5)):
W i;j ¼ W ki;j
Dt Dt k F iþ1;j F ki;j Gki;jþ1 Gki;j þ DtDki;j Dx Dy
ð4Þ
1 1 xW ki;j þ ð1 xÞ W kiþ1;j þ W ki1;j þ W ki;jþ1 þ W ki;j1 2 4 Dt Dt ð5Þ F i;j F i1;j ðGi;j Gi;j1 Þ þ DtDi;j þW i;j Dx Dy
¼ W kþ1 i;j
where Dt, Dx, Dy are the time step, the space step in x- and the space step in y-direction, respectively, and x is the diffusion factor. Finally, the non-rectangular computational domains are represented by a pseudo-computational rectangular domain which encloses the real domain. A detailed description of the FLOWR2D can be found in Tsakiris and Bellos (2014). An extensive grid convergence study of the model has been also performed for steady and dynamic evolution state applications and can be found elsewhere (Bellos and Tsakiris, 2016). Further, interesting applications of the model in urban environments can be found in Bellos and Tsakiris (2015a, 2015b). In order to calculate the effective rainfall depth which is transformed into direct runoff, the infiltration losses should be determined. Various empirical equations have been proposed and
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
used in the past for the infiltration simulation, such as the equations proposed by Kostiakov, Horton, Philip and several others physically-based such as the Green-Ampt model (Tayfur et al., 1993; Tayfur and Kavvas, 1998; Borah, 2011). In this study, both the empirical Kostiakov equation and the Green-Ampt model are used for the simulation of the infiltration process. The first one is chosen, due to its simplicity and the small number of parameters for determination, and the second one for its accuracy and broad acceptance. As known, according to Kostiakov, the infiltration rate f(t) can be calculated by the following equation:
f ðtÞ ¼ abt
b1
ð6Þ
where a and b are parameters estimated by experiments and related to the soil type, and t is the time elapsed from the beginning of the infiltration process. According to the Green-Ampt model, the infiltration rate is calculated as follows:
f ðtÞ ¼ K
wDh þ1 FðtÞ
ð7Þ
where K is the hydraulic conductivity of the soil, w is the wetting front soil suction head, Dh is the water content change, and F(t) is the cumulative infiltration determined by the following implicit equation:
FðtÞ FðtÞ ¼ Kt þ wDhln 1 þ wDh
ð8Þ
3. Description of the case-study The FLOW-R2D model is applied in a small catchment (about 5 km2) of the Halandri stream, in the north suburbs of Athens (Greece). Although the greater part of the catchment is natural, urban zones are identified in the downstream section of the catchment. No manmade structures are located in the hydrographic network of the catchment. In Fig. 1a, an aerial imagery of the catchment is shown with the location of the hydrometric station. Also the land use map based on the information of CORINE 2000 is presented in Fig. 1b. For simplicity purposes, the five zones of
333
CORINE 2000 are merged into three friction zones (Fig. 1c): the coniferous forest zone (friction zone A), the discontinuous urban fabric with the mineral extraction sites zone (friction zone B) and the transitional woodland-shrub with the sclerophyllous vegetation zone (friction zone C). The geology of the catchment consists of marble. The maximum, minimum and average surface elevation of the catchment are 1045 m, 325 m and 609 m asl, respectively. The maximum length of river reach of the watershed is 4.75 km. The average slope is about 32% and 12% for the catchment and the main river reach, respectively. Based on the above characteristics, the concentration time of the catchment is calculated 1.19 h by the Giandotti equation and 0.50 h by the Kirpich equation. Finally, various geomorphological factors are calculated for the catchment. The Gravelius factor is Kc = 1.43, the form factor F = 0.27, the catchment circularity Rc = 0.49 and the catchment elongation factor RL = 0.59. As known, the Gravelius factor is an indicator for the calculation of the compactness of the catchment, the form factor depends on the area and the length of the catchment, the circularity depends on the area and the perimeter of the catchment and the elongation factor depends on the area and the length of the catchment. On the 22nd of February 2013, a storm episode occurred in Athens creating significant casualties, especially in the northern suburbs of the city. In the case-study area, the precipitation depth reached 121.4 mm within 5 h, causing a flash flood and hitting the catchment of Halandri stream. The rainfall event was recorded at the meteorological station of Palea Penteli and the flood hydrograph at the hydrometric station of Nea Penteli. These data were recorded in the framework of DEUCALION Project (http:// deucalionproject.gr/). Based on these data, the runoff coefficient was estimated to be equal to 0.49. In Fig. 2, the hyetograph and the flood hydrograph recorded at the outlet of the catchment are presented. The storm episode of the case-study was divided in two parts: the first one is from the onset of the event at 5:30 up to 9:00 and the second from 9:00 to 11:30.
4. The hybrid method As mentioned above, the proposed hybrid method is implemented in three steps: (a) the Unit Hydrographs of the catchment
Fig. 1. (a) Aerial imagery of the Halandri catchment, (b) land use of the case-study according to CORINE 2000, (c) merged land use zones to the three friction zones.
334
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
Fig. 2. The hyetograph and hydrograph of the rainfall-runoff episode of the 22nd of February 2013.
are derived using a 2D hydrodynamic model. In the present case study, the FLOW-R2D model is used for their derivation. (b) The effective rainfall depth of the storm episode is calculated using an appropriate infiltration model, such as the Kostiakov equation and the Green-Ampt model. (c) The final flood hydrograph of the event is generated based on the principle of hydrograph superposition.
4.1. Unit hydrographs derivation In this section, the derivation of the Unit Hydrographs is presented, implementing the 2D FLOW-R2D model and using the source term of the Continuity Equation (Eq. (1)) for the simulation of the precipitation. Two types of Unit Hydrograph were derived for the rainfall durations: 15 min (UH-15) and 30 min (UH-30). As far as the friction is concerned, a roughness coefficient is calibrated for each friction zone. As easily understood, due to the high computational cost, no optimisation method could be implemented apart from the conventional trial and error method. For the calibration phase, various friction scenaria which have the form S1S2S3 are compared. The objective function is the Nash-Sutcliffe efficiency coefficient (Nash and Sutcliffe, 1970) between recorded data and the results of runoff simulation of the storm episode, obtained by the composition of the unit hydrographs:
Pn ðQ Q obs Þ2 NSE ¼ 1 Pi¼1 num 2 n i¼1 ðQ obs Q obs Þ
ð9Þ
where Qnum is the simulated and Qobs the observed flow, respectively. The first component of the friction scenario S1 represents the Manning coefficient for the friction zone A and takes values 1 and 2 for Manning friction coefficient values 0.050 and 0.070 s/m1/3, respectively. The second component of the friction scenario S2 represents the Manning coefficient for the friction zone B and takes values 1 and 2 for n values 0.020 and 0.035 s/m1/3, respectively. Finally, the third component of the friction scenario S3 represents the Manning coefficient for the friction zone C and takes values 1 and 2 for n values 0.025 and 0.050 s/m1/3, respectively. These
values are selected based on the related literature but increased slightly to take into account the sub-grid effects (Chaudhry, 2008). The computational cell size of the Digital Terrain Model (DTM) of the catchment is of 10 10 m. In total, the computational area consists of 50,811 cells. It should be noted that a preliminary investigation is performed in order to choose the appropriate grid size. It has been found that the coarser grids (e.g. 25 25 m) result in unrealistic outcomes, underestimating the main magnitudes of the phenomenon. On the contrary, finer grids (e.g. 5 5 m) increase considerably the computational burden. The time step equal to Dt = 0.004 s for the UH-15 and Dt = 0.005 s for the UH30, are chosen based on the authors experience. In real world applications, the diffusion factor commonly takes values around 0.90 (Bellos and Tsakiris, 2015b). However, the overland flow in the catchment is characterised by low water depths and no intense turbulence phenomena (apart from streams), and therefore the diffusion factor can take relatively higher values. In this work, the diffusion factor is set equal to x = 0.99. The wet/dry threshold is chosen hdry = 104 m, while a small water depth equal to 103 m over the entire catchment is assumed as the initial water depth. This thin film of water is used as initial condition. It can be noticed that using only the water depth threshold of the wet/dry modelling option, the time step should be further reduced to grasp the resulting high velocities which occur in the catchment. For the model simulation, solid boundaries are imposed to the catchment boundaries, in order to preserve the water mass balance. Since the simulation is referred to the effective rainfall, the infiltration process is not simulated except from a constant infiltration rate equal to the initial water depth for achieving 1 cm equivalent runoff depth. This rate is calculated f0 = 0.4 cm/h for the UH-15 derivation and f0 = 0.2 cm/h for the UH-30 derivation, respectively. In order to preserve the mass balance between the effective rainfall water volume and the runoff water volume, the results of the FLOW-R2D model are slightly modified through a correction factor (Table 1). The slight mismatch between these volumes appears due to the fact that the simulations are decided to be terminated when the contribution of the Unit Hydrographs spans up to a specific time limit, which was set to 3 h. This limit is imposed to save computational cost. However, the runoff process continues, but with very small rates. In Fig. 3, the two types of the Unit Hydrographs
335
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339 Table 1 Performance of the hybrid method against the observed data for the selected scenaria. Friction scenario
Unit hydrograph duration
Correction factor
Infiltration losses
NSE
MAE
RMSE
TPE (%)
PRE (%)
1.1.1 1.1.2 1.2.1 1.2.2 2.1.1 2.1.2 2.2.1 2.2.2
UH-15 UH-15 UH-15 UH-15 UH-15 UH-15 UH-15 UH-15
1.10 1.17 1.18 1.31 1.11 1.22 1.22 1.41
Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov
0.804 0.603 0.689 0.381 0.805 0.567 0.638 0.298
2.930 3.373 4.335 4.431 3.987 4.288 5.833 5.416
5.112 5.537 7.273 7.209 6.444 6.618 9.087 8.782
5.263 5.263 15.789 10.526 10.526 10.526 15.789 15.789
24.686 26.118 38.861 39.952 22.742 24.273 41.108 41.971
1.1.1 1.1.2 1.2.1 1.2.2 2.1.1 2.1.2 2.2.1 2.2.2
UH-15 UH-15 UH-15 UH-15 UH-15 UH-15 UH-15 UH-15
1.07 1.13 1.13 1.23 1.07 1.13 1.14 1.25
Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt
0.770 0.610 0.672 0.422 0.765 0.574 0.625 0.345
3.088 3.601 4.580 4.706 4.368 4.584 6.266 5.743
5.099 5.594 7.599 7.534 6.946 7.071 9.679 9.346
5.263 5.263 15.789 15.789 10.526 10.526 21.053 21.053
22.184 23.700 35.606 37.005 21.117 22.615 38.163 39.476
1.1.1 1.1.2 1.2.1 1.2.2 2.1.1 2.1.2 2.2.1 2.2.2
UH-30 UH-30 UH-30 UH-30 UH-30 UH-30 UH-30 UH-30
1.10 1.17 1.18 1.31 1.11 1.22 1.22 1.41
Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov
0.778 0.645 0.693 0.433 0.743 0.600 0.664 0.377
3.146 3.509 4.166 4.442 4.072 4.239 5.527 5.259
5.445 5.870 6.882 6.889 6.404 6.543 8.694 8.419
10.526 10.526 15.789 5.263 10.526 10.526 10.526 10.526
22.435 24.901 36.504 37.363 26.775 28.741 41.283 41.905
1.1.1 1.1.2 1.2.1 1.2.2 2.1.1 2.1.2 2.2.1 2.2.2
UH-30 UH-30 UH-30 UH-30 UH-30 UH-30 UH-30 UH-30
1.07 1.13 1.13 1.23 1.07 1.13 1.14 1.25
Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt
0.742 0.644 0.679 0.469 0.709 0.606 0.657 0.416
3.468 3.784 4.451 4.625 4.219 4.371 5.840 5.542
5.850 6.234 7.301 7.250 6.699 6.760 9.113 8.825
10.526 10.526 15.789 15.789 10.526 10.526 21.053 21.053
18.904 21.497 36.809 38.362 32.951 34.623 39.024 40.630
are shown, for all the friction scenaria derived from the FLOW-R2D model.
4.2. Effective rainfall determination Following the derivation of Unit Hydrographs, the effective rainfall is determined. As mentioned earlier, the Kostiakov equation and the Green-Ampt model were used for infiltration simulation. Since the catchment consists of marble, infiltration is expected to be high. Hence, the required parameters for each method are selected assuming that the infiltration process for marble is similar to that of sandy loam. Specifically, the Kostiakov equation parameters are estimated: a = 3.584 cm/h0.584 and b = 0.584 (Antonopoulos, 1999) and the Green-Ampt model
parameters are estimated: K = 1.09 cm/h, w = 2.67 cm and Dh = 0.283 (Tsakiris, 2013). In Fig. 4, the infiltration rate as calculated from the Kostiakov equation and the Green-Ampt model, are plotted parallelly to the rainfall intensity. It is observed that these methods converge, except of some discrepancies at the initial time period. 4.3. Flood hydrograph simulation The final step of the proposed method is the composition of the flood hydrograph based on the principle of hydrological superposition. In Figs. 5 and 6, the simulated hydrographs for each type of Unit Hydrograph, friction scenario and method used to simulate the infiltration process are shown, compared with the observed hydrograph. In Table 1, a comparison through the Nash efficiency
Fig. 3. The 15 min and the 30 min unit hydrographs derived from the FLOW-R2D for the selected friction scenaria.
336
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
Fig. 4. The rate of infiltration losses according to the Kostiakov equation and the Green-Ampt model plotted parallelly to the rainfall intensity.
Fig. 5. Observed data vs storm episode simulation results composing 15 min and 30 min unit hydrographs with the Kostiakov equation for simulating infiltration, for several scenaria.
Fig. 6. Observed data vs storm episode simulation results composing 15 min and 30 min unit hydrographs with Green-Ampt infiltration model, for several scenaria.
coefficient for all the scenaria, is presented. According to the NSE results, the best performance is achieved for the flood hydrograph created by the UH-15, the friction scenario 2.1.1, while the infiltration losses are determined by the Kostiakov equation. Apart from the Nash coefficient, additional indicators are calculated, for the method performance estimation against the observed data: the Mean Absolute Error (MAE), the Root Mean Square Error (RMSE), the time to peak error (TPE) and the peak rate error (PRE). The values of these indicators are also presented in Table 1.
An effort is also made to treat the parameters of the infiltration process simulation methods as black box parameters which should be calibrated in each case. Although the results of this procedure fit better to the observed data, the parameters take unrealistic values according to the related literature. For example using the UH-15 and the 2.1.1 friction scenario, the calibrated parameters for the Kostiakov equation are a = 211.970 cm/h0.183 and b = 0.183. Also for the Green-Ampt model the calibrated parameters take the following values: K = 0.003 cm/h, w = 2933.78 cm and Dh = 1.0. The
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
337
Fig. 7. Comparison between observed data and simulation results obtained by the hybrid method assuming the infiltration simulation methods models as black box models.
Fig. 8. Flood in the hydrographic network produced by FLOW-R2D model 30 min from the beginning of the rainfall event (UH-15, scenario 2.1.1).
calibration is made through the maximisation of the NSE values and the only physical constraint which is adopted is that Dh should be equal or less than one. The NSE values are optimised as 0.818 and 0.812 for the two hydrographs, respectively. The results obtained from this procedure against observed data are presented in Fig. 7. Finally, for illustration purposes, the flood in the hydrographic network of the catchment produced by the FLOW-R2D model is shown in Fig. 8. This snapshot refers to 30 min from the beginning of rainfall event, for the UH-15 and the 2.1.1 scenario. 5. The fully hydrodynamic method In this section, the entire storm episode is simulated using the FLOW-R2D model applied at a catchment scale. Details of the fully hydrodynamic model can be found elsewhere (Tsakiris and Bellos, 2014; Bellos and Tsakiris, 2015a). Due to the high computational cost, a single simulation run is performed using the output
obtained from the hybrid method. Specifically, the Kostiakov equation is used for the infiltration losses and the friction scenario 2.1.1 for the Manning coefficient values is adopted. All the input parameters of the model (diffusion factor and wet/dry threshold) remain unchanged in relation to the hybrid method. The time step is selected Dt = 0.004 s. Fig. 9 presents the comparison between the flood hydrographs obtained by the observed data, the fully hydrodynamic method and the best results from the proposed method. The performance of the fully hydrodynamic method towards measured data by the Nash performance coefficient is NSE = 0.860. As expected, the fully hydrodynamic method performs slightly better than the hybrid method but with much higher computational cost. In the present case-study, the computational time for the fully hydrodynamic method is three to four times greater compared to the hybrid method. Specifically, the hybrid method requires about 3–4 days CPU time for an i7 CPU/3.6 GHz/16 GB RAM, whereas the fully hydrodynamic method requires 10–12 days.
338
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
Fig. 9. Comparison between observed data and simulation results obtained by the fully hydrodynamic method and the hybrid method.
6. Discussion Although the rainfall depth is of the same magnitude at the two parts of the storm episode (rainfall peak of the first part is about 70% of the corresponding rainfall peak of the second part), a great difference can be noticed in the runoff records. The runoff peak which occurs during the first part is just about 30% of the corresponding runoff peak of the second part. This is explained by the fact that during the first part, the soil becomes fully saturated and in the second part, nearly all the rainfall depth is transformed into runoff. It can be seen that the Unit Hydrographs of 15 min rainfall duration give better results from the Unit Hydrographs of 30 min, in which significant oscillations in the simulated flood hydrograph are observed, most probably due to the time step. Besides, there are not significant differences in the results by the use of either the Kostiakov equation or of the Green-Ampt model for the infiltration process simulation, except during the first few time moments of the storm episode. It can be concluded that, unless there are available reliable data, the Kostiakov equation can be the appropriate choice due to its simplicity. Based on the Nash-Sutcliffe performance coefficient, the UH-15 with the 2.1.1 friction scenario performs better than all others. The above combination gives also better results investigating the capability of the proposed method to simulate the peak of the flood hydrograph. Specifically, the peak of the flood hydrograph simulated with the UH-15, the 2.1.1 friction scenario and the Kostiakov equation is 96% of the observed peak. The timing of the peak flow for both of the above cases is delayed by 20 min from the peak based on the observed data. As far as the calibration phase is concerned, friction seems to be also a critical factor for the flood hydrograph simulation. Indicatively, the simulated peak flow of the entire episode for the worst friction scenario 2.2.2 (according to the NSE values), is about 71% of the corresponding best friction scenario 2.1.1 and has 35 min time lag, for the UH-15 and for both methods of infiltration simulation. For the UH-30 and the worst friction scenario 2.2.2 is about 72% of the corresponding better friction scenario 1.1.1 and has 34 min time lag, for the UH-15 and for both methods of infiltration simulation. The most influential zone towards the final results is zone C, followed by zone B. Friction zone A seems to have no significant influence on the final results. This was rather expected, due to the fact that zone A represents a small percentage of the catchment area. An alternative approach for the calibration phase is to assume the parameters required for the infiltration simulation as black box parameters. This approach performs better if the simulated
flood hydrograph is compared to the observed data. On the other hand, the parameters take unrealistic values according to the related literature. Therefore, this approach is not recommended, unless the catchment case study is gauged and at least some flood events are recorded, in order to implement the calibration and the validation phase. A further aspect which should be examined in detail is the convergence of the proposed method as the computational grid becomes finer. It is to note that the great computational cost did not permit an extensive convergence analysis in the context of the specific case study. This analysis should probably be performed in a smaller catchment with milder land surface slopes, in which the computational domain will be represented by a smaller number of cells and greater time step. From the case studied there are strong indications that the proposed hybrid method could be successfully used in devising appropriate design tools for the flood hydrograph simulation of small catchments in various storm episodes, even if their duration is relatively long. Besides, in the context of this method, the calibration phase by the trial and error method is still applicable. Needless to say that the fully hydrodynamic method produces relatively better results than the hybrid method according to the NSE values. The peak of the flood hydrograph obtained from this procedure is 111% of the observed peak, while the time of the peak is just 4 min delayed from the respective observed peak. 7. Concluding remarks A novel hybrid method for the 2D flood simulation was presented in this study. The method combines the 2D hydrodynamic modelling and the unit hydrograph theory. The method was applied to a small catchment in a suburb of Athens for a single storm event. The results from the method were compared to the recorded data at a hydrometric station at the outlet of the catchment and the results derived by the application of the fully hydrodynamic simulation model, FLOW-R2D. For the case studied, the proposed hybrid method showed a satisfactory performance according to the Nash-Sutcliffe efficiency coefficient. Also, the obtained results are close to those produced by the fully hydrodynamic simulation, although the required computational time was considerably reduced. This finding is promising provided that it is further verified in other catchments and conditions. Then, the proposed hybrid method may be useful for devising effective hybrid tools for relatively accurate 2D flood simulations which can be attained at a reasonable computational cost.
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
Acknowledgments This paper is based partially on a previous short presentation made at the 9th World Congress of the European Water Resources Association, Istanbul, 24–27 June 2015. Acknowledgments are expressed to D. Koutsoyiannis and A. Efstratiadis for providing useful data on the case-study, to the National Cadastre & Mapping Agency of Greece for the provision of the DTM maps, to I. Papageorgaki for the handling of the DTM, and finally to the Computer Centre of NTUA for the access provided to the computer facilities. References Abderrezzak, K., Paquier, A., Mignot, E., 2008. Modelling flash flood propagation in urban areas using a two-dimensional numerical model. Nat. Hazards 50, 433– 461. Antonopoulos, V., 1999. Hydrology of Unsaturated Soil Zone. AUTH editions (in Greek). Bellos, V., 2012. Ways for flood hazard mapping in urbanised environments: a short literature review. Water Utility J. 4, 25–31. Bellos, V., Tsakiris, G., 2015a. Comparing various methods of building representation for 2D flood modelling in built-up areas. Water Resour. Manage. 29, 379–397. Bellos, V., Tsakiris, G., 2015b. 2D flood modelling: the case of Tous dam break. In: Proceedings of 36th World Congress of IAHR, the Hague, the Netherlands (e-proceedings). Bellos, V., Tsakiris, G., 2016. Grid coarsening and uncertainty in 2D hydrodynamic modelling. 4th European Congress of IAHR, Liege, Belgium (accepted for oral presentation). Borah, K.D., 2011. Hydrologic procedures of storm event watershed models: a comprehensive review and comparison. Hydrol. Process. 25 (22), 3472–3489. Chaudhry, M.H., 2008. Open Channel Flow, second ed. Springer. Costabile, P., Costanzo, C., Macchione, F., 2012a. Comparative analysis of overland flow using finite volume schemes. J. Hydroinformatics 14 (1), 122–135. Costabile, P., Costanzo, C., Macchione, F., Mercogliano, P., 2012b. Two-dimensional model for overland flow simulations: a case-study. Eur. Water 38, 13–23. Costabile, P., Costanzo, C., Macchione, F., 2013. A storm event watershed model for surface runoff based on 2D fully dynamic wave equations. Hydrol. Process. 27 (4), 554–569. Esteves, M., Faucher, X., Galle, S., Vauclin, M., 2000. Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed results. J. Hydrol. 228 (3–4), 265–282. Feng, K., Molz, G.J., 1997. A 2-D diffusion-based wetland flow model. J. Hydrol. 196 (1–4), 230–250.
339
Fiedler, F.R., Ramirez, J.A., 2000. A numerical method for simulating discontinuous shallow flow over an infiltrating surface. Int. J. Numer. Methods Flows 32, 219– 240. Gottardi, G., Venutelli, M., 2008. An accurate time integration method for simplified overland flow models. Adv. Water Resour. 31 (1), 173–180. Kalita, H.M., 2016. A new Total Variation Diminishing predictor corrector approach for two-dimensional shallow water flow. Water Resour. Manage. 30 (4), 1481– 1497. Kazezyilmaz-Alhan, C., Medina, M.A., 2007. Kinematic and diffusion waves: analytical and numerical solutions to overland and channel flow. J. Hydraul. Eng. 133 (2), 217–228. Kim, J., Warnock, A., Ivanov, V.Y., Katopodes, N.D., 2012. Coupled modeling of hydrologic and hydrodynamic processes including overland and channel flow. Adv. Water Resour. 37, 104–126. Liang, D., Falconer, R.A., Lin, B., 2007. Coupling surface and subsurface flows in a depth averaged flood wave model. J. Hydrol. 337, 147–158. Liang, J., 2010. Evaluation of Runoff Response to Moving Rainstorms. Dissertations (2009); Paper 73, Marquette University, Wisconsin, USA.. Liu, Q.Q., Chen, L., Li, J.C., Singh, V.P., 2004. Two-dimensional kinematic wave model of overland-flow. J. Hydrol. 291 (1–2), 28–41. McCormack, R.W., 1969. The effect of viscosity in hypervelocity impact cratering. AIAA Hypervelocity Impact Conference, Cincinnati, Paper 69-354. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models. Part I – A discussion of principles. J. Hydrol. 10 (3), 282–290. Neitsch, S.L., Arnold, J.G., Kiniry, J.R., Williams, J.R., King, K.W., 2002. Soil and Water Assessment Tool, Theoretical Documentation, Version 2000. Texas Water Resources Institute. Nguyen, P., Thorstensen, A., Sorooshian, S., Hsu, K., AghaKouchack, A., Sanders, B., Koren, V., Cui, Z., Smith, M., 2015. A high resolution coupled hydrologichydraulic model (HiResFlood-UCI) for flash flood modelling. J. Hydrol. http://dx. doi.org/10.1016/j.jhydrol.2015.10.047. Singh, J., Altinakar, M., Ding, Y., 2014. Numerical modeling of rainfall-generated overland flow using nonlinear Shallow-Water Equations. J. Hydrol. Eng. 20 (8), 04014089. http://dx.doi.org/10.1061/(ASCE)HE.1943-5584.0001124. Tayfur, G., Kavvas, M.L., Govindaraju, R.S., Storm, D.E., 1993. Applicability of St. Venant equations for two-dimensional overland flows over rough infiltrating surfaces. J. Hydraul. Eng. 119 (1), 51–63. Tayfur, G., Kavvas, M.L., 1998. Areally-averaged overland flow equations at hillslope scale. Hydrol. Sci. J. 43 (3), 361–378. Tsakiris, G., 2013. Water Resources: I. Engineering Hydrology. Symmetria editions (in Greek). Tsakiris, G., Bellos, V., 2014. A numerical model for two-dimensional flood routing in complex terrains. Water Resour. Manage. 28 (5), 1277–1291. Vacondio, R., Aureli, F., Ferrari, A., Mignosa, P., Dal Palù, A., 2016. Simulation of the January 2014 flood on the Secchia river using a fast and high-resolution 2D parallel shallow-water numerical scheme. Nat. Hazards 80, 103–125.
Contents lists available at ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
A hybrid method for flood simulation in small catchments combining hydrodynamic and hydrological techniques Vasilis Bellos ⇑, George Tsakiris Laboratory of Reclamation Works and Water Resources Management & Centre for the Assessment of Natural Hazards and Proactive Planning, School of Rural and Surveying Engineering, National Technical University of Athens, Athens, Greece
a r t i c l e
i n f o
Article history: Received 6 April 2016 Received in revised form 20 May 2016 Accepted 18 June 2016 Available online 20 June 2016 This manuscript was handled by K. Georgakakos, Editor-in-Chief, with the assistance of Jennifer Guohong Duan, Associate Editor Keywords: 2D flood modelling Small catchments FLOW-R2D model Unit hydrograph Kostiakov equation Green-Ampt model
a b s t r a c t The study presents a new hybrid method for the simulation of flood events in small catchments. It combines a physically-based two-dimensional hydrodynamic model and the hydrological unit hydrograph theory. Unit hydrographs are derived using the FLOW-R2D model which is based on the full form of two-dimensional Shallow Water Equations, solved by a modified McCormack numerical scheme. The method is tested at a small catchment in a suburb of Athens-Greece for a storm event which occurred in February 2013. The catchment is divided into three friction zones and unit hydrographs of 15 and 30 min are produced. The infiltration process is simulated by the empirical Kostiakov equation and the Green-Ampt model. The results from the implementation of the proposed hybrid method are compared with recorded data at the hydrometric station at the outlet of the catchment and the results derived from the fully hydrodynamic model FLOW-R2D. It is concluded that for the case studied, the proposed hybrid method produces results close to those of the fully hydrodynamic simulation at substantially shorter computational time. This finding, if further verified in a variety of case studies, can be useful in devising effective hybrid tools for the two-dimensional flood simulations, which are lead to accurate and considerably faster results than those achieved by the fully hydrodynamic simulations. Ó 2016 Elsevier B.V. All rights reserved.
1. Introduction It is widely accepted that the physically-based models cannot simulate the total complexity of the hydrological cycle using the full set of equations for each phenomenon (e.g. interception, infiltration, overland flow, evaporation, etc.), due to the large number of parameters which should be estimated, the bulk of data needed, and the prohibitive computational cost. Therefore, for the majority of the physically-based models only one or two principal processes of the hydrological cycle are simulated comprehensively, whereas all the other processes are simulated by simplified methods (i.e. SWAT software (Neitsch et al., 2002)). This study focuses on the hydrodynamic modelling for the simulation of the flood in a catchment created by a storm episode and especially on the two-dimensional (2D) modelling in which the full dynamic form of 2D Shallow Water Equations (2D-SWE) are used. It should be noted that in the majority of cases, the hydrodynamic simulation is applied only to the river channel or the floodplain (Abderrezzak et al., 2008; Bellos, 2012; Vacondio et al., 2016). In ⇑ Corresponding author. E-mail addresses: [email protected] (V. Bellos), [email protected] (G. Tsakiris). http://dx.doi.org/10.1016/j.jhydrol.2016.06.040 0022-1694/Ó 2016 Elsevier B.V. All rights reserved.
the following paragraphs, some of the efforts made for the 2D flood modelling at catchment scale, are briefly discussed. The first studies for flood flow simulation employing 2D distributed models, used simplified forms of mass and momentum conservation equations instead of the full dynamic form of 2D-SWE, due to the computational cost of the latter. Indicatively, interesting applications in this category were published by Tayfur et al. (1993), Feng and Molz (1997), Liu et al. (2004), Kazezyilmaz-Alhan and Medina (2007) and Gottardi and Venutelli (2008). However, the dynamics of the very shallow flows, such as the flow created by precipitation, (especially with low intensity rain) are quite complex and are characterised by mixed subcritical and supercritical conditions. According to some researchers, phenomena such as the backwater effect (Liang, 2010; Costabile et al., 2012b) or the river reach confluences, where momentum transfer is significant (Singh et al., 2014), was not able to be simulated satisfactorily by simplified approaches, such as the kinematic or diffusive wave methods. Therefore, the 2D Shallow Water Equations (2D-SWE) were used as an improvement of the simulation. Fiedler and Ramirez (2000) developed their algorithm solving the 2D-SWE, using the Finite Difference Method (FDM) and the McCormack numerical scheme. Infiltration is simulated through
332
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
the Green-Ampt method. In their study, numerical experiments were performed in a micro-scale simulating the rainfall-runoff process. Esteves et al. (2000) have also developed a numerical model solving the 2D-SWE with the same procedures (FDM and the McCormack numerical scheme and Green-Ampt method for the infiltration losses), which was validated using the experimental data derived from an experimental basin at a micro-scale, constructed in the Niger river (West Africa). Liang et al. (2007) combined 2D-SWE for overland flow with the Boussinesq equations for groundwater flow. Costabile et al. (2012a, 2013) have tested their model for various benchmark scenaria. They have also applied their model in a sub-basin of the Reno river (Italy), simulating an extreme rainfall episode. Kim et al. (2012) applied their 2D model, which incorporates precipitation, ice melting, infiltration and groundwater dynamics, in a real world episode which occurred in Peacheater Creek watershed of the Oklahoma state, USA. Singh et al. (2014) simulated a rainfall-runoff episode in Goodwin Creek watershed (Mississippi, USA). For the infiltration simulation the Green-Ampt method was used. Finally, the authors of this paper have proposed an in-house 2D fully hydrodynamic model, called FLOW-R2D, which is based on the solution of the Shallow Water Equations using a modified form of McCormack numerical scheme incorporating artificial viscosity (Tsakiris and Bellos, 2014). The model was tested recently in complex built-up areas (Bellos and Tsakiris, 2015a). A short description of this model is presented in the next section of this paper. Most of the recent applications of hydrodynamic models showed an increased accuracy in 2D flood simulations. However, although the speed of calculations has been considerably increased due to the technological innovations in computing facilities, the computational time for a fully hydrodynamic simulation (particularly in large watersheds) is still high and in many cases still unmanageable. This led researchers to redirect their efforts towards hybrid methods which exploit the merits of the hydrodynamic modelling but at the same time use also some simpler hydrological approaches aiming at reducing the computational time. A recent example in this direction is the hydrologicalhydraulic model proposed by Nguyen et al. (2015), called HiRes Flood-UCI. This model combines the hydrological model HLRDHM as a rainfall-runoff model, but instead of the classical hydrological routing scheme, it uses the 2D hydrodynamic BreZo model. The model has been successfully tested in a catchment of the Illinois river in USA using a decameter resolution. In the present study a novel hybrid method which combines both the hydrodynamic and the hydrological approaches is used, in order to simulate flood events in small catchments. The implementation of the method is performed in three steps: (a) the Unit Hydrographs of the catchment are derived using a 2D hydrodynamic model, (b) the effective rainfall depth of the storm episode is calculated using an appropriate infiltration model, (c) the final flood hydrograph of the event is generated based on the principle of hydrograph superposition. The performance of the proposed method is tested in a flash flood event which occurred in a small catchment of the Halandri stream, in a suburb of Athens-Greece, in February 2013. For the derivation of the Unit Hydrographs, the in-house 2D hydrodynamic model FLOW-R2D is used and for the infiltration process, both the Kostiakov equation and the Green-Ampt model, are used. The results from the proposed hybrid method in the above event are compared directly with the available recorded data and the results derived from the fully hydrodynamic simulation of the entire flood episode. 2. Theoretical background For the hydrodynamic part of the hybrid method, the numerical model FLOW-R2D is used, based on the solution of the twodimensional Shallow Water Equations (2D-SWE) by the Finite
Difference Method, in a cell-centred, non-staggered computational grid (Tsakiris and Bellos, 2014):
@W @F @G þ þ ¼D @t @x @y where
h uh 2 W ¼ uh ; F ¼ u2 h þ 0:5gh ; vh uv h rf D ¼ ghS0;x sb;x =q ghS0;y sb;y =q
ð1Þ vh G¼ uv h ; 2 v 2 h þ 0:5gh
ð2Þ
and h is the water depth, u is the horizontal component of flow velocity at x-direction, v is the horizontal component of flow velocity at y-direction, S0,x, S0,y are the bottom slopes for the x and y-directions, respectively, and q is the fluid density. Shear stresses sb,x and sb,y can be modelled by various methods such as the Manning equation which is used in the present study, as follows:
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
sb;x gn2 uð u2 þ v 2 Þ sb;y gn2 v ð u2 þ v 2 Þ ¼ ; ¼ q q R1=3 R1=3
ð3Þ
where n is the Manning coefficient and R is the hydraulic radius. The two terms, r and f in the Continuity Equation, represent the precipitation rate (source term) and the rate of losses, such as infiltration or drainage (sink term), respectively. It should be noted that infiltration and rainfall are not incorporated into the Momentum Equations of the model because their magnitude is very small in comparison with the other components. A modification of the explicit McCormack numerical scheme (McCormack, 1969) is used to solve the 2D-SWE. It is noted that this scheme has been also used by several researchers (i.e. Kalita, 2016). With this modification, artificial viscosity is added through a diffusion factor. The oscillatory errors are smoothed out while the shock capturing capability of the scheme is retained (the solution has still second order accuracy) and therefore discontinuities of the flow can be simulated. Wet/dry modelling is achieved through a water depth threshold which distinguishes the wet and the dry cells. In the dry cells the water depth and flow velocity are taken as zeros. In the following equations, the discretisation of 2D-SWE using the McCormack numerical scheme is obtained in two steps (predictor strep in Eq. (4) and corrector step in Eq. (5)):
W i;j ¼ W ki;j
Dt Dt k F iþ1;j F ki;j Gki;jþ1 Gki;j þ DtDki;j Dx Dy
ð4Þ
1 1 xW ki;j þ ð1 xÞ W kiþ1;j þ W ki1;j þ W ki;jþ1 þ W ki;j1 2 4 Dt Dt ð5Þ F i;j F i1;j ðGi;j Gi;j1 Þ þ DtDi;j þW i;j Dx Dy
¼ W kþ1 i;j
where Dt, Dx, Dy are the time step, the space step in x- and the space step in y-direction, respectively, and x is the diffusion factor. Finally, the non-rectangular computational domains are represented by a pseudo-computational rectangular domain which encloses the real domain. A detailed description of the FLOWR2D can be found in Tsakiris and Bellos (2014). An extensive grid convergence study of the model has been also performed for steady and dynamic evolution state applications and can be found elsewhere (Bellos and Tsakiris, 2016). Further, interesting applications of the model in urban environments can be found in Bellos and Tsakiris (2015a, 2015b). In order to calculate the effective rainfall depth which is transformed into direct runoff, the infiltration losses should be determined. Various empirical equations have been proposed and
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
used in the past for the infiltration simulation, such as the equations proposed by Kostiakov, Horton, Philip and several others physically-based such as the Green-Ampt model (Tayfur et al., 1993; Tayfur and Kavvas, 1998; Borah, 2011). In this study, both the empirical Kostiakov equation and the Green-Ampt model are used for the simulation of the infiltration process. The first one is chosen, due to its simplicity and the small number of parameters for determination, and the second one for its accuracy and broad acceptance. As known, according to Kostiakov, the infiltration rate f(t) can be calculated by the following equation:
f ðtÞ ¼ abt
b1
ð6Þ
where a and b are parameters estimated by experiments and related to the soil type, and t is the time elapsed from the beginning of the infiltration process. According to the Green-Ampt model, the infiltration rate is calculated as follows:
f ðtÞ ¼ K
wDh þ1 FðtÞ
ð7Þ
where K is the hydraulic conductivity of the soil, w is the wetting front soil suction head, Dh is the water content change, and F(t) is the cumulative infiltration determined by the following implicit equation:
FðtÞ FðtÞ ¼ Kt þ wDhln 1 þ wDh
ð8Þ
3. Description of the case-study The FLOW-R2D model is applied in a small catchment (about 5 km2) of the Halandri stream, in the north suburbs of Athens (Greece). Although the greater part of the catchment is natural, urban zones are identified in the downstream section of the catchment. No manmade structures are located in the hydrographic network of the catchment. In Fig. 1a, an aerial imagery of the catchment is shown with the location of the hydrometric station. Also the land use map based on the information of CORINE 2000 is presented in Fig. 1b. For simplicity purposes, the five zones of
333
CORINE 2000 are merged into three friction zones (Fig. 1c): the coniferous forest zone (friction zone A), the discontinuous urban fabric with the mineral extraction sites zone (friction zone B) and the transitional woodland-shrub with the sclerophyllous vegetation zone (friction zone C). The geology of the catchment consists of marble. The maximum, minimum and average surface elevation of the catchment are 1045 m, 325 m and 609 m asl, respectively. The maximum length of river reach of the watershed is 4.75 km. The average slope is about 32% and 12% for the catchment and the main river reach, respectively. Based on the above characteristics, the concentration time of the catchment is calculated 1.19 h by the Giandotti equation and 0.50 h by the Kirpich equation. Finally, various geomorphological factors are calculated for the catchment. The Gravelius factor is Kc = 1.43, the form factor F = 0.27, the catchment circularity Rc = 0.49 and the catchment elongation factor RL = 0.59. As known, the Gravelius factor is an indicator for the calculation of the compactness of the catchment, the form factor depends on the area and the length of the catchment, the circularity depends on the area and the perimeter of the catchment and the elongation factor depends on the area and the length of the catchment. On the 22nd of February 2013, a storm episode occurred in Athens creating significant casualties, especially in the northern suburbs of the city. In the case-study area, the precipitation depth reached 121.4 mm within 5 h, causing a flash flood and hitting the catchment of Halandri stream. The rainfall event was recorded at the meteorological station of Palea Penteli and the flood hydrograph at the hydrometric station of Nea Penteli. These data were recorded in the framework of DEUCALION Project (http:// deucalionproject.gr/). Based on these data, the runoff coefficient was estimated to be equal to 0.49. In Fig. 2, the hyetograph and the flood hydrograph recorded at the outlet of the catchment are presented. The storm episode of the case-study was divided in two parts: the first one is from the onset of the event at 5:30 up to 9:00 and the second from 9:00 to 11:30.
4. The hybrid method As mentioned above, the proposed hybrid method is implemented in three steps: (a) the Unit Hydrographs of the catchment
Fig. 1. (a) Aerial imagery of the Halandri catchment, (b) land use of the case-study according to CORINE 2000, (c) merged land use zones to the three friction zones.
334
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
Fig. 2. The hyetograph and hydrograph of the rainfall-runoff episode of the 22nd of February 2013.
are derived using a 2D hydrodynamic model. In the present case study, the FLOW-R2D model is used for their derivation. (b) The effective rainfall depth of the storm episode is calculated using an appropriate infiltration model, such as the Kostiakov equation and the Green-Ampt model. (c) The final flood hydrograph of the event is generated based on the principle of hydrograph superposition.
4.1. Unit hydrographs derivation In this section, the derivation of the Unit Hydrographs is presented, implementing the 2D FLOW-R2D model and using the source term of the Continuity Equation (Eq. (1)) for the simulation of the precipitation. Two types of Unit Hydrograph were derived for the rainfall durations: 15 min (UH-15) and 30 min (UH-30). As far as the friction is concerned, a roughness coefficient is calibrated for each friction zone. As easily understood, due to the high computational cost, no optimisation method could be implemented apart from the conventional trial and error method. For the calibration phase, various friction scenaria which have the form S1S2S3 are compared. The objective function is the Nash-Sutcliffe efficiency coefficient (Nash and Sutcliffe, 1970) between recorded data and the results of runoff simulation of the storm episode, obtained by the composition of the unit hydrographs:
Pn ðQ Q obs Þ2 NSE ¼ 1 Pi¼1 num 2 n i¼1 ðQ obs Q obs Þ
ð9Þ
where Qnum is the simulated and Qobs the observed flow, respectively. The first component of the friction scenario S1 represents the Manning coefficient for the friction zone A and takes values 1 and 2 for Manning friction coefficient values 0.050 and 0.070 s/m1/3, respectively. The second component of the friction scenario S2 represents the Manning coefficient for the friction zone B and takes values 1 and 2 for n values 0.020 and 0.035 s/m1/3, respectively. Finally, the third component of the friction scenario S3 represents the Manning coefficient for the friction zone C and takes values 1 and 2 for n values 0.025 and 0.050 s/m1/3, respectively. These
values are selected based on the related literature but increased slightly to take into account the sub-grid effects (Chaudhry, 2008). The computational cell size of the Digital Terrain Model (DTM) of the catchment is of 10 10 m. In total, the computational area consists of 50,811 cells. It should be noted that a preliminary investigation is performed in order to choose the appropriate grid size. It has been found that the coarser grids (e.g. 25 25 m) result in unrealistic outcomes, underestimating the main magnitudes of the phenomenon. On the contrary, finer grids (e.g. 5 5 m) increase considerably the computational burden. The time step equal to Dt = 0.004 s for the UH-15 and Dt = 0.005 s for the UH30, are chosen based on the authors experience. In real world applications, the diffusion factor commonly takes values around 0.90 (Bellos and Tsakiris, 2015b). However, the overland flow in the catchment is characterised by low water depths and no intense turbulence phenomena (apart from streams), and therefore the diffusion factor can take relatively higher values. In this work, the diffusion factor is set equal to x = 0.99. The wet/dry threshold is chosen hdry = 104 m, while a small water depth equal to 103 m over the entire catchment is assumed as the initial water depth. This thin film of water is used as initial condition. It can be noticed that using only the water depth threshold of the wet/dry modelling option, the time step should be further reduced to grasp the resulting high velocities which occur in the catchment. For the model simulation, solid boundaries are imposed to the catchment boundaries, in order to preserve the water mass balance. Since the simulation is referred to the effective rainfall, the infiltration process is not simulated except from a constant infiltration rate equal to the initial water depth for achieving 1 cm equivalent runoff depth. This rate is calculated f0 = 0.4 cm/h for the UH-15 derivation and f0 = 0.2 cm/h for the UH-30 derivation, respectively. In order to preserve the mass balance between the effective rainfall water volume and the runoff water volume, the results of the FLOW-R2D model are slightly modified through a correction factor (Table 1). The slight mismatch between these volumes appears due to the fact that the simulations are decided to be terminated when the contribution of the Unit Hydrographs spans up to a specific time limit, which was set to 3 h. This limit is imposed to save computational cost. However, the runoff process continues, but with very small rates. In Fig. 3, the two types of the Unit Hydrographs
335
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339 Table 1 Performance of the hybrid method against the observed data for the selected scenaria. Friction scenario
Unit hydrograph duration
Correction factor
Infiltration losses
NSE
MAE
RMSE
TPE (%)
PRE (%)
1.1.1 1.1.2 1.2.1 1.2.2 2.1.1 2.1.2 2.2.1 2.2.2
UH-15 UH-15 UH-15 UH-15 UH-15 UH-15 UH-15 UH-15
1.10 1.17 1.18 1.31 1.11 1.22 1.22 1.41
Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov
0.804 0.603 0.689 0.381 0.805 0.567 0.638 0.298
2.930 3.373 4.335 4.431 3.987 4.288 5.833 5.416
5.112 5.537 7.273 7.209 6.444 6.618 9.087 8.782
5.263 5.263 15.789 10.526 10.526 10.526 15.789 15.789
24.686 26.118 38.861 39.952 22.742 24.273 41.108 41.971
1.1.1 1.1.2 1.2.1 1.2.2 2.1.1 2.1.2 2.2.1 2.2.2
UH-15 UH-15 UH-15 UH-15 UH-15 UH-15 UH-15 UH-15
1.07 1.13 1.13 1.23 1.07 1.13 1.14 1.25
Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt
0.770 0.610 0.672 0.422 0.765 0.574 0.625 0.345
3.088 3.601 4.580 4.706 4.368 4.584 6.266 5.743
5.099 5.594 7.599 7.534 6.946 7.071 9.679 9.346
5.263 5.263 15.789 15.789 10.526 10.526 21.053 21.053
22.184 23.700 35.606 37.005 21.117 22.615 38.163 39.476
1.1.1 1.1.2 1.2.1 1.2.2 2.1.1 2.1.2 2.2.1 2.2.2
UH-30 UH-30 UH-30 UH-30 UH-30 UH-30 UH-30 UH-30
1.10 1.17 1.18 1.31 1.11 1.22 1.22 1.41
Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov Kostiakov
0.778 0.645 0.693 0.433 0.743 0.600 0.664 0.377
3.146 3.509 4.166 4.442 4.072 4.239 5.527 5.259
5.445 5.870 6.882 6.889 6.404 6.543 8.694 8.419
10.526 10.526 15.789 5.263 10.526 10.526 10.526 10.526
22.435 24.901 36.504 37.363 26.775 28.741 41.283 41.905
1.1.1 1.1.2 1.2.1 1.2.2 2.1.1 2.1.2 2.2.1 2.2.2
UH-30 UH-30 UH-30 UH-30 UH-30 UH-30 UH-30 UH-30
1.07 1.13 1.13 1.23 1.07 1.13 1.14 1.25
Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt Green-Ampt
0.742 0.644 0.679 0.469 0.709 0.606 0.657 0.416
3.468 3.784 4.451 4.625 4.219 4.371 5.840 5.542
5.850 6.234 7.301 7.250 6.699 6.760 9.113 8.825
10.526 10.526 15.789 15.789 10.526 10.526 21.053 21.053
18.904 21.497 36.809 38.362 32.951 34.623 39.024 40.630
are shown, for all the friction scenaria derived from the FLOW-R2D model.
4.2. Effective rainfall determination Following the derivation of Unit Hydrographs, the effective rainfall is determined. As mentioned earlier, the Kostiakov equation and the Green-Ampt model were used for infiltration simulation. Since the catchment consists of marble, infiltration is expected to be high. Hence, the required parameters for each method are selected assuming that the infiltration process for marble is similar to that of sandy loam. Specifically, the Kostiakov equation parameters are estimated: a = 3.584 cm/h0.584 and b = 0.584 (Antonopoulos, 1999) and the Green-Ampt model
parameters are estimated: K = 1.09 cm/h, w = 2.67 cm and Dh = 0.283 (Tsakiris, 2013). In Fig. 4, the infiltration rate as calculated from the Kostiakov equation and the Green-Ampt model, are plotted parallelly to the rainfall intensity. It is observed that these methods converge, except of some discrepancies at the initial time period. 4.3. Flood hydrograph simulation The final step of the proposed method is the composition of the flood hydrograph based on the principle of hydrological superposition. In Figs. 5 and 6, the simulated hydrographs for each type of Unit Hydrograph, friction scenario and method used to simulate the infiltration process are shown, compared with the observed hydrograph. In Table 1, a comparison through the Nash efficiency
Fig. 3. The 15 min and the 30 min unit hydrographs derived from the FLOW-R2D for the selected friction scenaria.
336
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
Fig. 4. The rate of infiltration losses according to the Kostiakov equation and the Green-Ampt model plotted parallelly to the rainfall intensity.
Fig. 5. Observed data vs storm episode simulation results composing 15 min and 30 min unit hydrographs with the Kostiakov equation for simulating infiltration, for several scenaria.
Fig. 6. Observed data vs storm episode simulation results composing 15 min and 30 min unit hydrographs with Green-Ampt infiltration model, for several scenaria.
coefficient for all the scenaria, is presented. According to the NSE results, the best performance is achieved for the flood hydrograph created by the UH-15, the friction scenario 2.1.1, while the infiltration losses are determined by the Kostiakov equation. Apart from the Nash coefficient, additional indicators are calculated, for the method performance estimation against the observed data: the Mean Absolute Error (MAE), the Root Mean Square Error (RMSE), the time to peak error (TPE) and the peak rate error (PRE). The values of these indicators are also presented in Table 1.
An effort is also made to treat the parameters of the infiltration process simulation methods as black box parameters which should be calibrated in each case. Although the results of this procedure fit better to the observed data, the parameters take unrealistic values according to the related literature. For example using the UH-15 and the 2.1.1 friction scenario, the calibrated parameters for the Kostiakov equation are a = 211.970 cm/h0.183 and b = 0.183. Also for the Green-Ampt model the calibrated parameters take the following values: K = 0.003 cm/h, w = 2933.78 cm and Dh = 1.0. The
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
337
Fig. 7. Comparison between observed data and simulation results obtained by the hybrid method assuming the infiltration simulation methods models as black box models.
Fig. 8. Flood in the hydrographic network produced by FLOW-R2D model 30 min from the beginning of the rainfall event (UH-15, scenario 2.1.1).
calibration is made through the maximisation of the NSE values and the only physical constraint which is adopted is that Dh should be equal or less than one. The NSE values are optimised as 0.818 and 0.812 for the two hydrographs, respectively. The results obtained from this procedure against observed data are presented in Fig. 7. Finally, for illustration purposes, the flood in the hydrographic network of the catchment produced by the FLOW-R2D model is shown in Fig. 8. This snapshot refers to 30 min from the beginning of rainfall event, for the UH-15 and the 2.1.1 scenario. 5. The fully hydrodynamic method In this section, the entire storm episode is simulated using the FLOW-R2D model applied at a catchment scale. Details of the fully hydrodynamic model can be found elsewhere (Tsakiris and Bellos, 2014; Bellos and Tsakiris, 2015a). Due to the high computational cost, a single simulation run is performed using the output
obtained from the hybrid method. Specifically, the Kostiakov equation is used for the infiltration losses and the friction scenario 2.1.1 for the Manning coefficient values is adopted. All the input parameters of the model (diffusion factor and wet/dry threshold) remain unchanged in relation to the hybrid method. The time step is selected Dt = 0.004 s. Fig. 9 presents the comparison between the flood hydrographs obtained by the observed data, the fully hydrodynamic method and the best results from the proposed method. The performance of the fully hydrodynamic method towards measured data by the Nash performance coefficient is NSE = 0.860. As expected, the fully hydrodynamic method performs slightly better than the hybrid method but with much higher computational cost. In the present case-study, the computational time for the fully hydrodynamic method is three to four times greater compared to the hybrid method. Specifically, the hybrid method requires about 3–4 days CPU time for an i7 CPU/3.6 GHz/16 GB RAM, whereas the fully hydrodynamic method requires 10–12 days.
338
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
Fig. 9. Comparison between observed data and simulation results obtained by the fully hydrodynamic method and the hybrid method.
6. Discussion Although the rainfall depth is of the same magnitude at the two parts of the storm episode (rainfall peak of the first part is about 70% of the corresponding rainfall peak of the second part), a great difference can be noticed in the runoff records. The runoff peak which occurs during the first part is just about 30% of the corresponding runoff peak of the second part. This is explained by the fact that during the first part, the soil becomes fully saturated and in the second part, nearly all the rainfall depth is transformed into runoff. It can be seen that the Unit Hydrographs of 15 min rainfall duration give better results from the Unit Hydrographs of 30 min, in which significant oscillations in the simulated flood hydrograph are observed, most probably due to the time step. Besides, there are not significant differences in the results by the use of either the Kostiakov equation or of the Green-Ampt model for the infiltration process simulation, except during the first few time moments of the storm episode. It can be concluded that, unless there are available reliable data, the Kostiakov equation can be the appropriate choice due to its simplicity. Based on the Nash-Sutcliffe performance coefficient, the UH-15 with the 2.1.1 friction scenario performs better than all others. The above combination gives also better results investigating the capability of the proposed method to simulate the peak of the flood hydrograph. Specifically, the peak of the flood hydrograph simulated with the UH-15, the 2.1.1 friction scenario and the Kostiakov equation is 96% of the observed peak. The timing of the peak flow for both of the above cases is delayed by 20 min from the peak based on the observed data. As far as the calibration phase is concerned, friction seems to be also a critical factor for the flood hydrograph simulation. Indicatively, the simulated peak flow of the entire episode for the worst friction scenario 2.2.2 (according to the NSE values), is about 71% of the corresponding best friction scenario 2.1.1 and has 35 min time lag, for the UH-15 and for both methods of infiltration simulation. For the UH-30 and the worst friction scenario 2.2.2 is about 72% of the corresponding better friction scenario 1.1.1 and has 34 min time lag, for the UH-15 and for both methods of infiltration simulation. The most influential zone towards the final results is zone C, followed by zone B. Friction zone A seems to have no significant influence on the final results. This was rather expected, due to the fact that zone A represents a small percentage of the catchment area. An alternative approach for the calibration phase is to assume the parameters required for the infiltration simulation as black box parameters. This approach performs better if the simulated
flood hydrograph is compared to the observed data. On the other hand, the parameters take unrealistic values according to the related literature. Therefore, this approach is not recommended, unless the catchment case study is gauged and at least some flood events are recorded, in order to implement the calibration and the validation phase. A further aspect which should be examined in detail is the convergence of the proposed method as the computational grid becomes finer. It is to note that the great computational cost did not permit an extensive convergence analysis in the context of the specific case study. This analysis should probably be performed in a smaller catchment with milder land surface slopes, in which the computational domain will be represented by a smaller number of cells and greater time step. From the case studied there are strong indications that the proposed hybrid method could be successfully used in devising appropriate design tools for the flood hydrograph simulation of small catchments in various storm episodes, even if their duration is relatively long. Besides, in the context of this method, the calibration phase by the trial and error method is still applicable. Needless to say that the fully hydrodynamic method produces relatively better results than the hybrid method according to the NSE values. The peak of the flood hydrograph obtained from this procedure is 111% of the observed peak, while the time of the peak is just 4 min delayed from the respective observed peak. 7. Concluding remarks A novel hybrid method for the 2D flood simulation was presented in this study. The method combines the 2D hydrodynamic modelling and the unit hydrograph theory. The method was applied to a small catchment in a suburb of Athens for a single storm event. The results from the method were compared to the recorded data at a hydrometric station at the outlet of the catchment and the results derived by the application of the fully hydrodynamic simulation model, FLOW-R2D. For the case studied, the proposed hybrid method showed a satisfactory performance according to the Nash-Sutcliffe efficiency coefficient. Also, the obtained results are close to those produced by the fully hydrodynamic simulation, although the required computational time was considerably reduced. This finding is promising provided that it is further verified in other catchments and conditions. Then, the proposed hybrid method may be useful for devising effective hybrid tools for relatively accurate 2D flood simulations which can be attained at a reasonable computational cost.
V. Bellos, G. Tsakiris / Journal of Hydrology 540 (2016) 331–339
Acknowledgments This paper is based partially on a previous short presentation made at the 9th World Congress of the European Water Resources Association, Istanbul, 24–27 June 2015. Acknowledgments are expressed to D. Koutsoyiannis and A. Efstratiadis for providing useful data on the case-study, to the National Cadastre & Mapping Agency of Greece for the provision of the DTM maps, to I. Papageorgaki for the handling of the DTM, and finally to the Computer Centre of NTUA for the access provided to the computer facilities. References Abderrezzak, K., Paquier, A., Mignot, E., 2008. Modelling flash flood propagation in urban areas using a two-dimensional numerical model. Nat. Hazards 50, 433– 461. Antonopoulos, V., 1999. Hydrology of Unsaturated Soil Zone. AUTH editions (in Greek). Bellos, V., 2012. Ways for flood hazard mapping in urbanised environments: a short literature review. Water Utility J. 4, 25–31. Bellos, V., Tsakiris, G., 2015a. Comparing various methods of building representation for 2D flood modelling in built-up areas. Water Resour. Manage. 29, 379–397. Bellos, V., Tsakiris, G., 2015b. 2D flood modelling: the case of Tous dam break. In: Proceedings of 36th World Congress of IAHR, the Hague, the Netherlands (e-proceedings). Bellos, V., Tsakiris, G., 2016. Grid coarsening and uncertainty in 2D hydrodynamic modelling. 4th European Congress of IAHR, Liege, Belgium (accepted for oral presentation). Borah, K.D., 2011. Hydrologic procedures of storm event watershed models: a comprehensive review and comparison. Hydrol. Process. 25 (22), 3472–3489. Chaudhry, M.H., 2008. Open Channel Flow, second ed. Springer. Costabile, P., Costanzo, C., Macchione, F., 2012a. Comparative analysis of overland flow using finite volume schemes. J. Hydroinformatics 14 (1), 122–135. Costabile, P., Costanzo, C., Macchione, F., Mercogliano, P., 2012b. Two-dimensional model for overland flow simulations: a case-study. Eur. Water 38, 13–23. Costabile, P., Costanzo, C., Macchione, F., 2013. A storm event watershed model for surface runoff based on 2D fully dynamic wave equations. Hydrol. Process. 27 (4), 554–569. Esteves, M., Faucher, X., Galle, S., Vauclin, M., 2000. Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed results. J. Hydrol. 228 (3–4), 265–282. Feng, K., Molz, G.J., 1997. A 2-D diffusion-based wetland flow model. J. Hydrol. 196 (1–4), 230–250.
339
Fiedler, F.R., Ramirez, J.A., 2000. A numerical method for simulating discontinuous shallow flow over an infiltrating surface. Int. J. Numer. Methods Flows 32, 219– 240. Gottardi, G., Venutelli, M., 2008. An accurate time integration method for simplified overland flow models. Adv. Water Resour. 31 (1), 173–180. Kalita, H.M., 2016. A new Total Variation Diminishing predictor corrector approach for two-dimensional shallow water flow. Water Resour. Manage. 30 (4), 1481– 1497. Kazezyilmaz-Alhan, C., Medina, M.A., 2007. Kinematic and diffusion waves: analytical and numerical solutions to overland and channel flow. J. Hydraul. Eng. 133 (2), 217–228. Kim, J., Warnock, A., Ivanov, V.Y., Katopodes, N.D., 2012. Coupled modeling of hydrologic and hydrodynamic processes including overland and channel flow. Adv. Water Resour. 37, 104–126. Liang, D., Falconer, R.A., Lin, B., 2007. Coupling surface and subsurface flows in a depth averaged flood wave model. J. Hydrol. 337, 147–158. Liang, J., 2010. Evaluation of Runoff Response to Moving Rainstorms. Dissertations (2009); Paper 73, Marquette University, Wisconsin, USA.