SEDIMENT MOVEMENT MODEL DEVELOPMENT Yves Zech, Professor1 and Sandra Soares-Frazão, PhD1,2 Benoit Spinewine1, Nicolas le Grelle1, Aronne Armanini3, Luigi Fraccarrollo3, Michele Larcher3, Rocco Fabrizi3, Matteo Giuliani3, André Paquier4, Kamal El Kadi4, Rui M. Ferreira5, João G.A.B. Leal5, António H. Cardoso5 and António B. Almeida5 Summary The present paper aims to present the issues and the scope of the IMPACT research project in the field of dam-break induced geomorphic flows, to give an overview of the experimental work carried out in the frame of the research program, to summarize the new developments in modeling, to outline the validation process, and to give some practical conclusions for the future of dam-break wave modeling. Introduction In a number of ancient and recent catastrophes, floods from dam or dike failures have induced severe soil movements in various forms: debris flows, mud flows, floating debris, and sediment-laden currents (Costa and Schuster, 1988). Other natural hazards also induce such phenomena: glacial-lake outburst floods and landslides resulting in an impulse wave in the dam reservoir or in the formation of natural dams subject to major failure risk. Figure 1. Entrained material from dam or dike failures (Capart, 2000) 1 Catholic University of Louvain (Université catholique de Louvain - UCL), Place du Levant 1, B-1348 Louvain-la-Neuve, Belgium, zech@gce.ucl.ac.be 2 National Fund for Scientific Research, Rue d’Egmont 5, B-1000 Bruxelles, Belgium, soares@gce.ucl.ac.be 3 University of Trento (Università degli Studi di Trento - UdT), Italy 4 CEMAGREF, France 5 Instituto Superior Técnico (IST), Lisbon, Portugal Paper - 6 98 Paper - 6 Fig. 1 presents some estimates of the volume of sediment material moved by such flows, gathered from published case studies (Capart, 2000; Capart et al. , 2001). In some cases, the volume of entrained material can reach the same order of magnitude (up to millions of cubic meters) as the initial volume of water released from the failed dam. Even when they involve comparatively small volumes of material, geomorphic interactions can lead to severe consequences because of localized changes or adverse secondary effects. In India, for instance, the Chandora river dam-break flow of 1991 stripped a 2 m thick layer of soil from the reaches immediately downstream of the dam (Kale et al., 1994). In the 1980 Pollalie Creek event, Oregon, the material entrained by a debris flow deposited in a downstream reach, forming a temporary dam that ultimately failed and caused severe flooding (Gallino and Pierson, 1985). Another cascade of events was that of the 1996 Biescas flood, Spain, where a series of flood-control dams failed (Benito et al., 1998). The problem with dam-break induced geomorphic flows is that they combine the difficulties of two types of flow: (1) alluvial flows, where the bed geometry evolves under the flow action, but with a sediment load small enough to play no dynamic role and (2) rapid transients involving such rapid changes and intense rates of transport that the granular component plays an active role in the flow dynamics, and that inertia exchanges between the bed and the flow become important (Capart, 2000). Research Issues and Scope The main goal of the “Sediment movement” IMPACT work package is, building upon the previously gathered information, to gain a more complete understanding of geomorphic flows and their consequences on the dam-break wave (Zech and Spinewine, 2002). Dam-break induced geomorphic flows generate intense erosion and solid transport, resulting in dramatic and rapid evolution of the valley geometry. In counterpart, this change in geometry strongly affects the wave behavior and thus the arrival time and the maximum water level, which are the main characteristics to evaluate for risk assessment and alert organization. Depending on the distance to the broken dam and on the time elapsed since the dam break, two types of behavior may be described and have to be understood and modeled. Near-field behavior In the near field, rapid and intense erosion accompanies the development of the dam- break wave. The flow exhibits strong free surface features: wave breaking occurs at the center (near the location of the dam), and a nearly vertical wall of water and debris overruns the sediment bed at the wave forefront (Capart, 2000), resulting in an intense transient debris flow (Fig. 2). However, at the front of the dam-break wave, the debris flow is surprisingly not so different as a uniform one. A first section is thus devoted to the characterization of the debris flow in uniform conditions. Figure 2. Near-field geomorphic flow (UCL) 99 Paper - 6 Behind the debris-flow front, the behavior seems completely different: inertial effects and bulking of the sediments may play a significant role. Surprisingly, such a difficult feature appears to be suitably modeled by a two-layer model based on the shallow-water assumptions and methods. The second section relates experiments, modeling and validation of this near- field behavior. Far-field behavior In the far field, the solid transport remains intense, but the dynamic role of the sediments decreases. On the other hand, dramatic geomorphic changes occur in the valley due to sediment de-bulking, bank erosion, and debris deposition (Fig. 3). The third section is devoted to experiments, modeling, and validation of the far-field behavior. Figure 3. Dam-break consequences in the far field Lake Ha!Ha! 1996 dam break (Brooks and Lawrence, 1999) Debris flow in uniform conditions Iverson (1997) reports some interesting information about various debris-flow events in USA, Peru, Colombia, and New Zealand. The main characteristics of this type of event are the involved volume, the run-out distance (sometimes hundredths of kilometers), the descent height (till 6000 m in the quoted examples), and the origin of the debris flow (mainly landslides and volcanic events). A debris-flow also occurs at the front of a dam-break wave, if the latter happens on mobile bed and/or banks. In this case, a high amount of sediments is generally mobilized, inducing a vertical velocity component able to form a kind of plug at the front of the wave. Experimental works (University of Trento UdT) To investigate the vertical structure of free-surface liquid-granular flows, it is of particular interest to be able to materialize steady uniform flow conditions. A re-circulating flume was designed and constructed for this purpose at the Università degli Studi di Trento, Italy. It consists in a tilting glass-walled channel linked with a conveyer belt, forming a closed loop for the circulation of both water and sediment (Fig. 4). 100 Paper - 6 From these experiments, it is possible to gather some information about the acting forces involved in such debris flow (Armanini et al., 2000). Also, the main characteristics of the debris flow may be measured, such as the distribution of the velocities and particle concentration in the normal upward direction. Both can be measured by Voronoï imaging methods, using the grains themselves as tracers (Capart et al., 1999). The concentration along the wall is deduced from 3-dimensional Voronoï cells built by use of stereoscopic imaging (Spinewine et al., 2003). The velocity of each particle may be decomposed into the sum of a mean velocity and of a random component, taking into account the relative motion of the particle compared to the mean value. It is thus possible to define a granular temperature Ts as the mean square value of the instantaneous deviation from the mean velocity (Ogawa, 1978). In analogy with thermodynamic temperature, the granular temperature plays similar roles in generating pressures and governing the internal transport rates of mass, momentum, and energy. Figure 4. Trento re-circulating flume – Photograph and plane view Modeling developments (UdT) Some physical similarities between rapid granular flows and gases has led to a great deal of work on adapting kinetic theories to granular materials, utilizing the idea of deriving a set of continuum equations (typically mass, momentum, and energy conservation) entirely from microscopic models of individual particle interactions. All of the models are based on the assumption that particles interact by instantaneous collisions, implying that only binary or two- particle collisions need to be considered. Particles are usually modeled in a simple way, ignoring surface friction. Furthermore, molecular chaos is generally assumed, implying that the random velocities of the particles are distributed independently. Jenkins and Hanes (1998) applied kinetic theories to a sheet flow in which the particles are supported by their collisional interactions rather than by the velocity fluctuations of the tur 101 Paper - 6 bulent fluid. The purpose of their analysis is the prediction of mean fluid velocity, particle concentration, and granular temperature profiles obtained as solutions of the balance equations of fluid and particle momentum and particle fluctuation energy. The flow of the mixture of particles and fluid is assumed to be, on average, steady and fully developed. The grains are taken to be identical spherical particles of diameter D composed of a material of mass density .s. The fluid is assumed to have a mass density .w . The constitutive relation for the particle pressure is taken to be the quasi-elastic approximation for a dense molecular gas proposed by Chapman and Cowling (1970): ss = Cs .s ... 1+ 4Cs 2( 21- - CCs )3 ... Ts (1) . s . where Cs is the grain concentration, the fraction (2-Cs)/2(1-Cs)3 is the radial distribution function at contact, describing the variation with concentration of the rate of collisions among the particles. In the same way, the constitutive relation for the particle shear stress is taken: 8 22 - Cs 21 . p. 52(1- Cs )3 .2 . du t= 21D. CT s ss s 5p 2(1- Cs )3 .. .1+ 12 ... 1+ 8 Cs (2 - Cs ).. . .. . dz (2) From experiments, it is possible to derive ss and ts by assuming that the buoyant weight of the grains is entirely supported by collisional granular contacts. In Fig. 5 comparison is made between so-derived experimental results and the theoretical relations in Eq. 1 and Eq. 2 (blue lines). A better fit is obtained by accounting for an added- mass effect by replacing .s in Eq. 1 and Eq. 2 by: . 1+ 2C .. sw .s '=.s .. 1+ 2(1- C ). .. (3) . ss . resulting in the red line in Fig. 5. More details can be found in Armanini et al. (2003). Figure 5: Particle pressure and shear stress: points represent experimental results, blue line theoretical kinetic relation, red line accounting for added mass effect 102 Paper - 6 Near-Field Geomorphic Flow Experimental approaches (Catholic University of Louvain UCL) Debris flow is only a part – in time and space – of a dam-break induced geomorphic flow. Other aspects due to the severe transient character of the flow are considered by means of idealized dam-break experiments. Typically, a horizontal bed composed of cohesionless sediments saturated with water extends on both sides of an idealized "dam" (Fig. 6). Figure 6. Scheme of a flat-bed dam-break experiment Upstream lies a motionless layer of pure water, having infinite extent and constant depth h0 above the sediment bed. An intense flow of water and eroded sediments is then released by the instantaneous dam collapse (Fig. 7). Figure 7. Idealized dam-break experiment (UCL) after 0.25 s (a), 0.50 s (b) and 1.00 s (c) In the experiments carried out at the Université catholique de Louvain whithin the frame of the IMPACT program, two materials have been used for representing the sediments: PVC pellets and sand, with rather uniform grain-size distribution. Two arrangements were tested: the flat-bed case with the same sediment level on both sides of the dam (see Fig. 6), and the stepped case where the upstream bed level is higher than the downstream bed level. Some of those experiment were proposed as benchmarks to the IMPACT partners for comparison with their numerical models. The measurement techniques were various: gauges, interface imaging by simple cameras, and particle tracking using tracers or the sediments themselves. 103 Paper - 6 Numerical modeling development (UCL) The near-field modeling generally relies on numerical methods, since analytical solutions (Fraccarollo and Capart, 2002), whilst clever, cannot take into consideration real-case geometry. Fig. 8 illustrates a simplified but fruitful approach to the problem (Spinewine, 2003; Spinewine and Zech, 2002a). Three zones are defined: the upper layer (hw) is clear water while the lower layers are composed of a mixture of water and sediments. In the original model (Ca- part, 2000), the concentration of sediment was assumed to be constant (Cs = Cb) and the upper part of this mixture (hs) was assumed to be in movement with the same uniform velocity as the clear-water layer (us = uw). According to those assumptions the shear stress was supposed as continuous along a vertical line. One of the main improvements brought to the model is to give new degrees of freedom to concentrations (Cs . Cb) and velocities (us . uw) between the three layers. In the frame of shallow-water approach, it is now possible to express the continuity of both the sediments and the mixture and also the momentum conservation with the additional assumption that the pressure is hydrostaticaly distributed in the moving layers, which implies that no vertical movement is taken into consideration: Figure 8. Assumptions for mathematical description of near-field flow .hw +.(hwuw )= -eb Cb - Cs (4a) .t .xC s .hs +.(hsus )= eb Cb (4b) .t .xC s .zb =-eb (4c) .t .(hwuw ) + . ...hwuw 2 + ghw 2 .. . + ghw .(zb + hs )= - t w - eb Cb - Cs ... uw if eb > 0 .. . (5a) .t .x . 2 ..x .w Cs . us if eb < 0 . .(hsus ) + . .. . hu2 + ghs 2 .. . + gh .. . .zb +.w .hw ... = t w -tb +.w Cb - Cs e . ... uw if eb > 0.. . (5b) .t .x . ss 2 . s ..x .s .x ..s .s .sCs b . us if eb < 0 . where eb is the erosion rate (negative is the case of deposition), resulting from the inequality between the shear stresses ts and tb on both faces of the bed interface: 104 Paper - 6 1 (t-tb ) (6) eb = s u .b s The shear stresses tw and ts are evaluated from the turbulent friction, while tb is related to the grain pressure by the soil cohesion and friction. The set of Eqs. 4-5 is solved by a second-order Godunov finite-volume scheme, where the fluxes are computed using the HLLC Riemann solver (Toro, 1997). Validation of the models (UCL, IST, UdT, Cemagref) Validation of the models for the near-field behavior was achieved through benchmarking (Spinewine and Zech, 2002b). The test consisted in the situation sketched in Fig. 6 with the following characteristic dimensions: a water layer of depth h0 = 0.10 m in the reservoir, and a fully saturated bed of thickness hs = 0.05 m. The bed material consisted of cylindrical PVC pellets with an equivalent diameter of 3.5 mm and a density of 1.54, deposited with a bulk concentration of about 60%. Fig. 9 presents a comparison between experimental observation at UCL and the model presented above. The first picture (Fig. 9a: time t = 0.2 s) clearly evidences the limitation of the model for the earlier stage of the dam-break: some features linked to the vertical movements are missed, like the splash effect on water and sediment. The erosion depth is slightly underestimated, partly due to a kind of piping effect under the rising gate, which is not included in the model. All those phenomena induce energy dissipation that is not accounted for in the model, what explains that the modeled front has some advance compared to the actual one. Looking at the second picture (Fig. 9b: time t = 0.6 s), it appears that some characters of the movement are really well modeled, such as the jump at the water surface, the scouring at the dam location, the moving layer thickness. The modeled front is yet ahead but this advance is the same as at the former time, which implies that the front celerity is correctly estimated. (a) (b) Figure 9. Comparison between experiments and numerical results (UCL) at times (a) t = 0.2 s and (b) t = 0.6 s The same test was run concurrently by the Impact teams to compare the characteristics of the various models. The model of the Technical University of Lisbon (IST) relies also on a three-layer idealization. Localized erosion / deposition processes are represented by vertical fluxes but not their impact on the thickness of the transport layer. The model features total (water and sediment) mass and momentum conservation laws, averaged over the flow depth .(zb +h)+.()hu =0 (7) .t .x 105 Paper - 6 . (.muh)+ .(.wuw 2 hw +.sus 2hs )+ 1 g . (.whw 2hw + 2.whwhs +.shs 2hs )= -g(.whw +.shs ). zb -tb (8) .t .x 2 .x .x and mass conservation equations of the transport layer and of the bed, respectively: .(Ch )+ .(C hu )= F (9) .t ss .x sss s (1-e0 ).zb = -Fs (10) .t where h = hs + hw, u = (ushs+ uwhw)/h represent the average velocity of the moving layers (whose thickness are hs and hw, respectively), tb is the bed shear stress, .m is the mean density of the layers such that .m h u = .w hw uw+ .s hs us, .s = .w(1+(s-1)Cs) is the transport layer density, e0 is the porosity, and Fs is the flux between the bed and the transport layer. In the IST model the dependent variables are h, u, zb and Cs. Closure equations are required for: hs, derived from the equation of conservation of granular kinetic energy; us, averaged from a power-law distribution; tb, quadratic dependence on the shear rate; and Fs, depending on the imbalance between capacity and actual transport. Further details can be found in Ferreira et al. (2003) and Leal et al. (2003). The model used by the University of Trento (Fraccarollo et al., 2003) considers constant concentration of sediment (Cs = Cb), and the upper part of this mixture (hs) is assumed to be in movement with the same uniform velocity as the clear-water layer (us = uw = u) in such a way that Eq. 4a-c may be combined in the following way: . . t (zb + hs + hw )+ . . x [(hs + hw )u]= 0 (11a) .(zb + hs )+ .[(hs + hw )u]= 0 (11b) .t .x and Eq. 5a-b are merged in the following form, where h = hs + hw and r = (s–1)Cs with s = .s/.w, the latter being the density supplement due to the presence of the sediment load. . [(h + rhs )u]+ .[(h + rhs )u2 + 12 gh2 + 12 r g hs 2 ]+ g (h + rhs ).zb =- tb (12) .t .x .x . w The Cemagref model RubarBE (El Kadi and Paquier, 2003) relies on the classical Saint- Venant equations extended to the whole cross section: .A +.Q = 0 (13) .t .x . . Qt + . . x ... . QA 2 ... . + gA . . zxw =-gSf (14) where A and Q are the section area and the discharge, zw is the water level, and Sf the friction slope. The conservation of bed material is expressed by the Exner equation, very similar to Eq. 10: (1-e ).As + .Qs = 0 (15) 0 .t .x where As is the bed material area and Qs the solid discharge. Only the water layer is taken into consideration, and the closure of the system is made by the solid discharge. 106 Paper - 6 -0.2 0 0.2 0.4 0.6 EXP CEM UdT UCL IST 0hz 0hgt The comparison of the various models with the experimental data is made in Fig. 10. Regarding the front celerity the results by Trento (UdT) take advantage of the calibration process, which involves these celerity as a calibration parameter. In contrast, their moving sediment layer is underestimated, due to the fact that the concentration of this layer is assumed to be the same as the bed material, which is not the case of the Louvain (UCL) and Lisbon (IST) models: in the reality, the concentration of this moving layer has to decrease to allow the movement of the particles. The erosion due to the front mobilization only appears in the Louvain and Cemagref (CEM) models. Even though Cemargref’s simple model cannot provide any results for the moving sediment layer, it still yields a valuable estimate for the water surface after the shock. The asymmetric treatment of erosion and deposition in Eq. 6 could explain the success for the UCL model in this regard. 0 4 8 12 16 20 Figure 10. Comparison between experimental and numerical results from the benchmark on dam-break wave over an initially flat erodible bed at x = 5 h0. For each set of results, the lower line corresponds to the fixed bed level, the middle line to the moving sediment layer and the upper line to the water surface Channel alteration in the far field The transition between near-field and far-field behavior is not absolutely clear. The debris- flow front resulting from the early stage of dam-break forms a kind of obstacle, which is progressively subject to piping and overtopping. That means that a sediment de-bulking occurs, and the solid transport evolves to a bed- and suspended-load transport with a particularly high concentration. The flow is highly transient and invades a part of the valley that was probably never inundated in the past. All the bank geotechnical equilibrium characteristics are ruined, in such a way that a dramatic channel metamorphosis may be expected. This corresponds to the so-called far-field behavior. A spectacular channel widening generally occurs due to bank scouring and collapse (see Fig. 3). This eroded material over-supplies the bed-load transport resulting in bed deposition and eventual generation of natural dams in the downstream reaches, which may rapidly collapse. 107 Paper - 6 Experimental approaches (Catholic University of Louvain UCL) Laboratory scale models of rivers give interesting information about geomorphic evolution but they are generally not used for sudden transients. Bank failure experiments are commonly carried out to study some fluvial mechanisms such as river meandering or braiding. Also channel-width adjustments during floods may be reproduced in laboratory (see e.g. Chang, 1992), but for cases where this evolution is rather progressive. The experiments carried out within the IMPACT project consist in a dam-break flow in an initially prismatic valley made of erodible material, as sketched in Fig. 11. Such experiments reproduce qualitatively well the features of fast transient geomorphic flows. The upstream part of the channel is fixed, i.e. neither the bed nor the banks can be eroded. The downstream part is made of uniform non cohesive material. A detailed description of the experiment can be found in le Grelle et al. (2004) Fixed bank Gate Sand bank Figure 11. Experimental set-up Figure 12. Bank erosion resulting from intermittent block failure The experience is launched by suddenly raising the gate. This releases a dam-break wave, which rapidly propagates down the channel and triggers a series of bank failures. The rapid erosive flow attacks the toe of the banks with the consequence that they become steeper near the bed and thus fail. Bank erosion then occurs in fact as a series of intermittent block failures (Fig. 12) that feed the flow with an important quantity of sediments. The channel enlargement due to bank failures is the most important in the immediate vicinity of the dam. The water depth there is greater, and the flow shows a two-dimensional expansion from the reservoir into the channel. After a relatively short time (about 10 s in the scale experiment), most of the geomorphic action has occurred. Only light bedload transport can be observed and the banks are no longer affected. Flow measurement is achieved using a laser sheet technique (le Grelle et al. 2004) that allows continuous measurement of the geomorphic evolution of a given cross-section during the flow. The overall principle of the method is to use a laser-light sheet to enlighten a given cross section and to film it during the whole duration of the experiment by means of a remote camera through the transparent side-wall of the channel. The trace of the imprinted laser line onto the digital images is then localized and projected back in 3D space using distinct projective transforms for the immerged and emerged portions. The results were found to be surprisingly reproducible, even though the bank erosion mechanism through intermittent block failures is quite stochastic. 108 Paper - 6 Numerical modeling development (UCL) The key issue in modeling geomorphic processes is to properly include bank failure mechanisms in the system. Indeed, such important geomorphic changes occur randomly and abruptly, and cannot be considered just as a continuous process such as bedload transport. Two different models were developed by UCL within the frame of the IMPACT project. First, a 2D extension of the model presented for the near field (Eq. 4-5) was developed, including a bank erosion mechanism. A detailed description of the method, summarized here, can be found in Spinewine et al. (2002) and Capart and Young (2002). The key idea is that by allowing separate water and fluid-like slurry layers to flow independently, the governing equations are fully equipped to deal with flow slides of bank material slumping into the water stream. Once failure is initiated, the post-failure flow can be captured just like any other pattern of water and sediment motion. A liquefaction criterion is needed to determine when and where portions of the banks are to be transformed from a solid-like to a fluid-like medium. Therefore, the following fundamental mechanism is assumed: activation of a block failure event occurs whenever and wherever the local slope exceeds a critical angle .c. An extended failure surface is then defined as a cone centered on the failure location and sloping outwards at residual angle .r < .c. Finally, sediment material above this cone is assumed to instantaneously liquefy upon failure. In order to account for the observed contrast between submerged and emerged regions, four distinct angles of repose are defined as indicated in Fig. 13: angles .cs and .rs apply to the submerged domain, and .ce and .re to the emerged domain. zz- s0 unstable metastable stable submerged emerged hw0 ( )( ) +( )xx y y - -2 2 1/2 00 Figure 13. Stability diagram for the 2D geostatic failure operator The second model selected for coupling with the above bank erosion mechanism is a one- dimensional scheme. It comprises a hydrodynamic finite-volume algorithm and a separate sediment transport routine (paper in preparation). The finite-volume scheme, developed with the aim of coping with complex topographies (Soares-Frazão and Zech, 2002), solves the hydrodynamic shallow-water equations, under the form of Eq. 13-14. The changes in cross-sectional geometry due to longitudinal sediment transport (bedload) over one computational time step are derived from the Exner continuity equation of the sediment phase (Eq. 15). In addition to sediment fluxes at the upstream and downstream faces of a cell, lateral sediment inflow resulting from bank failures must be considered. A failure is triggered by the submergence of a bank by a rise .h in water level that destabilises a prismatic portion of material as sketched in Fig. 14 that results in a lateral solid discharge qs. The final shape of the cross section shows a submerged slope of angle ae,s (angle of repose under the water level 109 Paper - 6 after erosion) while the emerged part gets the angle ae,e corresponding to the angle of repose of humid sand above the water level after the erosion process. Figure 14. Bank failure triggered by the sub-Figure 15. Deposition of the material eroded mergence of the bank from the banks The eroded material deposits into the channel as sketched in Fig. 15. The submerged portion deposits with an angle ad,s corresponding to the angle of repose under water while the emerged portion stabilizes at an angle ad,s (angle of repose above the water level after the deposition process). All those angles of repose are specific to the material used in the experiments and were measured by means of static and dynamic experiments. Finally, the numerical 1D model consists in solving in a de-coupled way the three different key steps of the process: (i) the hydrodynamic routing of the water, (ii) the longitudinal sediment transport and the resulting erosion and deposition, and (ii) the bank failures and the resulting morphological changes in the cross-section shape. Validation of the models (UCL, UT, Cemagref, IST) Validation of the models will be achieved through benchmarking at two different levels. A first benchmark concerns the idealized dam-break flow experiment presented in a previous section. The blind test was achieved by the involved partners and the comparison process is underway. The second level concerns the simulation of a real event, namely the Lake Ha!Ha! flood that occurred in the Saguenay region of Quebec in 1996 (Brooks and Lawrence, 1999). This second benchmark has just started and the blind modeling by the partners is in progress. Some preliminary comparisons for the first benchmark are presented in Fig. 16. The experimental measurements are compared to the results obtained by the 1D model developed by UCL. The overall agreement is good: the numerical model appears to follow quite accurately the progressive enlargement of the cross section. Conclusions The problem with dam-break induced geomorphic flows is that they combine several difficulties. They involve such rapid changes and intense rates of transport that the granular component plays an active role in the flow dynamics, and that inertia exchanges between the bed and the flow become important. Dam-break induced geomorphic flows generate intense erosion and solid transport, resulting in dramatic and rapid evolution of the valley geometry. In return, this change in geometry strongly affects the wave behavior and thus the arrival time and the maximum water level. 110 (a) (c) (b) (d) Paper - 6 Figure 16. Bank erosion benchmark: comparison between UCL 1D model and experiments Distance downstream from the dam: 0.50 m (a-b), 1.50 m (c-d) : initial situation, experiments, and numerical modeling In the near field, rapid and intense erosion accompanies the development of the dam- break wave, leading to an intense transient debris flow. The numerical models existing at this stage provide encouraging results. The jump at the water surface, the scouring at the dam location, and the moving layer thickness are fairly well represented. But the earlier stage of the dam-break flow is not so well modeled, since the vertical movements depart from the shallow- water assumptions. Finally, all those phenomena dissipate some energy, what is not represented in the models, what explains that the computed front is generally too fast at the beginning. For the far field behavior, the models at this stage can produce valuable results to compare with experimental data from idealized situations, but it is suspected that we are far from a completely integrated model able to accurately simulate a complex real case. A tentative answer to this could probably be given from the results of the second benchmark regarding the Lake Ha!Ha! test case, available after the last IMPACT meeting to be held in Zaragoza, Spain, in November 2004. Acknowledgements The authors wish to acknowledge the financial support offered by the European Commission for the IMPACT project under the fifth framework program (1998-2002), Environment and Sustainable Development thematic program, for which Karen Fabbri was the EC Project Officer. In addition, the authors acknowledge the financial support offered by the Fonds pour la Recherche dans l’Industrie et l’Agriculture, Belgium and by the Fonds National de la Recherche Scientifique, Belgium. 111 References Paper - 6 (1) Armanini A., Fraccarollo L., Guarino L., Martino R., Bin Y . (2000). Experimental analysis of the general features of uniform debris flow over a loose bed. Proceedings Second Int. Conf. on Debris-Flow Hazards Mitigation, Taipei, Taiwan, August 2000, 327-334. (2) Armanini A., Fabrizi R.B., Fraccarollo L., Giuliani M., Larcher M., Rosatti G. (2003) Dynamics of steady and unsteady granular debris flows with uniform material. In EC Contract EVG1-CT-2001-00037 IMPACT Investigation of Extreme Flood Processes and Uncertainty, Proceedings 3rd Project Workshop, Louvain-la-Neuve, Belgium 6-7 November 2003 (CD-ROM). (3) Brooks G.R, Lawrence D.E. (1999) The drainage of Lake Ha!Ha! reservoir and downstream geomorphic impacts along Ha!Ha! River, Saguenay area, Quebec, Canada. Geomorphology 28: 141-168. (4) Benito G., Grodek T., Enzel Y. (1998). The geomorphic and hydrologic impacts of the catastrophic failure of flood-control dams during the 1996-Biescas flood (Central Pyrenees, Spain). Z. Geomorph. 42, 417-437. (5) Capart H., Fraccarollo L., Guarino L., Armanini A., Zech Y. (1999). Recent developments in debris-flow characterisation by digital imaging measurements, CD-ROM proceedings of the 4th CADAM meeting Zaragoza, Spain, 367-380. (6) Capart H. (2000) Dam-break induced geomorphic flows. PhD thesis, Université catholique de Louvain, Louvainla-Neuve, Belgium (7) Capart H., Young D.L., Zech Y. (2001) Dam-break induced debris flow. “Particulate gravity currents”, special publication of the International Association of Sedimentologists (eds B Kneller, B McCaffrey, J Peakall, T Druitt), vol. 31, 2001, pp. 149-156. (8) Capart H., Young D.L. (2002) Two-layer shallow water computations of torrential geomorphic flows. Proceedings of River Flow 2002, Louvain-la-Neuve, Belgium, September 2002. (9) Chang H.H. (1992). Fluvial processes in river engineering. Krieger publishing company, Malabar, Florida. (10) Chapman S., Cowling, T.G. (1970). The Mathematical Theory of Non-Uniform Gases, 3rd Edn. Cambridge University Press. (11) Costa J.E., Schuster R.L. (1988) The formation and failure of natural dams. Bull. Geol. Soc. Am. 100: 1054-1068. (12) El Kadi Abderrezzak K. and Paquier A. (2003). IMPACT project: Dam-break waves over movable beds. In EC Contract EVG1-CT-2001-00037 IMPACT Investigation of Extreme Flood Processes and Uncertainty, Proceedings 3rd Project Workshop, Louvain-la-Neuve, Belgium 6-7 November 2003 (CD-ROM). (13) Ferreira R., Leal J., Cardoso A.H. and Almeida A.B. (2003). Sediment Transport by Dam- Break Flows. A Conceptual Framework Drawn from the Theories for Rapid Granular Flows. In EC Contract EVG1-CT-2001-00037 IMPACT Investigation of Extreme Flood Processes and Uncertainty, Proceedings 3rd Project Workshop, Louvain-la-Neuve, Belgium 6-7 November 2003 (CD-ROM). (14) Fraccarollo L., Capart H. (2002). Riemann wave description of erosional dam-break flows. Journal of Fluid Mechanics, 461, 183-228. (15) Fraccarollo L., Capart H., Zech Y. (2003). A Godunov method for the computation of erosional shallow water transient. Int. Journal of Numerical Methods in Fluids, 41:951-976. 112 Paper - 6 (16) Gallino G.L., Pierson T.C. (1985). Pollalie Creek debris flow and subsequent dam-break flood of 1980, East Fork Hood River Basin, Oregon. U. S. Geological Survey Water Supply Paper 2273, 22 pages. (17) Iverson R.M. (1997). The physics of debris flow. Reviews of Geophysics, 35, 3, 245-296. (18) Jenkins J.T., Hanes D.M. (1998). Collisional sheet flows of sediment driven by a turbulent fluid. Fluid Mechanics, vol. 370: 29-52, Cambridge University Press. (19) Kale V.S., Ely L.L., Enzel Y., Baker V.R. (1994). Geomorphic and hydrologic aspects of mon-soon floods on the Narmada and Tapi Rivers in central India. Geomorphology 10, 157-168. (20) Leal J.G.A.B., Ferreira R.M.L., Cardoso A.H., Almeida A.B. (2003). Overview of IST Group Results on the Sediment Benchmark. In EC Contract EVG1-CT-2001-00037 IMPACT Investigation of Extreme Flood Processes and Uncertainty, Proceedings 3rd Project Workshop, Louvain-la-Neuve, Belgium 6-7 November 2003 (CD-ROM). (21) le Grelle N., Spinewine B., Soares-Frazão S., Zech Y. (2004). Non-intrusive imaging measurements of the morphological evolution of a channel during a dam-break flow. Proceedings of River Flow 2004, Naples, Italy, June 2004 (22) Ogawa S. (1978). Multitemperature Theory of Granular Materials. Proc. U.S-Japan Seminar on Continuum Mechanical and Statistical Approaches in the Mechanics of Granular Materials, S. C. Cowin & M. Satake Eds., Gakujutsu Bunken Fukyukai, Tokyo, 208-217. (23) Soares-Frazão S., Zech Y. (2002) Dam-break in channels with 90° bend, Journal of Hydraulic Engineering, American Society of Civil Engineers (ASCE), Vol. 128, n°11, pp. 956-968 (24) Spinewine B., Capart H., Le Grelle N., Soares-Frazão S., Zech Y (2002). Experiments and computations of bankline retreat due to geomorphic dam-break floods. Proceedings of River Flow 2002, Louvain-la-Neuve, Belgium, September 2002 (25) Spinewine B., Zech Y. (2002a). Geomorphic dam-break floods: near-field and far-field perspectives. In EC Contract EVG1-CT-2001-00037 IMPACT Investigation of Extreme Flood Processes and Uncertainty, Proceedings 1st Project Workshop, Wallingford, UK, 16-17 May 2002 (CD-ROM). (26) Spinewine B., Zech Y. (2002b) Dam-break waves over movable beds: a “flat bed” test case. In EC Contract EVG1-CT-2001-00037 IMPACT Investigation of Extreme Flood Processes and Uncertainty, Proceedings 2nd Project Workshop, Mo-i-Rana, Norway, 1213 September 2002 (CD-ROM). (27) Spinewine B. (2003). Ecoulements transitoires sévères d’un milieu diphasique liquide / granulaire (eau / sédiments), en tenant compte du couplage rhéologique entre les deux phases et des composantes verticales du mouvement (in French). MSc thesis. Université catholique de Louvain, Louvain-la-Neuve, Belgium (28) Spinewine B., Capart H., Larcher M., Zech Y. (2003). Three-dimensional Voronoï imaging methods for the measurement of near-wall particulate flows. Experiments in Fluids, vol. 34, pp. 227-241 (29) Toro E.F. (1997) Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction. Springer-Verlag, Berlin, Germany (30) Zech Y., Spinewine B. (2002). Dam-break induced floods and sediment movement – State of the art and need for research. In EC Contract EVG1-CT-2001-00037 IMPACT Investigation of Extreme Flood Processes and Uncertainty, Proceedings 1st Project Workshop, Wallingford, UK, 16-17 May 2002 (CD-ROM). 113 Paper - 7 PROCESS UNCERTAINTY: ASSESSING AND COMBINING UNCERTAINTY BETWEEN MODELS Mark Morris, HR Wallingford UK Mohamed Hassan, HR Wallingford UK Francisco Alcrudo, UDZ, Spain Yves Zech, UCL, Belgium Karim Lakhal, Engineering School of Science and Technology of Lyon France Abstract The uncertainty associated with breach formation, flood propagation, and sediment movement is important in the risk management of these structures. Assessing risks involves identifying the hazards associated with each issue. This paper identifies the advances that have been made as a result of the IMPACT Project research. Guidance and relative importance of the uncertainty for each process is considered for the purposes of improved risk management of these structures. Scoping the Problem The wider picture Uncertainty is a general concept that reflects our lack of sureness about something, ranging from just short of complete sureness to an almost lack of conviction about an outcome! Two main sources of uncertainty include: 1. Natural variability: referring to the randomness observed in nature. 2. Knowledge uncertainty: referring to our state of knowledge of a system and our ability to measure and model it. Knowledge uncertainty may be further divided according to: 1. Statistical uncertainty: referring to the uncertainty resulting from the need to extrapolate a particular set of data. 2. Process model uncertainty: describing the uncertainty associated with using a process model based on incomplete knowledge of the process, data, or representation of reality. 3. Decision uncertainty: describing the strength of belief in the decision made and its robustness. This decision is likely to be linked to one or all of the above categories of uncertainty. Work on uncertainty1 within the IMPACT2 project focuses on knowledge uncertainty – and specifically process model uncertainty. 1 It should be noted that the word uncertainty means process model uncertainty for numerical models for the rest of this document. 2 IMPACT Project: Investigation of Extreme Flood Processes and Uncertainty. EC Contract No: EVG1-CT-2001-00037. www.impact-project.net 114 Paper - 7 Scope of work feasible within the IMPACT Project The objective of this work package is to identify and emphasise the uncertainty associated with the various components of the flood prediction process; namely breach formation, flood routing and sediment transport. The effect that uncertainty in each of these predictions has on the overall flood prediction process will be demonstrated through application to a real or virtual case study. Thus, the focus of work under IMPACT was to: a) Investigate uncertainty within modelling predictions for predicting breach formation, flood propagation and sediment transport b) Demonstrate how uncertainty within each of these modelling approaches may contribute towards overall uncertainty within the prediction of specific conditions (such as flood water level at a specific location) c) Consider the implications of uncertainty in specific flood conditions (such as water level, time of flood arrival etc.) for end users of the information (such as emergency planners). The scope of work under IMPACT does not allow for an investigation of uncertainty in the impact of flooding or in the assessment and management of flood risk. The assessment of modelling uncertainty provides essential information upon which a later assessment of the uncertainty in risk may be developed through further research. The IMPACT Approach The challenge of assessing overall modelling uncertainty is complicated by the need to assess uncertainty within two or more models, to somehow transfer a measure of uncertainty between these models and to develop a system that allows for the different complexities of the various models. Two basic approaches were adopted, namely sensitivity analysis and Monte Carlo analysis. However, whilst a breach formation model may be able to run hundreds or thousands of simulations within a period of hours, it is unrealistic to assume that a complex 2D flood propagation model can undertake a similar process without undertaking weeks or months of analysis. A compromise solution was adopted for IMPACT that combines sensitivity analysis, Monte Carlo simulation and expert judgement. Whilst this approach may provide an estimate of uncertainty which contains a degree of subjectivity (expert judgement) it also provides a mechanism that is achieved relatively simply and provides a quick indication of potential uncertainty. Sensitivity Analyses The uncertainty of various modelling parameters is examined here by first selecting some representative values for each parameter (e.g. the upper, most likely, and lower values). The modeller then runs the model using these values for each parameter. The output from the model for values other than the most likely value can then be compared with the output of the most likely value or a range is assigned to this model output in relation with the range of the input parameter. This comparison gives an indication of the uncertainty in output derived from each of the input parameters. For example, consider the situation whereby the broad crested weir equation is used in a numerical model and an assessment of the uncertainty contributed by the weir coefficient 115 Paper - 7 (Cd) is required. The following values may be considered as the upper, most likely, and lower values for Cd, 1.90,1.70, and 1.50. Say, for example, three runs have been undertaken and the values of the peak outflow for these runs were 187, 149 and 113 respectively. Based upon these results, the sensitivity results can be represented either as: “ A +12 % change in Cd produced a +25 % increase in the peak outflow and a –12 % in change Cd produced a –24 % decrease in the peak outflow” or as: “The uncertainty of calculating the peak outflow ranges from +25 % to –24 % if Cd values change by ± 12 %.” The advantages offered by this approach are that it is simple, quick and the relative importance of parameters can be identified. The disadvantages of this approach are that only a small number of values for each input parameter may be tested. The selection of the representative values of a parameter is, to some extent, a subjective process. Where uncertainty within a single model arising from a number of parameters is required, the following equation may be used to give one parameter measure of uncertainty (Runc): where: R is the parameter of interest, and x, y, and z are parameters upon which R depends. .R/.x reflects the relative importance of each of the input parameter x on the parameter of interest (R) (same for y and z). xunc reflects the uncertainty range of the input parameter x (same for y and z). This approach can be used for quick assessment of the uncertainty of parameters or when the use of the Monte Carlo analysis approach (described below) is difficult to apply due to, for example, excessively long model run time. Monte Carlo Analysis In this approach, an appropriate probability distribution is selected for each input parameter examined in the model (an example is given in Figure 1 for Cd). A number of runs are then undertaken by changing the values of all of the input parameters based upon their probability distribution. The values of the output parameter are then ranked and the distribution of results is plotted. Confidence limits may be assigned to this distribution (usually 5% and 95 % limits are selected). The range between these limits is then a quantified range for the uncertainty of the output parameter. Under this approach, the output inherently combines the uncertainty of the full range of input parameters. Regardless of whether one or n parameters are considered, no further analysis of output is required to find the overall range of uncertainty. This is advantageous, if the model can be run repeatedly within a reasonable time frame. 116 Paper - 7 0 1 2 3 4 5 6 1.4 1.5 1.6 1.7 1.8 1.9 2 Cd Likelihood Figure 1: Example triangular probability distribution of Cd Figure 2 shows an example of this approach where 1000 runs were undertaken, leading to the distribution shown. Taking the confidence limits as 5% and 95%, the range of uncertainty would be 55-180 m3/s with a likely peak outflow of 120 m3/s (at 50%). This translates to – 65 m3/s to + 60 m3/s uncertainty in the peak outflow generated from all of the selected input parameters. 0 20 40 60 80 10 0 12 0 14 0 16 0 18 0 20 0 020406080100120140160180200220240260280300More P eak o u tflo w ran ks Frequency Figure 2: Peak outflow distribution based on the Monte Carlo analysis approach The advantages of this approach are that a wider range of data is tested giving a better indication of uncertainty in comparison to a simple sensitivity analysis. In addition, a probability distribution is produced for the output parameters (e.g. Peak outflow). The main disadvantages of this approach are that the relative importance of each input parameter is not identified and a greater number of model runs is required in comparison to the simple sensitivity analysis approach (i.e. long run time). 117 Paper - 7 Balancing the Approach With the pros and cons of each approach in mind, the approach adopted by IMPACT was to: 1 Assess breach model uncertainty via sensitivity analysis and Monte Carlo simulation 2 Extract representative flood hydrographs from the breach model analyses representing upper, mid and lower scenarios for use in flood propagation 3 Assess flood propagation models through sensitivity analysis only 4 Either select flood propagation model parameters to match upper, mid and lower scenarios for running with upper, mid and lower scenario breach hydrographs – ending with three sets of model predictions or Select upper, mid and lower scenario parameters for application to each of the 3 breach hydrographs, resulting in 9 sets of model predictions, from which representative upper, mid and lower conditions may be extracted (see Figure 3) Figure 3: Linking uncertainty analysis between models Example Analysis of Uncertainty within a Breach Model The HR BREACH model has been used to develop and demonstrate the approach for uncertainty analysis. Data from the Norwegian Field Test #2 (homogeneous non-cohesive field), which was undertaken in Norway in 2002, was selected for analysis. This field test (Figure 4) was built mainly from non-cohesive materials (D50 ˜ 5 mm) with less than 5 % fines. The purpose of this test was to better understand breach formation and to identify the different failure mechanisms in homogeneous non-cohesive embankments failed by overtopping. More information on the field test programme may be found in Vaskinn (2004). Hassan (2002) gives an overview of the HR BREACH model whilst Morris(2005) gives an overall report on the IMPACT project. 118 Paper - 7 2 m 1.7 15m Clay Rock Figure 4: Geometry of Field Test #2 Sensitivity Analysis: To perform the sensitivity analysis, a range of input parameters for the HR BREACH model were selected. Their relative importance (impact) on the modelling results was then assessed, allowing prioritisation of parameters and hence selection of those most influencing the results. This in turn allows a reduction in the number of parameters involved in the Monte Carlo approach, hence allowing the analysis to be undertaken more quickly. The parameters initially selected can be divided into the following categories: 1. Hydraulic parameters: • The weir equation coefficient (Cd): a coefficient combining the effect of energy losses and approach velocity when calculating the flow through a weir. • Manning’s friction coefficient: a coefficient that represents the boundary roughness. 2. Sediment parameters: • Sediment median diameter (D50): representative size of sediment. • Sediment transport equation: an equation that is used to compute sediment transport rates. 3. Soil parameters • Soil density: a measure of how much mass is contained in a given unit volume. • Angle of friction: a parameter represents the internal friction between soil particles and similar to the angle of repose of the soil. • Cohesion: the force that holds together the molecules in a soil. 4. Model specific parameters: • Sediment flow factor: a factor to ‘adjust’ the sediment transport rates calculated by a sediment transport equation. • Breach width to depth ratio: a ratio used to distribute the sediment transport volume and hence update the breach shape. The following output parameters of the model were also selected to measure the relative importance of the above input parameters: 1. Peak outflow 2. Time to peak outflow 3. Final breach width 4. Final breach depth 119 Paper - 7 Table 1 below shows the base values used in the model to simulate the failure of the test case. These values were considered the ‘base’ for comparison as they represent the best estimate value for each parameter and they are either measured in the field, estimated, or default values of the model. Table 1 also shows the range values that were selected for each parameter to undertake the sensitivity analysis. This range of values was either based upon the variation of the measured data (e.g. D50) or judgement of what is a reasonable range for each parameter (i.e. Cd). 5 variations of each parameter were considered enough to present the valid range for this test case except for the sediment transport equation where only the three available equations in the model were used varying the sediment flow factor once. Therefore, 44 model runs were undertaken, in addition to the base run, varying only one parameter for each run and recording its effect on the selected output parameters. Table 1: Base value and variations of each parameter InputParameters Base Value Variations 1 2 3 4 5 CD 1.50 1.55 1.60 1.70 1.80 1.90 Mannings’n 0.020 0.025 0.030 0.035 0.040 0.045 D50 (mm) 4.65 2.00 3.00 4.00 5.00 6.00 Density (KN/m3) 21.15 19.00 20.00 21.00 22.00 23.00 Angle of Friction (o) 42 25 30 35 40 45 Sediment flow 1.0 0.5 1.5 2.0 2.5 3.0 factor Breach width to depth ratio 0.50 1.00 1.50 1.75 2.00 2.50 Sediment transport equation Yang Visser (1 flow factor) Chen (Sand) (1 flow factor) Visser (2 flow factor) Chen (Sand) (2 flow factor) Cohesion 0.9 0.0 2.0 2.0 7.0 10.0 (KN/m2) As an example of the output of the sensitivity analysis, Table 2 shows the minimum, maximum, and mean values of the peak outflow (obtained by varying the input parameters) and also the outcome of the base run for those parameters. It also shows the minimum, maximum, and the range of variation from the mean and the base values. Based upon the sensitivity analysis results, It was found that: 1-The sediment transport equation and the sediment flow factor are consistently showing the highest impact on all output parameters. 2-The angle of friction, however to a less extent than the two parameters above, is constantly showing a higher impact than other physical parameters. 3-Other physical parameters such as D50, Manning‘s n, and Cd and model parameters such as the breach width to depth ratio have a lesser effect than the above input parameters. 4-Cohesion has a low effect on the model output parameters. This is expected for a non- cohesive embankment where the value of cohesion is very low. 120 Paper - 7 5- The final breach depth is the least affected parameter by the variation of the input parameters, followed by, the time to peak with the peak outflow and the final breach width showing the most affected parameters by this variation. Table 2: Sensitivity of the peak outflow to the different input parameters Based upon the above and given that the sediment transport equation is not a suitable parameter for a Monte Carlo approach, the following input parameters were selected to undertake the Monte Carlo uncertainty analysis: 1-The sediment flow factor3 2-The angle of friction 3-The sediment median diameter (D50). This parameter was selected over the other parameters as it has an uncertainty range defined by the measured data and to test the effect of another physical parameter rather than a model parameter. The Monte Carlo Analysis: In order to perform a Monte Carlo analysis a probability distribution has to be assigned to each of the selected input parameters. For this exercise, a triangular probability distribution was chosen and assigned to each of the input parameters (see Figure 5). The basis for selecting the minimum, most likely, and maximum value for each distribution are either measured data or judgement of what is a reasonable range for each parameter for this test case. It was also assumed that there is no correlation between these parameters. Sediment Flow Factor Probability Distribution D50 Probability Distribution Angle of Friction Probability Distribution 0 0.1 0.2 0.3 0.4 0.5 0.6 2 4 6 D50 Value Likelihood 0 0.02 0.04 0.06 0.08 0.1 0.12 25 30 35 40 45 Angle of Friction Value Likelihood 0 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 2.5 3 Sediment Flow Factor Value Likelihood Figure 5: Selected probability distribution of input parameters 3 The sediment flow factor relates to the distribution of sediment movement around the wetted perimeter of the breach formation area (Mohamed (2002)). 121 Paper - 7 The question now arises as to how many runs should be undertaken to ensure that the results truly reflect the full range of possibilities. It was considered that 30-50 runs are needed for each parameter. This means that for 3 parameters a total number of (30)3 to (50)3 runs are required, which = 27000 – 125000 runs! This number of runs is clearly impractical in terms of time and effort. Reducing the number of runs from range from 30-50 per parameter to say 10 per parameter (tactical Monte Carlo) = 1000 runs, which is achievable. The adequacy of this number of runs can be checked by considering how the probability distribution of the output parameters is converging towards a consistent distribution with the number of runs. Figure 6 shows the probability distribution of the model peak outflows after undertaking the Monte Carlo analysis. Peak outflow 5% 50% 95% New0.00 10.00 20.00 30.00 40.00 50.00 60.00 Frequency 0-6570-7580-8590-95100-105110-115120-125130-135140-145150-155160-165170-175180-185190-195200-205210-215220-225230-235240-245250-255260-265270-275280-285290-295 Peak outflow ranks (m^3/s) Figure 6: Probability distribution of peak outflow The following conclusions were drawn from the results of the Monte Carlo analysis for this test case: 1-The probability distribution of each output parameter has converged towards a consistent distribution for all parameters. 2-The final breach depth is the least affected parameter by the variation of the model-input parameters. 3-The range of variation of the model output parameters is very similar to that obtained for the same parameters using the sensitivity analysis. 4-The probability distribution of all the monitored output parameters except the final breach depth converged to approximately a triangular distribution. 122 Defining the Uncertainty Range of the Model Further analysis has been done on the probability distribution of the peak outflow to define the uncertainty range of the HR BREACH model for these parameters and test condition. As shown in Figure 6, values that correspond to the 95%, 50%, and 5% bands in the probability distribution were used to represent the upper, medium, and lower bands for the model respectively. These bands are: 1- Upper band : 220-230 m3/s 2- Medium band : 150-160 m3/s 3- Lower band : 90 –100 m3/s A hydrograph within each of the above bands was then selected (Figure 3) to present the upper, medium, and lower limits of the model for the outflow hydrograph. The figure also shows that the base run hydrograph falls between the medium and lower bands in terms of the peak outflow value as well as the measured peak outflow value. 0 50 100 150 200 250 0 1000 2000 3000 4000 5000 6000 Time (s) Flow (m^3/s) Lower Meduim Upper Base Run Measured Peak Outflow Figure 3: Upper, medium, and lower hydrographs for the HR BREACH model in comparison with base run hydrograph and measured peak outflow value. From Figure 7 it can therefore be concluded that: - the base run peak discharge falls between the mid and lower uncertainty estimates - the range of uncertainty around peak discharge is approximately 90 – 210 m3/s – best and mid level estimate is ~ 125-155m3/s. It is also interesting to note that the base run estimate is closer to the actual field value than the mid range uncertainty estimate, falling within approximately 10% of observed value. Paper - 7 123 Paper - 7 Analysis of Uncertainty within Flood Propagation and Sediment Transport Models As highlighted earlier, the analysis of uncertainty within flood propagation and sediment transport models is not so simple. Within IMPACT, the analysis of uncertainty for the Tous Dambreak case study will be undertaken for both breach and flood propagation models to demonstrate application of the whole approach. However, this work is scheduled for completion during summer 2004 and hence cannot be reported here. Instead, the following two sections provide a short discussion of the issues relating to the prediction of uncertainty relating to flood propagation and sediment movement. Analysis of Uncertainty within Flood Propagation Models Uncertainty in flood propagation model results is a key issue because many important decisions depend upon the output of a flood propagation model: Risk maps, land use, evacuation time... depend heavily on the model flood forecast. Within flood propagation models the following sources of uncertainty can be listed: -Lack of knowledge as regards proper mathematical description of the flood. In fact it is certain that the mathematical models used (flat pond models, kinematic, dynamic, liberalised or non-linear shallow water models etc…) are not a good enough mathematical description of a flood. The extent to which one particular model deviates from reality is unknown a priori since some of the conditions that make it fail often depend upon the event itself. Unfortunately not much can be done in this respect. -Lack of knowledge concerning the description of the initial (for instance flood hydrograph) or boundary conditions. Among the latter the bathymetric description of the flooded area can be included, and is usually of paramount importance. However, with current topographic techniques (e.g. LIDAR) this is just a question of economics. -Uncertainties regarding physical parameters affecting the flood. The most typical is the bottom roughness and its distribution. Since flood propagation models combine a range of modelling parameters and formulae it is not a simple process to analyse and understand the effect that one (or more) input parameters could have on results under a variety of different conditions. The only way to obtain this relationship is by running the model with different values of every input parameter and to observe the effect on the outputs (i.e. Monte Carlo type analysis). However, for flood propagation models the problem is aggravated by the fact that state of the art propagation models usually require considerable computing power and can take hours or days to complete one simulation. These constraints typically prevent a Monte Carlo approach and encourage simpler sensitivity analyses to be undertaken instead. In the course of the IMPACT project, several attempts to determine the uncertainty of flood propagation models have been carried out depending upon the phase of the project. During model validation against laboratory data, many of the parameters that usually bear a certain amount of uncertainty are relatively well known. This is, for instance, the case with the Manning’s n or the topography that is described with high accuracy. Also, most of the characteristics of the flood are well controlled during the experiment, and usually subject to repeatability checks: Inflow hydrographs, timing, depth measurements etc… are known with great accuracy. Therefore the uncertainty analysis can concentrate on model parameters: Mesh influence, model numerical parameters and strategies have been tested. 124 Paper - 7 When coping with a real case situation some input data are not well known. This is the case when attempting to model a case study for which the inflow hydrograph is not certain and the bathymetry prior to the event is not exactly known. The strategy adopted to ascertain the influence of the inflow hydrograph on the flood characteristics and effects downstream entails making several runs of the flood propagation model with input hydrographs coming from the uncertainty analysis performed during breach formation work. A representative high, medium and low hydrograph as estimated from the breach models will be fed to the different flood propagation models. Water elevation at different locations will be monitored as a representative output variable to assess the influence of the inflow hydrograph. Despite the simplicity of the procedure, it is expected that it will provide a measure of the uncertainty associated with a representative input parameter. Analysis of Uncertainty within Sediment Movement Models Issues Floods from dam or dike failures induce severe soil movements in various forms. Other natural hazards also induce such phenomena: glacial-lake outburst floods and landslides resulting in an impulse wave in the dam reservoir or in the formation of natural dams subject to major failure risk. In some cases, the volume of entrained material can reach the same order of magnitude (up to millions of cubic meters) as the initial volume of water released from the failed dam. The risks associated to the sediment movement is thus substantial. Dam-break induced geomorphic flows generate intense erosion and solid transport, resulting in dramatic and rapid evolution of the valley geometry. In counterpart, this change in geometry strongly affects the wave behaviour, and thus the arrival time and the maximum water level, which are the main characteristics to evaluate for risk assessment and alert organisation. That means that the uncertainties affecting the sediment movement prediction may critically affect the whole prediction process. Uncertainties in the sediment movement modelling In the near field, rapid and intense erosion accompanies the development of the dam- break wave. Wave breaking occurs at the centre (near the location of the dam), and a nearly vertical wall of water and debris overruns the sediment bed at the wave forefront, resulting in an intense transient debris flow. Behind the debris-flow front, the behaviour is different: inertial effects and bulking of the sediments may play a significant role. Most of the processes involved in this kind of phenomenon are uncertain. The models are based on an idealisation of the dam break. The problem is represented in a vertical plane and the dam is supposed to instantaneously disappear without lateral effects. Only the valley- bed material is taken into account in the near-field solid transport, neglecting the material issued of the breaching itself. At the current stage of the models, the bed mobilisation modelling is not yet coupled with the breaching modelling. The models are promising for idealised situations but are still far to represent the real-life situations. 125 Paper - 7 Also the data needed for such a modelling are commonly difficult to get. The material constituting the reservoir bottom is not uniform, its thickness is not well known, this material seriously evolves with the time, above all in case of fine material. The material of the valley bed downstream from the dam is also heterogeneous: it consists of soils and rocks in an unpredictable arrangement. Measurement of this is tedious, difficult and expensive. For the far field, the point is to represent the valley evolution, with a succession of erosion and deposition supplied by the upstream solid transport and by the bank collapses. A part of the morphologic evolution may be modelled, above all locally, but for a reach of a few kilometres, there are so many stochastic phenomena involved that the cascade of events becomes unpredictable, forming a kind of uncertainty tree that is difficult to manage. Overall Progress and Conclusions A methodology has been developed combining sensitivity analysis, Monte Carlo analysis and expert judgement to allow assessment of modelling uncertainty and integration of uncertainty between different models. The approach adopted does not adopt a rigorous analysis of uncertainty, but includes subjective components as a trade off against developing a method that is practicable given modelling and user time constraints. The methodology has been applied to the HR BREACH model, using real field data to demonstrate a potential range in breach flood hydrograph results. Similar modelling for the flood propagation of these hydrographs has not yet been completed, but will be undertaken during summer 2004. The procedure will be demonstrated using the Tous Dam failure as a case study. The analysis will be extended from the dam downstream to the prediction of flood water levels at selected locations in a town some kilometres below the dam. Based upon the range of uncertainty predicted, implications may then be drawn for the use of such data by end users such as emergency planners. Analysis of uncertainty within the field of sediment movement is far less developed. The issues discussed above all suggest that our ability to model such phenomena is not yet mature enough to allow analysis in detail regarding the uncertainties. It is possible to carry out some sensitivity analyses in order to identify the most relevant among the uncertainty sources. However, the order of magnitude of these uncertainties is still too large to be of immediate value and application for end users. It may be concluded therefore that further detailed and fundamental research on this important topic is still needed. Future Direction of Work The analysis of uncertainty within dambreak or extreme flood predictions is clearly important if we are to adopt a risk based approach to flood risk management – including activities such as emergency planning, land use planning, defence optimisation etc. The scope of work within the IMPACT project was limited to an assessment of the issues and application of relatively simple techniques for assessing potential uncertainty. Clearly, there are multiple aspects to the assessment of uncertainty within the flood risk management concept, that require detailed and fundamental research work. One project where this issue will be considered in some considerable detail is the EC FLOODsite project (see www.floodsite.net ) 126 Paper - 7 which will run from March 2004 until February 2009. It is also understood that work currently being reported by the Dam Safety Interest Group (DSIG) should also significantly advance understanding in this area. Acknowledgements IMPACT is a research project supported by the European Commission under the Fifth Framework Programme and contributing to the implementation of the Generic Activity on “Natural and Technological Hazards” within the Energy, Environment & Sustainable Development programme. EC Contract: EVG1-CT-2001-00037. The financial support offered by DEFRA and the Environment Agency in the UK is also acknowledged. The IMPACT project team comprises Universität Der Bundeswehr München (Germany), Université Catholique de Louvain (Belgium), CEMAGREF (France), Università di Trento (Italy), Universidad de Zaragoza (Spain), Enel.Hydro (Italy), Sweco (formerly Statkraft Grøner AS) (Norway), Instituto Superior Technico (Portugal), Geo Group (Czech Republic), H-EURAqua (Hungary) and HR Wallingford Ltd (UK). References Mohamed, M.A.A., SAMUELS, P.G., MORRIS, M.W. and GHATAORA, G.S. (2002). Improving the accuracy of prediction of breach formation through embankment dams and flood embankments - Proc. of the Int. Conf. On Fluvial Hydraulics, Louvain-la-Neuve, Belgium, 3-6 Sept. 2002 - River Flow 2002, Bousmar & Zech (editors). Morris M W (2005). IMPACT Project Final Report. European Commission. FP5 Research Programme. Contract No. EVG1-CT-2001-00037. (Draft under development 2004). Vaskinn KA, Løvell A, Höeg K, Morris M, Hanson G and Hassan MAA (2004). Physical modelling of breach formation: Large scale field tests. Association of State Dam Safety Officials: Dam Safety 2004. Phoenix, Arizona Sept 2004. 127 Paper - 8 DETERMINATION OF MATERIAL RATE PARAMETERS FOR HEADCUT MIGRATION OF COMPACTED EARTHEN MATERIALS By 1Greg Hanson, and 2Kevin Cook Abstract. The timing and formation process of a dam embankment breach due to flood overtopping can dramatically impact the rate that water is released from a reservoir. This rate of water release directly impacts the hazard to life and property downstream of a breached dam. Therefore, dam embankment erosion and breaching from overtopping events is important to both engineers and planners alike, who must predict impacts on local communities and surrounding areas affected by flooding. Based on observations from overtopping tests the erosion process has been described as a four-stage process. A key erosion feature has been observed to be headcut formation and migration. Therefore, determination of the material parameter for predicting rate of headcut migration is important to modeling embankment erosion. An equation for predicting the material parameter based on results from a flume study is compared to results from embankment overtopping tests. Flume tests were conducted on 2 soil materials and embankment-overtopping tests were conducted on 3 soil materials. The flume tests and overtopping tests were compacted using similar compaction efforts. It was concluded that the headcut migration parameter was primarily dependent on compaction water content. A 4% change in compaction water content caused an order of magnitude change in the headcut migration parameter. 1 Greg Hanson, Research Hydraulic Engineer, USDA-ARS, 1301 N. Western St. Stillwater, OK 74075, Phone (405) 6244135 ext. 224, Fax (405) 624-4136, e-mail: greg.hanson@ars.usda.gov. 2 Kevin Cook, Civil Engineer, USDA-NRCS, State Office, Stillwater, OK 74075, Phone (405) 742-1257 e-mail: Kevin.Cook@ok.usda.gov. 128 INTRODUCTION Paper - 8 Interest in the occurrence and effects of overtopping of earth embankments by storm runoff has existed for years. Based on conclusions made by Ralston (1987) there are about 57,000 dams on the national dam inventory that have the potential for overtopping. Reservoirs overtop as a result of inflow exceeding the capacity of the reservoir storage and spillway outflow system, and since this risk can never be completely eliminated, the challenge is determining how these embankments will perform in advance of overtopping. One of the key factors in predicting embankment performance is determining the influence of soil materials on the processes and rate of erosion during overtopping and breach. Headcut migration and widening have been observed to be important erosion processes during embankment overtopping and breach (Ralston 1987, Al Qaser 1991, Hahn et al. 2000, Hanson 2003a). Observed Breach Morphology Ralston (1987), in his discussions of dam overtopping, distinguishes between cohesive and non-cohesive soils and their erosion characteristics. Overtopping of embankments with cohesive soils results in eventual degradation of the surface, formation of a discontinuity, and development of an overfall or headcut. The headcut advances progressively headward as the base of the headcut deepens and widens. Failure and breach occur when the headcut migrates through the upstream crest of the dam. The point at which the headcut migrates through the upstream crest has been named, “time of breach initiation,” ti. The point at which erosion reaches the toe of the upstream slope of the embankment has been named, “time of breach formation,” tf. Upstream headcut advance has been attributed to a combination of a) insufficient soil strength to stand vertically due to the height of the headcut face, stress relief cracking and induced hydrostatic pressure in the stress cracks and, b) loss of foundation support for the vertical face due to the waterfall flow plunging effect and its associated lateral and vertical scour. Ralston (1987) recognized that this type of erosion process was a three- dimensional process, in which not only upstream migration occurs but also lateral widening. The rate of widening has been observed to be a function of the headcut migration rate and both are important in determining the timing and amount of water discharge through the breach (Hanson et al. 2003a). Hanson et al. (2003a), based on observations and data recorded during seven overtopping tests, describe the erosion process of cohesive embankments during overtopping as a four stage process involving headcut development, headcut migration, and the three- dimensional aspects including widening: I. Flow initiates at t = t0. Initial overtopping flow results in sheet and rill erosion with one or more master rills developing into a series of cascading overfalls (fig. 1a). The cascading overfalls develop into a large headcut (fig. 1b and 1c). This stage ends with the formation of a large headcut at the downstream crest and the width of erosion approximately equal to the width of flow at the downstream crest at t = t1. II. The headcut migrates from the downstream to the upstream crest of the embankment. The erosion widening occurs due to mass wasting of material from the banks of the gully. This stage ends when the headcut reaches the upstream crest at t = t2 (fig. 1d). III. Lowering of the crest occurs during this stage and ends when downward erosion has virtually stopped at t = t3 (fig. 1e). The peak discharge and primary water surface lowering occurs during this stage for small reservoir. IV. During this stage breach widening occurs (fig. 1f). The peak discharge and primary water surface lowering occur during this stage (t3 < t < t4) rather than during stage III for large 129 Paper - 8 reservoirs. This stage may also be split into two stages, similar to observations by Visser (1998) for sand dike breaching, depending on the upstream head through the breach. Stages I and II (t = t2) encompass the time period up to breach initiation t = ti, and Stage III (t3 – t2) encompasses the time period referred to as breach formation t = tf. These stages as described are a generalization of the processes that were observed. a) Rills and cascade of small overfalls during Stage I at t = 7 min. b) Consolidation of small overfalls during Stage I at t = 13 min c) Headcut at downstream crest, transition from Stage I to Stage II at t = t1 = 16 min. d) Headcut at upstream crest, transition from Stage II to Stage III at t = t2 = 31 min, ti. e) Flow through breach during Stage III at t = 40 min. f) Transition from Stage III to Stage IV at t = t3 = 51 min, tf. Figure 1. Erosion processes during overtopping test (soil 1, embankment 1). (Hanson et al. 2003) 130 Paper - 8 Compaction Effects The headcut migration rate is a function of the soil material properties as well as the hydrodynamic forces and embankment geometry. The embankment materials are typically compacted cohesive soils. It has been observed that the nature and magnitude of compaction have a significant effect on the physical behavior of a soil (White and Gayed 1943, Powledge and Dodge 1985, Robinson and Hanson 1996, Hanson et al. 2003b). White and Gayed (1943) observed from overtopping tests on 0.3 m high embankments constructed in the laboratory that the rate of erosion of the cohesive embankments varied from test to test in such a complicated fashion that the tests could not be correlated numerically. They did observe, however, that the variations could be traced to the clay and water content, to which the erosion rates were very sensitive. Powledge and Dodge (1985) observed that increasing compaction from 95 to 102 percent of standard Proctor compaction resulted in reducing, by half, the erosion of small embankments in flume tests. Robinson and Hanson (1996) conducted large- scale flume studies on headcut migration of cohesive soils. Based on these studies, resistance to headcut migration was reported to increase over several orders of magnitude as compaction water content and compaction energy were increased (Hanson et al. 1998). Hanson et al. (2003b), based on results from seven embankment overtopping tests of 3 different soil materials, observed that since compaction efforts were similar, compaction water content played a major role in setting the rate of erosion of the embankments, including headcut migration and widening. Embankments are constructed of soil material and therefore are affected by these factors in an overtopping event. Therefore it is important to develop algorithms that incorporate hydrodynamic forces as well as soil properties in predicting the erosion processes occurring in an embankment failure during overtopping. Headcut Migration Prediction Predicting the rate of headcut migration has been observed to be one of the keys to predicting cohesive embankment failure during overtopping (Hanson et al. 2003a). Simple relationships for headcut migration prediction have almost universally focused on energy at the overfall as the driving mechanism (De Ploey, 1989, Temple 1992, Temple and Moore, 1997). One of the exceptions to the energy-based approach has been a stress-based approach proposed by Robinson and Hanson (1994). The energy-based formulations typically use some form of unit discharge q and headcut height H to describe hydraulic attack in terms of energy dissipation at the headcut. Temple (1992) proposed a simple model describing headcut migration dX/dt based on a material dependent coefficient C and a hydraulic attack parameter A such that: dX/dt = C(A) A = qaHb (1) (2) where: a and b = exponents. Temple and Moore (1997) used a value of 1/3 for both exponent values of a and b. At present there is no approach for determining the material dependent coefficient C other than based on observed migration rates dX/dt, and q and H. The objective of this paper is to develop a relationship for the material dependent coefficient C based on flume results and compare this relationship to embankment overtopping results. 131 Paper - 8 Figure 2. Schematic of flume set-up (Robinson and Hanson, 1996). TEST SETUP Flume Tests Headcut advance tests were performed in a 1.8-m wide and 29-m long flume with 2.4-m high sidewalls (fig. 2). The test section within the flume was constructed by placing soil in horizontal loose layers 0.15 to 0.20 m thick. A 0.86-m wide vibratory padfoot roller was used to compact each layer, and a hand-held pneumatic compactor was used to compact the soil against the flume walls. Compactive effort and water content were varied from test to test. Prior to testing, a near vertical overfall was preformed at the downstream end of the test section. Overfall heights varied from 0.9 m to 1.5 m. The surface of the fill was protected using carpet strips or a soil cement surface layer to minimize surface erosion and emphasize headcut migration. Following placement of soil in the flume and before testing, samples were taken from the downstream end of the placed soil. The dry unit weight gd and water content wc% were determined as an average of the values determined from undisturbed tube samples. Following soil sampling and headcut forming, the outlet basin was filled with water to the desired level, and flow was delivered to the test flume (fig. 3). Even though advance of the headcut was observed to often be in discrete steps due to mass failures of the soil material at the headcut face, the global rate of movement for a set of flow conditions and soil material properties appeared to be uniform. Therefore, advance rates for each test were determined based on linear regression of the observed headcut position versus time (fig 4). Headcut migration rates of two soils were examined in these flume experiments, Soil E and Soil F (table 1). Standard proctor tests on Soil E exhibited a maximum dry unit weight of 1.90 Mg/m3 at optimum water content of 12% (fig. 5) while Soil F exhibited a maximum dry unit weight of 1.96 Mg/m3 at optimum water content of 10.5%. A total of 46 tests were conducted using Soil E and Soil F, and in 6 test cases of Soil E and one test case of Soil F the compaction effort was similar to the compaction effort used in the embankment overtopping tests. The compaction water contents and dry unit weights for these seven tests are plotted on figure 5. 132 Headcut Migration (m) 6 5 4 3 2 1 0 Ar=0.83 m/h Paper - 8 012345 Time (h) Figure 3. Headcut migration test in large Figure 4. Typical observed headcut outdoor flume. migration during a flume test. Dry unit weight (Mg/m3 ) 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 2.05 2.10 Soil E Soil F Saturation = 100% Standard Laboratory Compaction Soil E Soil F 6 8 101214161820 Water Content % Figure 5. Compaction results for Soil E and Soil F in the laboratory and in the flume. 133 Paper - 8 Embankment Overtopping Tests Three large-scale embankments, two at 2.3 m and one at 1.5 m in height were constructed and tested. Two of the embankments, 2.3 m high (embankment 1) and 1.5 m high (embankment 2), had three test sections, 7.3 m and 4.9 m wide with 2 m and 1.5m trapezoidal notch overflow sections, respectively (fig. 6), and 3H:1V slopes on both the upstream and downstream sides. The three test sections for each embankment had three different soil materials: Soil 1, Soil 2, and Soil 3 (table 1). In order to test each soil material individually, the notches in the other test sections were filled on the upstream end with a soil plug. The height of the embankment at the notch crests was 1.83 m and 1.22 m for the 2.3 m and 1.5 m high embankments, respectively. The third embankment constructed and tested was 2.3 m high (embankment 3) and had a single 12 m test section of Soil 2 with an 8.2 m wide trapezoidal overflow notch. Embankment 3 was constructed and tested in the same location as embankment 1 (fig. 6). The soils were placed in loose lifts 0.15 m thick and compacted with 2 passes with vibration of a vibratory roller, resulting in a compactive effort similar to the seven flume tests previously described. Table 1. Properties of soils used in flume and embankment tests. Soil Parameters Soil E Soil F Soil 1 Soil 2 Soil 3 Gradation % Clay < 0.002 mm25 13 4 6 26 % Silt40 30 28 30 48 % Sand 35 57 68 64 26 Liquid Limit26 16 34 Plastic Limit9 13 17 Plasticity Index 15 3 NP NP 17 USCS CL SM SM SM CL Standard Compaction Values Maximum Dry Unit Weight (Mg/m3)1.90 1.96 1.84 1.86 1.79 Optimum Water Content % 12.0 10.5 9.0 10.5 14.0 1.7-m embankment 2.3-m embankment Reservoir Staff Gage 1.2-m modified Parshall Flume (Inflow) Supply Canal 1.2-m H Flume (Outflow) V – notch Weir (Outflow) Reservoir Drain 2 1&3 Figure 6. Schematic of embankment overtopping facilities. 134 Paper - 8 Chart recorders were utilized to record inflow and outflow hydrographs. An adjustable length (7 to 12 m) overhead rolling carriage with attached point gage was utilized to obtain bed profiles, cross sections, and water surface elevations during testing (fig. 7). Digital cameras were placed at fixed locations for photographic measurement of headcut location and headcut gully width (Hanson et al., 2002). Inflow to the reservoir during testing was supplied by a canal and measured at the test site with a modified Parshall flume for embankments 1 and 2 and a sharp crested weir for embankment 3. Maximum overtopping head attained prior to breach was 0.46 m, 0.30 m, and 0.30 m for the embankments 1, 2, and 3 respectively. The inflow discharge stabilized quickly during each test, and was then maintained at a relatively constant flow of about 1, 0.3, and 2 m3/s for embankments 1, 2, and 3 respectively. This relates to a unit discharge of approximately 0.37, 0.22, and 0.22 m3/s/m for embankments 1, 2, and 3 respectively. Staff gage measurements of the reservoir elevation were used to keep track of the volume of storage in the reservoir at any given time during the test. The outflow hydrograph was determined with a combination of methods. Flows were measured downstream of the reservoir with a H-flume and a V-notch weir. Reservoir elevation and storage records were also used for evaluating the breach outflow. Figure 7. Point gage and carriage. 135 Paper - 8 RESULTS Flume Tests The results from the seven flume tests that were compacted with similar compaction effort to the embankment studies indicated that the compaction water content played a significant role in setting the rate of headcut migration (table 2). The migration rate was observed to change 50 times from a compaction water content of 9.2% to 15.9%. The C Table 2. Relevant flume test measurements. Test Soil WC gd dX/dt qH C # (%) (Mg/m3) (m/s) (m2/s) (m) (s-2/3) 1 E 9.2 1.68 1.5x10-3 0.84 1.2 1.5x10-3 2 F 12.1 1.86 3.5x10-4 0.87 1.2 3.4x10-4 3 E 12.4 1.84 4.1x10-4 0.89 1.2 4.0x10-4 4 E 14.2 1.79 9.3x10-5 0.86 1.3 9.0x10-5 5 E 14.4 1.79 4.2x10-5 0.86 1.3 4.0x10-5 6 E 14.8 1.81 3.5x10-5 1.36 1.0 3.1x10-5 7 E 15.9 1.78 3.0x10-5 0.85 1.3 2.9x10-5 parameter was calculated from equation 1 using the flume test measurements of q, 10-2 H, and dX/dt (table 2). A plot of compaction water content wc% versus C (fig. 8) shows the correlation between these two parameters for these seven flume tests. The following equation can 10-3 be used to estimate C: C = 3000 (wc%)-6.5 (3) where C is in units (s-2/3). Embankment Overtopping Tests10-4 The results from the seven embankment overtopping tests indicate, as did the flume tests, that the compaction water content plays a significant role in influencing the rate of headcut migration (table 3). Stage II of the embankment 10-5 erosion for the seven tests was used to 8 9 10 11 12 13 14 15 16 17 18 determine and compare headcut migration Water Content % rates. Stage II is the phase of erosion in Figure 8. Compaction wc% versus C for flumewhich the headcut that has formed tests at a specific compaction effort. migrates from the downstream crest of the embankment to the upstream crest Table 3. Embankment overtopping results. (Hanson et al. 2003a). This stage Soil WC Emb. dX/dt q H C was chosen due to the fact that the (%) # gd (m/s) (m3/s) (m) (s-2/3)headcut height is typically equivalent (Mg/m3) to the embankment height at the 1 8.7 1 1.72 2.1x10-3 0.38 1.8 2.3x10-3 notch during this stage, discharge is 1 11.5 2 1.73 2.1x10-3 0.21 1.2 3.3x10-3 nearest its constant rate (typically 2 11.5 3 1.77 3.6x10-4 0.22 1.8 4.9x10-4 equal to inflow), and the headcut 2 12.1 1 1.73 1.9x10-4 0.34 1.8 2.2x10-4 migration rate is nearly constant (fig. 2 14.5 2 1.74 6.4x10-5 0.19 1.2 1.0x10-4 9). Values of C predicted from 3 16.4 1 1.65 3.9x10-5 0.39 1.8 4.4x10-5 equation 3 compare well to 3 17.8 2 1.67 1.1x10-5 0.21 1.2 1.7x10-5 C (s-2/3 ) Soil E Soil F C = 3000 (wc%)-6.5 Flume Data 136 Paper - 8 computed values of C from the embankment overtopping test results (fig. 10) with the exception of Soil 1 for test 4. These results indicate that a functional relationship between water content and the material parameter C of equation 1 for a given compactive effort can be used to make excellent predictions independent of material texture. These results also point out the consistency in results for headcut migration rates between 2-D tests in the flume and 3D tests of an embankment overtopping. It is significant to note, that even though the data set is small, the relationship appears to be independent of texture. 10-2 10 Embankment Data 6 8 10-5 10-4 10-3 10-2 Data from Camera Data from Point Gage dX/dt Stages I II II I IV C (s-2/3 ) Headcut Location (m) Soil 1 Soil 2 Soil 3 10-3 4 2 10-4 0 0 100 200 300 400 TIME (minutes) Figure 9. Headcut location versus time. (Hanson et al. 2003b) 10-5 Cpredicted= 3000(wc%)-6.5 (s-2/3) Figure 10. Predicted C versus computed C for embankment overtopping tests. CONCLUSIONS AND SUMMARY One of the key erosion processes that have been identified in embankment overtopping erosion and failure is headcut development and migration. The rate of headcut migration influences the timing of embankment breach and the discharge hydrograph associated with a breach. Development of a computer model to predict breach and the discharge hydrograph requires the identification of appropriate algorithms and definition of the input parameters. In the case of predicting breach this will require determination of a material parameter similar to the C value identified in equation 1. Evaluation of seven headcut migration tests conducted in a flume that were compacted at similar compactive efforts to seven embankment overtopping tests led to the development of equation 3 to predict the C coefficient based strictly on compaction water content independent of soil texture. These results show promise in developing a universal relationship for C to compaction water content and compaction effort. This type of development will be useful in evaluating existing embankments as well as planned embankments and the potential risk of breaching. 137 REFERENCES Paper - 8 Al-Qaser, G. N. 1991. Progressive failure of an overtopped embankment. Unpublished PhD Dissertation. Colorado State University. Fort Collins, CO. De Ploey, J. 1989. A model for headcut retreat in rills and gullies. CATENA Supplement 14, 81-86. Cremlingen, West Germany. Hahn, W., Hanson, G. J., and Cook, K. R. 2000. Breach morphology observations of embankment overtopping tests. Proceedings of the 2000 Joint Conference of Water Resources Eng., Planning, and Management, ASCE. Minneapolis, MN. 10 pp. Hanson, G. J., Cook, K.R., Hahn, W., and Britton, S. L. 2003a. Observed erosion processes during embankment overtopping tests. ASAE Paper No. 032066. ASAE St. Joseph, MI. Hanson, G. J., Cook, K.R., Hahn, W., and Britton, S.L. 2003b. Evaluating erosion widening and headcut migration rates for embankment overtopping tests. ASAE Paper No. 032067, ASAE St. Joseph, MI. Hanson, G. J., Robinson, K. M. and Cook, K. R. 1998. Effects of compaction on embankment resistance to headcut migration. Proceedings of ASDSO, 11-16 Oct. 1998, Las Vegas, NV. pp. 13-20. Powledge, G. R. and Dodge, R. A. 1985. Overtopping of small dams – An alternative for dam safety. Hydraulics and Hydrology in the Small Computer Age. Vol. 2. Proc. Hydr. Div. Specialty Conf., Orlando, FL. ASCE. pp. 1071-1076. Ralston, D. C. 1987. Mechanics of embankment erosion during overflow. Proceedings of the 1987 National Conference on Hydraulic Engineering, Hydraulics Division of ASCE. p.733-738. Robinson, K. M. and Hanson, G. J. 1994. A deterministic headcut advance model. Transactions of the ASAE 37:1437-1443. Robinson, K.M. and Hanson, G. J. 1996. Gully headcut advance. Transactions of the ASAE 39(1):33-38. Temple, D. M. 1992. Estimating flood damage to vegetated deep soil spillways. Applied Engineering in Agriculture 8(2):237-242. Temple, D. M., and Moore, J. S. 1997. Headcut advance prediction for earth spillways. Transactions of the ASAE 40(3):557-562. Visser, P. J. 1998. Breach Growth in Sand-Dikes. Communications on Hydraulic and Geotechnical Engineering. Report no. 98-1. Hydraulic and Geotechnical Engineering Division. Faculty of Civil Engineering and Geosciences. Delft University of Technology. White, C.M. and Gayed, Y. K. 1943. Hydraulic models of breached earthen banks. Institute of Civil Engineers. Vol. 3 of The Civil Engineer in War. p.181-200. 138 Paper - 9 SIMULATION OF CASCADING DAM BREAKS AND GIS-BASED CONSEQUENCE ASSESSMENT: A SWEDISH CASE STUDY Romanas Ascila, M. Sc. (CE) Principal Engineer, SwedPower AB, P.O. Box 527, SE-162 16, Stockholm, Sweden Phone: +46 8 739 66 91; Fax: +46 8 739 53 32; E-mail: romanas.ascila@swedpower.com Claes-Olof Brandesten, Lic. Eng. (CE) Head of Dam Safety Office, Vattenfall AB, SE-162 87, Stockholm, Sweden Phone: +46 8 739 70 95; Fax: +46 8 739 60 30; E-mail: claesolof.brandesten@vattenfall.com Abstract Dam owners are commonly liable for the consequences of failure of any of their dams. In Sweden, like in some other countries, this liability is also strict, meaning that the dam owner is responsible for the consequences whatever the reasons for incidents and failures may be. In some cases the consequences are potentially of such a large scale as to be uninsurable and may also be in excess of the dam owners financial capacity. Swedish legislation for public safety management makes specific reference to risk analysis and risk characterization for major hazards to people, property and the environment. Further, and increasingly, owners of hazardous facilities are required to identify the hazards, assess the risks, prepare a safety case demonstrating how the risks will be prevented or otherwise controlled, and set out a safety management system demonstrating how the safety case will be implemented and maintained. This paper presents how the integrated methodology with state-of-the-art techniques was introduced and practically tested to perform consequence analysis and assessment due to dam-break. This paper addresses different issues including methods for systematic data management, advanced hydraulic and hydrologic modeling, breach formation processes, geographical information systems (GIS) and developments in detailed consequence analysis. The study was performed on the dam cascade including eighth dams, which are situated on the river Lilla Lule Älv in Northern Sweden. 139 Paper - 9 1 Introduction Vattenfall as the main hydropower producer in Sweden and a major owner of high consequence dams must ensure and maintain dam safety at the highest possible level. The liability for the consequences due to dam incidents and possible dam-break requires that these consequences should be evaluated and understood in advance. Other activities that are planned to be implemented, such as strengthening of the downstream face, filter improvements or, particularly, heightening of the dam crest must be analyzed from the dam safety perspective. This study challenged in the building an integrated approach of having a seamless system for dam-break modeling, flood routing and consequence evaluation. 2 Description of the Cascade of Dams The dam cascade is situated on the river Lilla Lule Älv in the northern part of Sweden. This study covers eight dams of which six are hydropower dams and two are stop-dams (Figure 1). Generally, these dams can be identified by the power plant location and grouped into four dam sites: Parki, Randi, Akkats and Letsi. Short description of each dam group is presented below. Figure 1. Map of the Study Area 2.1 Parki The Parki Power Plant is situated between the lakes Parkijaure and Randijaure on the river Lilla Lule Älv. The power plant is supplied with a tubular turbine aggregate. The gross 140 Paper - 9 head is 14 m, with a maximum capacity of about 20 MW. The rock fill dam, which is common to the power plant and the regulation of the Lake Skalka, provides a storage capacity of 460 million m3. The power plants construction started during the autumn of 1967 and the plant was taken into service in 1970. The Parki Power Plant has a main dam with the crest length of 1 900 m, one stop-dam at Stainas with the crest length of 700 m, and three stop-dams at Parkijaure reservoir with the total crest length of about 650 m. 2.2 Randi The Randi Power Plant is situated between the lakes Randijaure and Vaikijaure on the river Lilla Lule Älv. The gross head is 25 m, with a maximum capacity of about 80 MW. The power plants construction started during the summer of 1973 and the plant was taken into service in 1976. The Randi Power Plant has an intake dam with the crest length of 100 m and a regulation dam with the crest length of 160 m. There is an intake tunnel of the length of 480 m. 2.3 Akkats The Akkats power plant is situated on the river Lilla Lule Älv about 4 km from the community of Jokkmokk. The plant is utilizing the gross head between the lake Vaikijaure and the storage of the Letsi power plant. The Gross head is 45.5 m, with a nominal turbine discharge of 385 m3/s and the capacity of about 146 MW. The dam, which is common to the power plant and the regulation of the lake Vaikijaure, will provide a storage capacity of 42 million m3. The main dam is constructed as an earth and rock-fill dam and is founded partly on bedrock, partly on moraine. The total length of the dam is 1515 m. A small stop-dam, which length is 240 m, is constructed to the right of the main dam. Construction was started during the autumn of 1969 and the plant was taken into service during the autumn of 1973. 2.4 Letsi The Letsi power plant was the first power plant in the river Lilla Lule Älv. The plant issituated about 16 km upstream of the confluence with the river Stora Lule Älv. The plant is utilizing the gross head between the lake Valjates and the storage of the Porsi power plant. The gross head is 136 m, with a nominal turbine discharge of 220 m3/s and a capacity of about 268 MW. The dam is constructed as an earth and rock-fill dam and is founded on bedrock. The length of the dam crest is 520 m. Construction was started in the end of 1960 and the plant was taken into operation during the spring of 1976. 3 Integrated Methodology 3.1 General Various initiatives have been taken earlier in integrating GIS and hydraulic applications in Sweden and elsewhere (Ascila, Brandesten 2002). The methodology developed for this project covers procedures for data collection and processing, hydraulic modeling, GIS integration and modeling, consequence analysis and dissemination of the results. 141 Paper - 9 The different components of work have been linked together in what can be denoted as a methodological framework (Figure 2), which describes the logic of the work and the linkages between work packages and models. All components used in the methodological framework compound an integrated scalable system which was implemented, tested and used on the Microsoft Windows® 2000 platform. Figure 2. Methodology Framework The methodology constitutes state-of-the-art, integrated and modular approach for analyzing consequences related to rivers and watercourses. The core of the framework is the GIS shell in which all geographical data is collected, stored, analyzed and presented. Coupled to the GIS are industry proven hydraulic models, for the purpose of solving issues and problems within the watercourse. According to the Figure 2, six Work Packages (WP) are defined: WP1 – Data Collection, WP2 – Site Description, WP3 – Data Processing, WP4 – Hydraulic Modeling, WP5 – Consequence Analysis and WP6 – Information Transfer. The subsequent chapters include brief introductions and explanations of the modules applied within the scheme; also brief descriptions of tools are given. 142 Paper - 9 4 Digital Elevation Modeling 4.1 Digital Elevation Model Over Land Areas A digital elevation model (DEM) of the study area, which extends from the beginning ofLake Parkijaure down to the confluence of the river Lilla Lule Älv and the river Stora Lule Älv, was generated (Figure 1). The aerial photos taken from the height of 9 200 m were the source of DEM generation. Digital photogrammetry techniques were applied for the processing of images and for the generation of final DEM. A DEM of 5x5m resolution, which has a vertical accuracy better than 0.5 m, was generated as a result of this work. Judgments of the data quality and corrections were being made routinely during the generation phase of the DEM. Data quality assessment and result inspections were of a primary interest for the high accuracy calculations such as hydraulic modeling and flood mapping in this study. A major problem existing with the elevation models generated from the stereo pair of aerial images is, that there are no possibilities to get the accurate values over the areas covered by forest, bushes and other features on the terrain. This is not suitable for the flood mapping applications. In order to obtain the correct land elevation values a semiautomatic method of forest removal was developed and applied. Generally the method can be described as identification and removal of the forest areas and interpolation with the help of known ground points. 4.2 Detailed DEM over dam sites Along with the generation of DEM for the whole study area, a set of detailed DEMs was produced. These elevation models were generated from aerial photos of the flight with the height of 1500 m. Aerial images taken from the low flying heights enable generation of the high-resolution elevation models (Figure 3). In this case a resolution of 0.5x0.5 m was achieved. Workflow of generation of these models was analogical as generation of the DEM for the whole area. Figure 3. 3D View of Detailed DEMs over the Parki Main Dam and Randi Regulation Dam 4.3 Bathymetry The methodology of generating elevation model using photogrammetrical technique is unable to produce bathymetry data over the river branches and the reservoirs. This data is 143 Paper - 9 very important in order to be able to develop the accurate hydraulic model for simulation of dam-breaks. In this project the bathymetry data was obtained by scanning and digitizing bathymetry maps (Figure 4). All maps were georectified to the common coordinate system to make possible integration with the elevation model over the land area. Figure 4. Scanned and Georectified Bathymetry Map With Satellite Image in the Background (left) and 3D View of the Digital Bathymetry Data (right) 4.4 Control of Dam Crest Levels with the Help of Detailed DEMs All dams analyzed in this study were built during the period of 1960-1973. Since the commissioning time some deformations may have occurred. In order to check the values known from the technical documentation a control with the help of the detailed elevation models was carried out. The spatial analysis performed with the help of GIS procedures indicates some declinations from the earlier known values of the dam crests (Figure 5). The new updated values were used as the input for the dam-break initiation levels in the hydraulic model. 144 Paper - 9 Figure 5. Control of Heights of the Parki Main Dam (left) and Stainas Stop-dam (right) 5 Hydrological Regime of Dam-break Scenarios 5.1 Scenarios and flows Swedish dam safety guidelines (RIDAS) issued by Swedenergy and the Guidelines for the calculation of design floods for dams provided by the Swedish Committee for Design Flood Determination (Flödeskommittén) has defined the requirements for the dam-break scenarios to be analyzed. According to these requirements the possible dam-break events must be analyzed both for the normal flow situations and for the high flow situations. These kinds of calculations serve to the purpose of consequence classification of the dams, consequence analysis due to dam-break as well as for the design flow calculations. Generally and in accordance with the guidelines, it is necessary to analyze (at least) three dam-break scenarios with the classified hydrological loads: • QP (Sunny Day Failure) – dam-break during normal operation situations. This scenario describes a situation when the dam-break occurs unexpectedly due to collapse of the dam body, piping or stability loss. Hydrological regime in the river corresponds to monthly averages of production discharges. The water levels in the reservoirs are maintained at the maximum water level. 145 Paper - 9 • Q100 (Risk Class 2) – dam-break in combination with the Risk Class 2 flood means that the dam-break occurs during a flood of 100-year return period. Causes for dam-break can be the same as above, and the initiation of the dam-break may start at the peak of the reservoir lever or at the peak of the flood. Causes for the dam-break can also be a combination of the flood together with failure of operation of the discharge facilities. The flows in the river must correspond to the floods of 100-year probability, which are calculated according to the methodology of calculating Risk Class 2 floods. According to the requirements for Consequence Class 1 (High Consequence Dam) and Consequence Class 2 (Low Consequence Dam) dam facilities, this flow should be possible to discharge at the max water level of the reservoir or, in some cases, slightly more than max water level. Initial water levels in the reservoirs are set up to max water level marker. The initial water levels in the rivers are calculated to be corresponding Risk Class 2 level. • Q10 000 (Risk Class 1) – dam-break in combination with the Risk Class 1 Flood means that the dam-break occurs in connection to the extreme flood situation which corresponds to the flow of 10 000 years probability. The causes for dam-break can be the same as both previous scenarios, while the flood will be much higher than Risk Class 2 flood. That means that it can exceed the maximum water level in the reservoir and freeboard leading to overtopping of the dam crest. In those cases when overtopping is not occurring the causes for dam-break are the same as in previous scenarios – peak of the flood or water level in the reservoir. The flows in the river system are based on the flows calculated according to the methodology for calculating of Risk Class 1 floods. The initial water levels in the rivers are calculated to be corresponding Risk Class 1 level. 5.2 Primary and secondary dam breaks or cascading dam-breaks Most of the Swedish rivers are regulated with many cascading dam facilities. In some cases, the dam-break upstream causes the secondary dam-breaks all the way downstream. The primary dam-breaks were analyzed in this study for the Parki and Randi dam sites. The secondary dam-breaks may follow as a consequence of a primary dam-break upstream. The primary dam-break by its nature can be defined as overtopping, piping or stability loss. The secondary dam-break is defined only as overtopping in this study. 5.3 Dam-break initiation and breaching In the beginning of dam-break calculations it is necessary to describe the situation at which the dam-break initiates. For the overtopping situation, this can be accomplished by defining the reservoir level (usually this corresponds to the lowest dam crest value, but in some cases even the level of the impervious dam core is used). For the piping situation either reservoir water level or the time moment can be used. It was possible to describe the breaching process by using the parametric approach i.e. a successful breach opening by time or by using the erosion based method. The erosion-based method was chosen in order to take into account the material properties of the dam. The method is based on the Engelund-Hansen formulation and is built in within the model. To define the maximum size of the breach opening limiting sections were identified. The erosion based approach was used both for simulation of breaching caused by overtopping and piping. 146 Paper - 9 5.4 Regulation strategies under the dam-break situation Different regulation strategies can lead to significant differences regarding the consequences due to dam-break. Realistic and conservative assumptions were applied when defining the regulation strategies under the different dam-break situations. Special emphasis was laid on existing instructions for regulation strategies, which are used in the real world situations when the dam-break incidents happen. These are the following: Main principles of regulation strategies: 1. Appropriate initial conditions in the river and in the reservoirs. This links to the descriptions of different hydrological load scenarios and flows. 2. It is assumed that the electricity production stops under the dam-break situation. 3. With the increasing inflows the spillway gates opens gradually until the full spillage is reached. 4. It is assumed that the gates cannot be operated when the dam-break occurs on the same facility. The level of the gates remains at the same position as before the dam- break initiation. 5. It is assumed that opening of gates of the downstream spillway is possible when the primary dam-break occurs. This allows the spilling of water from the downstream water storages and thus freeing some storage for the dam-break flood. 6 Linking of Hydraulic Model and GIS Hydraulic modeling and GIS technology become more powerful when they are coupled into an integrated system. This combination enables bi-directional data exchange; modeling is performed in the common geographic reference system, and once the system is developed, data update procedures become simple and effective. One of the most valuable and efficient procedures using this approach is the possibility of automatic extraction of the hydraulic parameters from digital elevation model. Among various parameters available, the most important for hydraulic modeling are river geomorphologic characteristics, representing the geometry of riverbed and floodplain. In this study all data was prepared and stored in the GIS shell. This enables the efficiency in maintaining and sharing data for the various purposes within this project. As indicated in the methodology, the integration of the hydraulic model into GIS was possible by using pre-and post-processing routines. They enabled bi-directional integrated communications between GIS and hydraulic model. The major advantages of this integration were the possibilities to automatically prepare the cross-section database and to import it into hydraulic model. Data generated by the hydraulic model was transferred to the GIS in order to perform further analysis. 7 Hydraulic Modeling 7.1 Hydraulic Model As indicated in the Methodological Framework (Figure 2), different hydraulic models can be applied for the 1-D flood routing, such as HEC-RAS, DAMBRK or MIKE 11, for the 1-D flood routing. It was decided to use MIKE 11 hydrodynamic model for this study. MIKE 11 is a 1-D model used for the simulation of hydrodynamic flow and sediment transport. More information about the model can be found in the references. 147 Paper - 9 The layout of the hydraulic model for river branches from the beginning of the Parki Reservoir and the confluence with the river Stora Lule Älv is represented in the Figure 6. The figure represents all waterways and components in the model, which were used for all simulations of the dam-break scenarios. Figure 6. Layout of the Hydraulic Model of the River Lilla Lule Älv Figure 7 represents the longitudinal profile of the main branch containing 4 dams. Included are also maximum and normal water levels of the first simulation scenario presented. 148 Paper - 9 Figure 7. Plot of Longitudinal Profile of the Dam Cascade 7.2 Model calibration Control calculations have been performed in order to calibrate the model to obtain the accurate water levels, which were known from earlier hydrological simulations. Model was calibrated by using the flow resistance values, which were associated to the land cover information available as a GIS layer. 8 Modeling Results and Discussion 8.1 General By utilizing the methods and techniques described in this paper, 16 dam-break scenarios were analyzed in this study. The content of some of the scenarios was obvious in the beginning, while the other scenarios have evolved after summarizing the results of previous ones. That means, that besides the analysis of current situation in the river system, a set of “what-if?” cases were analyzed. Due to limited extent of this paper, only the most significant scenarios are presented in this chapter. For the same reason only fragments of the inundation maps are presented in this chapter and these cover only the “hot spots” in the area. 8.2 Present Conditions with the Risk Class 1 Flood This scenario covers the analysis of the dams with the present properties and Risk Class 1 flood used as a hydrological load. By routing this flood trough the river system it was concluded that two dams in the system couldn’t pass the flood and were overtopped. 149 Paper - 9 These were the Stainas Stop-Dam and the Randi Intake Dam. The breaching of these dams resulted, however, in no secondary dam-breaks in the dam cascade. Further analysis of third party consequences indicated that these two dams should be categorized as low consequence dams that require a design flood capacity of the 100-year probability. Figure 8. Inundation due to Primary Dam-break in Stainas Stop-dam and Secondary Dam-break in Randi Intake Dam under the Risk Class 2 Hydrological Conditions 8.3 Heightening of the Stainas Stop-dam and Randi Dams This scenario evolves as a natural step after obtaining the results from the scenario described in the previous chapter, i.e. it has been necessary to model what would be the consequences if the dams that breached were heightened to the level capable to pass through the Risk Class 1 flood. After the describing of these measures in the model it was possible to simulate how the heightening of the dams would change the inundation and the consequences. The difference of the consequences is presented in the next chapter while the visual representation is provided in the Figure 9. 150 Paper - 9 Figure 9. Inundation due to Piping in the Parki Main Dam after Heightening of the Stainas Stop-dam and Both Randi Dams (Green Area) in Comparison with the Present Conditions Scenario (Red Area) It is obvious from the illustration that scenario described in this chapter gives lower consequences due to dam-break. Figure 10 presented below demonstrates the water levels compared with other dam-break scenario and normal flow situation. By using this approach itwas possible to find an optimal solution for the river Lilla Lule Älv regarding the measures to be implemented, that gives the lowest possible consequences due to dam-break. 151 Paper - 9 Water Level (m.a.s.l.) Comparison of the Water Levels of Three Dam-break Scenarios 13-8-1981 15-8-1981 17-8-1981 19-8-1981 21-8-1981 23-8-1981 25-8-1981 217.0 217.5 218.0 218.5 219.0 219.5 220.0 220.5 221.0 221.5 Water Level in Present Conditions Scenario Normal Water Level Water Level due to Piping in Parki Main Dam after heightening of the Stainas Stop- dam and Both Randi Dams Figure 10. Water Level Comparison 9 Consequence Evaluation Based on the information gathered from the hydraulic model in the GIS it is possible to determine the aerial extent of various land-types as well as individual objects that will be flooded as a consequence of dam-break. These consequences were analyzed within GIS application and some of the major scenarios are summarized in the Table 1 below: Table 1. Summary of the Consequences of Characteristic Dam-break Scenarios Inundation (ha) Dwelling Houses(Units Affected) Auxiliary Buildings(Units Affected) Transformer Stations (UnitsAffected) Railroads (km) Roads (km) Power Lines (km) Reference Scenarios Without Dam-break Scenario 10. Reference QP 0 0 0 0 0 0 0 Scenario 9. Reference Q100 (Risk Class 2 Flood) 622 15 29 0 0 26 10 Dam-break Scenarios with Risk Class 1 Flood (Q10 000) Scenario 1. Current Conditions 1700 75 103 2 24 202 212 Scenario 2. Current conditions and piping in Parki Main Dam 2833 134 169 2 19 528 306 152 Paper - 9 Scenario 5. Stainas Stop-dam and both Randi dams are heightened and piping in Parki Main Dam 2667 100 130 3 22 221 254 Scenario 7. Stainas Stop-dam and all Parki dams are heightened and piping in Randi Regulation Dam 1750 81 110 3 24 215 160 Dam-break Scenario with Risk Class 2 Flood (Q100) Scenario 11. Piping in Randi Intake Dam 553 12 26 0 0 32 8 10 Conclusions This study demonstrates the efficiency of the methods of applying an integrated approach by combining different technologies in order to understand how the regulated waterways are functioning. By having this understanding and by having the data necessary it is possible to make the decisions which are based on the material of high accuracy. This leads to better solutions when implementing dam safety measures or when dealing with other tasks such as emergency planning. All dam facilities are described in one model, which enables simulation of the whole river system i.e. the cascading dam-breaks, which leads to better understanding of possible consequences of such event. Results from the study presented in this paper demonstrated that hydrological upgrading is necessary for the Parki main dam. At present pre-studies on the hydrological upgrading, as well as general dam safety improvements are about to be initiated at both Parki and Randi facilities. It is planned to maintain the system for the future applications by adding new features and by updating data. 11 References 1. Parki – Technical Data (1968). Swedish State Power Board. 2. Randi kraftstation (1974). Vattenfall (in Swedish). 3. Parki – Technical Data (1971). Swedish State Power Board. 4. Parki – Technical Data (1971). Swedish State Power Board. 5. Consequence Analysis due to Dam-break: A Pilot Study for Messaure (2002). Report. SwedPower. 6. RIDAS - Swedish Dam Safety Guidelines (2002), Swedenergy AB. 7. Flödeskommittén -Riktlinjer för bestämning av dimensionerande flöden för dammanläggningar (1990). (In Swedish: Guidelines for the calculation of design floods for dams). Final report from the Swedish Committee for Design Flood Determination. Swedish State Power Board, Swedish Power Association and Swedish Meteorological and Hydrological Institute, Stockholm and Norrköping, Sweden. 8. MIKE 11 – Users Guide (2003). DHI Water & Environment. 9. MIKE 11 – Reference Manual (2003). DHI Water & Environment. 10. Flodvågsberäkning för Luleälven (1997). Vattenfall Hydropower AB. (In Swedish: Hydraulic calculations for river Lule Älv) 11. Flodvågsberäkning vid klass 2 flöde för 10 dammar i Luleälven (1997). Vattenfall Hydropower AB. (In Swedish: Hydraulic calculations with the Risk Class 2 flow for 10 dams in the river Lule Älv) 153 Paper - 10 TWO-DIMENSIONAL MODEL FOR EMBANKMENT DAM BREACH FORMATION AND FLOOD WAVE GENERATION By David C. Froehlich1 Introduction Earthen embankments that serve as dams to impound water or as levees to prevent rivers from overflowing sometimes fail catastrophically from the erosive action of water overtopping them. Gaps or breaches that form in the embankments allow water to flow through them without control, often producing floods that cause great damage or suffering (Figure 1). Characteristics of flood waves issuing from breached embankments depend largely on interactions between flow and the morphological development around the openings. A two-dimensional depth-averaged flow and model known as DaveF that allows breach development and the resulting flood wave to be simulated simultaneously is presented here. The governing partial differential equations are solved by means of a finite volume technique with explicit time discretization. This method is locally conservative, has built- in stability mechanisms such as upwinding, and allows for nonconforming meshes. The model can simulate transport processes that dominate rapidly varying flows in natural channels where depth-averaging is well- grounded. The model was used to simulate the controlled failure from overtopping flow of a large-scale experimental embankment six meters high composed of cohesive clay. Taking everything into account, good agreement was obtained between observed and calculated breach development. Field Test Embankment The experimental embankment was built in a narrow section of river channel downstream of the reservoir Røssvatnet near the city of Mo I Rana, Norway as part of a European Commission study known as the IMPACT Project. Normally there is no outflow from the reservoir and the downstream channel is dry. A 6-m high homogenous embankment about 36 m in length composed of cohesive silty clay (25% clay, 65% silt, and 10% sand) was constructed for the test (Figure 2). Characteristics of the embankment are summarized in Table 1. 1 Consulting Engineer, 303 Frenchmans Bluff Drive, Cary, North Carolina 27513-5662, Tel:919-468-8724, Email: dcfroehlich@aol.com. Figure 1. At least 20 people were killed when the Zeyzoun Dam in Syria failed on 4 Jun 2002. Several villages were flooded by depths up to four meters. 154 Paper - 10 To initiate breaching, a notch 0.5 m deep and about 2 m wide was cut through the center of embankment crests. Impoundments behind the embankments were filled to crest level, and upstream water-surface elevations were held nearly constant during the initial overtopping. Parameters for the test embankment are summarized in Table 1. Figure 2. Experimental embankment at test site near Mo I Rana, Norway just before start of overtopping Table 1. Experimental Embankment Characteristics Embankment parameter Value Height Crest elevation Crest width Crest length Embankment slope: Upstream Downstream Median grain size, D50 Porosity Dry unit weight, ?d Friction angle, f Cohesion Plasticity index Manning coefficient, n 6.0 m 670.81 m 2.0 m 36.0 m 2:12:1 0.010 mm 0.200 17,000 N 10o 25,000 N/m2 15 0.030 155 Paper - 10 Governing Equations The numerical model is based on conservation forms of the depth-averaged fluid, and momentum transport relations, which comprise a coupled system of nonlinear, hyperbolic, partial differential equations as follows: ¶ ¶ U t ×[(),( )] +Q =0 (1) +ÑFU GU where ï q ïï q ïï 0 ï ìüh 12 ïï ïï ï ïq21 ï ï ï ï qq ï ï ï¶z 1 ï U =íý í12 ýí 12 ý b q F =+gh G = Q =ígh +t ý (2) 1 bx ïï ïh 2 ïï h ï ï¶x r ïq2 ï qq ïï2 ïï ïîþ ï 12 ïï q2 +1 gh2 ïïgh ¶zb +1 tby ï î h þîh 2 þ î¶y r þ t = time, h = water depth, q1 and q2 = unit flow rates in the horizontal x and y Cartesian coordinate directions respectively, q =q12 +q22 =total unit flow rate, g = gravitational acceleration, zb = bed elevation, r = water mass density, and t bx and t by = bed shear stresses in the x and y directions respectively. Bed shear stresses are given by f bq1 q122 +q22 and t by =rcmf bq2 q122 +q22 (3) tbx =rcm hh where gn2 cf = y 2h1/3 (4) Is a dimensionless bed friction factor, n = Manning’s roughness coefficient, ? = units factor (1.0 for SI, 1.486 for U.S. Customary), and mb =1+çæ ¶ ¶ zb ÷ö 2 +çæ ¶zb ÷ö 2 (5) èx ø è¶y ø is a coefficient that accounts for a sloping bed. defined, although some guidance is available 156 Paper - 10 Embankment Erosion Erosion from embankment surfaces is accounted for using a simple empirical erosion rate formula. Transport of embankment soil is not considered, once eroded the sediment is considered to be removed from the vicinity of the embankment immediately and have no further effect. The rate of soil eroded from an embankment surface is given by E& = Kd (t -t ) (6) bc r s when tb >tc , where Kd = detachment rate constant that depends on original bed material properties, rs = sediment mass density, tb = bed shear stress acting in the flow direction, and tc = detachment threshold bed shear stress. Development of the bed is tracked by the mass conservation expression (1 -h) ¶ ¶ ztb =-E& (7) where h = bed material porosity. Embankment erosion calculations then require four parameters: Kd, tc , ?, and the Manning roughness coefficient of the soil surface. Values of Kd andtc for soils of different textural classes found from a series of onsite experiments carried out in narrow channels formed in 33 natural soils are presented in Table 2 (Flanagan and Livingston 1995). These coefficients have been found to provide estimates of embankment erosion and breach development that are in accordance with reason. Nevertheless, other appropriate sources of equivalent coefficients, or direct measures of soil erosion, can be used to obtain the needed coefficient values. Table 2. Detachment rate erosion coefficients and detachment threshold bed shear stresses for various soil textural classifications. Detachment Soil textural Detachment threshold bed classification rate constant, Kd shear stress, ct (kg/s/m2/Pa) (Pa) Clay loam 0.0048 4.7 Loam 0.0085 3.3 Sand 0.0250 2.1 Sandy loam 0.0100 2.5 Silt loam 0.0120 3.5 Clay 0.0089 2.9 Silty clay 0.0120 4.8 Silty clay loam 0.0053 3.2 157 Paper - 10 Model Formulation The depth-averaged surface-water flow and sediment transport equations are solved numerically using a two-dimensional, cell-centered, Godunov-type, finite volume scheme. Godunov-type methods for solving first-order hyperbolic systems of equations are based on solutions of initial value problems, known as Riemann problems, involving discontinuous neighboring states. These methods are found to be accurate and robust when used to solve partial differential equations describing fluid motion, being able to capture locations of shocks and contact surfaces. Values of the conserved variables are calculated for each of the volumes or cells. Cells can be any convex polygon, but are limited in DaveF to triangles and quadrilaterals. Bed slope source terms are taken into account by combining them with edge fluxes in a manner that leads to proper balance of forces. Discretization of integral forms of the conservation equations (1) by the finite volume method assures that the basic quantities, mass and momentum, will be conserved across discontinuities (Hirsch 1988). Applying the divergence theorem to (1) and integrating over an arbitrary cell Ei gives the basic finite volume equation ¶ ¶ òUdA +. [() ()] ×ndS +òQ dA (8) FU ,GU = 0 t Ei ¶Ei Ei where nº [,]nn º qq [cos ,sin] = outward unit normal vector to the boundary Si, ? = angle 12between the positive x direction and the vector n, and dA and dS are the area and boundary elements of the cell, respectively. Making use of the rotational invariance property of (1) (Toro 2001, p. 65), the normal flux component through a surface is given by n × [(),( )]FU GU =-1() (9) T FTU where T = T(?) = rotation matrix given as é10 0 ù T = êê 0 cosq sinq ú ú (10) ê0 -sinq cosq ú ëû Using the normal flux expression given by (9), the integral relation (8) becomes ¶ ¶ ò dA +. -1( )dS +òQ dA = 0 (11) U T FTU t Ei ¶Ei Ei Within each cell U is considered to be constant, and the flux across any edge is based on the states in the two adjacent cells. Letting U% º TU = variables transformed to the edge- normal/tangential directions, (11) becomes 158 Paper - 10 ¶-1 %% UdA+ (, ))dS+QdA TFUU =0 ¶ò . LR ò (12) t Ei ¶Ei Ei where the subscripts L and R denote cells to the left and right of an edge, respectively, when circumnavigating a cell in a counter-clockwise direction. Approximating the boundary integral in (12) by single point quadrature gives -1 %% -1 %%l (, ))dS =(, TFUU TFUU) . LR å LRj (13) ¶Ei j where lj = length of edge j. Godunov (1959) calculates the numerical flux by solving the local one-dimensional Riemann problem in the direction normal to the cell edge. Since exact solutions are comparatively time-consuming, many approximate Riemann solvers have been developed for fluid dynamics problems. An approximate Riemann solver suggested by Harten, Lax, and van Leer (1983), commonly known as the HLL solver, is used in to calculate edge fluxes. The HLL solver is straightforward to implement in comparison to some other methods, and it has proven robust in practice. The technique is founded on division of the Riemann problem solution space into three constant states separated by two waves traveling with celerities sL and sR. Based on this notion, numerical fluxes are approximated as follows: ì FU( % L ) forsL >0 FU U(,) = FU) fors %% . í ( % <0 (14) LR RR ï* %% FUU(,) fors 0 . LR L ££sR where * s ( %) -s ( %) +ss ( %% FUU( L , R ) = FU FU U -U) (15) %% R L L R LRRL sR -sL ** ** s =min{u%-gh , u%-gh } sR =max{u%R +gh R , u%+gh } (16) L LL u%* = 1 %% ( (uL +uR )-gh L +ghR )çæ hR -hL ÷ ö (17) 2 èhL +hR ø %%u h* = 21 (hL +hR )çç æ è 1 -12 gh uLR + -LghR ÷ö ø (18) and u%=ucosq +v sinq =velocity normal to the edge under consideration. 159 Paper - 10 Bed Slope Terms Centered approximation of bed slope source terms in (12) unfortunately leads to unbalanced forces for non-horizontal beds, which prevents retrieval of trivial solutions having horizontal water surfaces and motionless states when only no-flux boundary conditions are applied around the computational mesh. Bermúdez and Vázquez-Cendón (1994), Bermúdez et al. (1996) and Leveque (1998) address this predicament for other Godunov-type finite volume schemes. To provide a proper balance of pressure and gravitational forces using the HLL solver, bed slope terms are merged with edge flux vectors as follows: ì 0 ü ¢ % ïï 1 ïï ()%= FU() + 2 g h L + h )(zbR FU í ( R -zbL )ý (19) . 0 . where ¢()= modified normal flux vector. Source terms are modified accordingly as FU% ï 0 ï ï 1 ï Q¢= í t bx ý (20) ï r ï ï 1 ï ï t by ï î r þ The finite volume problem statement then becomes ¶-1 ¢ %%l ¢ UdA +. (, R ) + Q dA = 0 ¶tE ò ij TFUUL jE ò i (21) FU ¢() with ()%replaced by FU%in (14) and (15). Solution of the Discrete System The local one-dimensional problem given by (21) is solved using Strang splitting (1968) whereby the pure advection problem given by the homogeneous part, that is, ¶-1 %%l UdA + ¢(,) = 0 TFUU òå LRj (22) ¶tEi j is solved first by evaluating the area integral and applying forward Euler time integration to obtain 160 Paper - 10 adv n Dt æ-1 %% lö U = U -TFUU¢(,) (23) i iAi çå j LRj ÷ èø adv where Ui = the advection-only solution. This is followed by solution of the ordinary differential equation dUi adv )0 + QU¢( i = (24) dt that accounts for source terms due to bed friction and sediment erosion/deposition. Solving (24) with forward Euler time integration gives n+1 n Dt æ-1 %% lö adv Ai è j TFUU ¢ ø ( = (25) U = U -(,) -D tQU¢ )0 ii çå LR j ÷ i which is considered to be the standard splitting scheme (Toro 2001, p. 233). The Field Test Model A computational mesh consisting of a mixture of triangular and quadrilateral cells covering a 300 m length of river channel was used to simulate overtopping flow and breach development (Figure 3). Bed elevations defined by the mesh are shown in Figure 4. Erosive soil was modeled only in embankment areas, although eroded embankment soils deposited downstream could be re-entrained and transported further downstream. Initial conditions consisted of a level water-surface equal to the embankment crest elevations and motionless states in the upstream channel. Constant head boundary conditions maintaining the initial water-surface elevations were applied at the upstream end of the channel in each case, and free-outfall Figure 3. Finite volume mesh Figure 4. Color contours conditions (that is, consisting of triangles and showing bed elevations quadrilaterals covering a 300 m represented in the field test critical depth) was set section of river channel finite volume mesh. Water flows at the downstream end. containing the test embankment. from the bottom to top of the 161 Paper - 10 Field Test Simulations The computational mesh for the test embankment is shown in Figure 5 along with bed elevation isocolors. Cells along the embankment crest are mostly 1 m squares. A few smaller cells are used to define the sides of the notch cut through the crest at the embankment center. For erosion calculations, the embankment soil was considered to be a silt loam, and was assigned a detachment rate constant Kd = 0.0120 kg/s/m2, and a detachment threshold shear stress tc = 3.5 N/m2. Combined bed elevation isocolor and velocity vector plots for conditions at various times during the first four hours of the simulation are shown in Figures 6 through 13. Profiles along the dam crest for various simulation times given in Figure 14 show how the central portion of the breach develops. Transects of the embankment at the location of the initial notch are shown in Figure 15. From the figures depicting geomorphic development of the breach, it can be seen that only the downstream slope erodes for the first hour of the simulation. Erosion begins at the toe of the downstream slope, creating a sharp break in grade, forming a scarp or headcut that migrates upstream towards the crest. Experimental studies in homogeneous cohesive soils show that under certain conditions a headcut can maintain a vertical face (Leopold et al. 1964, p. 442), which is evident in at least one of the transects. The headcut reaches the crest of the embankment after about one hour of overtopping. As the upstream slope begins to erode downward, outflow from the breach increases. Maximum sustained outflow from the breach Figure 5. Finite volume mesh showing color contours of bed elevations at the test embankment. 162 Paper - 10 Figure 6. Test embankment at 0:00. Figure 7. Test embankment at 0:30. Figure 8. Test embankment at 1:00. Figure 9. Test embankment at 1:30. Figure 10. Test embankment at 2:00. Figure 11. Test embankment at 2:30. Figure 12. Test embankment at 3:00. Figure 13. Test embankment at 4:00. 163 Paper - 10 was reached after about three and one-half hours, at which time the embankment had eroded nearly to its base level, and the breach width, as shown in the dam crest profile plot, has nearly reached is maximum value. Photographs of the embankment at various times after the beginning of overtopping are shown in Figures 16 through 20. Simulated embankment centerline cross sections at various times are also shown in Figure 20 and can be compared to the photograph taken just after the crest was breached at about 1:45 h. Field Test #1 - Dam Crest Profile 372 0 10 20 30 40 Crest Distance (m) 0:00 1:00 1:30 2:00 3:00 6:00 Elevation (m) 370 368 366 364 Figure 14. Test embankment centerline profile for various simulation times showing geomorphic development of the breach. Embankment Transect at Notch 371 370 0 5 10 15 20 25 Distance (m) 0:00 0:30 1:00 1:30 2:00 3:00 4:00 5:00 Elevation (m) 369 368 367 366 365 364 Figure 15. Test embankment transect through the notch at various simulation times showing geomorphic development of the breach. 164 Paper - 10 Figure 16. Start of over overtopping flow. Figure 17. Beginning of erosion at the toe of the embankment as headcuts are formed. 165 Paper - 10 Figure 18. Erosion of the embankment after about one hour of overtopping. Figure 19. Breach after nearly reaching its maximum size. 166 Paper - 10 Figure 20. Embankment cross sections at various times and photograph of embankment as theembankment crest is breached, about time 1:45. Summary and Conclusions Erosion of a large-scale experimental earthen embankment from erosion by overtopping was simulated using a two-dimensional depth-averaged flow model that allows breach development and the resulting flood wave to be simulated simultaneously. The governing partial differential equations were solved by means of a Godunov-type finite volume technique using an approximate Riemann solver. The model can simulate all transport processes that dominate rapidly varying flows in natural channels where depth-averaging is well-grounded. The model also shows that erosion began at the toe of the downstream slope, creating distinct scarps or headcuts that migrated upstream towards the embankment crest, as observed during the test. The numerical simulations also predict the test embankment to form in the general shape of a trapezoid, first eroding downward to an erosion resistant base level, and then expanding laterally, also as observed. Taking everything into account, good agreement was obtained between the observed breach development and the calculated embankment erosion and flood wave generation. 167 References Paper - 10 Bermúdez, A., and Vázquez, M. E. (1994). “Upwind methods for hyperbolic conservation laws with source terms.” Computers and Fluids, 23, 1049-1071. Bermúdez, A., Dervieux, A., Désidéri, J. A., and Vázquez, M. E. (1995). “Upwind schemes for the two-dimensional shallow water equations with variable depth using unstructured meshes.” Computer Methods in Applied mechanics and Engineering, 155, 49-72. Flanagan, D. C., and Livingston, S. J. (1995). “WEPP user summary.” NSERL Report No. 11, National Soil Erosion Research Laboratory, U.S. Department of Agriculture, Agricultural Research Service, West Lafayette, Indiana. Froehlich, D. C. (1995). “Embankment dam breach parameters revisited.” Proceedings of the First International Conference on Water Resources Engineering, San Antonio, Texas, August 14-18, 1995, American Society of Civil Engineers, New York, NY, 887-891. Godunov, S. K. (1959). “A difference scheme for numerical computation of discontinuous solutions of hydrodynamics equations.” Matemsticheskly Sboraik, 47 (translated by U.S. Joint Publications research Service). Harten, A., Lax, P. D., and van Leer, B. (1983). “On upstream differencing and Godunov-type schemes for hyperbolic conservation laws.” SIAM Review, 25(1), 35-61. Hirsch, C. (1988) Numerical computation of internal and external flows, volume 1, fundamentals of numerical discretization. John Wiley and Sons, new York, New York. Leopold, L. B., Wolman, M. G., and Miller, J. P. (1964). Fluvial processes in geomorphology. W. H. Freeman and Company, San Francisco, California. LeVeque, R. J. (1998). “Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm.” Journal of Computational Physics, 146, 346-365. Pugh, C. A. (1985). “Hydraulic model studies of fuseplug embankments.” Report No. RECERC- 85-7, U.S. Department of the Interior, Bureau of Reclamation, Denver, Colorado. Roe, P. L. (1981). “Approximate Riemann solvers, parameter vectors, and difference schemes.” Journal of Computational Physics, 43, 357-372. Strang, G. (1968). “On the construction and comparison of finite difference schemes.” SIAM Journal on Numerical Analysis, 5(3), 506-517. Toro, E. F. (1997). Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer-Verlag, Berlin. Toro, E. F. (2001). Shock-capturing methods for free-surface shallow flows. John Wiley and Sons, Chichester, United Kingdom. 168