TWO-DIMENSIONAL FINITE ELEMENT RIVER MORPHOLOGY MODEL by JOSE A L F R E D O V A S Q U E Z B.Sc. (Eng.), University of Piura, 1992 M.Eng., University of British Columbia, 1999 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE F A C U L T Y OF G R A D U A T E STUDIES (Civil Engineering) THE UNIVERSITY OF BRITISH C O L U M B I A December 2005 © Jose Alfredo Vasquez, 2005 A B S T R A C T Alluvial streams are those that flow over a movable bed that they have formed by the transport of sediment, either suspended in the water current or dragged along the bottom. Bedload represents the fraction of that sediment that is transported along the channel bottom. Since bedload transport is associated with the removal (erosion) or deposit of sediment over the bed, it is the main transport mechanism responsible for the morphological changes in alluvial channels. Channel changes made by sediment deposition increase the risk of river flooding, block navigation canals, reduce the life span of reservoirs and modify the course of streams. While erosion can undermine the foundation of built structures (e.g. bridges, pipelines, bank revetments), scour riparian habitats or induce drawdown in the water table levels, just to mention a few effects. For those reasons, bedload transport in alluvial streams is a topic of high interest for both river engineers and fluvial morphologists. In this thesis, a two-dimensional numerical model based on the Finite Element method was developed for the morphodynamic simulation of scour and deposition in straight and curved alluvial channels. The model uses unstructured meshes formed by elements of triangular shape. The model was developed by coupling a solver for the bedload sediment continuity equation to existing fixed-bed hydrodynamic models to dynamically update the channel bed elevation as sediment is transport by the flow. The model was successfully applied to reproduce experimental observation of bed degradation and aggradation along straight flumes; as well as scour and deposition along curved channels. Data from a meandering river was also applied to test the model for full-scale applications. For the tests performed, the model proved stable and accurate. There are two versions of the model. A conventional vertically averaged (VA) model and a novel vertically averaged and moment (VAM) model. The V A M is a quasi-3D flow model that provides more detail of the vertical structure of the flow than a V A model for curved channels. The model has some unique capabilities for modeling rapidly varied flows and transcritical flow (supercritical flow and hydraulic jumps) over alluvial beds. u This is probably the first two-dimensional river morphology model based on unstructured triangular Finite Elements that has combined capabilities for bedload transport, secondary flow effect in bends, transverse bed slope and transcritical flow. The results of this thesis suggest that the model has the potential for future morphodynamic applications to dam-breaks, dam removal and steep mountain streams. 111 T A B L E OF C O N T E N T S A B S T R A C T ii T A B L E OF CONTENTS iv LIST OF T A B L E S viii LIST OF FIGURES ix LIST OF S Y M B O L S A N D ABBREVIATIONS xi A K N O W L E D G E M E N T S xv 1 INTRODUCTION 1 1.1 STATEMENT OF P R O B L E M A N D MOTIVATION 1 1.2 BRIEF HISTORICAL PERSPECTIVE 4 1.3 THESIS OBJECTIVES 9 1.4 LITERATURE REVIEW 12 1.4.1 Empirical and analytical models 12 1.4.2 Morphodynamic models in 1, 2 and 3 dimensions 18 1.4.3 Unstructured meshes and Finite Element hydrodynamic models 24 1.4.4 Curvature effects and VAMmodel 27 1.4.5 Transcritical flow morphodynamics 33 1.4.6 Relevant issues of literature review 35 1.5 S U M M A R Y 36 2 STRAIGHT ALLUVIAL CHANNELS 39 2.1 INTRODUCTION 39 iv 2.1.1 Bed aggradation and degradation 39 2.1.2 River2D Hydrodynamic Model 41 2.1.3 Bedload Sediment Model 44 2.1.4 Objectives 46 2.2 BED A G G R A D A T I O N TEST C A S E 47 2.2.1 Soni et al. experiment 47 2.2.2 Results 48 2.2.3 Discussion 49 2.3 BED DEGRADATION TEST C A S E 50 2.3.1 Suryanarayana experiment 50 2.3.2 Results 51 2.3.3 Discussion 52 2.4 KNICKPOINT MIGRATION TEST CASE 53 2.4.1 Brush and Wolman experiment 53 2.4.2 Results • 55 2.4.3 Discussion 58 2.5 S U M M A R Y 59 3 CURVED ALLUVIAL CHANNELS 61 3.1 INTRODUCTION 61 3.1.1 Bedload sediment model 62 3.1.2 Influence of secondary flow (helical motion) 64 3.1.3 Curvature and inertial adaptation 66 3.1.4 Influence of transverse bed slope 68 3.1.5 Objectives 69 3.2 E X P E R I M E N T A L D A T A A N D N U M E R I C A L M O D E L 70 3.3 RESULTS 73 3.4 DISCUSSION 76 3.5 S U M M A R Y 78 v 4 VERTICALLY-AVERAGED AND MOMENT OF MOMENTUM MODEL 79 4.1 INTRODUCTION.. 79 4.1.1 Objectives 82 4.2 THE N U M E R I C A L M O D E L 82 4.2.1 The VAM hydrodynamic model 82 4.2.2 The bedload sub-model 88 4.2.3 Bed slope effect 89 4.2.4 Secondary flow effect on the bed shear 89 4.3 TEST CASES A N D N U M E R I C A L M O D E L 91 4.3.1 Experimental data. -. 91 4.3.2 Numerical model 94 4.4 RESULTS 95 4.5 DISCUSSION 99 4.5.1 Prediction capabilities of the VAM model 99 4.5.2 Comparison with VA model 101 4.5.3 Future research 102 4.6 S U M M A R Y 104 5 SUMMARY AND CONCLUSIONS 105 5.1 S U M M A R Y 105 5.2 THESIS CONTRIBUTION 107 5.3 FUTURE R E S E A R C H 109 5.3.1 Convection-dominated applications 109 5.3.2 Sediment transport 110 5.3.3 Dam-break and dam removal Ill 5.3.4 VAM model and non-linear secondary flow correction 112 REFERENCES 114 vi APPENDIX A. Discretization of the sediment continuity equation 129 APPENDIX B. Running the computer models 144 APPENDIX C. Computer program: algorithms and functions 154 APPENDIX D. Curvature and inertial adaptation 167 APPENDIX E. Varying width applications 178 APPENDIX F. Convection-dominated applications: migrating bedforms 194 APPENDIX G. Qualitative simulations of dyke breach and dam break 213 vii LIST OF T A B L E S Table 1 Examples of some 2D and 3D numerical models for alluvial channel morphology. 20 Table 2. Main parameters of the experiments in alluvial bends performed at the Delft Hydraulics Lab (DHL), the Laboratory of Fluid Mechanics (LFM) curved channels (Struiksma et al. 1985) 71 Table 3. Main parameters of the DHL and L F M flumes (Struiksma et al. 1985) and the Waal River (Struiksma 1985) 92 viii LIST OF F I G U R E S Figure 1. Examples of (a) structured curvilinear mesh based on quadrilateral elements and (b) unstructured mesh based on triangular elements 11 Figure 2. Aggradation by overloading of the supply of solid discharge Aqs 17 Figure 3. Indexing of elements in (a) structured mesh and (b) unstructured mesh 24 Figure 4. Sketched of assumed linear vertical velocity profile for ID V A M model (flow moving in the x direction; not to scale) 30 Figure 5. Flow chart of morphological model. Bedload module is inside dashed box 46 Figure 6. Observed and computed longitudinal profiles at 15 and 40 minutes for Soni et al. (1980) bed aggradation experiment due to an overload of sediment equal to four times the equilibrium value 48 Figure 7. Observed and computed longitudinal profiles at 10 hours for bed degradation experiment due to sediment supply shut-off for two correction coefficient values of Cr = 1.0 and Cr =0.4 52 Figure 8. Refined Finite Element mesh around knickpoint and steep 10% slope reach (gray lines are bed contours every 0.001 m) 55 Figure 9. Initial conditions for knickpoint migration experiment (t = 0). Water surface elevation (WSE) was computed by the numerical model 56 Figure 10. Observed and computed longitudinal profiles at t = 2 hr 40 min for knickpoint migration experiment 57 Figure 11. Secondary circulation and deviation of bed shear stress along bends 65 Figure 12. Example of vector fields of depth-averaged velocity and bedload sediment transport, computed in a flat bed using the secondary flow correction given by equation (3-5) 67 Figure 13. Influence of the transverse slope S, on the direction a of sediment transport 68 Figure 14. Delft Hydraulics Lab (DHL) 140° curved flume and Laboratory of Fluid Mechanics (LFM) 180° curved flume. Straight reaches up/downstream of bends are 15 m for D H L flume and 10 m for L F M flume 72 ix Figure 15. Details of the computational meshes for the DHL and L F M curved flumes 72 Figure 16. Comparison between the computed and measured longitudinal bed profiles for D H L flume, experiments T l and T2. Profiles are 0.375 m from right and left banks 74 Figure 17. Comparison between the computed and measured longitudinal bed profiles for L F M flume, experiment T6. Profiles are 0.34 m from right and left banks 75 Figure 18. Redistribution of the transverse horizontal velocity and flattening of the vertical velocity profile caused by the secondary flow in bends 80 Figure 19. Linear vertical distributions of velocity assumed by the vertically-averaged and moment (VAM) model over the water depth h: u(z) in the x-direction and v(z) in the ^-direction 84 Figure 20. Simulation of the helical motion in a bend by the V A M model: (a) surface velocity moving towards outer bank, and (b) near-bed velocity moving towards inner bank 88 Figure 21. Study reach of the Waal River (after Struiksma 1985) 94 Figure 22. Detail of the computational mesh for the Waal River (lower bend) 95 Figure 23. Longitudinal bed profiles computed by V A M model for D H L experiment (VA results shown for comparison) 96 Figure 24. Longitudinal bed profiles computed by V A M model for L F M experiment (VA results shown for comparison) 97 Figure 25. Comparison between results computed by V A M model with 1969 soundings along profiles 50 m away from banks of the Waal River 98 Figure 26 Plan view of the computed water depths for Waal River showing the pools along outer banks and point bars along inner banks (initial depth was 5 m) 99 x LIST OF S Y M B O L S A N D A B B R E V I A T I O N S a = coefficient of empirical sediment transport equation. ac = centrifugal acceleration (m2/s) A = coefficient in secondary flow correction b = exponent of empirical sediment transport equation. B = channel width (m) c = total load concentration in volume C = Chezy friction coefficient (m1 /2/s) C* = C / yjg dimensionless Chezy coefficient Cr = correction coefficient for degradation tests D 5 0 = median grain size of bed material (m) D = constant E = constant e = base of natural logarithms fs = shape factor F = Froude number g - 9.8 m /s = gravitational acceleration Gs = 2.65 = specific gravity of sediment K = constant in Lisle et al. (2001) sediment continuity equation k = diffusion coefficient (m2/s) ks = effective roughness height (m) kj = shape factor for grain area &2 = shape factor for grain volume h = near-bed velocity coefficient xi h = depth of flow measured vertically (m) h0 = initial water depth (m) Lb = length of bend (m) rc = radius of curvature of streamlines (m) Rc = centerline radius of circular bend (m) Q = total inflow discharge (m3/s or L/s) qs = volumetric bedload sediment transport rate (m Is) qsx = qscosa = sediment transport rate in x-direction (m /s) qsy = qss'ma = sediment transport rate in ^-direction (m Is) qse = equilibrium sediment transport rate (m2/s) qsiN = upstream inflow of sediment transport (m2/s) qx = discharge per unit width in x-direction (m /s) . qy = discharge per unit width in y-direction (m2/s) S = water surface slope S0 = bed slope Sf = (U/Cflh = friction slope St = transverse bed slope in bends S,e = equilibrium transverse bed slope in bends t = time (s) T = total simulation time (s) U = - ^ / U Q + V Q = depth-averaged velocity (m/s) u = vertically -average velocity in x-direction (VA model) v = vertically -average velocity in _y-direction (VA model) UQ = vertically-averaged velocity in x-direction ( V A M model) xii vo = vertically-averaged velocity in ^-direction ( V A M model) ui = surface velocity in excess of tio (m/s) vi = surface velocity in excess of vo (m/s) Ub = uo - hu] = near-bed velocity in x-direction (m/s) Vb = vo - fovi = near-bed velocity in ^-direction (m/s) W = submerged weight of sediment grain (kg) x = horizontal co-ordinate, longitudinal (m) y — horizontal co-ordinate, transverse (m) Zb = bed elevation (m) a = direction angle of bedload transport B = inertial adaptation coefficient K = 0.4 = Von Karman constant S = direction angle of bed shear stress Ss = -tcm~x(Ah/rc) = deviation of bed shear stress relative to mean flow At = time step (s) Ax = average mesh size (m) X = sediment porosity p = water density (kg/m3) Tb = bed shear stress (N/m2) T* = dimensionless bed shear stress (Shield's parameter) D H L = Delft Hydraulics Lab 140° curved flume FD = Finite Difference FE = Finite Element FV = Finite Volume xiii G F E M = Galerkin Finite Element Method L F M = Laboratory of Fluid Mechanics 180° curved flume M O M = moment-of-momentum NET = non-equilibrium transport PDE = partial differential equation V A = vertically-averaged V A M = vertically-averaged and moment 2D two-dimensional 3D three-dimensional A K N O W L E D G E M E N T S Dr. Robert Millar, my thesis supervisor, provided the support, advice and encouragement to make this thesis possible, but above all that, he offered me his friendship. Dr. Peter Steffler gave not only the source code of the computer programs used in this research, but also invaluable assistance and guidance to develop the morphological numerical model. The reviews and comments of Dr. Michael Church, Dr. Gregory Lawrence and Dr. Sheng L i improved the final writing of the thesis, especially the thorough review of Dr. Church. Brief conversations with Edmond Y u helped me better understand the computer language. Cynthia Bluteau and Dianne Stewart made the initial tests of the model for the Waal River; their help is greatly appreciated. This thesis is dedicated to my wife Sofia and to my daughters Sofy and Winy. Their love, dedication and patience gave me the strength and encouragement to finish this work. xv 1 INTRODUCTION 1.1 S T A T E M E N T OF P R O B L E M A N D M O T I V A T I O N In 1998, during the last episode of El Nino in Northern Peru, the hyper-arid region of Piura experienced a total rainfall of about 2000 mm, more than an order of magnitude greater than the annual average of only 60 mm/year. This dramatic increase in rainfall produced a massive increase in the surface runoff and hence in the flow discharges of all the streams and rivers in the region. A l l streams, from small arroyos to large rivers, suffered major morphological changes, such as bed degradation, bank erosion and meander migration. Streams increasing their channels widths by more than 2 or 3 times were not uncommon. As a result, 68 bridges were damaged and 600 km of roads destroyed (Vasquez 2004). In Piura City alone, 2 of the 4 bridges suddenly collapsed, killing dozens of people. These extreme events are an interesting example of fluvial dynamics and at the same time my personal motivation for this thesis. I was an eye witness to these events. From the engineering perspective, knowledge of fluvial morphology is vital for the design of river structures (e.g. bridges, culverts, dykes, water intakes, dams, pipeline crossings, bank protections) and to assess the impact of certain river activities (e.g. sand and gravel mining, water diversion, channel dredging). Unfortunately, river engineers face a formidable problem. 1 Open channel flow in alluvial rivers has two movable boundaries, the water surface on the top and the riverbed on the bottom. The mutual interaction between water and the loose bottom boundary is extremely intricate. The movement of bottom sediment is produced by a complex three-dimensional turbulent flow, whose structure changes continuously in space and time. The properties of the sediment, such as size, shape, and imbrication, also vary in large extent and may dynamically change as water flows over the bed. At this time, a complete theoretical description of the interaction between flow and sediment transport has not been achieved. Most of the available .relationships for sediment transport are mainly empirical, lacking universal applicability. For a long time, river engineers could only resort to empirical methods, such as physical modeling, for dealing with complex river engineering problems (Jansen 1979). Physical or scale models are small reproductions of the actual river (prototype). Physical modeling is an excellent tool to study near-field phenomena, such as the interaction between rivers and local structures (e.g. bridges, pipelines, intakes). Since scale models are usually expensive and time consuming, their application to large scale or long term problems (e.g. general river erosion) is commonly prohibitive. One-dimensional numerical models emerged early as an alternative to simulate long river reaches, usually longer than 20 channel widths (USACE 1993), when the transverse detail was not deemed too important. As computer power increased, more sophisticated two and three-dimensional models also appeared. Numerical models usually have the advantage of higher speed and lower cost when compared with physical models. For water flow only, hydrodynamic models are fairly advanced today (e.g. Sinha et al. 1998, Lane 1998, Lane and Richards 1998, Lane et al. 1999, Meselhe and Sotiropoulos 2000, Booker et al. 2001, Ma et al. 2002, Lane et al. 2002, Morvan et al. 2002, Ferguson et al. 2003). As computer power increases, so do the capabilities of numerical modeling, which ensures that numerical models will become even more popular in the future. In the last 15 years, powerful three-dimensional (3D) sediment transport models have emerged that can provide impressive results that before were possible only using physical models (e.g. Shimizu et al. 1990, Olsen and Melaeen 1993, Olsen and Kjellesvig 1998, Olsen 2003). However, numerical models, which solve mathematical equations, are only as good as the 2 equations they solve. The sediment transport equations used by morphological models are still a major source of uncertainty and the main drawback of numerical modeling when compared to physical models. Numerical models can reduce the uncertainty of the sediment transport predictors if measured data are available for calibration (i.e. adjustment of parameters in equations to fit observed data); in such cases, good quantitative results usually can be obtained (e.g. Struiksma 1985, Struiksma et al. 1985, Gessler 1999, Lopez and Falcon 1999, Nagata et al. 2000, Busnelli et al. 2001, Darby et al. 2002, Defina 2003, Due et al. 2004). Even i f reliable data are not available for calibration, numerical models still can provide useful qualitative results. For example, they can be used to assess the morphological effects of different alternatives (e.g. length of bridge, orientation of groynes, depth of dredging, shape of guide bank, etc.). For many practical problems, such qualitative results are usually sufficient to select a suitable solution. Clearly, numerical river morphology models can be invaluable tools for river engineering. Although 3D morphology models provide high resolution and better description of the flow field, their practical application is still limited because of the exorbitant computation resources they demand. In fact, the movement of sediment load over the riverbed (bedload transport) happens over the two-dimensional (2D) bottom surface. The turbulent flow over that bed surface is indeed 3D; but in certain cases, for example when the channel width is much larger than the water depth (i.e. shallow flow), the water movement in the vertical direction can be ignored and the flow assumed to be 2D in the horizontal plane. 2D shallow water models, also known as depth-averaged or vertically-averaged (VA), are at least one order of magnitude faster than 3D models and can provide complete detail of the 2D riverbed topography. 2D models can be seen as a viable compromise between demanding, high-resolution 3D models and inexpensive, low-resolution ID models. Even in cases where 3D effects are not negligible, such as the helical flow along bends, 2D models can still be applied provided that the 3D flow effects over the riverbed are modeled somehow. For example using empirical secondary flow corrections (e.g. Struiksma et al l985, 3 Shimizu and Itakura 1989, DHI 1998, Nagata and Muramoto 2000, Kassem and Chaudry 2002) or additional moment of momentum equations (Ascanio and Kennedy 1983, Yeh and Kennedy 1993b, Jin and Steffler 1993, Ghamry and Steffler 2005). Only in some special circumstances, usually involving local structures, 3D effects cannot be resolved by 2D models, and physical models or full 3D models are required (e.g. local scour at bridge piers, submerged guide vanes, etc.). However, for most applications to natural rivers, 2D models are normally suitable. In summary, morphological changes in alluvial (mobile bed) rivers can have severe effects on structures located in or near rivers, in some cases with tragic consequences for human life. The morphological riverbed changes occur because of the sediment transport over the 2D bottom surface is driven by a 3D turbulent flow field. In many practical cases, the effects of the 3D flow field over the bed can be approximated by a simpler 2D flow model. In this thesis two versions of a 2D river morphology model were developed. The V A version is a conventional 2D model with a secondary flow correction for flow along bends. The V A M version replaces the secondary flow correction with a novel approach in which moment of momentum equations are used to model the 3D flow effects. Both versions compared very well with experimental data from various small physical models, as well as a full scale meandering river. 1.2 BRIEF HISTORICAL PERSPECTIVE The numerical river morphology model developed in this thesis, as is all present knowledge of fluvial hydraulics, is rooted in the discoveries and achievements of notable hydraulic engineers and mathematicians over the last 500 years.-As a tribute to them, a short historical review of fluvial hydraulics relevant to the present research is presented below. 4 Since the early development of human civilization, our life has been intimately linked to rivers. The greatest and oldest ancient civilizations flourished in the valleys of Mesopotamia, The Nile, The Indus and the Yangtze Kiang. As early as 4000 B.C. engineers were regulating rivers in China and building irrigation canals in Mesopotamia (Graf 1971, Simons and Senturk 1977, Julien 2002). Fascination and interest for rivers has been present in most civilizations, both ancient and modern. However, steady and continuous growth in the scientific knowledge of fluvial hydraulics did not begin until the Renaissance in Italy (Graf 1971). Leonardo da Vinci (1452-1519) is credited with making hydraulics a science (Graf 1971). He carried out experimental work in a flume with a glass wall (Simons and Senturk 1977) and wrote a monumental book, Del moto e misura dell'aqua, where in he re-discovered the law of water continuity first stated by Hero of Alexandria 1500 years before him (Julien 2002). Da Vinci correctly stated that, in a channel flow, surface velocity is higher than bottom velocity because of the greater resistance exerted by the solid bottom (Rouse and Ince 1957). Da Vinci also studied surface waves and hydraulic jumps (Graf 1971) and provided some of the first known diagrams of turbulent eddies (Levi 1995). Unfortunately, most of the writings of da Vinci were unknown to his contemporaries (Graf 1971). Delia Misura dell'acqua correnti, published in 1628 by Benedetto Castelli, a pupil of Galileo Galilei, seems to be the first treatise on hydraulics (Graf 1971). Castelli made the water continuity equation available and accepted by the engineering community of his time. In 1697, Domenico Guglielmini published Delia Natura de 'Fiumi, the first treatise exclusively devoted to fluvial hydraulics (Levi 1995), recognizing the variation of velocity with depth and the importance of channel slope and depth in relation to scour and deposition. Guglielmini rejected the established idea that rivers require a minimum bed slope for water to flow, he correctly understood that is the drop in water surface, not bed, which drives the movement of water (Levi 1995). Guglielmini stated that for scour to occur "it is required that the force that scours overcome the resistance of the earth or other material that forms the riverbed" (Graf 1971, Levi 1995). Guglielmini understood that "a river will not deepen its bed to infinity" by scour, but will "reach some kind of equilibrium" because "the bottoms and the widths of the beds are determined by Nature due to combination of action and resistance causes" (Levi 1995, Rouse and Ince 1957). Guglielmini was the first 5 to make statements about loose-boundary hydraulics and is considered by many as the father of fluvial hydraulics (Graf 1971). Later, the most notable developments in fluvial hydraulics were achieved in France. In 1786, Pierre du Buat published his seminal treatise Principes d'Hydraulique, the first fundamental book on hydraulics, condensing the current knowledge of his time (Simons and Senturk 1977). Du Buat conducted the first experiments on movable beds to determine the maximum flow velocity they can resist ("critical velocity" in modern terminology). He also vividly described the movement of dunes and the armouring of beds. Some early notions about alluvial bed equilibrium (stabilite) and how nonuniform flow conditions can affect that equilibrium, i.e. by scour and deposition, are stated in his book. Referring to meandering rivers he stated that "the water does not like straight lines". The distorted velocity distribution in curved channels and the importance of local scour were also recognized (Graf 1971). Du Buat's greatest contribution to the field of sediment transport is the shear-resistance concept (Graf 1971). He stated that the "improperly called friction" F is F = ARhS (1.1) where A is the area of the bed; Rh the hydraulic radius; and S the slope of the channel. This expression is similar to the modern version of the bed shear stress r t ' Tb ~ PgRhS Q 2) where p is water density and g is gravitational acceleration. Antoine Chezy, another French engineer contemporary to du Buat, measured the velocity of floating wax balls along the Seine River in 1769 (Newbury 1995), from which he developed his famous equation for the average flow velocity U of an open channel 6 U = CJ^S (1.3) C is a flow resistance factor, nowadays known as the Chezy roughness coefficient. In 1822, Louis Marie Henri Navier presented his paper Memoire sur les lois du mouvement des fluides applying for membership to the French Academie des Sciences (which he received). There, Navier extended the partial differential equations (PDE) of ideal fluid movement developed by Euler, by considering the hypothetical attraction and repulsion between adjacent molecules (Rouse and Ince 1957), based on the Newtonian model of molecules joined by springs (Levi 1995). Although the form of the equations is correct, Navier did not understand the concept of shear stresses in a fluid. Several other researchers re-derived Navier's equations based on different assumptions. In 1829 Simeon Denis Poisson derived them by adapting the equations of equilibrium of elastic solids to flow of compressible fluids. In 1843, Jean Claude Barre de Saint-Venant derived the equations considering the internal stresses in the fluid, a more general derivation valid also for turbulent flows, and considered more acceptable based on today's knowledge. Finally, the British mathematician George Gabriel Stokes derived the equations in the manner followed today, using the concept of viscosity (Levi 1995, Rouse and Ince 1957). Today, it is customary to call the system of PDE representing the water continuity and flow momentum, which describe fluid motion, as Navier-Stokes equations for three-dimensional flow, and Saint-Venant equations for the one or two-dimensional flow (depth-averaged). Although the Navier-Stokes equations completely describe the movement of water, they are not analytically solvable, numerical methods are required for that purpose. Moreover, Navier-Stokes equations do not take into account the effects of sediment transport on the river bed dynamics; for that, another equation for sediment continuity is needed. In 1925, the Austrian researcher Felix Maria Exner became interested in the development of dunes in rivers. During his research, he derived and employed a sediment continuity equation, known today as Exner's equation. Exner's fundamental contribution is that bed elevation Zb in 7 a given reach changes in time t only i f the volumetric rate of sediment transport per unit width qs entering the reach is different from the one leaving (Parker 2005) ( 1 - ^ = - ^ (1.4) dt ox where X is the porosity of the bed material and x the longitudinal direction. The sediment transport rate qs is a function of the properties of both flow (velocity, depth, shear stress) and sediment (grain size, density, etc.), which were poorly known in Exner's time (Parker 2005). However; Exner correctly reasoned that erosion occurs i f the flow velocity increases in the downstream direction, while deposition occurs i f velocity decreases (Graf 1971, Simons and Senturk 1977). The accurate determination of the rate of sediment transport qs moving over the riverbed to be used in equation (1.4) is still an unsolved problem. The idea that bedload transport is produced by the drag force of the fluid was first proposed by Paul Francois Dominique du Boys in 1879 (Leliavsky 1955, Rouse and Ince 1957, Graf 1971). Du Boys assumed that the sediment moved in layers according to a linear velocity distribution in the vertical and arrived at qs = Xxb\?b ~iTb)cA- The critical shear stress (Tb)cr is the value needed for sediment to start moving, while x is a characteristic sediment coefficient. Although du Boys' assumption was later proved incorrect (Rouse and Ince 1957, Graf 1971), several transport equations have a similar structure. Since du Boys' early contribution, several equations have been proposed based on different assumptions, as discussed in detail by Graf (1971). For example, in 1926 Schoklitsch proposed an equation based on unit discharge instead of bed shear stress. In 1934, Meyer-Peter proposed an empirical equation based on shear stresses assuming (zb)cr = 0.047. In 1942, Einstein avoided the use of a critical shear stress and proposed an equation based on the probability that instantaneous lift forces exceed the submerged weight of the particle. In 1966, Bagnold introduced an approach that assumes sediment transport is a function of the available stream power of the flow. Unfortunately, at this time no sediment transport equation 8 is universally accepted, leaving the user with the task of selecting the equation believed to be most suitable for a given application. The river morphology model developed in this thesis is based on a Finite Element numerical scheme that solves the two-dimensional Saint-Venant hydrodynamic equations plus Exner sediment continuity equation for bedload. The purpose of the model is to dynamically compute the bed evolution in alluvial channels caused by sediment scour and deposition. Several bedload sediment transport equations have been included in the model. 1.3 THESIS O B J E C T I V E S The main objective of this thesis is the development of a numerical river morphology model with the following features: [1] Two-dimensional (2D) in the horizontal plane (depth-averaged). [2] Solves Exner's equation for bedload transport. [3] Based on Finite Element Method and unstructured meshes (triangular elements). [4] Valid for both subcritical and supercritical flows (transcritical flow capabilities). [5] Capable of simulating scour and deposition in meandering rivers (secondary flow). [6] Capable of simulating bed scour and deposition caused by forced inflow of sediment. Another objective of the thesis is to assess the applicability of the 2D vertically-averaged and moment of momentum hydrodynamic model (VAM) developed by Ghamry (1999) to predict scour and deposition in meandering rivers. The 2 D - V A M is a quasi-3D flow model that has been successfully applied to simulate flow along fixed-bed curved channels (Ghamry and Steffler 2005), but it has not been applied to movable beds before. The hypothesis that the 9 V A M model does not require a semi-empirical secondary flow correction for bedload transport, as needed by conventional 2D models, was investigated and verified. As will be shown later in the literature review, at this time there is no 2D Finite Element (FE) river morphology model with all the 6 features listed above. 2D FE morphology models for meandering rivers are still uncommon, as most models are based on Finite Differences (FD) or Finite Volumes (FV) schemes. Existing 2D FE morphology models, such as RMA2-SED2D (King 2000) and TELEMAC-SUBIEF (Moulin and Slama 1999), are commonly intended for suspended sediment. For transcritical flow, most existing morphology models are one-dimensional (ID). A 2D V A M morphology model has not been developed previously. Some brief comments on the features of the model [numbers in brackets] and their relevance for river morphology follow below. The morphological simulation of meandering rivers [5], demands a 2D model [1] with bedload transport capabilities [2]. As du Buat noticed more than 300 years ago, the velocity distribution along bends has a strong lateral variability (e.g. Rozovskii 1957, Johanneson and Parker 1989b, Termini et al. 2001, Blanckaert and Graf 2004, Ghamry and Steffler 2005) that a ID model is unable to capture. Along alluvial bends, scour occurs along the outer banks and deposition along the inner banks, caused by the deviation of bedload transport direction induced by flow curvature (helical motion). Modeling of bedload transport [2] is vital to reproduce features such as point bars and pools found in meandering rivers. The use of unstructured or irregular meshes [3] to represent the planform geometry of rivers, as opposed to the structured meshes used by the vast majority of river models, is an advantage when dealing with complex geometries and strong flow gradients (Figure 1). Unstructured meshes are very flexible to accommodate practically any planform geometry, they can dramatically vary the resolution over the computational domain (Waddle and Steffler 2002). 10 Increasing the resolution with a finer mesh in areas of strong flow gradients, such as sharp hydraulic jumps, is a clear advantage when simulating transcritical flow [4]. Figure 1. Examples of (a) structured curvilinear mesh based on quadrilateral elements and (b) unstructured mesh based on triangular elements. 11 Popular 2D flow models (e.g. M I K E 21, R M A 2 , FESWMS) are intended for lowland rivers and are mainly applicable for estuarine or large riverine environments with low gradients (Papanicolaou et al. 2005). Morphology models for steep mountain rivers are mostly limited to ID flow (e.g. Lopez and Falcon 1999, Busnelli et al. 2001, Cui et al. 2003b, Papanicolaou et al. 2005), because of the challenges imposed by transcritical flow [4] in steep mountain rivers. Therefore, a single 2D river morphology model suitable for both lowland and mountain rivers (which implies [4]) is highly desirable for practical applications. In some cases, bed degradation (scour) or aggradation (deposition) is caused by an imposed or forced upstream inflow of sediment. A typical example is bed degradation below a dam caused by reduced or negligible inflow of sediment in the downstream reach. The model has the ability to take this effect into account [6]. 1.4 L I T E R A T U R E R E V I E W 1.4.1 Empirical and analytical models 1.4.1.1 Empirical models Until Exner developed the sediment continuity equation in the mid 1920's, morphodynamic modeling did not exist (Parker 2005). For centuries, river morphology was a descriptive science. A good example is the relationships between thalweg and planform geometry of meandering rivers developed by Louis Fargue in the late 1800's (Martin-Vide 2002, Hager 2003). Fargue was the first to use physical models of rivers in 1875 (Rouse and Ince 1957, Graf 1971, Hager 2003); his rules on the geometry of meandering rivers have been verified in some rivers, but they did not depend on any hydraulic parameter such as discharge or Froude number (Martin-Vide 2002, Hager 2003). One of the first attempts to provide quantitative expressions for river morphology was made by Robert Kennedy in 1895 from his analysis of irrigation canals in India (Rouse and Ince 1957, Graf 1971), with the so-called regime 12 (equilibrium) theory to determine the width, depth and slope of alluvial channels (Blench 1969, Graf 1971). Although such empirical expressions are still in use today to provide initial estimates of channel geometry (e.g. average depth), they do not provide enough detail for many practical applications (e.g. maximum scour depth in a bend). Moreover, empirical equations are usually derived from limited data that cover a narrow range of parameters (Blench 1969). According to Graf (1971) the application of empirical equations (e.g. regime equations) outside the range of the original data is "questionable and often dangerous". The motivation for empirical equations derived from observations of alluvial canals in the late 1800's and early 1900's can be easily understood from the words of Blench (1969): "Quantitative analysis and explanation of river behaviour... require formulas expressing the control of dynamic laws. Fluid dynamics cannot deduce these laws rigorously at present so they must be obtained, as closely as possible, from observation". Fortunately, in the last half of the past century considerable progress has been made in theoretical models based on sound physical laws of fluid and sediment motion, as discussed below. • 1.4.1.2 Analytical models The pioneer work of Exner on the morphology of alluvial beds deserves special attention. Exner was not only the first to derive the sediment continuity equation (Parker 2005), but he was also the first to show analytically both the convective (translation) and diffusive (dispersive) nature of his new equation. Exner proved that for alluvial channels of constant width it is possible to have migrating (convective) bedforms; his analytical model predicts a shape of migrating bedforms very similar to the observed shape of dunes (Leliavsky 1955, Graf 1971). But he also demonstrated that, for channels of varying width, the bedforms do not migrate (Leliavsky 1955, Graf 1971). These fixed bedforms caused by varying width have been observed in experiments reported by Tsujimoto (1987) and Pittaluga et al. (2000). Exner's pioneer ideas resemble, at a very crude level at least, the theory of "free" and "forced" bars of Seminara and Tubino (1989) and Tubino and Seminara (1990), which states that free migrating alternate bars become fixed when forced by continuous changes in channel width. It is 13 interesting to notice that none of these modern researchers (Tsujimoto 1987, Pittaluga et al. 2000, Seminara and Tubino 1989) makes any reference to Exner's previous model. It might be argued that Exner's model is not very accurate because of the simplifying assumptions he made, and that his model is intended for dunes. But he deserves the credit for being the first to provide sound analytical models for alluvial bed morphology. Neglecting friction and assuming that sediment transport depends linearly only on velocity by means of a erosion coefficient ag q s =aEU (1.5) (1 -A) Exner finds for a channel of constant width B, constant discharge Q and constant depth h (Leliavsky 1955, Graf 1971, Graf 1998) ^ = - c ^ ; c = -^4 (1.6) dt dx Bh2 Equation (1.6) is clearly convective; it represents the translation of a wave moving at celerity c. Equation (1.6) was used by Exner to model the translation of dunes (Leliavsky 1955, Graf 1971). Later, he included the effects of friction, a simplified momentum equation and other assumptions to arrive at (Leliavsky 1955) dt dx2 b dxdt dt2 D and E are constants. The second term in equation (1.7) is a diffusion term. The fact that Exner showed 80 years ago both the convective (equation 1.6) and diffusive (equation 1.7) nature of the sediment continuity equation is remarkable because even today there is an 14 ongoing debate on this issue (Lisle at al. 2001, Cao and Carling 2003, Cui et al. 2005, Cao and Carling 2005). According to Jansen et al. (1979) and Przedwojski et al. (1995) assuming that sediment transport is a simple power function of velocity qs=aUb (1.8) the following analytical models of bed movements in alluvial rivers can be proposed. Simple wave equation ( , - ^ + c * t . O , c-W*$>U (1.9) dt dx l-F Parabolic model dzh , d2zh ba b , k = — ^ (1.10) dt dx2 3 5 0 ( 1 - ^ ) Hyperbolic model dzb d2zb k d2zb n —t-k—f ^ = 0 (1.11) dt dx2 c dxdt y ' 15 where q =Uh = Q/B is the discharge per unit width or unit discharge; F is the Froude number; So is the initial bed slope of the channel; and k is the diffusion coefficient. Recently, Lisle at al. (2001) have proposed the following model dzb _ K 82zbfd (1.12) dt (1-/1) dx2 {dx K is an empirical constant that depends on the sediment transport equation used (Lisle et al. 2001 used Meyer-Peter Muller). The unspecified terms at the right hand side are unsteady flow terms which are small for F < 1 (Lisle et al. 2001). Notice the resemblance between equations (1.6) and (1.9) and also (1.7) with (1.11) and (1.12). Although some of Exner's assumption to derive his equations (1.6) and (1.7) might be criticized (e.g. equation 1.5), the form of the equations is correct and the roles of convection and diffusion are evident from the equations themselves. The parabolic model (eq. 1.10), also known as diffusion model, has been successfully applied to problems of alluvial riverbed aggradation and degradation by several researchers (Jansen et al. 1979, Soni et al. 1980, Gi l l 1980, Jain 1981, Gil l 1983a, b, Jaramillo and Jain 1984, Lu and Shen 1986, Gi l l 1987). Since equation (1.10) is parabolic, it is amenable to analytical solution for simple straight channels. For example, in the case of an aggrading channel due to overloading of the supply of sediment (Figure 2), subject to the initial and boundary conditions: zb(x,0) = 0 zb(0,t) = Ah(t) lim zb(x,t) = 0 16 The solution is (Graf 1998): zb(x,t) = Ah(t)-erfc f x ^ 2-Jkt. (1.13) where erfc is the standard error function. Equation (1.13) compares very well with the experimental data of Soni et al. (1980), proving that diffusion models can be suitable for some morphological applications, despite the criticism of Cao and Carling (2003, 2005). Figure 2. Aggradation by overloading of the supply of solid discharge Aqs In summary, analytical models show that the sediment continuity equation is in general a hyperbolic PDE with two main modes: convection and diffusion. Both modes of transport are present, but in certain cases, one may dominate over the other. Convection (advection) may dominate when migrating bedforms are present (e.g. dunes, alternate free bars, sediment waves, prograding deltas); while diffusion (dispersion) seems to be dominant in aggradation and degradation of straight channels and forced bars caused by varying width or channel curvature. 17 1.4.2 Morphodynamic models in 1,2 and 3 dimensions Analytical models provide valuable insight about river morphology. However, because they are based on simplifying assumptions about the flow and sediment equations, their application is usually limited to simple cases (e.g. low Froude number, constant discharge, and constant width) that may not be adequate for solving some practical engineering problems. This is because analytical solutions are available only for some types of relatively simple PDE's (e.g. parabolic model). For the general PDE's describing flow and sediment motion, only numerical solutions are possible. Typically, numerical models are classified as ID, 2D or 3D (Table 1). 1.4.2.1 ID models ID morphology models are based on the assumption that water and sediment follow only the longitudinal -streamwise- direction. The river is discretized in cross-sections; velocity, depth and bed elevations are computed as averages in those cross sections. The main advantage of these models is their fast computational speed, which makes them suitable for long-term simulations of relatively long reaches usually spanning several river widths (USACE 1993). Currently, several ID models exists with slightly different capabilities, such as HEC-6 (USACE 1991), F L U V I A L 12 (Chang 1988), M I K E 11 (DHI 2000), SEDROUT (Ferguson et al. 2001) and BRI-STARS (Molinas 2000). The main limitations of ID models is that they neglect the transverse or lateral changes across the channel, which are important in many practical case; e.g. scour and deposition along bends, alternate bars, braiding. Additionally, the sediment transport equations are highly non-linear functions of flow velocity or shear stress, and cross sectional averages of the flow field may provide unrealistic estimates of sediment transport (Ferguson 2003). Some ID models try to compensate for their inherent limitation using different strategies. F L U V I A L 12 uses simple secondary flow models and approximate equations to compute the 18 transverse bed slope in curved channels (Chang 1988). BRI-STARS (Molinas 2000) divides the modeled channel into several stream tubes (which is somehow similar to having several ID models in parallel) to obtain additional information across the channel. Both, BRI-STARS and F L U V I A L 12 can also compute the changes in channel width by bank erosion using extremal hypotheses (Chang 1988, Molinas 2000). M I K E 11 can be set as a network of ID channels, instead of a single ID thread, to provide additional information away from the main channel (Willems et al. 2002). 1.4.2.2 2D models 2D models in the horizontal plane depth-average the flow velocity, neglecting the vertical velocity component, which is suitable for shallow and wide rivers. 2D models have the advantage of providing information of the transverse flow field. However, in cases of highly 3D flow, such as the helical (spiral) motion in bends caused by the centrifugal acceleration (or centripetal depending on the frame of reference), 2D models require external secondary flow correction sub-models to simulate scour and deposition in bends and meandering rivers. . Several 2D models currently exist as shown in Table 1 and reviewed by Przedwojski et al. (1995) and more recently by Due et al. (2004). Many 2D models are commonly used by the engineering community as commercial packages. As discussed by Mosselman (2004) and Kassem and Chaudhry (2004), some of the information regarding commercial models is available in web sites and internal reports, as opposed to the peer-reviewed publications of research models. As a comparison cannot be made on the same ground, the following review of 2D modes has been divided into two categories: research and commercial. 19 Table 1. Examples of some 2D and 3D numerical models for alluvial channel morphology. Model Dim. Sche-me Bedload / suspended (Name) Application / Comments Koch and Folkstra (1981) 2D FD Yes / No Circular 180° and 320° bends. Struiksma(1985) 2D FD Yes / No The meandering Waal River. Struiksma et al. (1985) 2D FD Yes / No Circular lab channels. Shimizu and Itakura (1989) 2D FD Yes / No Meander loops; alternate and braided bars. Shimizu et al. (1990) 3D FD Yes / Yes Sine generated curve, natural bend. Spasojevic and Holly (1990) 2D FD Yes / Yes (MOBED2) Multiple size fractions DHI (1998) 2D FD Yes / Yes (MIKE 21C) Jamuna River. Moulin and Slama (1999) 2D FE No / Yes (TELEMAC-SUBIEF) River intake. Gessler at al. (1999) 3D FD Yes / Yes (CH3D-SED) Mississippi River. Jia and Wang (1999) 2D FE Yes / No (CCHE2D) Meandering channels. King (2000) 2D FE No / Yes (SED2D-WES) suspended sediment. Nagata et al. (2000) 2D FD Yes / No Bank erosion and migration Wu et al. (2000) 3D FV Yes / Yes (FAST3D) 180° bend, NET. Duan etal. (2001) 2D FE Yes /Yes (enCCHE2D) Meandering development. Brethour (2001) 3D FV No / Yes (FLOW3D) local scour of jet over sand Zhang and Wang (2001) 2D FE No / Yes Dyke break /transcritical flow, suspended NET. Froelich (2002) 2D FE Yes / No (FESWMS) unknown. Kassem and Chaudhry (2002) 2D FD Yes/No Circular bends. Darby et al. (2002) 2D FD Yes/No (RIPA) bank erosion and migration Kusakabe et al. (2003) 2D FD Yes / No Transcritical flow in expansion Olsen (2003) 3D FV Yes / Yes (SSIIM) Meandering from initial straight flume. Defina (2003) 2D FE Yes / No Alternate bars growth in straight flume. Due et al. (2004) 2D FV Yes / Yes Laboratory channels / NET Hervouet and Villaret (2004) 2D FE Yes / Yes (TELEMAC-SYSIPHE) groyne, trench. Wu (2004) 2D FV Yes / Yes East Fork River / NET Wu et al. (2005) 2D FV Yes / Yes Vegetation and NET in bends ID: One-dimensional; 2D: two-dimensional; 3D: three-dimensional. FE: Finite Elements; FV: Finite Volumes; FD: Finite Differences. NET: Non-equilibrium transport 20 Commercial 2D models M I K E 21C developed by the Danish Hydraulics Institute (DHI 1998) is probably one of the most sophisticated 2D morphology models available, according to Mosselman (2004). M I K E 21C is based on a structured curvilinear FD grid to simulate bedload and suspended sediment; as well as scour and deposition along bends. It is limited to low Froude numbers and sand bed rivers, typical of lowland regions. Recently, L i and Millar (2004) have incorporated an additional sub-model for gravel transport to simulate the Fraser River gravel reach. The 2D FE hydrodynamic model T E L E M A C - 2 D uses triangular elements and can be linked to independent sediment transport modules. Moulin and Slama (1999) used the module SUBIEF to simulate suspended sediment transport and the bed changes in a water intake. Recently, Hervouet and Villaret (2004) describe applications with the module SYSIPHE to simulate both bedload and suspended sediment in straight hypothetical flumes. From Moulin and Slama (1999) and Hervouet and Villaret (2004) it is unclear i f these modules can also be applied for alluvial bends or meandering rivers; i.e. i f they considered secondary flow and transverse bed slope corrections. The model SED2D (King 2000) can compute suspended sediment transport using the flow field computed by the popular FE element model RMA2. Since the model does not compute bedload, it cannot be applied for scour and deposition in bends. The popular hydrodynamic model FESWMS developed by the US Federal Highway Administration claims to have capabilities for sediment transport modeling in the manual of its latest release 3 (Froelich 2002). Unfortunately, there are no applications shown to verify that claim. 21 Research 2D models Some of the more interesting and cited research on alluvial bend morphology, both experimental and numerical, has been developed at Delft Hydraulics Institute in the Netherlands. Koch and Folkstra (1981) applied a 2D model to curved alluvial flumes of constant curvature (circular bends), with limited success because they used a rather simple flow model. The improved model of Struiksma et al. (1985) was more successful at reproducing the scour and deposition measured in curved flumes. Struiksma (1985) successfully reproduced the observed patterns of scour and deposition along a 10-km meandering reach of the Waal River (a branch of the Rhine River in The Netherlands). A l l these models have been based on structured FD grids. Some of the experimental data reported by these researchers (Koch and Folkstra 1981, Struiksma et al. 1985, Struiksma 1985) have been used to test the model developed in this thesis (see Chapters 3 and 4). Most of the research 2D morphological models are based on structured grids with quadrilateral elements (Fig. 1). Spasojevic and Holly (1990) developed the model MOBED2 with capabilities for multisize transport, applicable to gravel bed rivers. Nagata and Muramoto (2000) simulated the widening of a channel caused by bank erosion. Duan et al. (2001) used the Enhanced CCHE2D to simulate the evolution and migration of single thread meandering channels (CCHE2D uses quadrilateral FE, Jia and Wang 1999). Kassem and Chaudhry (2002) applied a boundary fitted FD model to simulate some of the alluvial bend experiments of Koch and Folkstra (1981) and Struiksma et al. (1985). Darby et al. (2002) simulated large amplitude meanders and bank erosion using a modified version of Delft's model RIPA. Applications using 2D FE models based on unstructured triangular elements (Fig. 1), as the one used in this thesis, are still uncommon. Defina (2003) reproduced the recent experiments of Lanzoni (2000) of alternate bar formation in straight channels from an initial flat bed. Termini (2002) simulated the bed topography in large meandering channels. However, since Termini's (2002) model considers only the convective acceleration of the flow, neglecting the 22 effects of both secondary flow and transverse slope on bedload transport, her model is probably only applicable to wide channels with mild curvature. Recently, research on non-equilibrium transport (NET) for unsteady flow modeling is receiving considerable attention (Nagata and Muramoto 2000, Zhang and Wang 2001, Due et al. 2004, Wu 2004, Cao et al. 2004, Wu et al. 2005). Conventional equilibrium bedload models assume that sediment transport rate corresponds to the transport capacity of the flow (Due et al. 2004), which may not be realistic for highly unsteady flows. 1.4.2.3 3D models 3D morphological models solve the full set of Navier-Stokes equations and therefore provide a high resolution of the flow field. Although the application of 3D models for river hydro-dynamics is now well stablished (e.g. Lane 1998, Sinha et al. 1998, Lane et al. 1999, Meselhe and Sotiropoulos 2000, Booker et al. 2001, Ma et al. 2002, Morvan et al. 2002, Lane at al. 2002, Ferguson et al. 2003); morphological applications are still uncommon (e.g. Shimizu et al. 1990, Gessler et al. 1999, Wu et al. 2000, Brethour 2001, Olsen 2003) because of the excessive computational time that they demand. Table 1 shows some of the 3D models reviewed here. Shimizu et al. (1990) successfully modelled'the bed evolution in sine-generated meanders using a 3D model, comparing also the results with the previous 2D model of Shimizu and Itakura (1989). Shimizu et al. (1990) found that 3D models produced better results than 2D models since the latter neglects the momentum transport by the secondary flow, necessary for the precise prediction of the flow field. Gessler et al. (1999) used the model CH3D-SED developed by the US Army Corps of Engineers to simulate the erosion and deposition in a reach of the lower Mississippi River. Wu et al. (2000) applied the model FAST3D to simulate scour and deposition along a 180° bend using a NET formulation. Brethour (2001) modified FLOW3D to simulate the transport of suspended sediment. One of the most impressive morphological applications is the 3D generation of a meandering channel from an initially 23 straight channel made by the model SSIIM (Olsen 2003). SSIIM uses a F V unstructured adaptive grid that changes as the channel planform evolves. Notably, SSIIM is not limited to single thread channels as previous 2D models (Nagata and Muramoto 2000, Duan et al. 2001, Darby et al. 2002); therefore, it can simulate complex processes such as meander cut-off and the initiation of braiding. 1.4.3 Unstructured meshes and Finite Element hydrodynamic models As shown in Table 1, 2D morphology models based on unstructured meshes and FE schemes are not as common as models based on FD or FV. Structured meshes based on quadrilateral elements have long dominated the field because of their relative computational simplicity. In structured meshes, the elements are commonly located along "rows" in both the longitudinal and transverse directions, as shown in Figure 1 .a. A simple indexing scheme (e.g. / in the longitudinal direction and j in the transverse direction) can be used to easily identify an individual element, and all the elements around it, as shown in Figure 3. For the structured mesh shown in Figure 3.a, given element (3,2) it is straightforward to find the upstream element (2,2). For unstructured meshes, a simple indexing is not possible. Elements do not follow any particular orientation. Two sets of data are needed to define unstructured meshes: a connectivity table of every element with its nodes (e.g. element (3) in Figure 3.b is attached to nodes 3, 5 and 7), plus a table with the spatial coordinates (x, y, z) of every node. (1,3) (2,3) (3,3) (1,2) (2,2) (3.2) (1,1) (2,1) (3,1) \6 7 / / {X) ( 2 ) \ 5 • (3) ( 5 ) \ ( 4 ) \ (a) Structured mesh (b) Unstructured mesh Figure 3. Indexing of elements in (a) structured mesh and (b) unstructured mesh. 24 FD models are restricted to structured meshes only. Numerical algorithms based on structured meshes are relatively simple, easier to program and faster to execute (Chung 2002), which explains the popularity of FD models. However, structured meshes lack flexibility for complex geometries. Additionally, a factor usually overlooked is that structured meshes may be quite complicated to generate for multiple-thread rivers. For some hydrodynamic simulations, the time required to generate a structured mesh might be greater than the computational time itself. Unstructured meshes based on triangular elements are extremely flexible to map any planform geometry (Fig. Lb), although at the expense of more complex algorithms and larger computational times. An advantage is that meshes based on triangular elements are very easy to generate. In contrast with FD schemes, FE models are almost exclusively based on unstructured meshes. F V schemes can be expressed on either structured or unstructured meshes, although unstructured F V are less common. There are some examples of FV models with unstructured triangular elements. Mohammadian et al. (2003) used an unstructured 2D FV morphology model to compute suspended sediment transport around a groyne by solving the convection-diffusion equation. Pena et al. (2002) used a 2D FV morphodynamic model based on vertex-centered triangular elements. The hydrodynamic model of Pena et al. (2002) has transcritical flow capabilities, verified with a dam-break test, although the morphological application reported was in subcritical flow. Below is a review of some FE hydrodynamic models with special emphasis on the research performed at the University of Alberta, which is the source of the FE hydrodynamic code used in this research. FE methods were developed in the mid 1950's for structural engineering applications (Chung 2002), a field where today they clearly dominate over any other method. Despite the great 25 success of conventional Galerkin FE methods (GFEM) when applied to solve the elliptic or parabolic PDE's of structural analysis; the initial applications of G F E M to solve the hyperbolic PDE's of Computational Fluid Dynamics (CFD) were not very successful. In the 1970's the general believe was that existing FD techniques were superior to FE for flow computations (Katapodes 1984a), supported by the formal demonstration of Dupont (1973) that G F E M do not achieve optimum accuracy for hyperbolic problems, as they do for parabolic and elliptic ones. G F E M are dispersive and nondissipative leading to spurious oscillation in the vicinity of hydraulic jumps (Katapodes 1984a). In comparison, the dissipative "upwinding" properties of some Finite Difference (FD) schemes help them damp those spurious oscillations (Appendix F). In conclusion, G F E M are inappropriate for hyperbolic or convection-dominated applications (Donea and Huerta 2003). To overcome the limitations with hyperbolic problems, the so-called Petrov-Galerkin methods were developed. The pioneer work of Dendy (1976) introduced a Dissipative Galerkin (DG) method based on new upwind weighted test functions. Katapodes (1984a) introduced a ID DG scheme for open channel flow that successfully damped the spurious oscillations around hydraulic jumps. Katapodes (1984b) applied a D G scheme to 2D hydraulic jumps with very good results. The D G scheme was also found to perform better than some FD schemes, such as Lax scheme for kinematic waves and Preissman scheme for dynamic waves (Katapodes 1984c). Hicks and Steffler (1992) improved the D G scheme by reconsidering the role of the characteristics in the determination of the upwind weighting and the use of the conservation form of the equations leading to a new Streamline-Upwind Petrov Galerkin (SUPG) scheme entitled the Characteristic-Dissipative-Galerkin (CDG) method. Since then, the C D G approach has been applied successfully to several flow conditions: 2D dam breaks over initially dry beds (Ghanem et al. 1995); ID flows over vertically curved beds (Khan and Steffler 1996a) and overfalls (Khan and Steffler 1996b); ID hydraulic jumps (Khan and Steffler 1996c); ID dam breaks in channels with varying width (Hicks et al. 1997); 2D flow in braided rivers (Christison et al. 1999); ID flow routing in the Peace River (Peters and Prowse 2001); 2D flow in the wandering Fraser River gravel reach (Yusuf 2001); 2D flow and recirculation zones in 26 lateral diversions (Vasquez 2005); among others. The University of Alberta freely distributes a 2D FE hydrodynamic model called River2D (Steffler and Blackburn 2002) based on the CDG approach (available at www.river2d.ca). River2D has unique capabilities for transcritical flow and wet/dry elements (Blackburn and Steffler 2001, Steffler and Blackburn 2002). For handling wetting and drying the model uses a simple groundwater model that is activated when water depth drops below a minimum value (Ghanem et al. 1995). Using this method allows elements to be partially wet or dry, as opposed to other models that turn-off the element completely when one or more nodes gets dry (e.g. R M A 2 , M I K E 21C, FESWMS). River2D can be considered as one of the most advanced 2D FE hydrodynamic models currently available. 1.4.4 Curvature effects and V A M model Pittaluga et al. (2000) stated that the ability of 2D morphological models to replicate the observed local topography of single thread channels is largely dependent on correcting the Saint Venant equations to account for secondary flows. Furthermore, they concluded that only through the inclusion of curvature effects can 2D models reproduce the transverse deformation of bed profile and flow characteristics. Similar conclusions were obtained previously by Kalkwijk and De Vriend (1980) and Struiskma et al. (1985). A characteristic feature of flow in bends is a strong 3D helical flow (also known as secondary current or transverse circulation) due to the difference of centrifugal forces between the upper and lower layers of the flow. This helical flow plays a significant role in the development of the transverse profile of bed and bottom shear stress. Because flow velocity usually is largest at or near the free water surface, centrifugal acceleration drives the fluid in the upper portion of the water column toward the outer bank. To fulfill continuity, fluid at the lower level moves towards the inner bank and may transport bed sediment. The exchange of momentum between 27 the primary (streamwise) and secondary (transverse) flows produces a redistribution of velocity, both across the width and over the vertical (Rozovskii 1957, Johannesson and Parker 1989b, Blanckaert and De Vriend 2003, Blanckaert and Graf 2004). Across the width, velocity increase along the outer bank of the bend and decreases along the inner bank; while at the same time the vertical velocity profile becomes more uniform (flatter). To simulate average bed variations in alluvial channel bends accurately, it is necessary to replicate the effects of spiral motion on the flow and bed-particle models (see Chapters 3 and 4 for details). A 3D mathematical simulation of natural meanders including the complexities mentioned previously is not yet practical because of exorbitant computational effort, complexity of numerical computations, and necessity of processing a large amount of computed results (Kassem and Chaudhry 1998). Therefore, approximations usually are made to obtain computational solutions. One such approximation is to use 2D depth-averaged water flow equations. Many investigators have attempted to predict the bed deformations in river bends analytically and numerically. A few analytical models have been developed; e.g., Kikkawa et al. (1976), Ascanio and Kennedy (1983), Odgaard (1986a,b) and Yeh and Kennedy (1993a,b). Such models are based on the balance of the dominant forces acting on a sediment particle moving along a radially inclined bed, usually applying the concept of conservation of moment-of-momentum (MOM). The forces are fluid drag and particle submerged weight. When these forces become equal, an equilibrium transverse bed slope is achieved. Due to the simplifying assumptions, • the analytical models are applicable only to relatively simple bend-flow conditions (e.g. weak curvature, constant radius of curvature). A number of 2D depth-averaged models have been developed for computing bed deformation in meandering channels (e.g. Struiksma et al.1985, Shimizu and Itakura 1989, Jia and Wang 1999, Nagata and Muramoto 2000, Kassem and Chaudry 2002), as discussed previously in section 1.4.2.2. Unfortunately, the vertical structure of the flow is lost in the depth-averaging process making it impossible for a 2D model to reproduce the effects of helical motion without further corrections. The common approach for depth-averaged morphological models is to employ the so-called secondary flow corrections (SFC) in an attempt to take into account the effects of secondary flow on the direction of bed shear stress. The angle Ss that the bed shear 28 stress vector makes with the depth-averaged flow direction is assumed to be a linear function of the flow curvature ratio h/rc, where rc is the local radius of curvature of the streamlines (Struiksma et all985) tanSs = A(h/rc) (1.14) A = 5 ~ 13, is a parameter that depends on friction. The linear SFC given by equation (1.14) neglects the momentum exchange between primary and secondary flows and it is strictly applicable only for weak curvatures (low values of h/rc). For strong curvatures, linear SFC tend to overestimate the effects of the secondary flow (Blanckaert and De Vriend 2003, Blanckaert and Graf 2004). An alternative to recover some information of the vertical structure of the velocity is to assume that the vertical profile follows a simple pre-determined profile such as linear or parabolic. Figure 4 shows an example of an assumed linear velocity profile in a ID flow, the depth-averaged velocity is u0 at mid-depth with velocity varying linearly from (u0 - u\) at the bottom to (u0 + u\) at the surface. Despite the non-physical slip velocity at the bottom, the linear profile assumption has the advantage that the ratio 2u//h can be used as an estimate of the vertical gradient of velocity for non-uniform flow (when flow accelerates/decelerates although depth-averaged velocity u0 may be the same, bed shear stress is different because of the change in velocity gradient, Elgamal 2002). The additional unknown u\ of the linear profile requires an additional equation. Instead of using empirical relationships (e.g. Rozovskii 1957), a physically-based equation derived from the conservation of moment of momentum (MOM) can be used (Jin and Steffler 1993). For 2D flow, exactly the same assumption can be made for the velocity in the (x, y) plane with corresponding depth-averaged velocities (u0, v0) and moment components (u\, v\). Although rather simple, such 2D vertically averaged model with moment equations (VAM) can reproduce at the simplest level the helical motion in bends. A short review of the history of the M O M approach and V A M model follows. 29 U0 Water surface Depth-averaged velocity Linear velocity profile Real velocity profile Bed Figure 4. Sketched of assumed linear vertical velocity profile for ID V A M model (flow moving in the x direction; not to scale). Ascanio and Kennedy (1983) applied the conservation of flux of M O M to obtain analytical relations for the vertical distributions of radial-plane velocity and radial shear stress, as well as the transverse bed profile, in a wide, erodible bed channel of constant centerline radius. The results were in good conformity with measured data. The reasoning of Ascanio and Kennedy (1983) was that in the spiral flow of a bend three principal radial forces are exerted on the flow: pressure, shear and centrifugal; and since two of them -shear and pressure- are unknown, the conservation of M O M must be introduced to close the system. They argue that the analysis is similar to that adopted in equations expressing balances of forces and moments to calculate supporting forces on a loaded, simply supported beam. One of the two unknown forces does not appear in the formulation of moments about one of the supports, and therefore the other can be calculated directly (Ascanio and Kennedy 1983). 30 Yeh and Kennedy (1993a) formulated a nonuniform, fixed-boundary, channel-bend flow as two coupled sets of equations, each comprising two simultaneous equations, derived from the cross-section-averaged equations (in circular cylindrical coordinates) of conservation of flux of M O M and the depth-averaged momentum and continuity equations. Bed shear stresses were related to the primary-flow shear stress and the primary and secondary velocities. The equations were solved numerically and the results found to be in satisfactory agreement with experimental data. Yeh and Kennedy stated that the moment formulation helped them elucidate the interplay among the secondary and primary flows, notably the observed flattening of primary velocity profiles and the radial (lateral) redistribution of the depth-averaged primary velocity. The fixed-boundary model was extended to erodible-bend flows (Yeh and Kennedy 1993b). Bed topography was computed from two-dimensional (2D) sediment continuity. The resulting two coupled sets of simultaneous equations were solved numerically. The model yielded accurate predictions of the velocities (primary and secondary), and 2D bed topography. The average radial bed slope and secondary rotational and translational velocities were observed to exhibit much larger variations along the channel than in fixed-boundary bends. Steffler and Jin (1993) argued that traditional depth-averaged equations that sacrifice details over the vertical dimension are valid only for horizontal length scales much larger than the flow depth. To incorporate more vertical detail into the Saint Venant equations, Steffler and Jin (1993) introduced the ID formulation of vertically averaged and moment (VAM) equations by taking moments around the mid-depth; assuming a linear vertical velocity profile (Fig. 4) and a parabolic non-hydrostatic pressure distribution. The ID V A M equations have been applied to model flow over curved beds (Khan and Steffler 1996a), overfalls (Khan and Steffler 1996b), and suspended sediment (Guo and Jin 1999). Jin and Steffler (1993) introduced a 2D depth-averaged model formulated in the Cartesian coordinate system to calculate the effects of the secondary flow on the depth-averaged velocities (velocity redistribution) in fixed-bed curved open channels. The mathematical model consisted of the depth-averaged continuity equation, the 2 momentum equations, and 2 M O M equations for closure purposes. The numerical analysis predicted satisfactory depth-averaged longitudinal and transverse velocities as well as reasonable secondary flows. The previous 31 studies of Ascanio and Kennedy (1983) and Yeh and Kennedy (1993a) used a cylindrical coordinate system that requires the radius of curvature as input, which is difficult to determine in complex natural rivers. In contrast, the model of Jin and Steffler (1993) implicitly accounts for the effects of curvature and does not require additional information about the radius of curvature. Ghamry (1999) extended to two dimensions the ID V A M formulation introduced by Steffler and Jin (1993). The V A M equations are derived by vertically averaging or integrating the 3D Reynolds equations after multiplying them by the vertical coordinate. The derivation is equivalent to the first moment about the mid-depth. Linear and quadratic distribution shapes are assumed for the horizontal velocity components, while quadratic distributions are considered for the vertical velocity and non-hydrostatic pressure, leading to a total of 10 equations. The derived equations could be represented as a quasi-3D model with more vertical detail than the conventional Saint Venant vertically-averaged (VA) 3-equation formulation (one continuity and two momentum equations). Ghamry and Steffler (2005) compared the V A M and V A models for flow in curved channels and verified that the V A M produced better predictions of the transverse velocity redistribution along the bends. Ghamry (1999) also proposed a simpler V A M 5-equation model assuming hydrostatic pressure and neglecting vertical velocity, i.e. it solves for h, u0, v 0, w, and v, (Fig. 4). Ghamry and Steffler (2003) demonstrated that the V A M 5-equation model produced results not significantly different than the V A M 10-equation model for flow in fixed bed curved channels. Ghamry and Steffler (2003) recommend the use of the V A M 5-equation model for cases involving secondary flows. The increment of computational time for the 5 V A M equation is about 40 to 50% more than the conventional V A 3-equation model, while the memory demand increases by 2.5 to 2.8 times. Other applications of the 2D V A M model are rapidly varied flows in a contraction-expansion and over a hemispherical bump (Ghamry and Steffler 2002); flow over bedforms (Elgamal 2002); and dispersion and mixing in rivers (Alberts et al. 2002). Elgamal (2002) demonstrated 32 that while the V A model failed to predict the longitudinal variation of bed shear stress over a field of dunes; the V A M model was able to reproduce very well the observed trends, because of the additional information of the vertical profile of both velocity and non-hydrostatic pressure. Kruger and Olsen (2001) developed their own 2D "extended" shallow water (ESW) equations derived from the ID model of Steffler and Jin (1993). This ESW model is similar to Ghamry's (1999) V A M 5-equation model. Kruger and Olsen (2001) compared a conventional V A model, the ESW and the 3D model SSIIM (Olsen 2003) to simulate an experiment of a standing shock wave in a channel contraction with Froude numbers of 4.0 and 6.0. The results of the V A model were poor; however, both the ESW and SSIIM led to good agreement with measured data. No morphological application of the 2D V A M or ESW hydrodynamic models has been found in the reviewed literature. 1.4.5 Transcritical flow morphodynamics As reported by Papanicolau et al. (2005) most morphodynamic models are intended for lowland river applications with small Froude numbers (subcritical flow). Morphodynamic models with capabilities for very steep (mountain) streams demand numerical schemes that can handle both subcritical and supercritical flows and hydraulic jumps (transcritical flow). Although there are some ID models intended for the morphodynamics of mountain streams, they do not always have built-in capabilities for transcritical flow (e.g. Lopez and Falcon 1999). Busnelli et al. (2001) performed a very interesting application of sedimentation upstream of a check dam (i.e. dams designed for debris flow control in steep mountain streams). They 33 applied a 1D FV model to reproduce some of the morphological experiments of Arminini and Larcher (2001). The experimental channel has a 5% slope that generates supercritical flow in the upstream reach. In the downstream end, a slit check dam forces subcritical flow and a hydraulic jump. A prograding sediment delta forms in the location of the hydraulic jump and advances downstream towards the dam. The numerical result compared well with measured data. Papanicolaou et al. (2005) applied the model 3ST1D for unsteady flow conditions occurring over transcritical flow regimes such as flows over step-pool sequences. 3 ST ID can handle transcritical flow and multisize sediment distribution. The model is ID and based on a FV shock-capturing numerical scheme. Dam removal is another active field of current research where transcritical modeling is important. The steep fronts of deposited sediment upstream from the dam (usually at the angle of repose) can generate supercritical flow and hydraulic jumps (Cui and Wilcox, in press). Cui et al. (in press, a, b) developed the ID model D R E A M for the morphodynamic simulation of sediment transport after dam removal, and the model was tested with laboratory and field data. The rapidly varied flows caused by dam break and dyke breach represent another increasingly active area of current research where transcritical morphodynamic capabilities are required (Capart and Young 1998, Ferreira et al. 2001, Leal et al. 2001, Zang and Wang 2001, Capart and Young 2002, Fraccarollo and Capart 2002, Spinewine and Zech 2002, Caleffi and Valiani 2002, Leal et al. 2003, Cao et al. 2004). Caleffi and Valiani (2002) applied two ID models to simulate the experiments of hydraulic jump formation by dam break wave over loose granular beds reported by Capart and Young (1998). The models were based on both FD and FV discretizations. The experiments of dam break were with initially dry and wet downstream channel. Both FD and FV models led to very similar results, that compared reasonably well with experimental data. Capart and Young (2002) developed a 2D F V model with 3 horizontal layers of pure water, slurry sediment-water mixture and saturated sediment. The model was applied by Spinewine et al. (2002) to simulate bankline retreat in a trapezoidal channel after a sudden dam break. 34 Kusakabe et al. (2003) simulated transcritical flow in a widening of a steep channel with a 2% slope, using a 2D FD model. The model predicted the scour hole downstream of the expansion, although its depth was under-estimated. Zang and Wang (2001) provide one of the few examples of 2D FE methods applied to transcritical flow morphology. The model solves the non-equilibrium suspended sediment transport equation applied to hypothetical cases of dyke breach and dam break. Bedload transport is not considered. It is unclear whether triangular or quadrilateral elements were used. 1.4.6 Relevant issues of literature review From the literature review, it is important to highlight some issues relevant for the objectives of the present research. 1. The sediment continuity or Exner equation is mathematically a hyperbolic differential equation that physically represent two main mechanisms: diffusion (dispersion) and convection (advection). In cases when convection is not important, the Exner equation can be reduced to a parabolic equation (diffusion model). 2. Diffusion seems to dominate in problems involving aggradation and degradation of alluvial channels caused by forced upstream sediment supply; as well as the generation of forced bars (point bars) and pools in curved channels. Convection dominates in cases when migrating bedforms are present, such as free alternate bars, dunes, prograding deltas, etc. 3. Conventional Galerkin Finite Element Methods are suitable for elliptic or parabolic equations, but not for hyperbolic ones (where convection dominates). For hyperbolic equations Petrov-Galerkin methods should be used. 35 4. Currently there are several advanced 2D river morphology models, but most of them are based on structured meshes and quadrilateral elements. Models based on triangular Finite Elements are still uncommon. 5. Morphodynamic models for transcritical flow are mainly limited to ID flow, although some examples of 2D models exist. It is unknown if any 2D FE model based on triangular elements with capabilities for bedload transport in transcritical flow presently exist. 6. Depth-averaged models with additional moment of momentum equations (VAM) have proved superior to conventional depth-averaged models for flow in curved open channels. However, 2D V A M models have never been applied for alluvial bend morphology. 1.5 S U M M A R Y For millennia, human civilization has flourished around rivers. During most of this time, our knowledge of rivers was based only on empirical observations and river engineering was no more than a trial-and-error endeavour. In the last five centuries, an unprecedented growth in our understanding of river morphology has emerged from both experimental and theoretical work. Although our knowledge is still incomplete and far from perfect, we can now reproduce by means of numerical models many of the observed features of alluvial rivers. The first morphological models assumed a one-dimensional approximation. It soon became evident that although ID models are suitable for long term and long reach simulations, they are not appropriate for accurately reproducing several important morphological features, notably scour and deposition along bends. 2D and 3D models emerged as an alternative. However, at present, 3D morphological models are still computationally too demanding for practical use, leaving 2D models as the tool of choice for river morphology. 36 Initially, 2D morphology models solved the flow and sediment equations using Finite Difference Method on rectangular orthogonal grids (meshes), which proved inadequate for the meandering planform of most rivers. A major advance was the introduction of curvilinear non-orthogonal meshes that could adapt exactly to the riverbanks. These boundary-fitted meshes are perfect for single thread channels with smooth banklines. However, their application becomes problematic for rivers with very complex geometries. That attracted interest to import the Finite Element Method widely used in structural modeling. The FE method has the great advantage of accommodating practically any conceivable planform geometry with meshes that can be easily generated by automatic algorithms. The resolution of these unstructured FE meshes can easily be increased in areas of interest, one of its main advantages. However, the initial attempts to apply FE to flow hydraulics in the 1970's were not very successful. The conventional Galerkin Finite Element Method, which proved extremely successful to solve the elliptic or parabolic equations of structural engineering, is not adequate for the hyperbolic equations of hydrodynamics, where convection is important. These initial difficulties plus the numerical complexity of the method, help explain the relative small share of FE in the pool of currently available 2D river morphology models. Fortunately, a new generation of Petrov-Galerkin FE methods overcame the initial problems with hyperbolic systems and have proved even more accurate than some FD methods for various hydrodynamic applications (Katapodes 1984a, Chung 2002). . The main objective of the present thesis is to develop a 2D river morphology model that exploits the advantages of unstructured triangular finite elements. The model was developed by programming and adding a solver for the bedload sediment continuity equation to existing state of the art computer codes of river hydrodynamics developed at the University of Alberta (Steffler and Blackburn 2002, Ghamry and Steffler 2005). Currently, the bedload model is limited to uniform sediment and intended for applications where convection is not dominant. The features of the model will be demonstrated in the following chapters of this thesis. 37 In Chapter 2, the model is applied to alluvial straight flumes. Two laboratory experiments of bed aggradation and degradation forced by overload and deficit of upstream sediment supply are tested. Additionally, an experiment of knickpoint migration in transcritical flow is also simulated. The latter test case demonstrates the unique transcritical capabilities of the model. In Chapter 3, two additional sub-models, a secondary flow correction and a transverse slope correction, are incorporated to simulate scour and deposition along curved alluvial channels. Three well-known laboratory experiments are used as test cases to assess the depth-averaged morphodynamic model. In Chapter 4, a novel approach for alluvial bend morphology is introduced. A depth-averaged model with additional moment-of-momentum equations is applied to alluvial curved channels without the need for a secondary flow correction. This new approach was successfully tested with both lab data and a full scale meandering river. Finally, the summary and conclusions of the entire thesis are discussed in Chapter 5. Some topics for further research are also suggested, such as enhancements of the model to deal with convection-dominated problems, multiple size sediment and dam-breaks over movable beds, among others. Several appendices contain additional information that complements this research. Topics included in the appendices are: details of the numerical discretization; flow curvature computation; applications to channels of varying width; problems with convection dominated applications; qualitative applications to dam breaks over movable beds; among others. 38 2 STRAIGHT ALLUVIAL CHANNELS 2.1 INTRODUCTION 2.1.1 Bed aggradation and degradation In a natural alluvial river an equilibrium or "regime" condition may exist between the How and sediment transport such that the riverbed elevation remains on average practically constant (Blench 1969). When this delicate equilibrium is altered, the riverbed responds by aggrading (sedimentation) or degrading (erosion) in an attempt to achieve a new equilibrium state. According to Graf (1998) riverbed aggradation may occur when the river is subject to a load of sediment higher than the flow is capable of transporting (e.g landslide); the water discharge is decreased (e.g. downstream from a water diversion); or the construction of an elevated fixed point in the channel (e.g. dams). Bed degradation may occur when upstream sediment supply is interrupted (e.g. below a dam); water discharge increases; or when a point downstream in the channel is lowered (e.g. upstream from a meander cut-off). Sometimes, both aggradation and degradation are observed in a given reach (Galay 1983). For example, a meander cutoff generates a sudden change in the bed slope known as "knickpoint"; upstream from the 39 knickpoint degradation occurs, while aggradation happens downstream (Brush and Wolman 1960). The description of bed aggradation and degradation mentioned above is simplified to one dimension. According to the principles of self-adjustment (Blench 1969) natural rivers have several degrees of freedom: width, depth, slope, roughness and sinuosity, that can change, with different time scales, to accommodate new flow and sediment conditions (USACE 1993). In laboratory flumes with fixed walls, the degrees of freedom are reduced to changes in water depth and bed slope only; allowing an easier study of the phenomenon. Under these simplified conditions.it is even possible to reduce the flow and sediment equations into a simple diffusion (parabolic) model amenable to analytical solution (Soni et al. 1980, Gi l l 1980, Jain 1981, Jaramillo and Jain 1984, Graf 1998), as was shown in section 1.4.1.2. Soni et al (1980) provide examples of experimental data of bed aggradation due to sediment overloading; while Suryanarayana (1969) provides data of bed degradation caused by sediment supply shut-off. These simple experiments are excellent benchmark tests for any morphological model. The high quality experimental data of Soni et al. (1980) have been extensively used by several researchers to tests their analytical or numerical models of bed aggradation (Gill 1980, Jain 1981, Gi l l 1983a, Jaramillo and Jain 1984, Park and Jain 1986, Gi l l 1987, Bhallamudi and Chaudhry 1991, Kassem and Chaudhry 1998, Cao et al. 2002). Bhallamudi and Chaudhry (1991) applied a one-dimensional (ID) Finite Difference (FD) scheme to solve the flow and sediment equations. Kassem and Chaudhry (1998) applied a 2D FD model to investigate the suitability of semi-coupled and fully coupled models. Fully coupled models are those that simultaneously solve the Saint-Venant (flow) and Exner (sediment) equations, while semi-coupled (or decoupled) models are those that solve the flow equations first, and then the sediment equations. Kassem and Chaudhry (1998) concluded that there is no significant difference between either approach when simulating the experiments of bed aggradation of Soni et al. (1980) and bed degradation of Suryanarayana (1969). Cao et al. (2002) using a ID FD model also investigated the suitability of coupled and decoupled models. Surprisingly, Cao et al. (2002) stated that for the Soni et al. (1980) experiment "decoupled modelling is unachievable and i l l posed mathematically" because the Froude number became temporarily 40 supercritical in the inflow cross section. Cao et al. (2002) made that argument despite the fact that several previous researchers have successfully replicated the Soni et al. (1980) experiment using decoupled models (Park and Jain 1986, Bhallamudi and Chaudhry 1991, Kassem and Chaudhry 1998). The experimental data of knickpoint migration reported by Brush and Wolman (1960) represent a very interesting and challenging problem because of the presence of transcritical flow, generated by a short oversteepened reach with 10% slope. Bhallamudi and Chaudhry (1991) simulated this experiment; however they did not observe a hydraulic jump in their computations (Bhallamudi, personal communication), likely because of a very coarse space discretization and a downstream boundary condition higher than reported in the experiment. The experimental data of Soni et al. (1980), Suryanarayana (1969) and Brush and Wolman (1960) were used to test the capabilities of the model developed in this thesis for straight alluvial channel morphology. The morphological model was derived by semi-coupling a bedload sediment transport module into the existing 2D hydrodynamic model River2D. Below is a description of the hydrodynamic model, followed by a description of the sediment transport module. 2.1.2 River2D Hydrodynamic Model River2D is a 2D vertically-averaged (depth-averaged) Finite Element model intended for use on natural streams and rivers and has special features for supercritical/subcritical flow transitions, ice covers, and variable wetted area. It has been developed at the University of Alberta by a research group led by Dr. Peter Steffler (freely available at www.River2D.ca). River2D is based on the 2D vertically averaged Saint Venant equations expressed in conservative form, which form a system of three equations representing the conservation of water mass and the two components of the momentum vector (Steffler and Blackburn 2002). 41 Conservation of water mass: dt dx dy Conservation of momentum equation in the x-direction: ^ + « + « + g ^ = ( S o x - S f x ) + ± -dt dx dy 2 dx V o x Jx' p Conservation of momentum equation in the ^-direction: dx dy dqy 8(uqy) d(vqy) g d h 2 , , ^ + ~dx- + ^ y - + 2 ^ y - = 8 h [ S ° > - S * } + dx dy (2.1) (2.2) (2.3) where h is the water depth, u and v are the vertically-averaged velocities in the x and y directions respectively, while qx and qy are the respective discharge intensities, related to the velocity components through qx =uh qy =vh (2.4) (2.5) Sox and Soy are the bed slopes in the x and y directions; 5^ and Sjy are the corresponding friction slopes. rxx, Txy, ryx, tyy axe, the components of the horizontal turbulent stress tensor, g is the gravitational acceleration and p is water density. The basic assumptions in equations (2.1) through (2.5) are: - The pressure distribution is hydrostatic, which limits the accuracy in areas of steep slopes and rapid changes in bed slopes. 42 The horizontal velocities are constant over the depth. Information on secondary flows and circulations is not available. Coriolis and wind forces are assumed negligible. The friction slope terms depend on the bed shear stresses which are assumed to depend on the magnitude and direction of the vertically averaged velocities. For example, in the x direction: AJx - ' Pgh ghd (2.6) Tbx is the bed shear stress in the x direction and C* is the dimensionless Chezy coefficient, which is related to the effective roughness height ks through C. =5.75 log 12 — V ks J A>£1 *, " 1 2 (2.7) e r i \ \ks J h_ ef_ k« < 12 (2.8) C* is related to Chezy's C coefficient through C.= C (2.9) The vertically-averaged turbulent shear stresses are modeled with a Boussinesq type eddy viscosity. For example: rdu dv^ — + — dy dx (2.10) 43 where vt is the eddy viscosity coefficient assumed by default as hylu2 + v2 n i n v, = 0.5 (2.11) The Finite Element Method used in River2D's hydrodynamic model is based on the Streamline Upwind Petrov-Galerkin weighted residual formulation (Hicks and Steffler 1992). In this technique, upstream biased test functions are used to ensure solution stability under the full range of flow conditions, including subcritical, supercritical, and transcritical flow. A fully conservative discretization is implemented which ensures that no fluid mass is lost or gained over the modeled domain. This also allows implementation of boundary conditions as natural flow or forced conditions (Steffler and Blackburn 2002). 2.1.3 Bedload Sediment Model By considering only bedload transport the 2D sediment continuity equation may be written as (Struiksma et al. 1995, Kassem and Chaudhry 2002, Cui et al. 2005) ( l _ A ) ^ + ^ k + -?^L = 0 (2.12) dt dx dy qsx and qsy are the components of the volumetric rate of bedload transport per unit length in the x and y directions (horizontal plane), X is the porosity of the bed material (a default value of 0.4 is assumed), t is time and Zb the bed elevation. The bedload transport rate is a function of the flow hydraulics and sediment properties. A very simple approach is to assume that sediment transport depends only on velocity according to an empirical power law relationship of the form (Kassem and Chaudhry 1998) A - l qsx = au(u2 + v 2 ) (2.13) 44 b-\ q =av(u2 +v 2 ) 2 (2.14) a and b are empirical parameters, b typically varies between 3 and 5. Equations (2.13) and (2.14) usually provide satisfactory results for uniform sand. Other sediment transport equations (e.g. Engelund-Hansen, Meyer-Peter Muller, Van Rijn) have also been coded, but they were not used in the applications shown in this chapter. Figure 5 shows the main steps of the semi-coupled morphological model. At every time step during the simulation, the hydrodynamic model (equations 2.1 through 2.11) computes water depth h and velocity components u and v. Those values are used in equations (2.13) and (2.4) to compute the sediment transport rates, which are applied in equation (2.12) to compute the changes in bed elevation Azb in a given time step At. A second order Runge-Kutta scheme is used to solve equation (2.12). The bed elevation is updated asz"bew =zfd + Azh . Next, the hydrodynamic model is invoked again to re-compute the flow field for the new bed topography. This procedure is repeated every time step until the final simulation time T is achieved. This approach, in which sediment equations are solved after the flow equations, is usually referred as decoupled or semi-coupled (Kassem and Chaudhry 1998). In contrast, in a fully coupled approach, equation (2.12) is solved simultaneously with the flow equations; Zb is treated as an additional unknown in a 4-by-4 system (h, u, v, Zb). The details of the numerical discretization are shown in Appendix A. Instructions of how to use the model are given in Appendix B. The logic of the program with its numerical algorithms and functions is explained in Appendix C. As mentioned in Chapter 1, the morphological model at its present stage is intended for morphological applications where convection is not dominant, because it uses a conventional Galerkin Finite Element Method (see Appendix A). This numerical scheme produces spurious oscillations i f applied for migrating bedforms, as detailed in Appendix F. 45 Input initial bed elevation and initial flow field. Solve flow hydrodynamics (equations 2.1 -2.11) ! Compute sediment transport rates ! (equation 2.13) Solve sediment continuity (equation 2.12) Output final bed elevation and final flow field. Figure 5. Flow chart of morphological model. Bedload module is inside dashed box. 2.1.4 Objectives The main objective in this Chapter is to test the capabilities of the model for simulating three flume experiments in straight alluvial channels, while at the same time assessing its numerical stability. The experiments selected provide real observed data to objectively assess the model's prediction capabilities. The three experiments are: (1) bed aggradation due to sediment 46 / Initial / Conditions > h, u, v Azb overloading (Soni et al. 1980), (2) bed degradation due to sediment feed shut-off (Suryanarayana 1969) and (3) knickpoint migration (Brush and Wolman 1960). The results demonstrate that the computed bed profiles agree reasonably well with measured data without any noticeable stability problem. 2.2 B E D A G G R A D A T I O N TEST C A S E 2.2.1 Soni et al. experiment Soni et al. (1980) performed experiments of bed aggradation due to sediment overloading in a straight alluvial flume. The experimental flume was 30 m long and 0.20 m wide covered by uniform sand with a median diameter of 0.32 mm. From 24 measurements of sediment transport rate for flow velocities between 0.2 and 0.6 m/s, they found that the value of parameters a = 0.00145 and b = 5.0 in equation (2.13) resulted in good agreement with measured transport rates. Starting from equilibrium conditions the inflow of sediment qsiN was increased up to several times the equilibrium value qse in order to induce bed aggradation. The measured bed profiles were averaged because of the presence of ripples and dunes, especially in the non-aggrading reach. About 5% of the injected sediment was deposited upstream of the section of sediment addition. The experiment with the highest inflow of sediment qsw = 4.0qse was selected for the numerical simulation. The initial conditions were: bed slope S0 = 0.036; depth h0 = 0.05 m and water discharge Q = 4 L/s. The flume was discretised using 336 triangular elements and 338 nodes. The average distance between nodes was set equal to the flume width Ax = 0.20 m. The roughness height was set to ks = 11 mm and remained constant during all the simulation. 47 2.2.2 Results Figure 6 shows the results of the computed bed and water surface profiles at / = 15 and t = 40 minutes for the Soni et al. (1980) experiment of bed aggradation. The numerical time step was set to At - 20 s. The Froude number computed in the upstream cross section became temporarily supercritical (F = 1.3) during the first minute, but returned to subcritical values afterwards. Figure 6. Observed and computed longitudinal profiles at 15 and 40 minutes for Soni et al. (1980) bed aggradation experiment due to an overload of sediment equal to four times the equilibrium value. 48 2.2.3 Discussion There is good agreement between the computed and measured results for the aggradation test case (Fig. 6). The model is able to capture reasonably well the transient changes in the evolution of the bed profiles. The difference between observed and computed values can be partially attributed to the varying bed roughness during the experiment due to the presence of bedforms (the model assumed constant roughness), the uncertainty in the sediment transport equation, and the small fraction of injected sediment deposited upstream of the first cross section (ignored by the numerical model). In a recent numerical simulation of this experiment, Cao at al. (2002) were forced to artificially reduce the friction factor to 80% relative to the steady-state case in order to obtain good agreement; such reduction was unnecessary in the present model. Cao et al. (2002) also reported that a decoupled model could not be applied to this case because it is not possible to define an upstream boundary condition for sediment transport when the flow becomes supercritical. The computed results (Fig. 6) contradict that conclusion, but agree with Park and Jain (1986), Bhallamudi and Chaudhry (1991) and Kassem and Chaudhry (1998) who successfully applied decoupled models for this problem. Good results when applying decoupled models were also found by Cui et al. (1996) when simulating aggradation of heterogeneous gravel. Cui et al. (1996) reported negligible differences between coupled and decoupled models to simulate bed aggradation despite having an upstream Froude number close to one. Park and Jain (1986) performed a numerical study of bed aggradation due to overloading using a ID Finite Difference approach. They stated that stability requirements impose a limit on the time step, based on the diffusion coefficient k (eq. 28 in Park and Jain 1986; see also Jansen et al. 1979 and Graf 1998) 49 (2.15) Applying equation (2.15) to the space discretization used in this research, Ax = 0.20 m, will lead to a maximum time step At = 1.8 s. However, the simulations remained stable with a time step about an order of magnitude higher: At = 20 s (Fig. 6). In contrast, all previous ID simulations performed by Park and Jain (1986), Bhallamudi and Chaudhry (1991), Kassem and Chaudhry (1998) and Cao et al. (2002) used time steps equal to or significantly smaller than the limit given by equation (2.15). It seems that the present application has the largest time step reported for the numerical simulation of the Soni et al. (1980) experiment. The probable reason for the enhanced stability of the model is the implicit time stepping algorithm used by River2D, which has fewer stability restrictions than conventional explicit time stepping algorithms commonly used by FD. 2.3 BED DEGRADATION TEST CASE 2.3.1 Suryanarayana experiment Suryanarayana (1969) performed experiments of bed degradation due to the shut-off of upstream sediment supply; similar to the effect of a dam construction in a river. The flume was 18.29 m (60 ft) long and 0.61 m (2 ft) wide covered by uniform sand with a median diameter of 0.45 mm. The parameters used in equations (2.13) and (2.14) were a = 1.786*10"4 and b = 3.878 (Kassem and Chaudhry 1998). Starting from equilibrium conditions, the upstream sediment supply was shut off leading to bed degradation. The initial flow conditions of the selected test case were: S0 = 0.007, h0 = 0.034 m and Q= 11.9 L/s. The results were given for 10 hours of flow time. 50 In the numerical model the flume was discretised using 70 triangular elements with 72 nodes. The average distance between nodes was set equal to the flume width Ax = 0.61 m. The roughness height was set equal to ks = 2.7 mm. The upstream boundary condition was set to constant and equal to qsfN = 0. Lu and Shen (1986) stated that since it takes a longer distance for clear water to pick up enough sediment from the channel bed to reach a high sediment transport capacity for a steeper slope condition, a reduction coefficient Cr = 0.4 must be introduced at the upstream end. An identical approach was also used by Kassem and Chaudhry (1998). For the present numerical model Cr was applied to the sediment fluxes across the sides of the most upstream elements attached to the inflow cross section. Two simulations were performed for Cr = 1.0 (no reduction) and Cr = 0.4. 2.3.2 Results Figure 7 shows the computed bed and water surface profiles at 10 hours after the beginning of the bed degradation experiment, using 2 correction coefficient values Cr = 1.0 and Cr = 0.4. The time step was At = 30 s. The Froude number remained subcritical during the entire simulation. The results of Figure 7 illustrate an example of extreme bed degradation. The water surface at / = 10 hours in the upstream end of the channel is indeed below the original bed elevation. 51 0.35 Observed Computed Cr = 0.4 Original bed Computed Cr= 1.0 6 8 1 0 1 2 Channel distance (m) Water surface Bed 1 4 1 6 1 8 Figure 7. Observed and computed longitudinal profiles at 10 hours for bed degradation experiment due to sediment supply shut-off for two correction coefficient values of Cr = 1.0 and Cr = 0.4. 2.3.3 Discussion There is a reasonable agreement between the measured and computed bed and water surface profiles as shown in Figure 7. The agreement was less accurate in the upstream reach where an adverse bed slope is present. Bedload sediment moving along this adverse slope is subject to a component of weight that acts against the direction of the flow, reducing sediment mobility. Such sloping effect is ignored by the simple equations (2.13) and (2.14) that depend only on flow velocity. This effect is somehow taken into account by the correction coefficient Cr applied to the upstream section (Lu and Shen 1986, Kassem and Chaudhry 1998). If the sediment flux is not reduced upstream (Cr = 1.0), the computed profiles fit almost perfectly 52 the observed profiles in the last two thirds of the flume (Channel distance > 6 m in Fig. 7); however, the water surface elevation is under-predicted in the upstream sections. When the correction Cr = 0.4 is applied, a better fit was obtained in the upstream reach, but to the detriment of the values in the downstream reach. Besides the adverse bed slope effect, another possible reason for the discrepancy in the upstream section is the application of the conventional equilibrium sediment transport equation (2.12). This equation assumes a local equilibrium condition in which flow transports sediment at its full capacity. In this case, that means that flow goes from imposed zero transport in the upstream inflow section to full capacity immediately downstream. Phillips and Sutherland (1989) argue that alluvial systems require a certain distance before the equilibrium state is recovered when upstream sediment boundary conditions are constrained or imposed. Phillips and Sutherland (1989) referred to this phenomenon as a spatial lag. Perhaps the application of a non-equilibrium sediment transport model, as proposed for highly unsteady flows (e.g. Nagata and Muramoto 2000, Zhang and Wang 2001, Due et al. 2004, Wu 2004, Cao et al. 2004, Wu et al. 2005) would be more appropriate for this case. However, this is a topic for further research, as discussed in Chapter 5. 2.4 K N I C K P O I N T M I G R A T I O N T E S T C A S E 2.4.1 Brush and Wolman experiment Knickpoints are points along longitudinal profiles of streams where the slope increases abruptly. When they are present in alluvial channels, in an attempt to reduce the bed slope the flow erodes the steep bed slope reach causing the knickpoint to migrate upstream. Brush and Wolman (1960) carried out experiments at the University of Maryland in a flume 15.85 m (52 ft) long and 1.22 m (4 ft) wide. Before starting each experimental run a trapezoidal channel 0.21 m wide and 0.03 m deep with rounded corners was moulded in noncohesive sand. A short 53 over-steepened reach with a length of 0.30 m (1 ft) and a fall of 0.03 m (0.1 ft), was located 10.8 m from the flume entrance. This oversteepened reach immediately downstream of the knickpoint had a 10% slope which was significantly steeper than the adjacent reaches. A total of 5 runs were performed with different slopes in the upstream and downstream reaches (between 0.125 and 0.88%). In every run, the slope of both water surface and bed below the knickpoint decreased with time, causing the position of the knickpoint to migrate upstream. As the knickpoint moved upstream the channel directly above it narrowed. At the lower end of the steep reach, sediment eroded from above deposited as a dune, which moved downstream causing the channel to locally widen. Run 1, which has more detailed information, was selected for the numerical simulation. For Run 1 the longitudinal bed slope changed at the knickpoint from 0.125% to 10% and back to 0.125% again (Fig. 9). The channel was moulded in sand 0.67 mm in diameter with initial water depth h0 = 0.0137 m and constant flow rate Q = 0.59 L/s. During Run 1 the width of the channel increased downstream of the knickpoint by bank erosion to about 0.25 m, while it got deeper and narrower upstream. At the beginning of Run 1, no sediment was moving in the reaches above and below the steep reach. The erosion and deposition resulted solely from the presence of the knickpoint. Sediment eroded above the knickpoint was deposited below. The head and toe of the steep reach moved upstream and downstream respectively. Sediment transport rate was not measured. In the numerical model the channel was assumed as rectangular with a fixed width of 0.21 m. Unlike the two previous test cases, a more irregular mesh was used in this case. The average size of the elements was about 0.20 m away from the steep reach and about 0.04 m around the knickpoint (Fig. 8). The mesh was refined close to the knickpoint to capture the strong flow gradients and the hydraulic jump. The mesh had 869 triangular elements and 563 nodes. The roughness height was set to ks = 0.3 mm. 54 Figure 8. Refined Finite Element mesh around knickpoint and steep 10% slope reach (parallel gray lines are bed contours every 0.001 m). The boundary conditions for the hydrodynamic model were constant discharge in the upstream inflow section, and constant water surface elevation in the downstream outflow section. Given that most sediment transport equations are derived for subcritical uniform flow conditions, there is considerable uncertainty in this case, where supercritical flow occurs. Bhallamudi and Chaudhry (1991) adopted a = (S» _ 1 7 and b = 4.2 in equations (2.13)-(2.14), Sf being the friction slope. The same values were applied in the present simulation using equation (2.6) for the friction slope. 2.4.2 Results Figure 9 shows the initial water surface profile computed before the beginning of the morphological simulation. Away from the steep reach, the flow remained subcritical with a Froude number of approximately 0.5, which corresponds to a velocity of about 0.2 m/s and a depth of 0.014 m, in agreement with observed values (Table 3 in Brush and Wolman I960). As the flow approaches the knickpoint, it accelerates and becomes supercritical. The maximum velocity of about 0.7 m/s is reached over the steep reach, the Froude number being as high as 3.5 with a water depth of only 0.004 m. 55 Downstream from the toe of the steep reach, the flow returns to subcritical through a hydraulic jump. The sharp front of the jump was captured by the model without any noticeable spurious oscillation, as shown in Figure 9. 0 .69 -j 0 .68 -0 . 6 7 -g c Q 0 .66 -w > 0 .65 <D UJ •o 0 .64 CD 0 . 6 3 -0.62 -0.61 -10 S0 =0.125% Knickpoint Hydraulic jump So = 10%v So = 0.125% 11 Water surface Bed 12 Channel distance (m) Figure 9. Initial conditions for knickpoint migration experiment (t = 0). Water surface profile was computed by the numerical model. A very small time step was used At = 0.01s to prevent large bed changes during the initial period of the simulation The results showed that both the computed bed and water surface profiles were rather smooth, without any noticeable numerical oscillation. A comparison between the computed profiles and the bed profiles measured by Brush and Wolman (1960) at / = 160 min is shown in Figure 10. By this time, the intensity of the hydraulic jump has dramatically decreased since the maximum value of the Froude number was only slightly larger than 1.0. The steep bed slope has reduced to only 0.44% from its original value of 10%. The computed bed profile was smoother than the observed profile and 56 tended to follow the average trend of the measurements. The new upstream location of the knickpoint agreed well with the observed value. 0 . 6 9 0) 111 -a 0 . 6 4 CD 0 .63 0 .62 0.61 - - • Original Computed • Measured Knickpoint Water surface Bed 6 8 10 Channel distance (m) 12 14 16 Figure 10. Observed and computed longitudinal profiles at t = 2 hr 40 min for knickpoint migration experiment. Bhallamudi and Chaudhry (1991) also simulated this experiment using a ID model. However, they did not observe a hydraulic jump in their numerical simulations (Bhallamudi, personal communication 2004), likely because the space discretization used (Ax = 0.30 m) was too large to capture the sharp front of the jump, and also because their downstream water level was higher (ha = 0.03 m). 57 2.4.3 Discussion The morphodynamic simulation of the knickpoint migration experiment is a challenging problem because of the steep 10% bed slope that gives rise to a hydraulic jump and transcritical flow conditions. ID morphological models with transcritical flow capabilities exist (e.g. Busnelli et al. 2001, Fraccarollo & Capart 2002, Papanicolaou et al. 2005, Cui at al.a,b, in press), but they are still rare in 2D (Zhang and Wang 2001, Cappart and Young 2002, Kusakabe 2003). The shock-capturing techniques incorporated in River2D, added to the flexibility of the Finite Element method to increase mesh resolution around areas with strong gradients (Fig. 8), made it possible to simulate the sharp hydraulic jump front (Fig. 9) without any numerical oscillation. This test case would be intractable for most 2D commercial models currently available (e.g. RMA2, M I K E 21C, FESWMS). Busnelli at al. (2001) argue that hydraulic jumps tend to disappear rapidly in alluvial beds under constant water discharge. The experiments reported by Belial et al. (2003) show that the amplitude of the hydraulic jump is drastically reduced by a movable sediment front. In the final computed bed profiles shown in Figure 10, the jump is barely noticeable. The agreement between measured and computed bed profiles is reasonable (Fig. 10), considering the large uncertainty in the sediment transport equation and the fact that both bank erosion and deposition have been completely ignored by assuming a fixed width. Brush and Wolman (1960) reported that the channel became wider downstream of the knickpoint and narrower upstream; the maximum changes in width were about ±20% of the original width. Recently, Cantelli et al. (2004) and Wong et al. (2004) have demonstrated the importance of width changes in the morphological evolution of steep channels after dam removal, which can be considered as a knickpoint-like problem. Ignoring those effects may be responsible for some of the discrepancy observed in Figure 10. 58 A very interesting potential application of the model is the simulation of sediment transport after dam removal, which is a topic of high current interest (Cantelli et al. 2004, Wong et al. 2004, Cui and Wilcox, in press, Cui at al. in press a,b). The hydraulic conditions after dam removal shown in Cui and Wilcox (in press) are practically identical to those depicted in Figure 9, with the over-steepened reach having a slope close the angle of repose of sediment. The computed bed profiles shown in Figure 10 also have some similarities with the measured bed profiles after dam removal reported by Cantelli et al. (2004). 2.5 S U M M A R Y This chapter describes the first applications of the 2D river morphology model applied to straight alluvial channels. The model was developed by semi-coupling a solver for the bedload sediment continuity equation to the existing River2D hydrodynamic model, making it possible to dynamically update the channel bed elevation as sediment is eroded or deposited by the flow. Thus, the model can reproduce the transient evolution of the bed as it changes through time. The sediment transport equations are not fully coupled to the flow equations; i.e., at every time step the sediment continuity equation is solved after the flow field has been computed, in a fashion usually referred as decoupled or semi-coupled. The results of the model challenge conclusions of several previous researchers who have heavily criticized decoupled models stating that they are unstable, not capable of simulating rapidly changing boundary conditions or supercritical flow (Lyn 1987, Correia et al. 1992, Cao et al. 2002). None of the computed results showed any of those problems, even when relatively coarse meshes and large time steps were used. In agreement with Kassem and Chaudhry (1998), it appears that the stability problems reported by previous researchers might be due to the numerical approach they employed to solve the mathematical equations, rather than being a fundamental limitation of decoupled models. 59 The model has demonstrated its capabilities to simulate successfully bed level changes in straight alluvial channels, even in cases with rapidly varying boundary conditions, supercritical flow and hydraulic jumps. The model remained stable and provided satisfactory agreement for three experimental test cases: (1) bed aggradation, (2) bed degradation and (3) knickpoint migration. The model uses a flexible unstructured mesh based on triangular elements, which can be easily refined around areas where more detail is sought, for example in areas with strong gradients, such as hydraulic jumps. Another feature of the model is that the upstream boundary condition can be set as free or with an imposed sediment load. The latter condition allows simulating cases of aggradation or degradation caused by changes in the upstream sediment supply, e.g. below a dam. The transcritical flow capabilities of the model make it possible to deal with problems involving supercritical flow and hydraulic jumps, which are rather uncommon for most 2D models (see section 1.4.5). Potential applications of this capability include the simulation of bed changes after dam removal (Cui and Wilcox, in press) or dam-breaks (Appendix G). The model is currently limited to bedload dominated flows (suspended sediment not important or negligible) and uniform sediment. However, since flow and sediment equations are not fully coupled, it is feasible to easily incorporate in the future additional modules for suspended sediment or non-uniform sediment (gravel beds). Applications of the model to curved alluvial channels subject to secondary flow conditions are presented in the next chapter. 60 3 CURVED ALLUVIAL CHANNELS 3.1 I N T R O D U C T I O N Bends are very common features of natural rivers and their importance in river engineering practices is well recognized (see section 1.2). Flow of water and sediment in erodible bends is extremely complex for the presence of a three-dimensional (3D) spiral motion, which generates erosion on the outer bank and deposition of sediment on the inner bank. Erosion can destroy or undermine structures located in the river, while deposition (point bars) may block water intakes or navigation canals. Thus, prediction of erosion and deposition patterns in erodible rivers is a very important and challenging task for river engineers, who usually resort to 2D numerical models for that purpose. Unlike flow in straight alluvial channels, where sediment and mean (depth-averaged) velocity have the same direction, bedload sediment moving along bends follows a path different from 61 that of the mean flow (Fig. 12). There are two reasons for this. First, the centrifugal acceleration induced by curvature generates a secondary motion in the transverse direction (perpendicular to the primary flow) that causes the bed shear stress to deflect towards the inner bank (Fig. 11). Second, the scour on the outer bank and deposition on the inner bank generate an important transverse slope. A sediment grain moving along this transverse slope is subject to an additional component of weight in the direction perpendicular to the flow (Fig. 13). The combined effect of the secondary motion and the transverse slope cause bedload direction to be different from the mean flow velocity direction computed by depth-averaged models. In order to apply the depth-averaged morphology model introduced in Chapter 2 to compute bed changes in curved alluvial channels, it had to be enhanced by adding two additional sub-models. These sub-models take into account the effects of the secondary flow on the direction of the bed shear stress and the influence of the lateral bed slope on the direction of bedload sediment transport, as discussed below. 3.1.1 Bedload sediment model By considering only bedload transport, the 2D sediment continuity equation may be written as (Struiksma et al. 1985, Struiksma 1989, Kassem and Chaudhry 2002) ( 1 _ A ) ^ + ^ + ^ L = 0 (3.1) dt dx dy qsx and qsy are the components in the x andy directions (horizontal plane) of the volumetric rate of bedload transport per unit length qs, X is the porosity of the bed material (a default value of 0.4 is assumed), t is time and zb the bed elevation. The components of the sediment transport rate depend on the sediment transport direction a 62 asx=asC(>sa (3.2a) qsy=qs since (3.2b) The rate of sediment transport qs can be computed from many of the equations available. In this study, the Engelund-Hansen sediment transport equation was adopted (Kassem and Chaudhry 2002) qs =0.05 — jg(Gs-\)D350 (T*)5N (3.3) g C is the Chezy roughness coefficient, Gs is the specific gravity of the sediment, D5o the median grain diameter, g is the gravitational acceleration and T* is the dimensionless shear stress (Shield's parameter) u2+v2 — 2 (3-4) u and v are the depth-averaged components of the velocity vector in the x and y directions respectively. Notice that in a straight flume, where the direction of bedload and mean flow coincide, the direction of bedload in equation (3.2) is simply \&na = (v/u). However, that is not the case in curved flumes as explained below. 63 3.1.2 Influence of secondary flow (helical motion) The difference between flow in a straight and in a curved channel is the presence of centripetal acceleration in the latter (or centrifugal depending on the reference frame). Assuming a simple ID streamline, the centripetal acceleration is ac = u2/rc, where rc is the radius of curvature of the streamline. Since the flow velocity distribution in the vertical is zero at the bottom and maximum close to the water surface, ac reaches a maximum value close the water surface. Therefore, in a cross section transverse to the bend the surface layers of water are pushed in the outward direction. To fulfill mass continuity (zero mass flux in the transverse direction), water in the lower layers (close to the bed) moves in the inward direction (Fig. 11). When this secondary flow circulation combines with the primary flow motion along the longitudinal direction, it generates a 3D helical or spiral motion, characteristic of flow in bends. The inward motion near the bed transports sediment from the outer bank, where scour occurs, towards the inner bank, where deposition occurs. Therefore, the magnitude of the bed changes in a bend is linked to the intensity of the secondary motion. The secondary flow causes the direction 8 of the bed shear stress to deviate from the direction arctan(v/w) of the mean depth-averaged flow velocity (Fig. 11). A 2D depth-averaged model cannot simulate this 3D effect. Therefore, semi-empirical secondary flow corrections must be introduced assuming a locally fully developed curved spiral flow (Rozovskii 1957, Struiksma etal. 1985) ( arctan - arctan A — I rc) The term Ss = -arctan(Ah/rc) is the deviation of the bed shear stress relative to the streamlines, caused by the secondary flow (Fig. 11). A is a parameter of order 10 and can be estimated as a function of the von Karman constant xrand Chezy's roughness coefficient. 64 (3.6) mean velocity x bed shear stress streamline y Figure 11. Secondary circulation and deviation of bed shear stress along bends. According to equation (3.5) the intensity of the secondary motion depends on the flow curvature parameter h/rc. Higher values of h/rc will tend to increase both the scour and deposition along the bend. Unfortunately, equation (3.5) is strictly applicable to weakly or moderately curved bends only; for strongly curved bends, it tends to overestimate the effects of the secondary motion (Blanckaert and Graf 2004). 65 3.1.3 Curvature and inertial adaptation The local radius of curvature of the streamlines is defined as (Larson and Hostetler 1990) - 3 U ^ = 7 ^ 7 7 (3-7) Assuming steady state (duldt = 0, dv/dt = 0), equation (3.7) can be approximated as (see Appendix D for details) (u2+v2)3'2 du idu 2 dv dv -uv v — + u vuv — dx dy dx dy (3.8) Equation (3.8) has the advantage of being applicable to channels of varying curvature, i.e. it is not limited to bends of constant radius of curvature (circular bends). Starting from a flat bed, the curvature will slowly adapt to the final bed topography. Equation (3.8) can be applied at every time step during the simulation to dynamically compute the curvature; however, that process is time consuming and might lead to numerical instabilities i f the velocity gradients vary suddenly. An alternative approach is to apply equation (3.8) only at the beginning of the simulation to compute the initial curvature. Since that initial curvature lags behind the final value, an inertial adaptation equation is required (Rozovskii 1957, Struiksma et al. 1985) C g 4u2 +v2 dx r , \ \ r c j + V-dy f 1 ^ \Tc J J + — = \u x ac \u\ (3-9) 66 The calibration parameter varies normally between 0.4 and 2.0. Equation (3.9) reduces to equation (3.8) i f the spatial gradients of curvature vanish (e.g. fully developed bend flow). Therefore, equation (3.9) affects mostly the areas where curvature changes, such as the entrance and exit of a bend; its main effect is to shift the computed profiles towards the downstream direction (see Appendix D). Although the radius of curvature remains constant during the simulation, depth and velocity do not; thus, the bed shear stress direction S changes dynamically in time according to equation (3.5). An example of the resulting vector field in a flat bed curved channel can be seen in Figure 12. Mean flow Velocity Sediment transport / / / / / / / / , / / / / / / / / / / / / / / / / / A . / / / / / / / / / / / / / / / / ' / / / / / J / / / ! / / / / / i / / / / / / / / / y y y / „ / / / / Figure 12. Example of vector fields of depth-averaged velocity and bedload sediment transport, computed in a flat bed using the secondary flow correction given by equation (3.5). 67 3.1.4 Influence of transverse bed slope Consider a grain of sediment of submerged weight W moving along a straight flume with a transverse slope S, corresponding to lateral bed inclination angle <f> (Fig. 13). The grain will be subject to the forces generated by the bed shear stress in the direction of the mean flow u and the component of the weight in the down slope direction Wtan (j>. The resultant of both forces forms an angle a with the x-axis, which is the direction the grain will finally move (i.e., the direction of bedload transport qs). Figure 13. Influence of the transverse slope St on the direction or of sediment transport. The force caused by the bed shear stress is kjD2 where kj is a parameter that depends on the shape of the grain (&/ = TC/4 for a single spherical particle). Similarly, the submerged weight is W = pg(Gs-l)k2D3, where is another shape parameter (fo - n/6 for a single spherical particle). The angle a of the resultant force is then pg(Gs -\)k2D' A 1 dzb tancc-^- feV s ' 2 tan^ = b-k,D\ ^ f,r* dy (3.10) 68 where p is the density of water and fs = k]/k2 (fs = 3/2 for a single spherical particle) is a shape factor that normally ranges between 1 and 2 (Struiksma et al. 1985). For a 2D case with secondary flow, equation (3.10) becomes (Koch and Folkstra 1981, Struiksma et al. 1985, Talmonetal. 1995) sin<5 tan a = 1 dzb cos<5 -fsr* dx (3.11) In an alluvial bend starting from a flat bed, sediment will be transported by the secondary flow as shown in Figure 12, scouring the outer bank and creating a transverse slope similar to Figure 13. As the lateral slope grows, transverse sediment transport will start to decrease up to the point when an equilibrium condition is reached. At this equilibrium state, the up-slope bed shear stress forces are balanced by the down-slope gravitational forces; then the maximum values for both scour and deposition are achieved. 3.1.5 Objectives The main objective in this chapter is to assess the capabilities of the model to simulate the bed elevation changes in alluvial curved channels. Three laboratory experiments were selected for comparison; two of them have a weak curvature, while the other has a strong curvature. As expected, the computed bed profiles agreed better with the measurements of the weakly curved channels. 69 3.2 E X P E R I M E N T A L D A T A A N D N U M E R I C A L M O D E L Three experiments of bed deformation in alluvial channels (Struiksma et al. 1985) were selected to assess the numerical model. These experiments were performed in weakly and strongly curved alluvial channels composed of uniform sand with a specific gravity Gs = 2.65. Table 2 summarizes the main characteristics of the experiments. Two experiments, called T l and T2, were performed in the Delft Hydraulics Laboratory (DHL) 140° curved flume, which is a relatively long {LJB - 29.32) and weakly curved channel (h/Rc = 1/150 or 1/120), as shown in Figure 14. The other experiment, T6, was performed at the Laboratory of Fluid Mechanics (LFM) 180° curved flume, which is a channel with a relatively shorter (LJB = 7.85) and narrower bend (B/h = 8.5) with a much stronger curvature (h/Rc = 1/21), as shown in Figure 14. These experiments were selected because they provide the opportunity to assess the influence of channel curvature on the performance of the numerical model. Several researchers (e.g., Koch and Flokstra 1981, De Vriend and Struiksma 1983, Struiksma et al. 1985, Struiksma 1989, Johannesson and Parker 1989a and Kassem and Chaudhry 2002) have used some of these experiments for their analytical or numerical models. Most of those models were based on curvilinear coordinates and structured grids. In contrast, the present model uses unstructured Finite Element meshes made of triangular elements in a Cartesian coordinate system. The meshes of both flumes had the computational nodes spaced roughly 0.20 m apart. The D H L flume mesh was made of 6051 elements connected in 3456 nodes, while the L F M flume mesh was made of 3831 elements with 2150 nodes (Figure 15). The boundary conditions for the hydrodynamic model were constant discharge in the upstream inflow section, and constant water surface elevation in the downstream outflow section. To avoid sediment flux across the sidewalls, the direction of sediment transport was set equal to the direction of the mean flow at the lateral sidewalls: tana = v/u (see Fig. 12). The bed elevation of both the inflow and outflow sections was kept constant (dziJdt = 0). 70 To minimize the influence of boundary conditions, the length of the straight reaches upstream and downstream of the bends were artificially increased in the numerical model. The adopted lengths were 15 m for the DHL flume and 10 m for the L F M flume (Figure 14). Table 2. Main parameters of the experiments in alluvial bends performed at the Delft Hydraulics Lab (DHL), the Laboratory of Fluid Mechanics (LFM) curved channels (Struiksma etal. 1985). Parameter Symbol D H L Flume T l D H L Flume T2 L F M Flume T6 Discharge Q (mJ/s) 0.047 0.061 0.170 Flume width B(m) 1.5 1.5 1.7 Water depth h(m) 0.08 0.10 0.20 Flow velocity U(m/s) 0.39 0.41 0.50 Water slope 0.236 0.203 0.180 Chezy coefficient C (m1 /2/s) 28.4 28.8 26.4 Median grain diameter D50 (mm) 0.45 0.45 0.78 Sediment transport qs (m2/s) 7.2 x 10"6 6.9 x 10"6 13 x 10"6 Shield's parameter r * ( - ) 0.26 0.27 0.28 Bend radius Rc(m) 12 12 4.25 Bend length Lc(m) 29.32 29.32 13.35 Relative bend length Lc/B (-) 19.55 19.55 7.85 Width-depth ratio B/h (-) 18.75 15.00 8.50 Depth-radius ratio h/Rc(-) 1/150 1/120 1/21 71 Figure 14. Delft Hydraulics Lab (DHL) 140° curved flume and Laboratory of Fluid Mechanics (LFM) 180° curved flume. Straight reaches up/downstream of bends are 15 m for D H L flume and 10 m for L F M flume. Figure 15. Details of the computational meshes for the DHL and L F M curved flumes. 72 3.3 R E S U L T S Following the experimental conditions, all the simulations started from an initial flat bed and continued until bed equilibrium was reached (i.e. negligible bed changes observed, dzbldt « 0). A constant shape factor value offs = 2.0 was used for the three experiments, as suggested by Kassem and Chaudhry (2002). However, the values of the inertial adaptation parameter were different for the three experiments: /?= 2.0 for T l , J3= 1.0 for T2 and /?= 0.4 for T6. Figure 16 shows the comparison between the measured and computed equilibrium bed profdes for D H L flume, experiments T l and T2. Similar results for flume L F M , experiment T6, are shown in Figure 17. The longitudinal bed profiles shown are located some distance away from the right and left banks: 0.375 m (0.255) for the DHL flume and 0.34 m (0.205) for the L F M flume. The relative depth in the ordinate of Figures 16 and 17 represents the equilibrium water depth relative to the initial depth. Values of relative depth larger than 1.0 represent scour, while values smaller than 1.0 represent deposition. Along the bend, scour is observed in the outer (concave) bank and deposition in the inner (convex) bank. However, that pattern switches in the straight reach downstream from the bend, where decaying alternate bars develop. The numerical model was able to successfully reproduce that alternate bar pattern. The best overall results were obtained for experiment T2. The agreement between the measured and computed profiles is excellent, the magnitude and location of the bed oscillations is well reproduced. This experiment was also more numerically stable with relatively large time step At - 120 s. At the beginning of the simulation the maximum bed shear deviation angle was Ss = 5°, proving that curvature effects were not very intense. 73 Figure 16. Comparison between the computed and measured longitudinal bed profiles for DHL flume, experiments T l and T2. Profiles are 0.375 m from right and left banks. 74 Distance / Width 0 2 4 6 8 10 0 i ' ' ' ' 1 Figure 17. Comparison between the computed and measured longitudinal bed profiles for L F M flume, experiment T6. Profiles are 0.34 m from right and left banks. During the numerical simulation of experiment T l , a bar was observed migrating from downstream of the first point bar and towards the bend exit. This migrating bar finally ended up merging with the bar located at the bend exit. For T l , a smaller time step At = 100 s was used to prevent instabilities. At the beginning of the simulation the maximum bed shear deviation angle was Ss = 4° Experiment T6 was simulated using a time step At = 60 s. The initial bed shear stress deviation angle reached a high maximum value 4 = 37° (see Fig. 12), indicating that the effects of the secondary flow are very strong in this case, as expected from the corresponding high values of the flow curvature (h/Rc = 1/21). It is also observed in Figure 17 that some degradation is predicted downstream of the bend, which has not been observed in the experiments. 75 In general, the time step At was selected such that the maximum bed change Azt in any given time step is smaller than 5% of the original water depth. That usually prevents the model from becoming unstable because of large changes in the bed elevation during the computations. 3.4 DISCUSSION The model seems to reproduce the mean features of the bed elevation changes along the three alluvial bends simulated. Notably, the overshooting of the bed profile downstream of the bend entrance, where maximum scour and deposition occur (Struiksma et al. 1985, Parker and Johannesson 1989, Struiksma 1989) is correctly predicted. Although it is unknown if a migrating bar was observed during the experiment T l , there are experimental observations of bars migrating along bends. Whiting and Dietrich (1993) reported migrating bars in a large amplitude meandering bend with a width to depth ratio of B/h = 17. However, the bars remained static when B/h = 15. Those results seem to agree with our numerical findings, as we observed a migrating bar for experiment T l (B/h = 19), but not for experiments T2 (B/h = 15) or T6 (B/h = 9). In general, the inclination of the transverse bed slope S, oscillates along the bend. From an initial flat bed upstream (St = 0), the transverse slope tends towards the bend equilibrium transverse slope (St —> Ste). However, depending on the flow and sediment conditions, the local value of the transverse slope may overshoot (St> Ste) downstream of the entrance (Struiksma et al. 1985, Parker and Johannesson 1989, Struiksma 1989), causing a large point bar on the inner bank and a deep pool on the outer bank, which can be observed looking at the measured profiles in Figures 16 and 17. Around the center of the bend, preaches a minimum value, to increase again around the bend exit. A l l those features were captured in the model, especially for the D H L flume (Fig. 16). For the L F M flume, the minimum value of S, around the mid-bend was not clearly reproduced (Fig. 17). 76 Quantitatively, the magnitude of the maximum scour and deposition seems to be accurately predicted for experiment T2 and T l . For T6, the model slightly over-predicts the height of the point bar, while it under-predicts the maximum scour. The results suggest that, for strongly curved channels, sediment transport is underestimated along the outer bank (smaller scour) and overestimated along the inner bank (higher point bar), which is expected since the simple secondary flow model (equation 3.5) neglects the feedback mechanism between the primary and secondary flows. The secondary flow not only re-circulates mass, but also redistributes momentum from the primary flow. Since the top, fast moving layers of the flow move outwards, velocity along the outer bank increases, while it decreases along the inner bank. 2D depth-averaged models cannot reproduce such momentum redistribution. An alternative way to enhance the capabilities of a 2D vertically-averaged model for flow in bends is the application of the moment of momentum (MOM) approach (Ascanio and Kennedy 1982, Jin and Steffler 1993, Ghamry and Steffler 2003, Ghamry and Steffler 2005)..Ghamry and Steffler (2005) have demonstrated that solving additional equations for the M O M allows a 2D depth-averaged model to simulate better the velocity redistribution of the primary velocity, but also the tendency of the secondary motion over the vertical profile. The application of a vertically-averaged and moment of momentum model is described in the next chapter. The results of the model are not perfect but probably accurate enough for most practical engineering applications. The L F M flume represents an extreme case of high curvature. In natural rivers, the radius of curvature is usually several orders of magnitude larger than the flow depth (Yalin 1992), leading to very small values of h/Rc, even smaller than those of the D H L flume. Therefore, the morphological model should provide good results for applications to real meandering rivers. 77 3.5 S U M M A R Y The bedload sediment transport model introduced in the previous chapter was extended to simulate the bed level changes in alluvial bends by incorporating effects of the secondary flow and transverse slope on the direction of bedload. Unlike most previous analytical or numerical formulations (e.g., Koch and Flokstra 1981, De Vriend and Struiksma 1983, Struiksma et al. 1985, Struiksma 1989, Johannesson and Parker 1989 and Kassem and Chaudhry 2002), the new model is based on a flexible unstructured mesh, which has the advantage of adapting to practically any planform geometry. Since the model can compute the local radius of curvature of the streamlines from the flow field, it can be applied to real meandering rivers of varying curvature (simpler models usually assume the meandering river as a series of connected circular segments of constant curvature). The model has successfully reproduced the equilibrium bed profiles of three experiments in alluvial bends. The main features of the bed profiles, such as the overshooting downstream of the bend entrance and the alternate bar pattern downstream were correctly predicted. The model provides better results for weakly curved bends, where the effects of secondary currents are less intense. However, this is a limitation of any 2D vertically-averaged model with a simple secondary flow sub-model that corrects the direction of the bed shear stress, but not its magnitude, which is also affected by the redistribution of momentum caused by the secondary flow. An alternative formulation using the vertically-averaged and moment of momentum model (which does not require an additional secondary flow correction) is presented in the next chapter. As the model is currently limited to uniform sediment, it cannot simulate bend sorting. That can constitute a topic for further research. 78 4 VERTICALLY-AVERAGED AND MOMENT OF MOMENTUM MODEL 4.1 I N T R O D U C T I O N Natural rivers are seldom straight along reaches farther than a few channel widths. Meandering is a clear example of the tendency of rivers to follow curved paths. The main difference between flow in straight and curved channels is the presence of a centripetal force in the latter. The centripetal forces along river bends induce a secondary flow in the transverse direction to the primary longitudinal flow, an effect also referred as helical or spiral motion. This secondary flow causes redistribution of velocities, boundary shear stresses and sediment transport along the bends, shaping the morphology of alluvial meandering rivers. The helical motion is in essence a three-dimensional (3D) flow that modifies the velocity distribution in both horizontal and vertical planes. In a fully developed bend flow, the horizontal velocity is higher along the outer bank than along the inner bank (i.e. a transverse gradient of horizontal velocity develops); while the vertical velocity profile becomes more uniform (i.e. the profile flattens) when compared to a straight channel. The momentum 79 advected by the secondary flow in the transverse direction (Fig. 18) causes these redistributions of velocity. transverse velocity A h outer inner bank bank Figure 18. Redistribution of the transverse horizontal velocity and flattening of the vertical velocity profile caused by the secondary flow in bends. Another effect of the secondary flow is on both the magnitude and direction of the bed shear stresses. Along bends, the bed shear stress deviates from the main longitudinal direction of the primary flow, pointing towards the inner bank of the bend. At the same time, the magnitude of the bed shear stress is locally amplified compared to its value along straight channels. The 3D structure of the helical motion can be resolved numerically via a full 3D flow model. Although 3D models are presently available, they are computationally too expensive for most practical river morphology applications. For that reason, two-dimensional vertically-averaged (VA) models still dominate this application field. However, V A models have a serious drawback: the vertical structure of the flow velocity is lost in the depth-averaging process. Therefore, V A models intended for river morphology must rely on so-called secondary flow 80 corrections in an attempt to take into account the effects of the helical motion on the direction of the bed shear stress, and hence on sediment transport. Traditional secondary flow corrections have several limitations: They assume fully developed bend flow, neglect the feedback between the primary and secondary flows, require information about the radius of curvature of the streamlines and need to account for the inertial adaptation of the secondary flow. These simple models can estimate the deviation in the direction of the bed shear stress, but not the amplification of its magnitude. However, they have been successfully applied in many cases (e.g. Struiksma et al. 1985, Struiksma 1985, Nagata et al. 2000, Kassem and Chaudhry 2002, Darby et al. 2002), but seem to be limited to small or moderate curvatures, when the secondary flow is not very strong (see Chapter 3). Recently, Blackaert and de Vriend (2003) and Blanckaert and Graf (2004) have introduced a nonlinear secondary flow correction for V A models, applicable even for strong curvatures, that incorporates the secondary flow effects on both magnitude and direction of the bed shear stresses. This nonlinear model seems promising, but it has not yet been tested on alluvial beds. Furthermore, it requires an additional semi-heuristic model for the inertial adaptation of the flow and seems to be limited to circular bends. An alternative approach to recover some of the information of the vertical structure of the flow is the application of the conservation of moment of momentum (Ascanio and Kennedy 1983, Johannesson and Parker 1989b, Yeh and Kennedy 1993a, Yeh and Kennedy 1993b, Jin and Steffler 1993, Ghamry 1999, Ghamry and Steffler 2002, Ghamry and Steffler 2005). Jin and Steffler (1993) introduced a V A model with two extra moment of momentum equations to allow the calculation of the effect of secondary flow on the depth-averaged velocity field. Ghamry (1999) developed the 2D vertically-averaged and moment (VAM) model by adding extra equations for the conservation of moment of momentum. This V A M model can be seen as a quasi-3D model because it can approximate, at the simplest level, some of the 3D structure of the helical motion which are beyond reach for conventional V A models. 81 Recently, Ghamry and Steffler (2005) demonstrated that a V A M model can predict the lateral redistribution of velocity in curved channels and the vertical velocity profile (Figure 18) significantly better than the conventional V A model. They recommended using the V A M for modeling strongly curved open channels. However, all previous applications of the V A M in curved channels implied non-erodible flat beds without sediment transport. The capabilities for bed morphology simulations with varying topography have not been tested before. 4.1.1 Objectives The main objective of the study described in this chapter is to verify i f the V A M model can be applied to simulate the morphological bed evolution in curved alluvial channels without a secondary flow correction. Initially, the model was tested using flume data from two well-known experiments (Struiksma et. al 1985) in curved alluvial channels, one weakly curved and the other strongly curved, to test the influence of curvature in model performance. Finally, data from The Waal River (Struiksma 1985) were used to test the model for a real scale application. The results showed that the V A M implicitly accounted for the effects of the secondary circulation and was able to successfully predict the equilibrium bed profiles for the three test cases. The V A M model required about 50% more computational time than a conventional V A model. 4.2 T H E N U M E R I C A L M O D E L 4.2.1 The V A M hydrodynamic model The V A M hydrodynamic model is derived by vertically averaging or integrating the 3D Reynolds equations after multiplying them by the vertical coordinate. Ghamry (1999) developed a general 2D model of 10 equations that can assume linear or quadratic vertical 82 distribution shapes for the velocity components in the three spatial coordinates plus non-hydrostatic pressure distribution. The derived equations represent a quasi-3D model with more vertical detail than conventional Saint Venant V A 3-equation formulation (one continuity and two momentum equations, see equations 2.1 through 2.3). The general V A M model can be solved for 3, 5, 8, 9 or 10 equations depending on the desired level of vertical detail (Ghamry 1999, Ghamry and Steffler 2002, Ghamry and Steffler 2005). Ghamry and Steffler (2002) compared the V A M 10-equation model with the simpler V A M 5-equation model for flow in curved fixed bed channels and concluded that there is no significant difference between the two of them. Therefore, the V A M 5-equation model was adopted for this study. The term " V A M " in this chapter is used exclusively to refer to the V A M 5-equation model (Fig. 19). Adopting the x-y plane in the horizontal and z direction forming the vertical coordinate, the V A M assumes the linear distribution of velocities u(z) and v(z) shown in Figure 19. u(z) = (u0- Uj) + 2 Uj(z-zb)/h (4.1) v(z) = (v 0 -v,) +2v1(z-zb)/h (4.2) Where zb is the bed elevation, h is the water depth, u(z) is velocity in the x-direction, v(z) is velocity in the ^-direction, u0 and v0 are the depth averaged velocity components, ui and v/ can be interpreted as the velocities at the water surface (z = zb +h) in excess of the means u0 and v0 respectively. A finite slip velocity (u0 - u,, v0 - v7) is assumed at the bed. Having applied the previous distributions, the following equations are left: The vertically averaged continuity equation: ^ + * k + ^ = 0 (4.3) dt dx dy 83 t z Figure 19. Linear vertical distributions of velocity assumed by the vertically-averaged and moment (VAM) model over the water depth h: u{z) in the x-direction and viz) in the y-direction. The vertically averaged momentum equation in the x-direction: dax , d dt dx V " J + -qxq J 1 + -3 d(hux ) d(huxvx) v dx dy + gh — (h + zb) dx 84 1 dhar 1 dhrxv ] p dx p by *» + - T x _ =0 (4.4) The vertically averaged momentum equation in the >>-direction: dt dx ylx f „2\ + -dy + • <9x <9y l ^ _ _ L ^ + _Lr = 0 p dy p dx p y'h (4.5) The moment of momentum equation in the x-direction: du. rqxu^ dt dx + v, 5y v « y 5w, 3 T ~ d ) 7 + 2 4crv 5z 4r^ dz 4r V hp dx hp dy hp hp (4.6) The moment of momentum equation in the v-direction: dt d_ + w, Sx + /z 5x 2 4<rr dz 4ryx dz 4r. + • hp dy hp dx hp hp (4.7) = 0 Where qx= u0h is the flow discharge in x-direction per unit width, qy = v0h is the flow discharge in the ^-direction per unit width, z is the mid channel depth, f is the vertically averaged total turbulent shear stress, T is the bed shear stress, a is the total turbulent normal stress, p is water density and g the gravitation acceleration. The above set of equations (4.3) through (4.7) is called the (5-equation) V A M model. The vertically averaged total turbulent shear and normal stresses appearing in equations (4.4)-(4.7) 85 are approximated according to the Boussinesq model (assuming laminar stresses are negligible) as follows: d (" ^ dx qx (4.8) ° y =*yy=~n \*yydz = 2P^ ^ V h J (4.9) T = T xy " yx = — r dz Pv> v h j dy + • dx V h J (4.10) T — T xz zx zh+h h J xz PV2 'du^ \dz j = pv2 \ n J (4.11) T = T yz zy = - \ryzdz = pi \dz j PVz („ \ \ n J (4.12) Where vh is the vertically averaged turbulent exchange coefficient or eddy viscosity in the horizontal direction (x-y plane); and vz is the vertically averaged turbulent eddy viscosity in the vertical direction. For simplicity, the case of bed-dominated turbulence is assumed and values of the order of vh = 0.5u,h and vz - 0.07 u,h are used, u* being the shear velocity defined as: \f- \ 2 + (T V (4.13) V r J V r J 86 The bed shear stresses, appearing in equations (4.4) - (4.7) are approximated according to: T (4.14) r (4.15) * Where C* is the dimensionless Chezy Coefficient and is related to the effective roughness height, ks, through: The closed system defines a model for 2D or a quasi-3D flow in open channels. This model solves for the dependent variables h, qx, qy, ut, and v,. This results in a 5-equation by 5-unknown (5x5) model (the V A M 5-equation model). If ui and v/ are forced to zero and their corresponding equations (4.6) and (4.7) are eliminated, then a 3-equation by 3-unknown (3x3) model is left, equations (4.3) - (4.5), which is the traditional St. Venant or V A model (similar to equations 2.1 - 2.3). The main advantage of the V A M model for the simulation of flow in bends is its ability to simulate the helical flow motion in the bend without resorting to an empirical secondary flow correction (as used in Chapter 3 for the V A model). The helical (spiral) motion in a bend is the result of the presence of a primary flow and a secondary flow. The primary flow is the main flow moving in the downstream direction driven mostly by gravity. The secondary flow is a (curvature induced) transverse circulation in which the surface water moves towards the outer bank of the bend, while the near-bed water moves in the opposite inward direction (Figure 20). C. =5.751og 12 — (4.16) 87 The near-bed velocity (and corresponding bed shear stress) transport bedload sediment in the transverse direction from the outer bank, where erosion occurs, towards the inner bank, where deposition occurs. The deviation of the near-bed velocity and bed shear stress directions from the direction of the primary flow is one of the more important effects of curvature on sediment transport, and hence on erosion and deposition along the bend. a. surface velocity b. near-bed velocity Figure 20 Simulation of the helical motion in a bend by the V A M model: (a) surface velocity moving towards outer bank, and (b) near-bed velocity moving towards inner bank. 4.2.2 The bedload sub-model The sediment transport sub-model solves the sediment continuity or Exner equation for bed load ( l - A ) ^ + ^ + ^ = 0 (4.17) dt dx dy 88 Zb is the bed elevation (Fig. 19), X is the porosity of the bed material (X = 0.4 was adopted for all simulations), qsx = qscosa and qsy = qssma are the components of the volumetric rate of bedload transport per unit width in the x- and ^-direction, respectively. Angle a is the direction of bed load transport in the horizontal x-y plane. 4.2.3 Bed slope effect When the bed is flat, the direction a of bedload transport coincides with the direction 8 of the bed shear stress. However, in sloping beds the gravity force induces a deviation between a and 8 which can be taken into account by (Struiksma et al. 1985) fs is a shape factor that ranges between 1 and 2; dz^dx and dzbldy are the bed slopes in the x-and y- directions respectively. 4.2.4 Secondary flow effect on the bed shear Up to this point, equations (4.17) and (4.18) are common to both the V A and V A M formulations. The difference between them arises in the direction 8 of the bed shear stress. The V A model computes the water depth h and depth-averaged velocities components u and v in the x and y directions, from which 8 can be computed by applying an equation derived for a locally fully developed curved spiral flow (Rozovskii 1957, Struiksma et al. 1985) tanor = cos 8 -sin<5 - 1 dzb / > * dy 1 dzb fsr* dx (4.18) 89 SVA = arctanl -arctan A — I rc ) (4.19) The term tan" l(Ah/rc) is the deviation of the bed shear stress due to the spiral flow effect; rc is the local radius of curvature of the streamlines. The parameter A is a function of von Karman constant /rand roughness C* I — (4.20) The V A M model computes h, UQ and VQ, plus uj and v/ (Fig. 19) via the moment of momentum approach (Ghamry and Steffler 2002, 2003). The bed shear stress direction can be computed using equation (4.21) without information of the radius of curvature, but introducing the calibration parameter kb. 'VAM = arctan rv0 -kbvx ^ (4.21) The assumption behind equation (4.21) is that the components Tbx and Tby of the bed shear stress have the same direction as the near bed velocity components: tan<5= (rby I rbx) = (vb I ub). From a sensitivity analysis (Vasquez et al. 2005) the value kb = 0.8 was adopted. Along straight reaches the depth-averaged velocity components (u,v) computed by the V A model are identical to the depth-averaged components (uo,vo) computed by the V A M model. However, along curved channels, because of velocity redistribution generated by the secondary circulation (Figure 18), those values are no longer identical. Therefore, the magnitude of the depth-averaged bed shear stresses along bends differs from one model to the other. The V A M 90 model uses equations (4.14) and (4.15) to compute the bed shear stress; the V A uses identical equations but replacing (uo,vd) by (u,v) (see equation 2.6). Attempts in V A M model to use the near-bed velocity information (ub,Vb) instead of the depth-averaged values (uo,vo) to compute bed shear stresses and sediment transport rates did not provide any improvement. Therefore, equation (4.21) was applied for the direction of the bed shear stress, but equation (4.22) was applied for the dimensionless bed shear stress (Shield's parameter) « * - , " ° + " ° (4.22) C ' ( G , - 1 ) D 5 0 C'(G,-\)D. 50 1/2 C = C*g is the Chezy roughness coefficient, GS is the specific gravity of the sediment, D5Q is the median grain diameter and U the depth-averaged velocity. 4.3 T E S T C A S E S A N D N U M E R I C A L M O D E L 4.3.1 Experimental data 4.3.1.1 Laboratory flumes Three test cases were selected for numerical simulation. Two of them were experimental curved flumes (Struiksma et al. 1985), while the other was a real meandering river (Struiksma 1985). Their main characteristics are summarized on Table 3. 91 Table 3. Main parameters of the DHL and L F M flumes (Struiksma et al. 1985) and the Waal River (Struiksma 1985). Parameter Symbol D H L Flume L F M Flume Waal River Discharge Qim'/s) 0.061 0.170 1350 Flume width B(m) 1.5 1.7 250 Water depth h(m) 0.10 0.20 5.0 Flow velocity U(m/s) 0.41 0.50 1.08 Water slope S(%) 0.203 0.180 0.012 Chezy coefficient C (m1 /2/s) 28.8 26.4 44 Median grain diameter D50 (mm) 0.45 0.78 1.3 Sediment transport qs (m2/s) 6.9 x 10"6 13 x 10'6 4.3 x 10"4 Shield's parameter r * ( - ) 0.27 0.28 0.28 Bend radius Rc(m) 12 4.25 >1288 Width-depth ratio B/h (-) 15.00 8.50 50.0 Depth-radius ratio h/Rc(-) 1/120 1/21 < 1/260 One of the experiments was performed at the Delft Hydraulics Laboratory (DHL) 140° curved flume, which is a weakly curved channel, as shown by its low flow curvature h / Rc = 1/120 (see Fig. 14). The sediment transport equation adopted for this channel was the Engelund-Hansen equation qs = g(Gs-l)D350 (t*f2 * 0.00057*75 (4.23) 92 The other experiment was performed at the Laboratory of Fluid Mechanics (LFM) 180° curved flume, which is a channel with a stronger flow curvature h / Rc = 1/21 (see Fig. 14). Initially the Engelund-Hansen equation was also applied for this flume, but it did not provide good results (Engelund-Hansen equation 4.23 is proportional to mean velocity raised to power 5). Instead, an empirical equation proportional to the depth-averaged velocity raised to power 4 was adopted for this flume 4.3.1.2 Meandering river The Waal River is the main branch of the Rhine River in the Netherlands. When the Rhine River flows into the Netherlands, it splits into two main branches, the Nederrijn to the North and the Waal to the South. The Waal River is one of the busiest navigation arteries in the world (Jansen 1979). Struiksma (1985) computed using a 2D V A model the equilibrium bed topography of the Waal River. The study reach covers a meandering reach of about 10 km divided in 10 arcs (Fig. 21). The length of the arcs varies between 0.6 and 2 km, with radii of curvature larger than 1.3 km. Two data sets of bed soundings 50 m away from both banks were available for comparison with the numerical model. The sediment transport was computed using a modified Meyer-Peter Muller equation (Struiksma 1985) qs = 0.000205c/4 (m2/s) (4.24) qs = 0.00078(0.65*72 -0.047) 3/2 (m2/s) (4.25) 93 The simulation was carried out with a steady flow of 1350 m3/s (Table 3). Identical conditions were assumed for the simulations with the V A M model. Figure 21 Study reach of the Waal River (after Struiksma 1985). 4.3.2 Numerical model The computational meshes for the D H L and L F M flumes were the same previously employed in the morphological simulation applying the V A model (see Chapter 3, Fig. 14 and 15). The computational mesh for the Waal River was made of nodes roughly 50 m apart, giving a total of 2266 nodes and 3622 triangular elements. Figure 22 shows a detail of that mesh around a bend. The boundary conditions for the hydrodynamic model were constant discharge in the upstream inflow section and constant water surface elevation in the downstream outflow section. In the upstream inflow section u; and v/ were approximated to the values for fully developed turbulent flow (logarithmic velocity profile) ui = UQ I 2.5 and v/ = vo I 2.5. To avoid sediment \ 94 flux across the sidewalls, the direction of sediment transport was set equal to the direction of the mean flow at the lateral sidewalls: tana = vo I Uo- The bed elevation of the inflow section was kept constant. Figure 22. Detail of the computational mesh for the Waal River (lower bend). The value of the shape factor was set to fs = 2.0 for the flume experiments and fs = 1.25 for the Waal River. The time step for the Waal River was 3 h, and the final simulation time 100 days. 4.4 R E S U L T S Figure 23 shows the measured longitudinal bed profile along the DHL flume; as well as the profiles computed by the V A M model. The results of a V A model from a previous application (see chapter 3) are also shown for comparison. The results of the V A M agree well with measured data; it correctly predicted the "overshoot" of the bed profile immediately downstream of the bend entrance. In addition, the alternating bed profiles downstream of the bend (switching from erosion to deposition and vice versa) were correctly reproduced. The V A M model accurately predicted the beginning and end of the erosion and deposition bed profiles. In general, the difference between the results of both V A and V A M models was small. 95 However, the V A model seems to reproduce better the length of the pool and point bar near the bend entrance. Distance / Width 1.4 Figure 23. Longitudinal bed profiles computed by V A M model for D H L experiment (VA results shown for comparison). Figure 24 shows the measured longitudinal bed profile along the L F M flume, as well as the profiles computed by the V A M models (again, previous V A results are shown for comparison). The results of V A M model agree reasonably well with the measured data. The V A M predicted the overshoot downstream of the bend entrance, but also the second oscillation of bed profile around the bend exit, with a minimum transverse bed slope around the middle of the bend, although the predicted wave length of the oscillation was not very accurate. 96 Distance/Width 4 6 8 10 Figure 24. Longitudinal bed profiles computed by V A M model for L F M experiment (VA results shown for comparison). The results for the Waal River are shown in Figures 25 and 26. Figure 25 shows the comparison between the computed profiles along both right and left margins, with the bed sounding made in 1969. The agreement with measured data is very good. The location of the areas of scour (pools) and deposition (point bars) along the banks is well reproduced as shown also in Figure 26. Starting from an initial water depth of 5 m everywhere, the model predicts equilibrium depths varying between 2.7 and 7.6 m. 97 Figure 25. Comparison between results computed by V A M model with 1969 soundings along profiles 50 m away from banks of the Waal River. 98 Depth WAAL RIVER SIMULATED BY VAM MODEL Qin = 1350.0000 Figure 26 Plan view of the computed water depths for Waal River showing the pools along outer banks and point bars along inner banks (initial depth was 5 m). 4.5 DISCUSSION 4.5.1 Prediction capabilities of the V A M model The computed bed profiles (Figure 23 through 26) demonstrate that the V A M model is capable of simulating the bed deformation in alluvial curved channels subject to helical motion without the need of an external secondary flow correction. The additional information of the near bed velocity provided by solving two additional moment of momentum equations allows the V A M model to estimate the direction of the bed shear stress, which is vital to predict the scour and deposition patterns along bends. 99 Figures 23 and 24 show that the mean features of the bed profile along the bends, such as the overshoot downstream of the bed entrance and the alternate bars downstream of the bend exit, are correctly simulated by the V A M model. Figures 25 and 26 prove that the V A M model is not only suitable for circular laboratory flumes, but also for real scale meandering rivers with channels of varying curvature. The model accurately predicted the locations of areas of scour (pools) and deposition (point bars) along the banks of the river. In contrast with conventional V A models with secondary flow corrections, the V A M model does not require explicit information about the radius of curvature of the streamlines. The effects of curvature on the direction of bed shear stresses are implicitly accounted for by computed velocity field near the bed. This computed flow field also incorporates the effects of the inertial adaptation of the secondary flow; no additional equation for this effect is needed. It can be argued that most of the relevant physics of the curved flow is embedded in the V A M equations. Only an empirical coefficient kb is needed to correct the near-bed velocity estimate, probably because of the crude approximation provided by the pre-assumed linear profiles of vertical velocity. Similar to parameter A in the secondary flow equation (4.19), the coefficient kt is probably dependent on bed roughness. For example, Elgamal (2002) developed an empirical equation for kb as a function of relative roughness h/ks applicable for flow over large ID dunes (with flow separation and reattachment in the vertical plane). Unfortunately, the experimental data used in this study do not provide information on the direction of bed shear stress to develop a similar expression. Nevertheless, the adopted value kb = 0.8 turned out to be adequate for both flume and river applications. It seems that the prediction capabilities of the V A M deteriorate with increasing curvature. The computed results for the weakly curved D H L flume and Waal River are better than those of the strongly curved L F M flume. 100 4.5.2 Comparison with V A model The difference between V A and V A M models is small for weak curvatures, as shown by the results of the D H L flume (Figure 23). However, the initial expectation that the V A M should provide significantly better results than the V A model for highly curved channels was not entirely verified, as seen from the results of the L F M flume (Figure 24). It is not completely clear which model performs better for the L F M flume. Both models underestimate the depth of the scour hole (pool) downstream from the bend entrance. The V A M model captures better the oscillation of the bed profile, but it underestimates its wavelength and amplitude near the bend exit. The V A model predicts larger bed changes than the V A M model and correctly predicts the amplitude of scour and deposition at the bend exit, but it erroneously predicts scour along the inner bank, downstream from the bend exit (distance/width around 9 in Figure 24). Blanckaert and de Vriend (2003) state that secondary flow corrections used by V A models, which neglect the momentum transfer between the primary and secondary flows, tend to overestimate the effects of the secondary circulation when curvature becomes strong. That might be the explanation for the higher values predicted by the V A model for the L F M flume. The momentum transfer between the primary and secondary flow (Figure 18) results in an outward increase of the primary velocity (Rozovskii 1957, Johannesson and Parker 1989b, Yeh and Kennedy 1993a). The V A M equations contain momentum convection terms b\huivi)l dx in equation (4.4) and d(huivi)l dy in equation (4.5) that represent the momentum exchange between primary and secondary flows, responsible for the velocity redistribution (Ghamry 1999). Since those terms become important for strong curvatures, it might explain why the bed profiles computed by the V A and V A M models differ when the channel curvature increases, as shown when comparing Figures 23 and 24. By comparing the results of the V A and V A M models (Figures 23 and 24), it seems that the V A model tends to predict a slightly higher point bar and slightly shallower pool, when 101 compared with the V A M model. The difference is too small to be conclusive, but it seems to follow the expected tendency. Given that the V A model neglects the transverse redistribution of velocity (Figure 18), it over-predicts the velocity magnitude along the inner bank (higher point bar), while under-predicts the velocity magnitude along the outer bank (shallower pool). In the present bedload transport application, the potential of the V A M model has only been partially exploited; only the information of the near-bed velocity was used to calculate the direction of the bed shear stress. However, the information of the vertical velocity structure should be relevant for other applications involving suspended sediment transport (Guo and Jin 1999) or transport of pollutants in the main water body (Albers et al. 2002). The surface velocity might be relevant for applications involving floating debris or ice. For more complex applications involving several physical processes (e.g. bedload and suspended sediments plus pollutants), the additional 50% computational expense of the V A M model would be justified by the enhanced prediction capability of the model gained from more detailed information of the flow structure. 4.5.3 Future research The V A M sediment transport model introduced here is currently limited to bedload only. Although suspended sediment transport has not been tested yet, it should be expect that the additional information of the vertical velocity profile provided by the V A M model should be beneficial to improve the estimates of suspended sediment concentration in the water column. Guo and Jin (1999) provide an example of the application of a type of ID V A M model for suspended sediment transport. The bedload sub-model could also be upgraded to handle multiple sediment fractions to simulate effects such as grain sorting or bed armouring. The parameter kb has a strong influence on the computed bed profiles (Vasquez et al. 2005). The value kb = 0.8 adopted in this study led to good results, but it cannot be claimed to be a 102 universal value. Research is needed to study the correlation of kb with relevant hydraulic variables, such as roughness. Experimental measurements of both direction and magnitude of bed shear stress along bends would be very beneficial for that purpose. Although the V A M model with linear distributions of vertical velocity has proven successful for alluvial bend morphology, it is possible to replace it by a V A M model with parabolic velocity distributions. That change would probably require a different correction coefficient kb for the near-bed velocity direction. Whether or not the parabolic distribution is better than the linear one for morphological simulations is unclear. Ghamry and Steffler (2002) have shown that the computed flow field is not very sensitive to the shape of the velocity distribution although they found that, near the walls of curved channels, the VAM-parabolic simulated the secondary flow much better than the VAM-linear. Research is needed to determine whether switching to the VAM-parabolic is beneficial at all. The 2D V A M model has been compared in the present study with a simpler 2D V A model, with successful results. A similar comparison should also be made with a more complex full 3D model. It might be possible that in certain cases the bend curvature is so strong that the 3D effects cannot be properly captured by the simple linear (or parabolic) vertical velocity profile assumed by the 2D V A M model (Figure 19), leaving the full 3D model as the only alternative. That might happen for example, i f more than one strong re-circulation zone is present in the vertical plane. If such limitation for the applicability of the 2D V A M model exists, it needs to be investigated. 103 4.6 S U M M A R Y The 3D information of the vertical structure of the flow in bends (helical motion), which is lost in conventional 2D vertically-averaged (VA) models due to the depth-averaging process, can be partially recovered by assuming a simple linear vertical distribution of horizontal velocity and solving two additional equations of the moment of momentum. This 2D vertically-averaged and moment (VAM) model can reproduce the main effects of the flow in bends in a quasi-3D fashion. The direction of the bed shear stress can be computed by the V A M model without relying on any additional secondary flow correction as needed by V A models. In the present study, the 2D V A M hydrodynamic model developed by Ghamry and Steffler (2003) was semi-coupled with a bedload transport module to compute erosion and deposition along erodible bends composed of uniform sediment (sand). Both the V A and V A M models led to good agreement with experimental data. For weakly curved channels, both models gave similar results. However, when curvature becomes strong the results of the two models start to differ from each other. This is because the moment of momentum equations of the V A M model can capture the momentum transfer between the primary (streamwise) and secondary (transverse) flows along the bend, which is neglected by the V A model. However, a clear superiority of the V A M model over the V A model could not be demonstrated from the results obtained. 104 5 SUMMARY AND CONCLUSIONS 5.1 S U M M A R Y Alluvial rivers are those that shape their own channel with the sediment they are actively transporting. Their morphology is determined by the mutual interaction between the flow of water and the sediment transported by that flow. Three conservation laws can describe the physics of such a morphodynamic process mathematically: conservation of water mass, conservation of flow momentum and conservation of sediment mass. The analytical solution of those equations is feasible only for a very few cases under highly restrictive simplifying assumptions. For most practical applications, those equations must be solved numerically. Traditionally, the Finite Difference (FD) method has been applied successfully to solve the flow and sediment equations (Koch and Flokstra 1980, Struiksma 1985, Kassem and Chaudhry 1998, Kassem and Chaudhry 2002, Darby et al. 2002). However, the main disadvantage of the FD approach is that it applies a structured (regular) mesh, which can be limited for modeling complex geometries. As an alternative, the Finite Element methods provide the possibility of using flexible unstructured meshes; this is the method used in the present research. 105 A morphological model, in two different versions, was developed in this thesis based on two-dimensional Finite Element hydrodynamic models. One of the versions is a conventional depth-averaged or vertically-averaged (VA) model (Steffler and Blackburn 2002), which was applied in Chapter 2 to straight alluvial channels and in Chapter 3 to curved alluvial channels. The other version, uses the vertically-averaged and moment of momentum (VAM) hydro-dynamic model (Ghamry and Steffler 2003), which was applied to curved flumes and a meandering river in Chapter 4. In Chapter 2, the V A model was initially applied to simulate bed aggradation/degradation in straight flumes caused by imposed upstream sediment supply and by transcritical flow (supercritical flow and hydraulic jump). The morphological simulation in transcritical regime is probably one of the first performed with a Finite Element model and indicates the potential for applying the model to mountain river morphology. In Chapter 3, two additional sub-models were added to the V A model to simulate scour and deposition in bends: secondary flow correction and bed slope correction. The model successfully replicated the patterns of scour and deposition along bends, such as point bars and scour pools, and the alternate bars downstream from the bends. In Chapter 4, the V A M model with a bed slope correction was applied to simulate scour and deposition in laboratory bends and a meandering river. No secondary flow correction was applied as the additional moment of momentum equations supplied the necessary information about the secondary flow. This is the first application of a river morphology model based on V A M hydrodynamic equations. 106 5.2 THESIS C O N T R I B U T I O N a) As a direct result of the present study, a Finite Element 2D depth-averaged model for river morphology with the following features was developed: - flexible unstructured FE mesh adaptable to complex geometries, with the possibility of refining around areas of interest or where flow gradients are intense (e.g. hydraulic jumps); - the inflow of sediment can be set as free or forced, the latter being useful for cases involving aggradation or degradation for overload or deficit in the supply of upstream sediment; - one of the first 2D FE river morphology models with transcritical flow capabilities (subcritical - supercritical); - secondary flow effects in bends are taken into account by a linear secondary flow correction sub-model. The curvature of streamlines and its inertial adaptation are also computed by the model. Therefore, it is applicable to meandering rivers of varying curvature; - possibility of selecting different sediment transport equations from the available list (e.g. Engelund-Hansen, Meyer-Peter Muller, Van Rijn) or incorporating additional ones; - the effect of the transverse bed slope on the direction of bed load is incorporated into the model. Currently, the main limitations of the model are: - the model is suitable for sediment transport problems not dominated by convection, such as the ones described in Chapter 2, 3 and 4. The model should not be applied for convection-dominated applications when migrating bedforms appear (see Appendix F); 107 - bedload is computed only for the median grain size, which is valid for uniform sediment. Multiple grain sizes and phenomena such as bed sorting or armouring (e.g. gravel beds) are not simulated at present; however, there is no conceptual obstacle for incorporating non-uniform sediment in the future; - suspended sediment is ignored, i.e., the model is valid only for bedload dominated rivers; - at present, the model assumes all finite elements as completely wet. Partially wet elements have not yet been considered; - all the bed is considered erodible. There is no provision for non-erodible layers, such as rock outcrops; - some problems appear in the sidewalls with sudden changes in direction (see Appendix E). - bank erosion is not considered. The limitations listed above reflect the early stage of development of the model. Future enhancement could overcome those limitations. b) The thesis provides additional empirical evidence that decoupled sediment transport models are not unstable as usually claimed in literature, even in cases involving rapidly changing boundary conditions and transcritical flow. Surprisingly, the time steps used in the present application for aggradation in straight channels are larger than others found in the literature (e.g. Park and Jain, Kassem and Chaudhry 1998, Cao and Egashira 2002). c) It has been demonstrated that the vertically-averaged and moment of momentum model (VAM) can be successfully applied to simulate the morphodynamics of meandering rivers without additional external sub-models for the effects of the secondary flow. This constitutes the first application of the 2D V A M model for river morphology. 108 5.3 F U T U R E R E S E A R C H 5.3.1 Convection-dominated applications This thesis is a first contribution to a potentially large field of future research. One important focus for research could be aimed to broaden the range of applications of the V A model. One of the main limitations of the model is that it cannot be applied for convection-dominated problems, such as bedform migration. Appendix A details the numerical discretization used by the model, based on the conventional Galerkin Finite Element Method (GFEM), which is known to be suitable for parabolic equations (diffusion), but not hyperbolic equations where convection is important (Dupont 1973, Katapodes 1984a, Hirsch 1987, Donea and Huerta 2003). G F E M are "centered" schemes, which give the same weight to information from both upstream and downstream. Centered schemes are appropriate for diffusive problem; but not when a wave or front is moving through the domain. Appendix F shows that, for convective problems, numerical schemes require "upwinding" properties; i.e. more weight should be given to upwind or upstream information when computing moving waves. A solution to this^ problem could be replacing the G F E M by a Petrov-Galerkin method, such as the Streamline Upwind Petrov Galerkin (SUPG) method currently used by the hydrodynamic model (Ghanem et al. 1995, Steffler and Blackburn 2002); this is the preferred option because it is consistent with the hydrodynamic model, but it is not the only one. In Appendix A it was demonstrated that a Finite Volume (FV) discretization based on unstructured triangular elements is also applicable for this model; therefore, another alternative could be to use one of the many available FV upwinding schemes. Once upwinding is introduced into the model, it could be applied to several interesting problems such as deposition upstream of check dams (Arminini and Larcher 2001, Busnelli et 109 al. 2001); migration of free alternate bars (Lanzoni 2000, Defina 2003); sediment waves (Cui et al. 2003a,b); among others. Some of those applications are discussed in Appendix F. 5.3.2 Sediment transport The model is currently limited to bedload-dominated channels with uniform sediment. For lowland sand bed rivers that may not be an important limitation (e.g. Waal River, Chapter 4); but that is not the case for mountain rivers covered with heterogeneous bed material (Papanicolau et al. 2005). Multiple size transport in a gravel bed could be incorporated using one or more of the existing approaches (e.g. Parker 1991, Kleinhans and Van Rijn 2002, Wilcock and Crowe 2003). Such a model would likely require some changes in the code as explained in Appendix C. Since the bedload and hydrodynamic models are semi-coupled (independent), the complexity of the sediment transport equations is not limited as in fully coupled models. Therefore, it is feasible to include complex sediment models (e.g. active layer, armouring). It is recommended to use the experimental data of Cui et al. (1996) to perform initial tests of a gravel transport model. To upgrade the model for suspended sediment transport would require solving the convection-diffusion equation for suspended sediment with a mechanism for the exchange between bedload and suspended material. Changes in the algorithms are also necessary to include the effect of partially wet elements in the sediment fluxes, as explained in Appendix C. The sediment continuity equation solved by the model assumes an equilibrium condition; i.e., the flow transports sediment at its full capacity, meaning that sediment transport depends only on local hydraulics. This assumption is valid to estimate the final equilibrium morphology of 110 an alluvial channel, but it might not be valid in cases such as highly unsteady flows (Wang 1999) or constrained upstream sediment supply (Phillips and Sutherland 1989), which requires a certain adaptation distance to achieve equilibrium conditions. In such cases non-equilibrium bedload transport should be considered. Recent examples of this type of application in 2D models are given by Nagata and Muramoto (2000), Wu (2004), Due et al. (2004) and Wu et al. (2005). Research can be conducted to determine i f non-equilibrium transport provides a real improvement over the current equilibrium approach. Potential test cases are degradation caused by sediment supply shut-off previously studied in chapter 2 (Figure 8) for straight channels and the experimental data of Yen and Lee (1995) of unsteady flow in bends. 5.3.3 Dam-break and dam removal Based on both safety and environmental concerns, dam decommissioning is becoming increasingly common in many parts of the world. The morphological impacts of the sediment released after dam removal is a topic of high current interest (Cantelli et al. 2004, Wong and Cantelli 2004, Cui and Wilcox, in press, Cui et al., in press). The challenge of this application is the potential presence of transcritical flow over the front of the deposited sediment (Cui and Wilcox in press), very similar to the pattern observed in the knickpoint migration experiments of Brush and Wolman (1960) and successfully simulated in Chapter 2. Most existing models intended for dam removal are currently ID (e.g. Cui et al., in press a, b), which may limit the analysis for dams where a 2D representation would be more realistic. The V A morphology model seems very promising for this application. Encouraged in part by the potential threat that global warming might have on the future safety of dams, there is renewed interest in studying dam-break and dyke breach phenomena (e.g. Fracarollo and Capart 1998, Ferreira et al. 2001, Zang and Wang 2001, Leal et al. 2001, Capart and Young 2002, Caleffi and Valiani 2002, Spinewine and Zech 2002, Fracarollo and Capart 2002, Leal et al. 2003, Spinewine et al. 2004). However, most numerical models are ID. In 111 Appendix G, some qualitative runs show that the model has the potential for applications in highly unsteady flows caused by dam break and dyke breach. This application will require the implementation of upwinding as discussed in section 5.3.1 for the presence of sharp moving fronts of sediment; probably also non-equilibrium transport (section 5.3.2) would be beneficial given the highly unsteady nature of the flow (Cao et al. 2004). 5.3.4 V A M model and non-linear secondary flow correction Although the V A M model with linear distributions of vertical velocity has proven successful for alluvial bend morphology, it is possible to replace it by a V A M model with parabolic velocity distributions. That change would probably require a different correction coefficient kb for the near-bed velocity direction (kb = 0.8 was used for the linear distribution). Whether the parabolic distribution is better than the linear one for morphological simulations is unclear. Ghamry and Steffler (2002) have shown that the computed flow field is not very sensitive to the shape of the velocity distribution; although they found that near the walls of curved channels the VAM-parabolic simulated the secondary flow much better than the VAM-linear. Research is needed to determine whether switching to the VAM-parabolic is beneficial at all. The V A M morphology model has proven successful for river morphology of meandering rivers. However, there are two main obstacles for this kind to model to become widely used: the 50% increase in computation time relative to the simpler V A model and the large number of existing V A models already in use. To overcome these problems, a VAM-type approach could be incorporated into existing 2D depth-averaged models to replace or enhance the linear secondary flow correction currently in use. Two preliminary tests, explained in Appendix D, have been successful in this objective: - The V A M model was run only once to compute the initial radii of curvature of the streamlines rc, which were later fed into the conventional secondary flow correction of 112 the V A morphology model (equation 3.5) to simulate alluvial bends with successful results. This procedure did not increase the computational time because the morphological simulation was performed with the V A model. The V A M replaced only the computations of radius of curvature and inertial adaptation done by the standard V A model (equations 3.8 and 3.9). The practical advantage of this approach is that the empirical calibration coefficient (3 for the inertial adaptation is not needed. - The V A M was run only once to compute the momentum velocity components (w/, v/), which were assumed as constants and fed into the V A model to compute the direction of the bed shear (using an equation similar to 4.21). The computed bed topography was in good agreement with measured data. Additionally, the information of (w/, v/) was also used to compute the bed shear stress redistribution caused by the secondary flow, which the V A model neglects. The results were encouraging. Another possibility (not tested yet) is to use the V A M to generate a non-linear secondary flow correction equation for the V A model (i.e. proportional to some power of the flow curvature h/rc) applicable for both weakly and strongly curved channels. The idea behind this research would be to investigate i f the V A M could be used as a pre-processor to compute the secondary flow features that can be later input into an existing depth-averaged morphology model. Optionally, the V A M could be invoked to update the curvature as needed. 113 REFERENCES Albers, C , Steffler, P .M. and Katapodis, C. (2002). Depth averaged and moment equation method for simulating vertical shear dispersion mixing in rivers. Proceedings of the International Conference on Fluvial Hydraulics River Flow 2002, Louvain-la-Neuve, Belgium, Vol . 1, pp. 91100. Arminini, A . and Larcher, M . (2001). Rational criterion for designing opening of slit-check dam. Journal of Hydraulic Engineering, 127(2): 91804. Ascanio, M . F. and Kennedy, J. F. (1983). Flow in alluvial-river curves. Journal of Fluid Mechanics, 133, 16. Belial, M . Spinewine, B., Savary, C. and Zech, Y . (2003). Morphological evolution of steep-sloped river beds in the presence of a hydraulic jump: experimental study. Proceedings of the XXXIAHR Congress, Thessaloniki, Greece, Theme C, Vol . 1, pp. 131140. Bhallamudi, S.M. and Chaudhry, M . H . (1991). Numerical modeling of aggradation and degradation in alluvial channels. Journal of Hydraulic Engineering, 117(9): 1145-1164. Blackburn, J. and Steffler, P .M. (2001). River 2D: Two-dimensional depth averaged model of river hydrodynamics and fish habitat. Tutorial, University of Alberta, Edmonton, Canada. Blanckaert, K. , Glasson, L. , Altinakar, M . , Jagers, H.R.A. and Sloff, C. J. (2003). A quasi-3D model for flow in sharp open-channel bends. Proceedings of the XXX IAHR Congress, Thessaloniki, Greece, Theme C, Vol . 2, pp. 317-324. Blanckaert, K . and Vriend, H. J. De (2003). Nonlinear modeling of mean flow redistribution in curved open channels. Water Resources Research 39(12): 6.1-6.14. Blanckaert, K. and Graf, W.H. (2004). Momentum transport in sharp open channel bends. Journal of Hydraulic Engineering, 130(3): 186-198. Blench, T. (1969). Mobile-bed fluviology: A regime theory treatment of canals and rivers for engineers and hydrologists. The University of Alberta Press, Edmonton. Blondeaux, P. and Seminara, G. (1985). A unified bar-bend theory of river meanders. Journal of Fluid Mechanics, 157: 449-470. Booker, D. J., Sear, D. A . and Payne, A . J. (2001). Modelling three-dimensional flow structures and patterns of boundary shear stress in a natural pool-riffle sequence. Earth Surface Processes and Landforms, 26(5): 551576. 114 Bradford, S.F. and Katapodes, N.D. (2000). The anti-dissipative, non-monotone behavior of Petrov-Galerkin upwinding. International Journal for Numerical Methods in Fluids, 33:581608. Brethour, J. M . (2001). Transient 3D model for lifting, transporting and depositing solid material. Proceedings of the 3rd International Environmental Hydraulics, Tempe, Arizona. Brush, L . M . and Wolman, M.G. (1960). Knickpoint behaviour in noncohesive material: A laboratory study. Bulletin of the Geological Society of America, 71: 59-73. Busnelli, M . M . , Stelling, G.S. and Larcher, M . (2001). Numerical morphological modeling of open-check dams. Journal of Hydraulic Engineering, 127(2): 105-114. Caleffi, V . and Valiani, A . (2002). A mathematical model for dam-break over movable bed. Proceedings of the International Conference on Fluvial Hydraulics River Flow 2002, Louvain-la-Neuve, Belgium, Vol . 2, pp. 10011012. Callander, R.A. (1969). Instability and river channels. Journal of Fluid Mechanics, 36: 465-480. Cantelli, A. , Paola C. and Parker, G. (2004). Experiments on upstream-migrating erosional narrowing and widening of an incisional channel caused by dam removal. Water Resources Research, 40(3).doi: 10.1029/2003WR002940, 2004. Capart, H. and Young, D.L. (1998). Formation of a jump by the dam-break wave over granular bed. Journal of Fluid Mechanics, 372:165-187. Capart, H . and Young, D.L. (2002). Two-layer shallow water computations of torrential geomorphic flows. Proceedings of the International Conference on Fluvial Hydraulics River Flow 2002, Louvain-la-Neuve, Belgium, Vol . 2, pp.10011012. Cao, Z., Day, R., and Egashira, S. (2002). Coupled and decoupled numerical modeling of flow and morphological evolution in alluvial rivers. Journal of Hydraulic Engineering, 128(3): 306-321. Cao, Z. and Carling, P.A. (2003). On the evolution of bed material waves in alluvial rivers. Earth Surface Processes and Landforms, 28:437-441. Cao, Z., Pender, G., Wallis, S. and Carling, P.A. (2004). Computational dam-break hydraulics over erodible sediment bed . Journal of Hydraulic Engineering, 130(7): 689-703. Cao, Z. and Carling, P.A. (2005). Further perspectives on the evolution of bed material waves in alluvial rivers. Earth Surface Processes and Landforms, 30:115-120. 115 Carey, F.C and Oden, J.T. (1986). Finite Elements: Fluid Mechanics. Volume VI. Prentice-Hall, New Jersey. Chang, H.H. (1988). Fluvial Processes in River Engineering. John Wiley and Sons. United States. Christison, K.J . , Steffler, P. M . and Katopodis, C. (1999). Two-Dimensional Hydrodynamic Modeling of a Braided River - A Case Study. Proceedings of the Third International Conference on Eco-hydraulics, Salt Lake City, US. Chung, T.J. (2002). Computational Fluid Dynamics, Cambridge University Press. Correia, L.R.P., Krishnappan, B.G. and Graf, W.H. (1992). Fully coupled unsteady mobile boundary flow model. Journal of Hydraulic Engineering, 118(3): 476-494. Cui, Y.T. , Parker, G., and Paola, C. (1996). Numerical simulation of aggradation and down-stream fining. Journal of Hydraulic Research, 34(2): 185-204. Cui, Y . , Parker, G., Lisle, T.E., Gott, J., Hansler-Ball, M.E. , Pizzuto, J.E., Allmendinger, N.E. and Reed, J .M. (2003a). Sediment pulses in mountain rivers: 1. Experiments. Water Resources Research, 39(9), 1239, doi:10.1029/2002WR001803. Cui, Y . , Parker, G., Pizzuto, J.E. and Lisle, T.E. (2003b). Sediment pulses in mountain rivers: 2. Comparison between experiments and numerical predictions. Water Resources Research, 39(9), 1240, doi:10.1029/2002WR001805. Cui, Y . , Parker, G., Lisle, T.E., Pizzuto, J.E. and Dodd, A . M . (2005). More on the evolution of bed material waves in alluvial rivers. Earth Surface Processes andLandforms, 30:107-114. Cui, Y . and Wilcox, A . (in press). Development and application of numerical modeling of sediment transport associated with dam removal. Sedimentation Engineering, ASCE Manual 54, Volume II. Cui Y . , Parker G., Braudrick C, Dietrich WE and Cluer B. (in press a). Dam Removal Express Assessment Models (DREAM). Part 1: Model development and validation. Journal of Hydraulic Research. Cui Y , Braudrick C, Dietrich WE, Cluer B, Parker G. (in press b). Dam Removal Express Assessment Models (DREAM), part 2: Sensitivity tests/sample runs. Journal of Hydraulic Research. Darby, S.E., Alabyan, A . M . and Wiel, M.J.V. de (2002). Numerical simulation of bank erosion and channel migration in meandering rivers. Water Resources Research, 38:2.2.21. 116 Defina, A . (2003). Numerical experiments on bar growth. Water Resources Research. 39(4), 1092, doi:10.1029/2002WR001455. Dendy, J. E. (1974). Two Methods of Galerkin-Type Achieving Optimum L 2 rates of Convergence for First-Order Hyperbolics. SIAM Journal of Numerical Analysis, 11: 637-653. DHI (1998). MIKE21C, user guide and scientific documentation, Danish Hydraulic Institute, Horsholm, Denmark. DHI (2000). MIKE11, A microcomputer based modeling system for rivers and channels, user and reference manual, Danish Hydraulic Institute, Horsholm, Denmark. Donea, J. and Huerta, A . (2003). Finite Element Methods for Flow Problems, John Wiley and Sons, London. Duan, J.G., Wang, S.S.Y., and Jia, Y.F . (2001). The applications of the enhanced CCHE2D model to study the alluvial channel migration processes. Journal of Hydraulic Research, 39(5): 469-480. Due, B . M . , Wenka, T. and Rodi, W. (2004) Numerical modeling of bed deformation in Laboratory Channels. Journal of Hydraulic Engineering, 130(9): 892604. Dupont, T. (1973). Galerkin Methods for First-Order Hyperbolics: An Example. SIAM Journal of Numerical Analysis, 10:890-899 Elgamal, M . H . (2002). Applications of the moment approach to low-flow regime bedforms. PhD thesis, University of Alberta, Edmonton, Canada. Federici, B. and Seminara, G. (2003). On the convective nature of bar instability. Journal of Fluid Mechanics, All: 125-145. Ferguson, R. I. (2003). The missing dimension: effects on lateral variations on 1-D calculations of fluvial sediment transport. Geomorphology, 56:14. Ferguson, R. I., Church, M . , and Weatherly, H . (2001). Fluvial aggradation in Vedder River: Testing a one-dimensional sedimentation model. Water Resources Research, 37(12): 3333347. Ferguson, R. I., Parsons, D.R., Lane, S.N. and Hardy, R.J. (2003). Flow in meander bends with recirculation in the inner bank. Water Resources Research, 39(11),1322. Ferreira, R .M.L. , Leal, J.G.A.B. and Cardoso, A . H . (2001). On the structure of the solution of the dam-break problem over a cohesionless erodible bed. Proceedings of the XXIX IAHR Congress, Beijing, China, Theme C, Vol . 1, pp. 319-326. 117 Fraccarollo, L. Capart, H . (2002). Riemann wave description of erosional dam-break flows. Journal of Fluid Mechanics, 461: 181228. Fredsoe, J. (1979). Meandering and braiding of rivers. Journal of Fluid Mechanics, 84: 609-624. Froelich, D.C. (2002). User's manual for FESWMS Flo2DH: Two dimensional depth-averaged flow and sediment transport model, Release 3. Federal Highway Administration. Publication No. FHWA-RD-03-053. Virginia, USA. Gessler, D., Hall B., Spasojeviv, M . , Holly, F., Pourtaheri, H . and Raphelt, N . (1999). Application of 3D mobile bed, hydrodynamic model. Journal of Hydraulic Engineering, 125(7):737-749. Galay, V.J . 1983. Causes of river bed degradation. Water Resources Research, 19(5): 1057-1090. Ghamry, H . K. (1999). Two dimensional vertically averaged and moment equations for shallow free surface flow. PhD Thesis, University of Alberta, Edmonton. Ghamry, H.K. and Steffler, P .M. (2002). Effect of applying different distribution shapes for velocities and pressure on simulation of curved open channels. Journal of Hydraulic Engineering, 128(11): 969-982. Ghamry, H.K. and Steffler, P .M. (2003). Two dimensional vertically averaged and moment equations for rapidly varied flows. Journal of Hydraulic Research, 40(5): 579-587. Ghamry, H.K. and Steffler, P .M. (2005). Two-dimensional depth averaged modeling of flow in curved open channels. Journal of Hydraulic Research, 43(1): 4225. Ghanem, A. , Steffler, P., Hicks, F. and Katopodis, C. (1995a). Two-dimensional modeling of flow in aquatic habitats. Water Resources Engineering Report No. 95-S1. University of Alberta, Canada. Ghanem, A. , Steffler, P., Hicks, F. and Katopodis, C. (1995b). Dry area treatment for two-dimensional finite element shallow flow modeling. Proceedings of the 12th Canadian Hydrotechnical Conference, Ottawa, Canada. Gi l l , M . A . (1980). Discussion on "Aggradation in streams due to overloading" by Soni, J.P and Garde, R.J., and Ranga Raju, K . G . Journal of Hydraulic Engineering, 106(HY11): 1955-1958. 118 Gil l , M . A . (1983a). Diffusion model for aggrading channels. Journal of Hydraulic Research, 21(5): 355-367. Gi l l , M . A . (1983b). Diffusion model for degrading channels. Journal of Hydraulic Research, 21(5): 369-378. Gi l l , M . A . (1987). Nonlinear solution for aggradation and degradation channels. Journal of Hydraulic Research, 25(5): 537-546. Graf, W. H . (1971). Hydraulics of sediment transport, McGraw-Hill, New York. Graf, W. H . (1998). Fluvial Hydraulics. Flow and Transport Processes in Channels of Simple Geometry, John Wiley & Sons, West Sussex. Guo, Q.C. and Jin, Y .C . (1999). Modeling sediment transport using depth averaged and moment equations. Journal of Hydraulic Engineering, 125(12): 1265269. Hager, W. (2003). Fargue, founder of experimental river engineering. Journal of Hydraulic Research. 41(3): 227-233. Hervouet, J .M. and Villaret, C. (2004) Integrated approaches for modeling sediment transport. Proceedings of the 2 n d International Conference on Fluvial Hydraulics River Flow 2004, Naples, Italy, Vol . 1, pp. 495-500. Hicks, F.E. and Steffler, P .M. (1992). Characteristic dissipative Galerkin scheme for open-channel flow. Journal of Hydraulic Engineering, 118(2): 337-352. Hicks, F.E., Steffler, P .M. and Yasmin, N . (1997). One-dimensional dam-break solution for variable widths channels. Journal of Hydraulic Engineering, 123(5): 462168. Hirsch, C. (1987). Numerical computation of internal and external flows. Volume 1: Fundamentals of numerical discretization, John Wiley and Sons, Brussels, Belgium. Ikeda, S., Parker, G. and Sawai, K . (1981). Bend theory of river meanders. Journal of Fluid Mechanics, 112: 361377. Jain, S.C. (1981). River bed aggradation due to overloading. Journal of Hydraulic Engineering, 107(HY1): 120-124. Jansen, P.Ph., Bendegom, L. van, Berg, J. van den, Vries, M . de and Zanen, A . (1979). Principles of River Engineering: the non-tidal alluvial river. Pitman, London. Jaramillo, W.F. and Jain, S.C. (1984). Aggradation and degradation of alluvial-channel beds. Journal of Hydraulic Engineering, 110(8): 1075085. Jin, Y . C . and Steffler, P .M. (1993). Predicting Flow in Curved Open Channels by Depth-Averaged Method. Journal of Hydraulic Engineering, 119(1): 109-124. 119 Johannesson, H . and Parker, G. (1989a). Linear theory of river meanders. River Meandering, A G U Water Resources Monograph, Washington, Vol . 1, pp. 18213. Johannesson, H . and Parker, G. (1989b). Velocity redistribution in meandering rivers. Journal of Hydraulic Engineering, 115(8): 1019-1039. Julien, P.Y. (2002). River Mechanics. Cambridge Univeristy Press, London. Kalkwijk, J.P.T. and Booij, R. (1986). Adaptation of secondary flow in nearly-horizontal flow. Journal of Hydraulic Research, 24(1): 19-37. Kalkwijk, J.P.T. and Vriend, H.J. De (1980). Computation of the flow in Shallow River Bends. Journal of Hydraulic Research, 18(4):327-342 Koch, F.G. and Flokstra, C. (1981). Bed level computations for curved alluvial channels. Proceedings of the XLX Congress of the IARH, New Delhi, India, Vol . 2, pp. 357-388. Kassem, A . A . and Chaudhry, M . H . (1998). Comparison of Coupled and Semicoupled Numerical Models for Alluvial Channels. Journal of Hydraulic Engineering, 124(8): 792502. Kassem, A . A . and Chaudhry, M . H. (2002). Numerical Modeling of Bed Evolution in Channel Bends. Journal of Hydraulic Engineering, 128(5): 507-514. Kassem, A . A . and Chaudhry, M . H. (2004). Closure of "Numerical Modeling of Bed Evolution in Channel Bends". Journal of Hydraulic Engineering, 130(1): 82-83. Katapodes, N.D. (1984a). A dissipative Galerkin scheme for open channel flow. Journal of Hydraulic Engineering, 110(4):450-466. Katapodes, N.D. (1984b). Two-dimensional surges and shocks in open channels. Journal of Hydraulic Engineering, 110(6):792512. Katapodes, N.D. (1984c). Fourier analysis of dissipative F E M channel flow model. Journal of Hydraulic Engineering, 110(7):927-944. Kleinhans, M.G. , and Rijn, L. C. van (2002). Stochastic prediction of sediment transport in sand-gravel bed rivers. Journal of Hydraulic Engineering, 128(4): 41825. Khan, A . A . and Steffler, P .M. (1996a). Vertically averaged and moment equations model for flow over curved beds. Journal of Hydraulic Engineering, 122(l):3-9. Khan, A . A . and Steffler, P .M. (1996b). Modeling overfalls using vertically averaged and moment equations. Journal of Hydraulic Engineering, 122(7):397-402. Khan, A . A . and Steffler, P .M. (1996c). Physically based hydraulic jump model for depth-averaged computations. Journal of Hydraulic Engineering, 122(10):540-548. 120 Kikkawa, H. , Ikeda, S. and Kitagawa, A . (1976). Flow and bed topography in curved open channels. Journal ofHydraulic Engineering, 102(HY9): 1327-1342 King, I. (2000). User Guide to SED2D-WES Version 4.5. US Army Waterways Experimental Station. Vicksburg, US. Kruger, S. and Olsen, N.R.B. (2001). Shock-wave computations in channel contractions. Proceedings of the XXUC IAHR Congress, Beijing, China, Theme D, Vol . 2, pp. 881-887. Kusakabe, S., Michue, M . , Hinokidani, O. and Fujita, M . (2003). A numerical simulation of flow pattern and bed variation on widening steep slope channels. Proceedings of the XXX IAHR Congress, Thessaloniki, Greece, Theme D, Vol . 1, pp. 335-342. Lane, S.N. (1998). Hydraulic modeling in hydrology and geomorphology: A review of high resolution approaches. Hydrological processes, 12:1131150. Lane, S.N., and Richards, K.S. (1998). High resolution two-dimensional spatial modeling of flow processes in a multithread channel. Hydrological processes, 12:1279-1298. Lane, S.N., Bradbrook, K.F. , Richards, K.S., Biron, P.A. and Roy, A . G . (1999). The application of computational fluid dynamics to natural river channels: three-dimensional versus two-dimensional approaches. Geomorphology, 29:20. Lane, S.N., Hardy, R.J., Elliot, L. and Ingham, D.B. (2002). High resolution numerical modeling of three-dimensional flows over complex river bed topography. Hydrological processes, 16:2262272. Lanzoni, S. (2000). Experiments on bar formation in a straight flume 1. Uniform sediment. Water Resources Research, 36(11): 3337-3349. Larson, R.E. and Hostetler, R.P. (1990). Calculus with analytical geometry. DC Heath and Company. Toronto. Leal, J.G.A.B., Ferreira, R.M.L. , Franco, A . B . and Cardoso, A . H . (2001). Dam-break waves over movable beds experimental study. Proceedings of the XXIX IAHR Congress, Beijing, China, Theme C, Vol . 1, pp. 23639. Leal, J.G.A.B., Ferreira, R.M.L. , Franco, A . B . and Cardoso, A . H . (2002). Dam-break waves on movable beds. Proceedings of the International Conference on Fluvial Hydraulics River Flow 2002, Vol . 2, pp. 981 -990. Leal, J.G.A.B., Ferreira, R.M.L. and Cardoso, A . H . (2003). Dam-break wave propagation over a cohesionless erodible bed. Proceedings of the XXX IAHR Congress, Thessaloniki, Greece, Theme C, Vol . 2, pp. 26268. 121 Leliavski, S. (1955). An introduction to fluvial hydraulics. Constable & Company Ltd., London. Levi, E. (1995). The science of water: the foundation of modern hydraulics. Translation of El agua segun la ciencia. A S C E Press, New York. L i , S. and Millar R.G. (2004). A computational distributed gravel budget for the lower Fraser River, British Columbia. Proceedings of the 2 n d International Conference on Fluvial Hydraulics River Flow 2004, Naples, Italy,Vol. 1, pp. 787-791. Lisle, T.E., Cui, Y . , Parker, G., Pizzuto, J.E. and Dodd, A . M . (2001). The dominance of dispersion in the evolution of bed material waves in gravel-bed rivers. Earth Surface Processes and Landforms, 26:1409-1420. Lomax, H. , Pulliam, T.H., and Zingg, D.W. (2001). Fundamentals of Computational Fluid Dynamics, Springer-Verlag. Lopez, J.L. and Falcon, M . A . (1999). Calculation of bed changes in mountain streams. Journal of Hydraulic Engineering, 125(3): 261270. Lyn, D.A. (1987). Unsteady Sediment-Transport Modeling. Journal of Hydraulic Engineering, 113(1): 15. Lu, J.Y., and Shen, H.W. (1986). Analysis and comparison of degradation models. Journal of Hydraulic Engineering, 112(4): 28299. Ma, L. , Ashworth, P.J., Best, J.L., Elliot, L., Ingham, D.B., Whitcombe, L.J. (2002). Computational fluid dynamics and the physical modeling of an upland river. Geomorphology, 44: 375-391. Martin-Vide, J.P. (2002). Ingenieria de rios (River Engineering). Ediciones UPC, Barcelona, Spain. Meselhe, E.A. and Sotiropoulos, F. (2000). Three-dimensional numerical model for open-channels with free-surface variations. Journal of Hydraulic Research, 38(2): 115-121 Mohammadian, A. , Azad, F.L and Tajrishi, M . (2003). Two dimensional numerical simulation of flow and geo-morphological processes near headlands by using unstructured grid. Proceedings of the XXX IAHR Congress, Thessaloniki, Greece, Theme D, Vol . 1, pp. 615-622. Molinas, A . (2000). User's manual for BRI-STARS. Federal Highway Administration, Publication No. FHWA-RD-99-190 122 Morvan, H. , Pender, G., Wright, N.G. , and Ervine, D.A. (2002). Three-Dimensional Hydrodynamics of Meandering Compound Channels. Journal of Hydraulic Engineering, 128(7): 672382. Mosselman, E. (2004). Discussion of "Numerical Modeling of Bed Evolution in Channel Bends" by Kassem and Chaudhry (2002). Journal of Hydraulic Engineering, 130(1): 82. Moulin, C. and Slama, E.B. (1998) The two-dimensional transport module SUBIEF. Applications to sediment transport and water quality processes. Hydrological Processes, 12: 11811195. Nagata, N . and Muramoto, Y . (2000). Numerical analysis of river channel process with bank erosion. Journal of Hydraulic Engineering, 126(4): 241252. Newbury, R. (1995). Rivers and the art of stream restoration. Natural and Anthropogenic Influences in Fluvial Geomorphology. Geophysical Monograph 89, Vol 1. A G U Press. N H C (2003). Bedford Channel hydraulic study. Northwest hydraulic consultants. North Vancouver, Canada. Odgaard, A.J . (1986a). Meander Flow Model I: Development. Journal of Hydraulic Engineering, 112(12): 1117-1136. Odgaard, A.J . (1986b). Meander Flow Model II: Applications. Journal of Hydraulic Engineering, 112(12): 1137-1150. Ollivier-Gooch, C. (2001). Mech 510: Computational Methods in Transport Phenomena I. Course notes, University of British Columbia. Vancouver, Canada. Olsen, N.R.B. (2000). CFD Algorithms for Hydraulic Engineering, Department of Hydraulic and Environmental Engineering. The Norwegian University of Science and Technology, Trondheim. Olsen, N.R.B. (2003). Three-dimensional CFD modeling of self-forming meandering channel. Journal of Hydraulic Engineering, 129(5): 366-372. Olsen, N . R. B. and Melaaen, M . C. (1993). Three-dimensional numerical modeling of scour around cylinders. Journal of Hydraulic Engineering, 119(9): 1048-1054. Olsen, N . R. B. and Kjellesvig, H . M . (1998). Three-dimensional numerical flow modeling for estimation of maximum local scour depth. Journal of Hydraulic Research, 36(4):570-579. Owen, D.R.J and Hinton (1980) A simple guide to Finite Elements, Pineridge Press Ltd. 123 Papanicolaou, A . N . , Bdour, A. , and Wicklein, E., (2005). One-dimensional hydrodynamic sediment transport model applicable to steep mountain rivers. Journal of Hydraulic Research, 42(4): 357-375. Park, I. and Jain, S.C. (1986). River-bed profiles with imposed sediment load. Journal of Hydraulic Engineering, 112(4): 267-280. Parker, G. (1976). On the cause and characteristic scales of meandering and braiding rivers. Journal of Fluid Mechanics, 76: 457-480. Parker, G. (2005) ID sediment transport morphodynamic, on-line e-book accessed on July 2005 at www.ce.umn.edu/~parker/morphodynamics_e-book.htm Parker, G., Klingeman, P. C , and McLean, D. (1982). Bedload and size distribution in paved gravel-bed streams. Journal of the Hydraulics Division, 108(HY4): 542271. Parker, G. and Johannesson, H . (1989). Observations on several recent theories of resonance and overdeepening in meandering channels. River Meandering, A G U Water Resources Monograph, Washington DC, 379-415. Pena, E., Fe, J. and Puertas, J. (2002). A 2D numerical model using finite volume method for sediment transport in rivers. Proceedings of the International Conference on Fluvial Hydraulics River Flow 2002, Louvain-la-Neuve, Belgium, Vol.1, pp. 691698. Peters, D.L and Prowse, T.D. (2001). Regulation effects on the lower Peace River, Canada. Hydrological Processes, 15: 3183194. Phillips, B.C. and Sutherland, A.J . (1989). Spatial lag effects in bed load sediment transport. Journal of Hydraulic Research, 27(1): 115-133. Pittaluga, M . B., Federici, B., Paola, C , Seminara, G., and Tubino, M . (2000). The morphodynamics of braiding rivers: Experimental and theoretical results on unit processes. Proceedings of Gravel Bed Rivers V, Wellington, New Zeland, Vol . 1, pp. 141181. Przedwojski, B., Blazejewski, R. and Pilarczyk, K.W. (1995). River Training Techniques: Fundamentals, Design and Applications. Balkema, Rotterdam. Rouse, H. and Ince, S. (1957). History of hydraulics. Iowa Institute of Hydraulic Research. University of Iowa. Rozovskii, J.L. (1957). Flow of water in bends of open channel. Academic of Science of the Ukrainian SSR, Kiev. 124 Seminara, G. and Tubino, M . (1989) Alternate bars and meandering: free, forced and mixed interactions. River Meandering, A G U Water Resources Monograph, Washington DC, 267-320. Schmautz, M . (2004). Gravel river widening by bank erosion in a straight stretch of river-investigations on a numerical model. Proceedings of the 2 n d International Conference on Fluvial Hydraulics River Flow 2004, Naples, Italy, Vol . 1, pp. 99-107. Schumm, S.A., Mosley, M.P. and Weaver, W.E. (1987). Experimental Fluvial Geomorphology. John Wiley & Sons. New York. Spasojevic, M . and Holly, F . M . (1990). 2-D bed evolution in natural watercourses-New simulation approach. Journal of Waterway, Port, Coastal, and Ocean Engineering, 116(4):425-443. Shimizu, Y . and Itakura, T. (1989). Calculation of bed variation in alluvial channels. Journal of Hydraulic Engineering, 115(3): 367-384. Shimizu, Y . , Yamaguchi, H. and Itakura, T. (1990). Three-dimensional computation of flow and bed deformation. Journal of Hydraulic Engineering, 116(9): 1090-1108. Sinha, S.K., Sotiropoulos, F. and Odgaard, A.J . (1998). Three-dimensional numerical model for flow through natural rivers. Journal of Hydraulic Engineering, 124(1): 1124. Simons, D.B. and Senturk, F. (1977). Sediment transport technology. Water Resources Publications, Fort Collins, Colorado. Spinewine, B., Capart, H. , Grelle, N . le, Frazao, S.S. and Zech, Y . (2002). Two-layer shallow water computations of torrential geomorphic flows. Proceedings of the International Conference on Fluvial Hydraulics River Flow 2002, Louvain-la-Neuve, Belgium, Vol . 2, pp. 10011012. Spinewine, B. and Zech, Y . (2002). Geomorphic dam-break floods: near-field and far-field perspectives. 1st IMPACT Workshop, HR Wallingford. Spinewine, B., Delobbe, A. , Elslander, L. and Zech, Y . (2004). Experimental investigation of the breach processes in sand dikes. Proceedings of the 2 n d International Conference on Fluvial Hydraulics River Flow 2004, Naples, Italy, Vol . 2, pp. 983-991. Soni, J.P., Garde, R.J. and Raju, K.G.R. (1980). Aggradation in streams due to overloading. Journal of the Hydraulics Division, 106(HY1): 117-131. 125 Steffler, P .M. , and Blackburn, J. (2002). River2D: Two-dimensional depth averaged model of river hydrodynamics and fish habitat. Introduction to depth averaged modeling and user's manual, University of Alberta, Edmonton, Canada. Steffler, P .M. and Jin, Y .C . (1993). Depth averaged and moment equations for moderately shallow free surface flow. Journal of Hydraulic Research, 31(1): 5-17. Struiksma, N . (1985). Prediction of 2-D bed topography in rivers. Journal of Hydraulic Engineering, 111(8): 1169-1182. Struiksma, N . (1989). Analysis of a 2-D bed topography model for rivers. River Meandering, A G U Water Resources Monograph, Washington, DC, 151180. Struiksma, N . , Olesen, K.W., Flokstra, C , Vriend, H.J. De (1985). Bed deformation in curved alluvial channels. Journal of Hydraulic Research, 23: 57-79. Suryanarayana, B. (1969). Mechanics of degradation and aggradation in a laboratory flume. PhD thesis, Colorado State University, Fort Collins, United States. Talmon, A . M . , Mierlo, M.C. Van and Struiksma, N . (1995). Laboratory measurements of the direction of sediment transport on transverse alluvial-bed slopes. Journal of Hydraulic Research. 33: 495-517 Termini, D. (2002). Numerical simulation of bed topography in large meadering channels. Proceedings of the International Conference on Fluvial Hydraulics River Flow 2002, Louvain-la-Neuve, Vol . 2, pp. 721729. Termini, D., Cristo, C. Di and Greco, M . (2001). Numerical simulation of experiments in a meandering flume. XXIXIAHR Congress Proceedings, Beijing, China, Theme E, Vol . 1, pp. 869-875. Tsujimoto, T. (1987). Non-uniform bed load transport and equilibrium bed profile. Proceedings of the XXII IAHR Congress, Lausanne, Switzerland, Fluvial Hydraulics, Vol . 1, pp.177-182. Tubino, M . and Seminara, G. (1990). Free-forced interactions in developing meanders and suppression of free bars. Journal of Fluid Mechanics, 214: 13159. U S A C E (1991). HEC-6, scour and deposition in rivers and reservoirs, User's manual. US Army Corps of Engineers, Davies, California. U S A C E (1993). River Hydraulics, Engineering and Design. US Army Corps of Engineers, Manual No. 1110-5416. Vasquez, J.A., (2004). Bridges: worse than a war. www.pepevasquez.com/bridges, web page last updated January 2004. 126 Vasquez, J.A., (2005). Numerical simulation of flow diversions. Proceedings of the 17 Canadian Hydrotechnical Conference, Edmonton, Canada. Vol . 1, pp. 845-853. Vasquez, J.A., Millar, R.G. and Steffler, P .M. (2005). Vertically-averaged and moment of momentum model for alluvial bend morphology. Proceedings of the 4th IAHR Symposium on River, Coastal and Estuarine Morphodynamics, Urbana, Illinois, Vol . 2, pp. 711-724. Vriend, H.J. De and Struiksma, N . (1983). Flow and bed deformation in river bends. Proceedings of Rivers '83, ASCE, Vol . 1, pp. 810-828. Whiting, P.J. and Dietrich, W.E. (1993). Experimental studies of bed topography and flow patterns in large-amplitude meanders: 2. Mechanisms. Water Resources Research, 29:3615-3622. Wilcock, P.R. and Crowe, J.C. (2003). Surface-based transport model for mixed-size sediment. Journal of hydraulic engineering, 129(2), 120-128 Willems, P., Vaes, G., Popa, D., Timbe, L. and Berlamont, J. (2002). Quasi 2D river flood modeling. Proceedings of the International Conference on Fluvial Hydraulics River Flow 2002, Louvain-la-Neuve, Vol . 2, pp. 721729. Yalin, M.S. (1992). River Mechanics. Pergamon Press, Oxford, New York. Yeh, K .C . , and Kennedy, J.F. (1993a). Moment Model of Nonuniform Channel-Bed Flow. I: Fixed Beds. Journal of Hydraulic Engineering, 119(7): 776-795. Yeh, K.C. , and Kennedy, J.F. (1993b). Moment Model of Nonuniform Channel-Bed Flow. II: Erodible Beds. Journal of Hydraulic Engineering, 119(7): 796-815. Yen, C.L. and Lee, W. T. (1995) Bed topography and sediment sorting in channel bend with unsteady flow. Journal of Hydraulic Engineering. 121(8): 591-599. Yusuf, F. (2001). Application of a two-dimensional hydrodynamic model to the Fraser River Gravel Reach. Master's thesis. University of British Columbia. Waddle, T. and Steffler, P .M. (2002). R2D_mesh: Mesh generation program River2D Two-dimensional depth averaged finite element, introduction to mesh generation and user's manual, University of Alberta, Edmonton, Canada. Wang, Z .Y. (1999) Experimental study on scour rate and river bed inertia. Journal of Hydraulic Research, 37(1), 17-37. Wang, S.S.Y. and Wu W. (2004) River sedimentation and morphology modeling-The state of the art and future developments. Proceedings of the 9th International Symposium on River Sedimentation, Yichang, China, Vol . 1, pp. 71-94. 127 Wong, M . , Cantelli, A. , Paola, C. and Parker, G. (2004). Erosional narrowing after dam removal: theory and numerical model. Proceedings of the A S C E World Water and Environmental Resources 2004 Congress, Salt Lake City, Vol . 1, pp. 10. Wu, W. (2004). Depth-averaged two-dimensional numerical modeling of unsteady flow and nonuniform sediment transport in open channels. Journal of Hydraulic Engineering. 130(10): 591-599. Wu, W., Rodi, W. and Wenka, T. (2000). 3D numerical modeling of flow and sediment transport in open channels. Journal of Hydraulic Engineering. 126(1): 185. Wu, W., Shields, F.D., Benett, S.J. and Wang, S.S.Y. (2005) A depth-averaged two-dimensional model for flow, sediment transport, and bed topography in curved channels with riparian vegetation. Water Resources Research, 41, W03015, doi:10.1029/ 2004WR003730. Zang, X . Z . and Wang, G.Q. (2001). Flow analysis and scour hole computation of dyke-breach. Proceedings of the XXLX IAHR Congress Proceedings, Beijing, China, Theme D, Vol . 2, pp. 270-277. 128 APPENDIX A DISCRETIZATION OF THE SEDIMENT CONTINUITY EQUATION SEDIMENT CONTINUITY EQUATION (EXNER'S EQUATION) dzb , dQsx , dasy ( l - A ) - ^ + - ^ + ^ - = 0 (A. l ) dt dx dy Equation (A. l ) is a partial differential equation (PDE) representing the conservation of sediment mass inside an infinitesimal control volume; qsx and qsy being the cartesian components of volumetric sediment transport rate per unit width qs (m3/s/m) in the x and y directions (horizontal plane), zb bed elevation, t time and X the porosity. It is assumed that qs can be computed using empirical equations from the previously computed flow field and that 129 porosity is constant (usually between 0.3 and 0.4); leaving zb as the only unknown. Using a depth-averaged hydrodynamic model to compute the flow field, equation (A. l ) can be applied to compute the evolution of bed levels through time. N U M E R I C A L D I S C R E T I Z A T I O N Equation (A. l ) is not solved numerically over the entire domain (the alluvial bed) but, instead, over a finite number of discrete points distributed over that domain. The purpose of such a discretization process is to transform the continuous PDE into a finite number of algebraic equations that can be readily solved by a computer. The unknowns in these equations are the bed elevation changes at discrete computational points, called nodes. The nodes are contained in small "elements" or "cells" which form a computational mesh. The mesh can be structured (regular) or unstructured (irregular). For unstructured meshes with triangular elements, there are two possibilities for the numerical discretization of the PDE, using either Finite Elements (FE) or Finite Volumes (FV). The main difference between the two methods is that the FV solution is represented by control volume averages, while the FE solution is represented by local shape functions (Ollivier-Gooch 2002). In other words, FV provide a stepwise approximation to the solution (average value within the element), while the FE provide a piecewise approximation to the same solution, as shown in Figure A - l for a one-dimensional (ID) case. X Figure A - l . Example of a bed elevation profile computed by both FE and FV. 130 For a 2D triangular mesh, the computational nodes are located at the triangle's vertices for FE and inside the element for cell-centered F V (Figure A-2). If the FV method is applied to solve equation A . 1 for the bed elevation change Azb, the result will be the average bed change Azeb inside the element. If the FE method is used instead, the result will be the bed changes Azbi, tszb2, Azb3 computed at the three vertices of the element. Azbi Azb2 FE element FV element Figure A -2 . Computational nodes (black circles) in both FE and FV formulations The main advantages of the FV formulation are that it can be easily derived from physical conservation laws. That makes the formulation not only simpler, because it is amenable to physical interpretation, but also conservative by ensuring that no mass is artificially gained or lost over the domain (some numerical discretizations cannot guarantee mass conservation). Being simpler means also being easier to code and debug, which is a practical advantage when programming the algorithm. However, for the present application that uses a FE hydrodynamic model to compute the flow field, the information of bed elevation at the vertices of the triangular elements is required. Therefore a mixed approach was applied. FV is used to compute the average bed changes inside all the elements and FE basis functions are used to project those values into the nodes. The word "node" hereafter is used exclusively to refer to FE computational nodes (triangle's vertices, Fig. A-2). 131 In the following sections the FV discretization is applied to compute the average bed changes in the elements. Later, the FE discretization is applied to derive a projection for that average element value into the nodes. FV DISCRETIZATION In general the FV method requires integrating the PDE in the control volume C Y (Lomax et al. 1999, Chung 2002) ( l - A ) f ^ d V + [ f ^ + ^ V = 0 dt The last equation arises from defining^ = qsxi +qsyj. Using the Gauss Theorem and assuming a fixed mesh (mesh is not moving), (1 -X) jf [yzbdV + jAqs-MA = 0 (A.3) In which n is the unit outward vector normal to the boundary and A the area of the element. The F V Method is not concerned with the details of the solution within a control volume; instead it is the average bed elevation inside the control volume zb =y^zbdV that is actually computed. Then, equation (A.3) can be written as follows 132 dt hdA (A.4) This equation states that the average bed elevation zb in the control volume changes at a rate determined by the net flux of sediment across the boundaries of the control volume. For a two-dimensional case, the equation can be further simplified to A (A.5) Where Qs, the flux integral, represents the net sediment flux across the three boundaries of a triangular element. Referring to Figure A-3, the flux integral for a given element e comprises three normal fluxes across the three sides 2, 7 and 11 (numbered in counter clockwise direction), 0 7 fevj) 3 2 (x2y2) Q Q.2 1 (xi,yi) ii Figure A-3. Normal sediment fluxes in a triangular element 133 (A.6) Then, equation (A.5) can be solved explicitly for every element as, 1 ^61-2+82-3+63-3-1 y\-x)Ae 1 ]Q: (A.7) Equation (A.7) has a very simple physical interpretation. In a given time step, the volume of solid material Vs leaving the element across its boundaries is Vs = QsAt. Taking into account porosity, the actual total volumetric change in the element is Vs I ( 1 - ^ ) - Therefore, the average bed level change in a particular element is found by dividing the total volume by the area of the element: Az^ = Q*At l(Ae'(1- A)). Since the flux of sediment leaving (entering) through a boundary is exactly the same entering (leaving) the next element, sediment continuity is enforced through the domain; i.e., this discretization method is conservative, no mass is lost or gained. These two features of the FV method, being conservative and amenable to physical interpretation, are some of the attractive characteristic of this discretization method. The flux integral can be computed as Qs = q„(2jh + qn(ijh + qn(ujhu where qn(i-j) is the normal outward sediment rate across side i-j and /,•_,• is the length of side i-j. The normal sediment rate can be expressed as a function of average sediment transport rate qs across the sides in both the x and j ' directions, being (x,,.y,) the coordinates of node / (Figure A-3). Qn(i-j) ~ asx(i-j) asy{i-j) ( V J \ (A.8) Since the hydrodynamic model is based on a FE formulation, velocities and flow depth (from where sediment transport rate qs can be computed) are known at the element's nodes. To 134 compute the average values along each element boundary, a linear variation of qs is assumed along the element side from one node to the other: = (qs^+qsiJ^)/2. Thus, the components of the sediment flux become Ql r ?«i 2 = 2 J -3 = V 2 J •1 = 2 J ' iyl-y\)-(^1-^3)-rqSyi + q ^ Qsy2 + Qsyl (x2 X,) (x3 x 2 ) ^ 3 + ^ 1 (Xj X 3 ) (A.9a) (A.9b) (A.9c) Finally, the net sediment flux becomes 0 = 2 ^ 2 - ^ 3 y3-y\ yi-yit qsx\ Qsx2 Qsxl I 2 1X3 X 3 X | X, Qsy\ °sy2 (A.10) The row vectors containing the node coordinates are constants which depend only on the geometry of the mesh. The column vectors containing the sediment transport rate at the nodes are updated every time step after the flow field is re-computed for the new bed topography. Replacing (A. 10) into (A.7) allows one to explicitly compute the average bed elevation change for every element. Although this FV formulation is straightforward, it cannot be applied alone for the present case. As mentioned before, the FE hydrodynamic model requires information of the nodal bed elevations, not the average value as given by equation (A.7). FE basis functions are used to project the value of the average element bed change into the three nodes, as explained below. 135 FINITE E L E M E N T S H A P E F U N C T I O N S For the FE formulation, once the values of the variables at the 3 nodes are known, the values inside the element are approximated by local shape functions Nj(x,y) (also known as trial, basis or interpolation function; Hirsch 1987). For example, zeb(x, y) = fjNf(x,y)zbl = N\z\ + Ne2ze2 + Ne3ze3 = {Nf}T{zebl} (A . l 1) Figure A -4 . Linear shape functions for 2D triangular elements The shape functions Nt have the property of being exactly Nt =1 at node i and decreasing linearly to zero towards the two remaining nodes. The triangles use linear shape functions as shown in Figure A - 4 . Thus, the derivates in equation (A.l) can be approximated as 136 dt = N? dz\ dz\ + N2 dt dt + Ne3 dt dt dq dx dN? dx dx 1sx2 +• dN dx dN? dx iltxi) dq dN? dy dy Isyl + • dNj dy 'Clsy2 dy dN? ~dy~ K , } (A. 12) Replacing the approximate values given by equation (A. 12) directly into equation (A. l ) will lead to a residual (discretization error). To minimize that residual the original PDE is multiplied by weighting (test) functions Wj such that the integral over the domain A is zero (Chung 2002). rs (\-X)\wi-^dA+ f dt )AWJ 8qsx , ^ + - sy dx dy \dA = 0 (A. 13) In the conventional Galerkin Finite Element Method (GFEM), the weighting function are set as identical to the linear shape function Wj = Nj. Multiplying equation (A. l ) by the shape function Nj (j - 1,2,3), integrating over the area and applying integration by parts to the second term leads to the so-called "weak" (or variational) FE formulation (1-/1) f N,.^dA = f lA 1 Fit JA dt dN, dx 8N, dy i sy dA-^{qs-n)Nj ds (A. 14) 137 Notice that i f the weighting functions inside the element were not linear but constants Wj= 1, equation (A .H) becomes similar to (A.3). If linear approximations are used in conjunction with triangular elements, the conventional G F E M is equivalent to a FV formulation (Lohner 2001). FV methods can be seen as G F E M by taking N, as polynomial and Wt = 1 inside the element and Wt = 0 outside (Hirsch 1987, Lohner 2001). The last term in equation (A.H) represents the flux of sediment across the boundaries of the elements. For internal (inter-element) boundaries, i f a continuity condition is imposed such that the outflux of sediment across an element boundary must equal the influx across the same boundary into the adjacent element, then all the inter-element boundary integrals will exactly cancel. This leaves only the global boundary. For no flow boundaries (e.g. sidewalls), the integral is zero as normal sediment flux is zero. Only in cases of imposed inflow of sediment, does the integral need to be evaluated along the inflow boundary. In all other cases, it can be simply ignored. A simple way to relate nodal with averaged bed changes is by placing (A. 12) into (A.H) with Wj = Nj and integrating using the following analytical equation (Owen and Hinton 1980) (A. 15) The left hand side of (A.2) becomes, inside an element, = ( 1 - A ) [ M , J (A. 16) Developing for j = 1,2,3 the following local square mass matrix [M0]e is found 138 [M}jf = \A{Ni}{NJ}TdA = £ N^dA \a N2N] dA £ JV 3 TV, <£4 J^TV, 7V 3^4 lN2N2dA [NJdA w v r = - 1 (A. 17) The integrals were computed from (A. 15) as f NNdA = JA ' J Ae/6 i = j Ael\2 i* j (A. 18) Similarly, the second integral of equation (A .H) can be expressed as a function of the element convection matrices [Cx]e and [Cy]e (Donea and Huerta, 2003) 'cw; A \~dx~, V J J><>' dx dA {qesJ-[Cx)e{gesxi} \~dx~ , \ J YsydA = 8Nej dx dA y J syi J (A. 19) The derivatives of the linear shape functions are constant inside an element and equal to 139 dNf 1 dx 2 A1 em I , dNe3 \ dN? 1 5y 2Ae (x 2 x 3 ) , 5x 2 v f dNe7 - 1 5y 2v4e ( x 3 X | ) ; dx 2Ae dm - 1 dy 2Ae ( x , - x 2 ) (A.20) The integral of the shape function, from (A. 15), is simply lNjdA = Ae/3, j = 1,2,3 (A.21) The convection matrices become y-L-y* v 3 ~ v i y\-yi y2 - ^ 3 ^3 - v i ^1 - ^ 2 ^ 2 - ^ 3 ^ 3 - ^ 1 y\-yi [cvr=-X^ X^ X^ X^ Xj X2 X~2 X-^ X^ Xj Xj X2 2^ "^1 *^ 2 (A.22) The second integral is then ^.«L4 = [ C J { ^ } + [ C J { ^ } = 140 y2-y3 y*-y\ y\-y2 y2-y^ y^-y\ y\-y2 y2-y^ y*-y\ y\-y2 qSx\ Qsxl l ? « 3 J ^2 "^3 *^2 *^3 *^3 "^ 2 X 3 X 3 X] asy2 asy3 (A.23) By comparing equation (A. 10) with equation (A.23) the relationship with the F V flux integral now becomes evident: + • sy •«l dx dy \N:dA = -3 3 Q. Q. Q. (A.24) By placing equations (A. 16), (A. 17) and (A.24) into (A 13), the discretization of the sediment continuity equation into an algebraic system at the element level is achieved: (\-A)Ae 6At I 1 2 1 2 'tell 1 3' 1 2 1 1 2 1 .2 1 2 1 teeb3 A4, Azeb2 2At (1 - A)Ae Q. Q. Q. (A.25) Finally, applying equation (A.7), the local system becomes "1 1 2 1 2 teebx Aze 1 2 1 1 2 teeb2 • = 2-Aze • (A.26) 1 _2 1 2 1 teeb3 Aze 141 Equation (A.26) [M]e{Azb}e = {RHS}e = {2Azb}erelates the local changes of bed elevation in the three nodes Azbi, Azb2, Azb3 of the triangular element with the average element bed elevation change Az^ inside the element, as computed by the FV approach. The assembled global system for all the nodes is then [M]{Azb} = {RHS} (A.27) S U M M A R Y OF C O M P U T A T I O N A L P R O C E D U R E The mass matrix [M] and convection matrices [Cx] and [Cy] are dependent only on the mesh geometry; they do not change in time. The local element matrices (A. 17) and (A.22) are computed and then assembled into global matrices using conventional matrix assembly rules of structural analysis (Owen and Hinton 1980). This process is performed only once before the morphodynamic simulation starts. During the morphodynamic simulation, the following procedure is performed every time step to compute the bed elevations at the nodes of the unstructured mesh: 1. From the values of velocity and water depth computed at the nodes by the hydrodynamic model and considering the effects of secondary flow and transverse bed slope, the components of sediment transport rate at every node, {qsx}and {q^jare found from a selected sediment transport equation. 2. Equations (A. 10) and (A.7) are applied explicitly to compute the average bed elevation changes {Szeh} in all the elements in the mesh. Those values are used to assemble the 142 {RHS} vector. The z'-th term of that vector is found by adding the average bed changes of all the elements attached to node /: rhst = 2 ^ ( A z | ) . 3. Equation (A.27) is solved to find the nodal values of the bed elevation. Those values are fed into the hydrodynamic model to update the nodal values of the flow field. L I M I T A T I O N S The numerical discretization using G F E M or conventional FV as explained here is applicable for the bed level deformation caused by changes in the upstream sediment supply (Chapter 2) or secondary flow in bends (Chapters 3 and 4). However, it is not suitable for convection-dominated bed changes (e.g. migrating bedforms), as explained later in Appendix F. 143 APPENDIX B RUNNING THE COMPUTER MODELS V A A N D V A M M O D E L S There are two main versions of the computer model, Vertically-averaged (VA) and Vertically-averaged and Moment (VAM), both written in C language using Microsoft Visual Studio.net™. The V A M version (Chapter 4) runs in a DOS window under a command prompt (Figure B- l ) . It was made by adding the file bedload.c with its respective functions into the CDG2D source code developed by Ghamry (1999) during his PhD research. This is a more research-oriented version, lacking its own graphical user interface (GUI). The visualization of results is made through the external program R2D_VAM10.exe. The V A version (Chapters 2 and 3) runs under Windows environment; the computational engine is fully integrated with a user-friendly GUI (Figure B-2). This version was made by adding the module morphology.c into the program River2D (Steffler and Blackburn 2002). This version is also known as River2D-MORphology or River2D-MOR. 144 nt\R2Dsed\D CDG2D Depth Averaged Hydrodynamic Simulation Development Iters i o n . Intended for research and evaluation puri Copyright: ft. Ghanem, H. Ghamry, P. S t e f f l e r and J . <Pepe) Uasquez Input command l e t t e r and integer code, for example: i 1 = input a mesh <.cdg f i l e ) , o 1. = output a mesh <.cdg f i l e ) , s 0 » run to steady state <fixed simulation time a c t u a l l y ) , s i - run to steady state <variable simulation time a c t u a l l y ) , b 1 = run bedload sediment transport , d -1 = run transient s imulat ion, q 0 • qui t a p p l i c a t i o n . Input Data F i l e Name filename.cdg Figure B - l . Command prompt of the V A M morphological model _ <3 X P i e "it V < . f ) p i * H ' t Bed Elevation MMM2.043 • l 769 I 495 " 1.221 |p0.946 . 0.672 ' II I 24 |l||-0.151 B-0.425 ^^ -0 699 Time 100d0h0m0 0s E X A M P L E O F R I V E R 2 D - M O R GUI Run Morphology Qin = 1350.0000 Qout = 1349 9131 x - 4182.158485, y - 3O2O.069S95 CAP NUM | Figure B-2. GUI of V A River2D-MORphology model 145 The source codes of both V A and V A M hydrodynamic models were made available by Dr. Peter Steffler from the University of Alberta, who also provided invaluable assistance and advice regarding the programming of the morphology module and its linkage to the hydrodynamic model. R U N N I N G T H E V A M M O D E L Getting started, file input il Figure B - l shows the initial window when the V A M program R2Dsed.exe is invoked. In the command prompt the user must type the command il to input the file with extension *.cdg. The C D G file contains the information of the mesh (elements, nodes) and the boundary conditions (water inflow discharge Q and downstream water level HBc)- After the input file is loaded a new window appears asking for Enter Actual Number of Variables = <3 or 5> Only two options are possible: 3 or 5 variables. The 3-variable option solves using the conventional depth-averaged or V A procedure; while 5 variables will solve using additional moment of momentum equations or V A M model. Using 5 variables will require about 50% more computational time than 3 variables, and more than twice the memory. After the number or variables is entered, the program gives the Actual number of unknowns, which is the number of nodes multiplied by the number of variables. After that, it goes back to the command prompt. 146 Flow simulation si Next, at the command prompt si should be entered to call the hydrodynamic model. The screen will show the message: Time = 0.0 Input first delta t, Max. Delta t, goal uchange, and model final time Time represents the current time of the simulations and it is set to zero before the model is run for the first time. The hydrodynamic model uses a variable time step for steady-state simulations. It starts with a small first delta t and as the solution converges, the time step increases up to a maximum specified by Max. Delta t. When the current time equals or exceeds the model final time the simulation stops. The variable uchange controls the change in velocity between one time step and the next, a typical value is 0.05. The rest of the variables are problem dependent. It is a good idea to start with a very small first delta t of 0.1 or 0.01 s; a typical value for Max. Delta t is 100 s and for model final time 1000. Figure B-3 shows the output screen at the end of a hydrodynamic simulation. The number of iterations, current time t, time step dt and net outflow are shown. As mentioned before, when / is larger than model final time, the simulation stops. The time step has increased steadily until it reached its maximum value of 100, that is a good indication that the model is converging. The net outflow is the difference between upstream inflow and downstream outflow, it should be a small value compare to the flow in the channel. Figure B-3 shows net outflow is small and constant, meaning steady state has been reached and the solution has converged. 147 28, t = 29.041265s dt = S.442858s net outflow = 0.0709564 29, t = 34.56?083s dt = 5.525818s net outflow = 0.0699175 30, t = 40.749497s dt = 6.182414s net outflow = 0.0478879 31 , t = 48.464968s dt = 7.715471s net outflow = 0.027312 32, t = 59.451126s dt - 18.986157s net outflow = i.0104997 33, t = 75.930362s dt = 16.479236s net outflow = 0.000729857 34, t = 100.649216s dt = 24.718854s net outflow = -0.0016524 35, t = 137.727497s dt = 37.078281s net outflow = -0.000796012 36, t = 193.344919s dt - 55.617422s net outflow = -0.000149396 37, t = 276.771052s dt = 83.426133s net outflow = ~3.18907e-005 38, t = 376.771052s dt = 100.i00i00s net outflow = 3.00099e -005 39, t = 476,771052® dt = 100.000000s net out f lo w = ~-3 .145 3 6e -005 48, t = 576.771052s dt = 100.000000S net outflow = ~3.16399e -005 41, t = 676.771B52s dt = 100.000000s net outflow = -3.16441e -005 42, t = 776.771052s dt = 100.000000S net outflow = -3.16424e -005 43, t = 876.771052s dt = 100.000000s net outflow = -3.l6422e -005 44, t - 976.771652s dt = 100.000000s net outflow = -3.16422e -005 45, \ „ t = 1076.771852s dt = 100.000000s net outflow = -3.16422e -005 T ime = 1B76.7?1§52 Input f i r s t d e l t a t . Max. D e l t a t;, g o a l uchange, and model f i n a l t i n e 10 100 1.05 Figure B-3. Flow simulation command 57. If the solution has not converged, si can be invoked again, as shown at the bottom of Figure B-3. It should be remembered that the new model final time must be larger than current Time for the hydrodynamic model to run. Output file command ol At any time after the input file has been loaded, an output file can be generated using the command ol which works similar to command il. The name of the output file should be carefully selected because this command will overwrite any file with similar name without a warning message. Bed load simulation command bl After a converged hydrodynamic solution is achieved the bed load sediment transport model can be invoked to perform computations of bed elevation changes. At this stage, all the 148 sediment parameters (e.g. grain size, coefficients of sediment transport equations, etc.) must be changed from within the program debugger; external users can input only the time step and the number of time steps; the total simulation time is the product of those two variables. During every time step the program will provide the Maximum relative bed change which is the highest value of the ratio Azb/H between bed elevation change Az* and a representative water depth H, usually the downstream water level HBC- The ratio Azb/H is expressed as a percentage. To avoid stability problems, the time step should be selected in such a way that the Maximum relative bed change remains small, usually smaller than 5%. At the beginning of the simulation, when bed changes are larger, the time step should be small. Since the bl command can be used many times consecutively, it may be advisable to run with small dt at the beginning and increase this value later in the simulation i f the Maximum relative bed change decreases. Quitting application command qO Once the simulation has finished, the program can be exited by using the command qO. R U N N I N G R I V E R 2 D - M O R P H O L O G Y M O D E L Since River2D-MOR is based on the hydrodynamic model River2D (www.River2D.ca), the interested reader is referred to Blackburn and Steffler (2001) and Steffler and Blackburn (2002) for details on how to run this model. Here only a brief description is presented of the additional features exclusively related to morphological modeling. 149 The input file to River2D-MOR must be a C D G file with a converged steady-state solution to provide the necessary initial conditions to the morphological model (nodal velocities and depths). The procedure to manage files and to obtain a steady-state solution is explained in detail in Blackburn and Steffler (2001) and it is not repeated here. The present version of River2D-MOR has been written by overwriting some of the functions available in River2D. Run Morphology replaces Run Transient (Fig. B-6); while the vector plots of sediment transport replace the vector plots of unit discharge (Fig. B-5). Therefore, River2D-MOR cannot be applied to simulate transient flow in a fixed bed. Currently work is underway in the GUI to include options for both fixed and movable beds in transient conditions. Also, work is underway to include an additional dialog box for sediment transport options; at the moment, all variables such as sediment size, sediment transport equation and others can only be changed from inside the source code (any change requires recompiling the program). To run the morphological model, from the main menu select Flow / Run Morphology. The dialog box Run Morphology, shown in Figure B-2, will pop-up asking the user for data to run the model. The user must enter the value of the Present time (initial time, usually 0), Final time, Time step increment At and Goal At. A l l units are in seconds. Time step increment At is a initial time step that will grow up to Goal At should be set to the same value. Once the data are input, click the button Run in the lower part of the dialog box. The options available in River2D to make animations (AVI videos), extract profiles and output C D G files at given time intervals, continue to be available for the morphodynamic simulations. In the Run Morphology dialog box press the Output Options button and follow the steps explained in Blackburn and Steffler (2001). The standard output option to extract water surface profiles has been replaced by bed profiles. 150 Visualization One of the most attractive features of River2D is that it allows one to see the evolution of any variable in real-time; while the computations are actually taking place (most commercial programs require an additional post-processor for visualization, like the V A M model does). This capability is very useful to see the morphodynamic evolution of the bed through time. Additional variables have been added for visualization: bed shear stress deviation angle (Fig. B-4), sediment transport rate (Fig. B-5) and bed changes (Fig. B-6). Q lfm_steady.cdg River?D Figure B-4. Deviation angle of the bed shear stress computed by the secondary flow correction model. 151 - X i i U S a *!D|Ai*|Eli 3| ml f| Bedload qs (cm2/s) VECTOR AND SCALAR FIELD OF SEDIMENT TRANSPORT 0 323 I — — " ™ 0 . 2 9 4 [jfiSR — 0 2 6 5 Qin = 0.1700 Qout = 0.1700 \ v > s \ \ v A \ \ \ \ \ \ \ \ \ \ \ v \ X \ \ \ A A \ \ V \ \ \ \ \ \ \ A A A \ \ \ X \ \ \ \ V \ A \ \ V \ Time 1 .OOOs Ptot ! c None a Sedirtent transport Scale (AnowLenglh/UnitVectat 1 Gnd Spacing Mnimum Depth 11-nt \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ . \ \ V-S'.v \ \ v \ \ U \ U \ \ \ X \ \ \ \ \ \ v \ \ \ \ M i > \ 1 I \ M i i i M 1 11 I I i i I I I I M I N i l l I I w,',V/\ CAP NUM Figure B-5. Vector and scalar field of sediment transport rate. « g m yon t Bed Change ^ _ 0 1 0 3 0.083 0.041 0.030 I.I 001 -0.063 H-0OB4 -0.105 Time 32m 0.000s BED CHANGES RELATIVE TO INITIAL TOPOGRAPHY Qin = 0.1700 Qout = 0.1700 X - -2.515394, y - 0.B62O9I Figure B-6. Colour shading plot of bed changes. The bed shear deviation angle (Fig. B-4) is computed by the secondary flow correction sub-model. This angle is linearly proportional to the flow curvature (h/rc), which is the ratio between water depth and radius of curvature of the streamlines. It represents the deviation of the bed shear stress from the direction of the depth-average flow velocity. Bedload sediment transport rate can be visualized as both vector and scalar fields (Fig. B-5). The direction of bedload is computed taking into account both the effects of the secondary flow and the transverse bed slope. The Cartesian components of bedload qsx and qsy can also be visualized as scalar fields on the screen. Figure B-6 shows a snapshot of the bed evolution during a morphodynamic simulation. The bed changes are computed by subtracting the current bed elevation from the initial bed elevation. Positive values represent deposition while negative values represent scour. At any time during the simulation, data from any of the variables shown on the screen can be extracted into a Comma Delimited File (CSV) for further analysis, for example using a spreadsheet like Microsoft Excel™. The data extracted can be the computed variables at all the nodes, or values interpolated into 2D grids or ID sections. Sections can be straight with uniformly distributed points, or along any arbitrary path i f the coordinates of points along that path are specified in an input CSV file. The latter option was used to extract the bed profiles along curved channels (e.g. Fig. 23, 24 and 25). The commands to perform these data extraction operations become available when selecting Display from the menu, as explained in detail by Steffler and Blackburn (2002). 153 APPENDIX C C O M P U T E R P R O G R A M : ALGORITHMS AND FUNCTIONS T H E T W O P R O G R A M S The computer codes for the morphological modules of both the V A model (River2D-MOR) and V A M model use the same algorithms and functions. The main difference between the two programs is how the main functions are called and the secondary flow computation. The main function in River2D-MOR is BedMorphoiogy () , which is contained in the file morphology.c and it is called by the transient flow function transientThreadProc () located in the file River2DDoc.cpp. The pseudo-code of transientThreadProc () is shown in Figure C - l . Before beginning the time loop, the global matrices [Cx], [Cy] and [M] (see Appendix A) are assembled and also the radius of curvature of the streamlines computed. Within the time loop, at every time t the morphology model is called to compute the bed changes, and immediately after the transient flow model is called to update the flow field, until the final time tflnai is achieved. 154 Assemble global matrices [Cx], [Cy] and [M] Compute curvature While t < tfmai Call BedMorphology() Call Transient() t <-1 +At End time loop Free memory from matrices End Figure C - l . Pseudo-code of transientThreadProc () In the V A M model, the main function is bedload () which is contained in the file bedload.c. bedload() has a structure similar to transientThreadProc!). The global matrices are assembled first and then a time loop is performed. In the time loop the bed levels are computed first and then the flow field is updated by calling the hydrodynamic model. The V A M does not require a secondary flow correction; therefore the radius of curvature of the flow is not computed. Since the functions and routines are the same for both V A and V A M models they will be explained in the context of the model River2D-MOR only. M A I N F U N C T I O N BedMorphology () Figure C-2 shows the pseudo-code of the main function BedMorphology {), which is called at every time step during the morphodynamic simulation (see Figure C-l ) . 155 The basic task of this function is to solve the system (see Appendix A , equation A.27). [M]{Azb} = {RHS} (C.l) Where {Azb} is the unknown bed elevation change at the nodes of the mesh and {RHS} is a vector computed from known sediment transport rates at the nodes. The mass matrix [M] remains constant in time, but the vector {RHS} is updated as the flow field changes. Begin Compute bed slopes Az^/Ax, Azb/Ay Compute secondary flow correction Compute {qsx}n and {q^Y Compute {RHS} n+'/2 Solve [M]{Az 6} n + ' / 2 ={RHS} n+Yl Compute {z 6} n + ' / 2 = {zb}n + {Azb}n+'/2 Compute {^} n + ' / 2 and {qsy}n+'/2 Compute {RHS}n+1 Solve [M]{Azb}n+l ={RHS}n+l Compute {zb}"+l = {zb}n + {AzA}" + 1 Check maximum bed slopes End Figure C-2. Pseudo-code of BedMorphology () 156 First, the bed slopes along the two Cartesian direction Az^/Ax, Az^/Ay and the secondary flow correction are computed. This information is required for the direction a of the bedload sediment transport, from which the two components qsx and qsy can be computed. Next, the vector {RHS} is computed using a FV formulation (equations A . 10 and A.7). The program advances in time using a second order Runge-Kutta scheme. The bed elevation is first computed at mid-time n+Vi (predictor). The {RHS} vector evaluated at n+Vi is used later to compute bed elevation at time n+\ (corrector). If selected by the user, an algorithm to check i f the bed slopes are not larger than the angle of repose can be called. This algorithm can check and reduce bed slopes in places where they are too steep. By default, this option is turned off. Computing bed slopes The direct evaluation of the bed slopes at the nodes is impossible because the shape functions exhibit a discontinuity in slope there (see Fig. A - l and Fig. A-4). They are obtained using the so-called "recovery of derivatives" procedure (Lohner 2001), by which nodal bed slopes are approximated by solving the system (C.2a) (C.2b) 157 Computing sediment transport rate vectors {qsx} and {qsy} Sediment transport at the nodes is computed by the function s e d i m e n t t r a n s p o r t ( ) . Input to this function are the nodal values of bed slopes and the bed shear angles 5 computed by the secondary flow correction algorithm. Also, other variables such as sediment median size dso and sediment shape factor fs are also read by this function. The pseudo-code is shown in Figure Compute mean velocity U = -yju2 + Compute dimensionless shear stress r* tana = (sin S-(Azb/Ay)/(f/))/(cos5-(Azb/Ay) /(//)) Compute qs qsx = ^cosa ; qsx = ^s ina If node is on no-flow boundary make qs parallel to mean flow qsx = qs(u/U ) ; qsx = qs (v/U ) Else qsx = qsy=0 C-3. For every node: Read ks, h, u, v, Az^/Ax, AzblAy, 5 If node is wet (h > hmi„): Next node Figure C-3. Pseudo-code of function sediment_transport(). 158 Function sediment_transport( ) reads the updated hydrodynamic flow values of water depth h and velocity components (u, v) at every node. Sediment transport rate is computed only at wet nodes, where water depth h is positive and greater than a minimum user-defined value hm\n. In nodes where h < hmm sediment transport is set to zero. The sediment transport angle, a, corrected by both secondary flow and transverse bed slope, is applied only to interior nodes. For nodes located along fixed sidewalls, it is assumed that sediment transport follows the direction of the mean flow (tangential to the walls). This is done to avoid having sediment leaving or entering through the sidewalls. The dimensionless shear stress r* is computed from the mean velocity U =yjul + v o u s m 8 Shield's equation (C.3). Gs = 2.65 is the specific gravity of the sediment and d5o is the median size of the bed sediment. T * = — (C.3) Ch(Gs-l)d50 Four options are currently available to compute sediment transport rate. Using the equations of Engelund-Hansen, Meyer-Peter Muller or Van Rijn, plus a simple power equation of the type qx =aUb, with coefficients a and b specified by the user. Computing vector {RHS} Function compute_RHS () computes {RHS} using information of sediment transport components {qsx} and {qsy}, as well as information about porosity, time interval and variables i e q u i l and q s i n . This function is called twice during the Runge-Kutta time stepping (Figure 159 C-2). In the first sweep (predictor) the time interval is set to At/2 and in the second sweep (corrector) the full value At is used. For every element: Read three nodes nh i = 1, 2, 3 Read coordinates xit yh i= 1,2,3 Compute Qx_2, Q2_2, g3_, using equation (A.9) If i e q u i l — 1: If side i-j on the inflow cross section: h-j=[(xj-xd2 + (yj-yd2f5 Qt-j = qsM'li-j. End if End i f Qs = 0,-2 +02-3 +03-1 rhsm = rhsni - 1/3(0,) At /( l -X), i = 1, 2, 3 next element Figure C-4. Pseudo-code of function compute_RHS(). i e q u i l is a decision variable for the upstream boundary condition of sediment supply. The default value is i e q u i l = 0, meaning that sediment supply is in equilibrium condition; i.e., the flow transports sediment at its full capacity. If i e q u i l = 1, the inflow of sediment is forced to be equal to qsw; which was the option used for aggradation/degradation test cases caused by overload/deficit of upstream sediment supply shown in Chapter 2. For example, for a case of 160 degradation below a dam with zero or negligible upstream sediment supply the variables i e q u i l =1 and qs\N = 0 should be used. The {RHS} vector is computed using the Finite Volume formulation as explained in Appendix A (equations A.6 through A. 10). The advantage of FV is that setting boundary conditions is straightforward. For example, i f the normal inflow of sediment upstream is forced to a value qSTM ( i e q u i l = 1), the model checks if the side i-j of the element is located on the inflow cross section. If true, the normal sediment flux through that side becomes <_Vj = qs\Uli-j, where /,_, is the length of side i-j. Computing nodal bed changes Once the vector {RHS} is known, nodal bed changes can be computed by solving the linear system of equations (C.l). This system has to be solved two times at every time step for the two sweeps of the Runge-Kutta scheme. Additionally, a similar system of equations (C.2) has to be solved also at every time step to compute the bed slopes. This large system of NnodesxNnodes, where Nnodes is the number of nodes, has to be solved implicitly 4 times every time step, representing the most computationally demanding part of the model. To alleviate the computational expense of solving this complete linear system, an explicit approximate approach has been implemented by lumping matrix [M] into the diagonal mass matrix [ML] (Donea and Huerta 2003, Lohner 2001, Defina 2003). The diagonal value my in [ML] is the sum of all the numerical values of the row / in original mass matrix [M]: (Nnodes m Lij ~ 2> y * = J (C.3) 7=1 0 161 Using the lumped matrix is less accurate but more stable (Donea and Huerta 2003). It is easier and faster to compute since the approach becomes fully explicit and does not require solving the full system (Carey and Oden 1986). Thus, the bed elevation change in node / in a given time step is found by the function compute_dz_node () as (Defina 2003) A , S i = ^ (C.4) mUj A similar explicit approach is also used to compute the bed slopes (equation C.2). Checking maximum bed slopes Function check_max_siope () is intended to ensure that bed slopes remain equal to or lower than a maximum bed slope defined by the angle of repose <j) of the bed material: Azb I Ax < tan<fi and Azb I Ay < tan <f>, following Olsen's (2000) sand slide algorithm. First, all the elements are checked to verify i f any of them has a side A-B steeper than tan^ 5. If that were the case, the function sand_s i ide () is called. As ZA and zg are the bed elevations of nodes A and B respectively, bed elevation corrections AZA and AZB are added to those nodes to reduce the bed slope to tan^ (Fig. C-5) (zB+AzB)-(zA+AzA)=ta^ ( C 5 ) I A-B 162 Physically, this can be understood as a "sand slide" of a portion of the bed that is too steep. The higher point B decreases its elevation, while the elevation of the lower point A is raised as shown in Figure C-5. Figure C-5. Reduction of bed slope along side A-B of an element. Figure C-6. Two group of elements that share node A and node B. 163 In order to ensure mass conservation, the volume of sediment before and after AZA and AZB are applied must remain the same. YAA and YAB being the areas of the two groups (or patches) of elements that share nodes A and node B (Fig. C-6), mass conservation can be expressed as: YAA ZA + YABZB = YAA (ZA+AZA) + YAB (ZB+AZB); leading to Once AzA is known from (C.6), AzB is computed from (C.5). Finally, bed elevations at nodes A and B are updated. This algorithm has not been thoroughly tested and by default it is not applied. force_BC(dz_node) It was observed that sometimes the elevation of one node in the inflow cross section overshots, especially in simulations involving aggradation. That leads to instabilities that make the model diverge. To avoid that problem the function forceBco forces the inflow cross section to remain horizontal. This function simply computes the average bed elevation change of the inflow section and applies that value to all the nodes in that section, forcejc () is called after the nodal elevations are computed. This function is valid only for a flat inflow cross-section with evenly spaced nodes fully underwater. For applications in a river, where the inflow cross section may be compound with dry nodes, this function is likely to require testing and changes. (ZB-ZA) + 1A-B W (C.6) 164 POTENTIAL ENHANCEMENTS TO THE MODEL Multiple size sediment transport Implementing more complex sediment transport equations in the future will require some changes in the function sediment_transport( ). For multiple size fractions, it is likely that vector {Qs}Nnodesx\ would need to be expanded into a matrix [qs]NnodesxNfraci where Nfract is the number of sediment size fractions. Also, since r* (equation C.3) is dependent on sediment size, angle a would probably need to be computed for every size fraction too. Additionally, function computeRHS o would likely require changes in order to compute the normal sediment fluxes for every size fraction, plus an algorithm for the exchange of sediment between surface and sub-surface layers. Partially dry elements Another potential enhancement in function computeRHS () could be the treatment of partially dry elements. The current algorithm assumes the entire element is wet. Figure C-7 illustrates an example of an element with 2 wet nodes and 1 dry node (node 3). Sediment transport only happens in the wet portion below the waterline (shaded area). However, the present algorithm will consider that sediment is eroded or deposited in the whole element (over-estimation of effective area), being possible for a dry element to change its bed elevation. Also, the model will linearly interpolate sediment transport assuming qS3 = 0 (dashed lines), when in reality sediment transport is zero at the waterline (shaded area), thus over-predicting the real sediment flux across the element. Further investigation is needed on this subject to decide a suitable approach to deal with partially dry elements. 165 qsx2 Figure C-l. Sediment fluxes (x-direction) and wet area in a partially dry element. Solving for the mass matrix Donea and Huerta (2003) report that for purely convective transport, the lumped mass matrix [ML] is more stable (larger time steps) but it exhibits some oscillations and a phase lag not observed when the full matrix [M] is used. For traveling waves, the accuracy of the lumped matrix approach is seriously compromised. Carey and Oden (1986) argue that, for the unsteady Navier-Stokes equation, using the lumped matrix leads to poor accuracy (first order) and a reduction of time step for stability reasons. For solving convective problems (e.g. migrating bedforms) it might be worthwhile to solve for the full mass matrix. Lohner (2001) proposes a simple procedure: [ML]({Azb}k+]-{Azb}k) = {RHS}k-[M]{Azb}k , *=0,. . . , nUer, {Az»}o= 0 (C.7) Lohner (2001) states that usually no more than 3 iterations (niter < 3) are employed. Notice that the current approach used by the model (equation C.4) computes {Azb} i only. 166 APPENDIX D C U R V A T U R E AND INERTIAL ADAPTATION C U R V A T U R E Knowledge of the local radius of curvature of streamlines rc is necessary to apply the secondary flow correction of depth-averaged models (see Chapter 3). For circular bends, rc could be assumed to be equal to the constant radius or curvature of the bend. However, for channels of varying curvature (e.g. natural meandering channels), rc can be computed from the definition of curvature found in standard calculus textbooks (e.g. Larson and Hostetler 1990) \u\ rc=Tzr±r1 (D.l) Where u =ui + vj = (dxf dt)i +(dy/ dt) j is the velocity vector in the x-y coordinate system, t is time and 167 is the centripetal acceleration as shown in Figure D - l . The components of the acceleration vector are du _ du ^ du dx du dy dt dt dx dt dy dt (D.3a) dv _ dv dv dx dv dy dt dt dx dt dy dt (D.3b) Assuming steady state (duldt = 0, dvldt = 0), equation (D.3) can be approximated as du du du — = u \- v — dt dx dy (D.4a) dv dv dv — = u — + V — dt dx dy (D.4b) The vector product of velocity and acceleration in the x-y plane is \uxaj = i j k u v 0 du I dt dv/dt 0 (D.5) 168 By replacing equation (D.4) into (D.5) and developing, reduce equation (D.l) into rc =" (u2+v2)3'2 du 2 du 2 dv dv uv v hu — + uv-dx dy dx dy (D.6) — Streamline u Figure D - l . Radius of curvature, velocity and centripetal acceleration in a point of a streamline. Equation (D.6) has a structure similar to equation (51) in Kalkwijk and Booij (1986). However equation (51) in Kalkwijk and Booij (1986) should not be used because it contains some typographical mistakes in the signs of several terms. INERTIAL ADAPTATION The radius of curvature computed by equation (D.6) assumes that the flow field adapts instantaneously to changes in curvature. However, in regions where the curvature of streamline varies, the secondary flow will adapt gradually. The flow field reaches equilibrium with the curvature only in infinitely long bends (Blanckaert & Vriend 2003). If equation (D.6) is 169 applied, the computed bed profiles turn to be slightly shifted upstream, as shown in Figure D-2. The radius of curvature given by (D.6) needs to be adjusted for the inertial effects using a semi-heuristic inertial adaptation equation (Rozovskii 1957, Struiksma et al. 1985, Darby 2002) Vw2 +v2 ( d ( 0 d ( 0 \ tl — + v— dx dy ) 1 lux a. • —13 \u\ (D.7) where Kf=P Ch (D.8) is the secondary flow adaptation length, which is a function of water depth h and roughness C. P is calibration parameter of order 1 and g the gravitational acceleration. Equation (D.7) is applied as an iterative correction to the radius computed by equation (D.6) 1 1 'c,k+\ 'c,k dx \rc,kJ + v-dy \.rc,kJ k 1,2,... yijter (D.9) The initial radius rc \ is computed according to (D.6). It was found that only one correction inner =2) was necessary. The inertial adaptation affects mainly areas of varying curvature. Its effect is to shift the computed profiles downstream as shown in Figure D-3. The amount of shifting depends on the empirical parameter /3. 170 Figure D-2. Resulting bed profiles for both (a) D H L and (b) L F M flumes, computed without inertial adaptation (j3 = 0). 171 0.30 Distance / Width Figure D-3. Effect of inertial adaptation algorithm on computed curvature of streamlines for L F M flume. V A M C U R V A T U R E The disadvantage of the secondary flow correction explained above is that inertial adaptation with its empirical parameter /? is needed. The experiments simulated in Chapter 3 for a vertically-averaged (VA) model showed a variability of p between 0.4 and 2.0. In contrast, the vertically-averaged and moment of momentum (VAM) model introduced in Chapter 4 does not require explicit information regarding the radius of curvature or its inertial adaptation. However, the V A M requires about 50% more computational time than a V A model. It would be desirable to find an approach to compute the secondary flow correction for a V A model with some of the advantages of the V A M model, but without an increase in computational time. One possibility is to use the V A M model to compute the curvature for the V A model. Another is to use the moment velocity component (wi,vi). Both approaches are explained below. 172 Radius of curvature from V A M model The direction 8 of the bed shear stress is computed by the V A model using 'VA = arctan — -arctan A — I rc) (D.10) The parameter A is a function of von Karman's constant K and roughness (Struiksma et al. 1985) 1- (D. l l ) For the V A M model, 8 is computed as SVAM = arctan (D.12) where (u0> v0) are the V A M depth-averaged velocity components; {u\,v\) are the moment velocity components in excess of the mean at the surface; and coefficient ^=0.8. Making 6 = arctan(w,v) and assuming 8VA =8VAM allows one to find the radius of curvature as ~ = 4i; t a n(0-<W) r„ Ah (D.13) Equation (D.13) is computed only once, before the morphological simulation starts, to compute the value rc, which is used in equation (D.10) by the conventional V A model. There is practically no increase in computational time since the V A model is being used for the morphological simulation. The additional time needed to compute the V A M flow field is in part compensated because equation (D.13) replaces equations (D.6) and (D.7). Besides, the 173 computational time needed to compute the initial steady flow fields is negligible compared to the total time needed for the morphological simulation. A n important advantage of using equation (D.13) is that the inertial adaptation equation, and hence calibration parameter /3, are not needed. Figure D-4 shows the results of applying this approach for the strongly curved L F M flume. It is evident that the model reproduces very well the beginning and end of the bed profiles, even for strong curvatures. This is not possible with V A models i f inertial adaptation is ignored (Fig. D-2). Figure D-5 shows similar results for the Waal River. 174 c o '•S3 re I a) 0) WAAL RIVER: Right margin Computed *-«^r; j sounding '' 4 5 6 7 Distance (km) 3.5 2.5 + 1.5 0.5 -0.5 -1.5 + -2.5 -3.5 8 9 10 C o '•IS (0 > <D •o a> CO WAAL RIVER: Left margin ' \ Computed / \ sounding -3 4 5 6 Distance (km) 3.5 2.5 1.5 -0.5 -1.5 -2.5 -3.5 9 10 Figure D-5. Bed profiles for Waal River computed by V A model using curvature computed by V A M model (inertial adaptation not necessary). 175 V A M velocity components (ui,vi) Another possibility tested for the secondary flow correction of the V A model was the use of the («i,vi) velocities computed by V A M model, ( 8VA = arctanl v - khv u - kbux (D.14) The results for the L F M are shown in Figure D-6. 1.8 2 Distance/Width 4 6 • BEND 10 Measured Figure D-6. Bed profiles for L F M flume computed by V A model using (u\,v\) velocities computed by V A M model (inertial adaptation not necessary). 176 Similar to the previous approach, the V A M flow field was computed only once at the beginning of the simulation to provide the information of (u],v\). Equation (D.14) is similar to equation (D.12), but with two important differences: the depth-averaged velocity field (u,v) in equation (D.14) is computed by the V A model, and the ( « i , V i ) field is not updated, it remains "frozen" during the morphological simulation. By comparing Figures D-4 and D-6 it can be seen that both approaches provide reasonable results. Using the radius of curvature computed from the V A M model (Fig. D-4) provides an excellent prediction of the beginning of the bed profile, but a rather constant transverse slope in the downstream half of the bend. When the («i,vi) field is used instead, the profiles are slightly shifted downstream, but the oscillation of the transverse bed slope, increasing again around the bend exit, is more clearly defined (Fig. D-6). Summary It has been shown that the flow field computed by the V A M model can be used to provide useful information for the secondary flow correction of the traditional V A model. The advantage of using such information is that the semi-heuristic inertial adaptation equation and its calibration parameter /? become unnecessary. In addition, i f equation (D.14) is used, the need for computing the radius of curvature is eliminated altogether. 177 APPENDIX E VARYING WIDTH APPLICATIONS P R O B L E M S IN C O R N E R S As demonstrated by the morphological applications shown in Chapters 2, 3 and 4, the proposed FE river morphology model is suitable to simulate bed level changes in straight and curved channels not dominated by convection. However, all those applications are for channels of constant width. Although the model can take into account the effects of varying width on bed morphology, there are still some stability problems in places where velocity direction changes rapidly, such as at corners along sidewalls. Because of the "slip" (zero stress) condition used by the hydrodynamic model along no-flow boundaries, a tangential velocity is computed at the boundary (Ghanem et al. 1995). In cases of sudden changes in the direction of the velocity along the boundaries (referred to here as "corners"), the nodal value of velocity in the corner could overshoot (Fig. E- l ) . Typically, in 178 convex comers the local velocity can be considerably larger than the average value in the channel, while in concave corners, the opposite effect is observed. These high/low nodal velocities at corners do not affect the stability of the hydrodynamic model but, unfortunately, they do have a strong effect on bed level computations. Bed level changes depend on the spatial gradient of sediment transport which, in turn, depends on a high power of velocity. Coiners can produce extremely high gradients in sediment transport rate which translate into acute local changes in bed elevation, as shown in Figure E-2. These nodal bed changes at the corners can be orders of magnitude larger than in the rest of the channel. Since these acute local bed perturbations tend grow and propagate, they can turn the model unstable. Some attempts to solve the problem, such as forcing bedload to zero along the sidewalls were unsuccessful. It is unknown if changing the boundary condition to no-slip in the hydrodynamic model could be more realistic. Ghanem et al. (1995) stated that a no-slip condition will negatively affect the solution close to the boundary, unless enough small elements are used. Molls and Chaudhry (1995) demonstrated that, for flow along curved channels, the slip condition provides results which are significantly more realistic than a no-slip condition. The influence of a no-slip condition on the bedload model is a topic for further research. A P P L I C A T I O N S Although some stability problems remain, three applications of channels with varying width are shown in this Appendix to demonstrate that the model is able to capture the influence of channel narrowing and widening on the bed levels. Only the first application, a fictitious channel with gradual transitions, was free from stability problems. The other two applications were stabilized using smoothing algorithms that are explained later. 179 Velocity UPSTREAM PORTION OF URBAN REACH OF PIURA RIVER Qin = 1000.000 Clout = 1000.000 Figure E - l . Example of high/low nodal velocities at convex /concave corners. Figure E-2. Typical pattern of local scour and deposition around convex corners. 180 C O N T R A C T I O N A N D E X P A N S I O N IN A S T R A I G H T F L U M E This is a qualitative application using a gradual change in the width of a straight flume by means of 3 m long transitions (Figure E-3) to avoid the problems with corners mentioned before. The channel is 31 m long and 10 m downstream from the inflow section (x = 10 m) a gradual contraction reduces the width from 2.5 to 2 m. The narrower section of the channel extends 5 m, to return again to its original width of 2.5 through a gradual expansion. The inflow is Q = 0.25 m3/s and downstream depth 0.2 m. The channel is initially flat with a constant bed elevation at zb = 10 m. The sediment transport rate qs is assumed as a power of mean velocity U: qs = 0.00145L/5 (identical to that of the Soni et al. (1980) experiment, see Chapter 2). Contraction Expansion 3 m 5 m 3 m / \ 2.5 m i A ~ 2 m V x = 10 m Figure E-3. Plan view of channel with gradual contraction and expansion. The purpose of this qualitative case study was to assess the ability of the model to predict the bed level changes caused exclusively by gradual changes in the width (with negligible secondary flow effects). Since this is a fictitious case, there are no measured data for comparison; however, the computed bed profiles agree reasonable well with what is expected from observations in natural alluvial channels, as discussed later. 181 The evolution in time of the bed profiles can be seen in Figure E-4 up to a final time of 180 minutes. The initial and final longitudinal profiles of water surface and velocity are shown in Figures E-5 and E-6 respectively. Initially (/ = 0), the contraction generates a large backwater effect upstream (Figure E-5). The reduction in channel width generates intense flow acceleration along the gradual contraction. Similarly, as the channel recovers its original width along the gradual expansion, the flow rapidly decelerates (Figure E-6). The immediate response of the alluvial bed is scour along the contraction, because here the flow acceleration generates a large positive sediment transport gradient dqjdx while, along the expansion, the opposite happens; a large negative value of sediment transport gradient -dqjdx produces an area of intense sediment deposition (Figure E-4). 10.06 10.04 l — i — i — i — i — r t = 0 — t = 20 min - t = 90 min 1 = 120 min 1 = 180 min 10 15 20 25 Distance along centerline (m) 30 35 Figure E-4. Temporal evolution of bed elevation. 182 10.225 10.220 CONTRACTION EXPANSION 10.195 10 15 20 25 Distance along channel centerline (m) Figure E-5. Smoothing of water surface slope caused by scour and deposition. CONTRACTION EXPANSION 0 10 15 20 25 30 35 Distance along channel centerline (m) Figure E-6. Smoothing of the longitudinal mean velocity. Deposition along the channel expansion is temporary. The sediment initially deposited comes from the upstream portion of the channel subject to intense scouring. The initial bed deposit forms a dune-like bedform that gradually migrates and disappears in the downstream direction. The bed elevation has changed as i f it were trying to achieve a straight water surface slope and constant velocity (Fig. E-5 and Fig. E-6), which is precisely the behaviour observed in natural sand bed rivers (Chang 1988, Martin-Vide 2002). Martin-Vide (2002) states that natural alluvial rivers tend to maintain a uniform water surface slope by changing their bed elevations i f necessary. Chang (1988) uses an extremal river regime hypothesis to argue that rivers try to minimize their stream power expenditure per unit length, which is equivalent to a linear (straight-line) water surface profile along the channel; he uses the example of San Dieguito River passing through a bridge contraction to demonstrate this effect using the model F L U V I A L 12. The apparent tendency of sand bed rivers to modify their bed elevation zb in order to achieve quasi-uniform conditions can be easily explained from the ID sediment continuity equation: dzh 1 dq — L = — (E.l) dt (I-A.) dx Bed changes in time dzbldt will continue as long as sediment transport rate qs changes in the downstream direction: dqJdx^O. Since qs is a function of velocity U, scour happens in areas of flow acceleration (dU/dx > 0, e.g. contraction) and deposition in areas of flow deceleration (dUldx < 0, e.g. expansion), as observed in Figure E-4 for t = 20 min. As scour and deposition proceed, velocity becomes more and more uniform (dUldx « 0, Fig. E-6), the gradient of sediment transport drops (dqjdx « 0), so bed level changes practically stop (dzb/dt « 0, Fig. E-4). The bed has achieved an equilibrium condition. Since water surface slope is a function of flow velocity (S ~ U2), uniform velocity (dUldx « 0) means constant slope (Fig. E-5), in agreement with the observations of Chang (1988) and Martin-Vide (2002). 184 B E D F O R D C H A N N E L (FRASER R I V E R ) Bedford Channel is a small branch of the Fraser River, located at Fort Langley, British Columbia, Canada. This side channel is located between Fort Langley on the left bank and McMillan Island on the right bank. It has a length of about 4.2 km measured from Tavistock Point downstream to Endsleigh Point upstream (Fig. E-7). Its width varies between 110 and 190 m, being generally wider in the upstream part. Bedford Channel is used mainly for recreational and fishing purposes. Because of the wider section, the upstream reach of the channel is subject to intense deposition that hampers navigability during low flows in the Fraser River. The average bed elevation in this area is -3 m for the Fraser River and -2.1 m for Bedford Channel. A qualitative numerical simulation was perform to determine if the current bed profile could be reproduced by the model starting from an initial flat bed. The modelling area extended from about 0.9 km m downstream of Tavistock Point to about 3.5 km upstream from Endsleigh Point (Fig. E-8). The mesh was formed by 4477 nodes and 8009 elements, with most of the elements located inside Bedford Channel, especially around Tavistock and Endleigh Points (Fig. E-8). The initial bed elevation was assumed as -3.0 m. The flow rate was assumed constant and close to the mean annual flood in the Fraser River, Q = 10,000 m3/s. The downstream water surface elevation was assumed fixed at 4.0 m. Roughness height was assumed constant ks = 0.5 m. Sediment transport rate was computed using the Meyer-Peter Muller equation for a sediment size of 1.65 mm (this value is not the real sediment size in the river). The model ran with an initial time step of 100 s, which was gradually increased up to 3600 s. The final simulation time adopted was 115 days. 185 a. Figure E-7. Bedford Channel in the Fraser River, British Columbia, Canada, (a) Aerial photograph (after nhc 2003) and (b) plan view of the centerline. 186 Distance BEDFORD CHANNEL Qin = 10000.0000 Figure E-8. Computational mesh for Bedford Channel. Initial bed level Computed - - Observed Pipeline 1000 2000 3000 4000 Distance from Tavistock Point (m) 5000 Figure E-9. Comparison between measured and computed bed profiles along the centerline of Bedford Channel. 187 The initial runs were unstable because of the bed oscillations generated around the corners (similar to Fig. E - l and Fig. E-2). To prevent these bed oscillations from growing, a bed smoothing algorithms was applied every time step. The smoothing algorithm averages the bed elevations zb over the entire computational domain, dampening spikes in the computed bed elevations. The smoothing algorithm works by using the nodal bed elevation to compute the average bed elevation inside the element, then the element elevation is used to recompute the nodal bed elevations. This averaging process between nodes and elements is repeated several times, eventually the bed is smoothed and all spikes removed. The computed bed levels along Bedford Channel were smooth but consistently lower than observed. This was because the starting bed level of -3 m is lower than the average bed level of Bedford Channel of -2.1 m. Therefore, the computed bed profile was raised +0.9 m. The results are shown in Figure E-9. Despite all the assumptions made for this application (some of them not very realistic), the agreement between the computed and observed profile as shown in Figure E-9 is remarkable, especially in the upstream reach. In the downstream reach the model predicts bed levels lower than currently observed. However, the model is not considering the presence of a buried gas pipeline (Figure E-9) which has a riprap protection that prevents the natural extend of bed scour. The bed-smoothing algorithm is useful in dampening spikes, but unfortunately it also smears the bed, reducing the transverse bed slopes. Preliminary tests of this algorithm in the curved L F M flume (Fig. 14) prevented the point bar and pool from adequately developing. Therefore, it is not recommended for further use. An algorithm for smoothing the bed changes Azb was used instead for the next application. Even with the smoothing algorithm, the time step had to be severely restricted (At = 100 s at the beginning and At = 1 h by the end), the model required several days to run. This computational time is excessive i f compared with the 30 minutes of CPU time required by the Waal River application (Fig. 21). Because of the large CPU time no more tests were attempted for Bedford Channel. 188 CONTINUOUS SINUSOIDAL EXPANSION AND CONTRACTION Tsujimoto (1987) developed an analytical model to predict the bed deformation in an idealized alluvial channel with continuous expansion and contraction, produced by changing the channel width B(x) according to a sinusoidal function B(x) = B0-Bl sm(kx) ; kx = 2n/L (E.2) An example of the experimental transverse bed profiles in a channel with varying width according to equation (E.2) is shown in Figure E-10, along with the theoretical predictions of Tsujimoto's (1987) analytical model. The profiles show rich bed variability. The transverse concavity of the bed profiles switches from downward across the narrowest section (central bar at expansion, kx = TC/2) to upward across the widest section (contraction, kx = 3TC/2). The main characteristics of the profiles are predicted by Tsujimoto's analytical model; however, across the expansion two small lateral bars are present (Fig. E-10) which were not predicted by his model. The only experimental data provided by Tsujimoto (1987) are related to channel plan geometry: B0 = 0.19 m, B\ = 0.04 m and L = 0.40 m. No more details of the experiment are given such as sediment size, inflow, downstream water level, bed slope, roughness or sediment transport rate, preventing a quantitative comparison with any particular experiment. The following values were adopted for the numerical model: Q = 0.0056 m3/s, h = 0.07 m, S0 = 0.001, ks = 0.0014. The computational mesh can be seen in Figure E - l 1, the nodes are roughly 1 cm apart. The initial runs exhibited stability problems close to the walls, typical of the problems with "corners" explained before (Figures E - l and E-2). Also, it was noticed that the bed shear deviation angle computed by the secondary flow correction was rather high. To stabilize the model, a smoothing algorithm was applied to the bed changes Azb. The coefficient A of the 189 secondary flow correction (equation 3.5) was arbitrary assumed as A = 2 (its value computed from equation (3.6) should be around A « 10 for this case) to reduce the deviation angle. b a 4 Figure E-10. Transverse bed profiles in channel with sinusoidal varying width. The widest section is located at kx = 3n/2, and the narrowest section at Ax = n/2 (after Tsujimoto 1987). The results of the computed bed elevations are shown in Figures E - l 1 and E-12. The time step was At = Is and the final time 1000 s. Although some oscillations were still present (notice rugged bed elevation in Figure E - l 1), the general features of the bed topography agree with the experiments. By comparing Figures E-10 and E-12, it is evident that the numerical model was able to reproduce the changes in the concavity of the transverse bed profile. The model predicted the central bar in the contraction; as well as the two lateral bars in the expansion, which the analytical model of Tsujimoto ignored (Fig. E-10). It was noticed that the direction a of bedload transport over the lateral bars oscillated from one time step to the next. That oscillation was eliminated when either the secondary flow correction or the transverse bed slope correction was turned off. It seems that the compound effect of both corrections is causing the problem. One probable reason might be that around the peak of the 190 bar, the sign (orientation) of transverse bed slopes changes, causing a sudden change in a. This is a topic for further research. i5w%?5S"SS"S&1 W W Bed Elevation CONTINUOUS EXPANSION AND CONTRACTION Qin = 0.0056 Time 16m 40.000s Figure E-l 1. Plan views of FE mesh and smoothed computation of bed elevations for channel with sinusoidal varying width. 191 0.015 E —•—Expansion —•— Contraction -0.15 0 0.15 Transversal distance (m) Figure E- l2 . Transverse bed profiles computed by numerical model at expansion (kx = 3TC/2) and contraction (kx = n/2). SUMMARY In this Appendix, some additional qualitative tests have been performed to illustrate the potential of the model to simulate the influence of varying channel width on the bed topography of alluvial channels. The results are promising, as demonstrated by the qualitative comparison between model predictions and observations. Particularly interesting are the results in a channel with continuous expansion and contraction. The rich topographic variability, with oscillating longitudinal bed profile, central and lateral bars, and alternating concavity in the transverse bed profiles are all captured by the model. It seems that the model has the necessary ingredients to reproduce the most important bed features of alluvial channels in the two horizontal dimensions (longitudinal and transverse). 192 However, some stability problems related to changes in the flow direction along the side walls are still preventing the application of the model to its full extent. These instabilities have been somewhat controlled by smoothing algorithms with some degree of success, but they do not eliminate the cause of the problem, leading to excessive computational times or rugged bed profiles. More research is still needed to find a suitable solution to this problem. 193 APPENDIX F CONVECTION-DOMINATED APPLICATIONS: MIGRATING B E D F O R M S SEDIMENT CONVECTION AND DIFFUSION The sediment continuity (Exner) equation is a hyperbolic partial differential equation (PDE) that physically can represent two mechanisms of sediment movement: convection and diffusion (Jansen et al. 1979, Przedwojski et al. 1995), also referred as translation and dispersion by Lisle et al. (2001) and explained in Figure F - l . In some cases, one of the mechanisms dominates over the other. For instance, for bed aggradation and degradation due to changes in the upstream inflow of sediment (see Chapter 2), several researchers have reduced the Exner equation into a diffusion equation, which is a type of parabolic PDE (Jansen et al. 1979, Soni et al. 1980, Gi l l 1980, Jain 1981, Gil l 1983a, b, Jaramillo and Jain 1984, Lu and Shen 1986, Gill 1987). 194 A, Dispersing wave ^ B. Dispersing wave with merging delta C. Translating wave D. Translating and dispersing wave Figure F - l . Idealized sediment waves according to Lisle et al. ( 2 0 0 1 ) . The one dimensional sediment continuity equation dt dx (F.l) where zb is bed elevation, qs sediment transport rate, X porosity, / time and x downstream coordinate, can be simplified according to Lisle et al. ( 2 0 0 1 ) as dt (\-X) dx' dx dx + . (F.2) K is an empirical constant that depends on the sediment transport equation used (Lisle et al. 2001 used Meyer-Peter Muller), F is the Froude number and h is water depth. The unspecified 195 terms at the right hand side (RHS) are unsteady flow terms which are small for F < 1 (Lisle et al. 2005). The RHS of equation (F.2) has two terms that express the influence of diffusion and convection. The first term (d2zb I dx2) represents the influence of diffusion (dispersion) on the bed changes, and it is similar to several previous diffusion models used for aggradation and degradation (Jansen et al. 1979, Soni et al. 1980, Gi l l 1980, Jain 1981, Gil l 1983a, b, Jaramillo and Jain 1984, Lu and Shen 1986, Gi l l 1987). The second term that contains the Froude number is related to convection (translation). According to equation (F.2) aggradation (dzb I dt > 0) is related to bed profiles concave upwards (82zb / dx2> 0); while degradation (dzb I dt <0) to bed profiles with downward concavity (d2zb I dx2 < 0), which agrees with the profiles observed and computed in Chapter 2 for bed aggradation (Fig. 6) and degradation (Fig. 8). Moreover, for the knickpoint migration experiment also discussed in Chapter 2, the initial downward concavity of the bed profile at the knickpoint should lead to scour (degradation) at that point according to equation (F.2), as in fact observed in the experiments and verified by the numerical model (Fig. 2-7). It seems that all the experiments in straight flumes simulated in Chapter 2 were diffusive, with convection playing a minimum or negligible role. For the curved channels subject to secondary flow simulated in Chapters 3 and 4, it is not straightforward to demonstrate which mechanism, diffusion or convection, plays the major role. The bed level changes in bends are characterized by the presence of bars, which are large bed forms that scale with the channel width. Bars can be classified as free or forced (Seminara and Tubino 1989, Tubino and Seminara 1990). Free bars arise from an instability of the alluvial bottom of a channel. The bars develop in an alternating sequences of riffles and pools separated by diagonal fronts, commonly referred as alternate bars. Alternate free bars migrate in the downstream direction at a speed smaller than the flow velocity. Forced bars arise from forcing effects, such as channel curvature (bends) or variations of channel width. In alluvial 196 bends, curvature induces a secondary flow which causes regions of depositions at the inner bends and scours at the outer bends with a structure similar to that of alternate bars except that they are steady and do not migrate. Based on the translational nature of free bars, it might be hypothesized that for free bars convection should be important, but not for forced bars. Federici and Seminara (2003) demonstrated that the instabilities that give rise to free alternate bars have a convective nature. Therefore, it is likely that convection also plays a minor role in the bend alluvial simulations discussed in Chapters 3 and 4. N U M E R I C A L I M P L I C A T I O N S The physical mechanisms of diffusion and convection have a strong impact on the mathematical nature of the sediment continuity equation (i.e. hyperbolic or parabolic PDE) and hence on the numerical scheme used to solve the equation. It can be concluded, based on the morphological applications shown in Chapters 2, 3 and 4, that the Finite Element river morphology model developed in this research is suitable to simulate bed level changes in straight and curved channels not dominated by convection. However, as explained in Appendices A and C, the numerical discretization based on the conventional Galerkin Finite Element Method (GFEM) using a lumped mass matrix should not be applied in problems where convection dominates (Donea and Huerta 2003, Carey and Oden 1988, Hirsch 1987, Chung 2002). Such problems may arise in practice when large bedforms migrate in the downstream direction (e.g. alternate free bars, sediment pulses, prograding deltas). It will be shown in this Appendix that the model breaks down for convection-dominated applications. Some suggestions to overcome this problem are also briefly discussed at the end, as a topic for further research. 197 To illustrate this limitation of the model, a simple analysis using Finite Differences to solve simplified versions of the one-dimensional sediment continuity equation will be introduced in the next section. Later, some simulations using River2D-Morphology for migrating bedforms and the problems that arise will be briefly discussed. C E N T E R E D A N D U P W I N D FINITE D I F F E R E N C E S (FD) Diffusive problem For aggradation and degradation due to overload or deficit of upstream sediment supply it is commonly assumed that the Exner equation can be approximated by a diffusion model or parabolic PDE (Jansen et al. 1979, Soni et al. 1980, Gil l 1980, Jain 1981, Gi l l 1983a, b, Jaramillo and Jain 1984, Lu and Shen 1986, Gi l l 1987). Assuming sediment transport as a function of flow velocity U: qs = al/ ^ = k^- , k= (F.3) dt dx1 35 0 (1-A) where k is the diffusion coefficient and So the initial bed slope. Notice that this equation is similar to (F.2) i f only the first term in the RHS remains. Equation (F.3) will be discretized in space and time (Fig. F-2) using an explicit Finite Difference (FD) second order central (SOC) scheme "+l ^ ^ ^ ' ^ ~2Z"m + Z * V , ) ^ ( F ' 4 ) 198 The SOC is referred to as central difference because information from both upstream ( M ) and downstream (/+1) is used to compute the bed elevations; it is also second order accurate in space. The scheme is explicit in time because the bed elevation at the new time level n+\ is computed from known values at the previous time step n. Time n +1 n Ax Ax -< • At i-1 i i +1 Space Figure F-2. Space and time discretization using Finite Differences. Equation (F.4) was applied to a hypothetical straight alluvial flume of L = 700 m with parameters a = 0.0004, b = 5 (Engelund-Hansen equation), U= 2 m/s, X = 0.4 and So = 0.0001. which led to a diffusion coefficient k = 355.6 m2/s The boundary condition for Zb(x,f) were assumed for simplicity as * ' ( < W = 0.001 m/s. * ^ - 0 m/s dt dt fort>0 (F.5) 199 Park and Jain (1984) suggested that, for stability reason, the time step At and space interval Ax should be limited by A , < ^ (F.5) 2k The time and space intervals were set to At = 1 s and Ax = 20 m. The time step is 1.8 times larger than the value given by equation (F.5), but the model remained stable. Figure F-3 shows the time evolution of the computed aggradation bed profiles along the upstream part of the channel. 0.60 0.30 - f -0 50 100 150 200 x(m) Figure F-3. Bed profiles computed every 5 s using a FD second order central (SOC) scheme for a hypothetical diffusive aggradation problem. Figure F-3 demonstrates that a SOC scheme produces stable and wiggle-free solution for diffusive problems. 200 Convective problem Assuming again sediment transport as a function of flow velocity only: qs = al/, the sediment transport gradient becomes ^ = abUb~x ™-dx dx (F.6) For a straight channel with constant unit discharge q = Uh dq dx * =-abUb-^^ = - a b U l ' 9 h h dx h dx (F.7) Water depth h is the difference between the water surface elevation H and bed elevation zb, or h = H-zb. The water surface slope is then 5 = dH/dx, reducing (F.7) to dx h > V dx j = -ab Ub( F2 dzb C2 Ig dx (F.8) In equation (F.8) S has also been expressed as a function of the Froude number F, Chezy roughness C and gravitational acceleration g. Assuming that the Froude number is small or the water surface slope is negligible S = dH/dx « 0, the Exner equation reduces to: dzh abUb dz, dt (\-X)h dx (F.9) or 201 dt dx (F.10) Which is a hyperbolic PDE known as the wave equation with celerity c = bqs I ((\-A)h). Jansen et al. (1979) propose a more general wave equation taking into account also the water continuity and momentum equations, steady-state and negligible hydraulic resistance dzh c dzh dt (l-F2) dx = 0 (FM) It can be seen that for small values of F (the assumption made to simplify equation F.8), equation (F.l 1) reduces to (F.10). Equation (F.10) wil l be discretized using two FD schemes: SOC, expressed in equation (F.12) and First Order Upwind (FOU) expressed in equation (F.l3) Zbi ~ Z b i + C ~ ~ Ax C 7 n _ 7" Zb(i+\) Zb(i-\) V (F.12) ZbiX - Z b i +C~^(Zb(>) ~Zb(i-\)) (F.l 3) Equation (F.l3) is referred as "upwind" or "upstream" because it uses information from the upstream node (z'-l); but not from the downstream node (z'+l) as the SOC does (Fig. F-2). 202 Both SOC and FOU schemes were applied to a hypothetical test case of sand wave propagation in a straight channel, loosely inspired in the qualitative test case of Cui et al. (2003b). Initially, a large parabolic bedform 2 km long and 1 m high is centered at x = 0, in a channel with an average slope So = 0.0001. The hydraulic parameters were assumed as U = 2 m/s, h - 5 m, which gives a low value of the Froude number F = 0.28 (for that value {/(l-F2) = 1.08, meaning that equations (F.10) and (F.l 1) are practically identical). The sediment parameters were assumed as a = 0.004, b = 5 and A = 0.4; which gives a celerity c = 0.021 m/s. The Courant stability condition imposes a limit in the feasible combination of time step At and space interval Ax in relation with celerity c c — <1 (F.14) Ax The time step was set to At = 2000 s and the space interval to Ax = 100 m, satisfying equation (F.14). The computed bed profiles can be seen in Figure F-4. From the results shown in Figure F-4 it is evident that the SOC scheme is unstable, while the FOU is stable. The SOC scheme produces strong numerical oscillations which eventually turn the model completely unstable (only the first time steps are shown in Fig. F-4). The FOU scheme remains stable, although the tails of the parabola are becoming progressively smoother. The purpose of this exercise is to demonstrate that, although central schemes are stable for diffusion equations (Fig. F-3), they cannot be applied to solve convective problems, which require numerical schemes with some "upwinding" properties. The explanation is that in a wave information propagates in a front moving in the downstream direction at a certain speed. When the central scheme is applied, information from the downstream side, where the front might not yet have arrived, is used in the computation of the front, which is physically incorrect. Only information from upstream or upwind of the front should be used. 203 SECOND ORDER CENTRAL -2000 -1000 0 1000 2000 3000 x(m) FIRST ORDER UPWIND 1.4 i x(m) Figure F-4. Computed bed profile every 10,000 s using central and upwind FD schemes for a convective problem (migrating bedform). 204 A P P L I C A T I O N S OF R I V E R 2 D - M O R P H O L O G Y The fact that convective problems require upwind numerical schemes is quite relevant for the present research because the G F E M discretization used by River2D-Morphology lacks upwinding properties (Appendix A). The following examples show the instability problems that arise when River2D-Morphology is applied to solve the Exner equation in cases when convection dominates. Sand Wave This application is similar to the parabolic sand wave simulated before using FD. The results using River2D-MOR are shown in Figure F-5. Distance (m) Figure F-5. Time evolution of computed bed profiles every 10 days using River2D-MOR for a hypothetical migrating sand wave. 205 The results shown in Figure F-5 are not directly comparable with Figure F-4 because River2D-M O R solves the full system of equations derived from sediment continuity plus water continuity and momentum instead of the highly simplified equation (F.10) solved by the FD schemes shown in Figure F-4. However, it is evident that the G F E M method used by River2D-M O R is unstable for this application. As the sediment waves advances, its leading front becomes increasingly steep, becoming more similar to a shock wave. After a certain point, the front breaks down and strong oscillations take over the solution. Free alternate bars This qualitative application is loosely based on the experimental run PI505 on bar formation in a straight flume performed by Lanzoni (2000) and numerically simulated by Defina (2003) using a two-dimensional FE model. The experimental flume was 55 m long and 1.5 m wide covered with uniform sediment with median size Dso = 0.48 mm. The flow rate was 0.03 m3/s and the normal depth 0.044 m. The longitudinal slope was 0.452% and measured sediment transport rate 94.5 litres per hour, which was mainly transported as bedload. The experiment started with an initial plane bed, which was first altered by the presence of two-dimensional ripples. Later, free alternate bars developed in the second half of the flume. Defina (2003) used a much larger flume length of 117.5 m for her numerical experiments to better capture the alternate bars in the downstream reach of the flume. She introduced an initial bottom disturbance (small bump) a few millimetres high to trigger the development of the bar train. For the present simulation, the length of 117.5 m was also adopted. The flume was discretized using nodes spaced roughly at 0.2 m, making a total of 3255 nodes and 5617 triangular elements. The roughness height was assumed as 0.0086 m. An empirical sediment transport equation with a = 0.0118 and b = 5 was adopted to match the observed sediment transport rate for the average flow velocity. Porosity was set to X = 0.4 and time step At = 1 s. 206 It was found that alternate bars develop even from an initially perfectly flat bed, with no upstream perturbation. A typical result is shown in Figure F-6. The computed alternate bars were observed to slowly migrate in the downstream direction. Figure F-6. Early development (t = 908 s) of alternate bars from an initial plane bed computed by River2D-MOR (distorted scale). As demonstrated by the results shown in Figure F-6, the model is able to capture the main features of the free alternate bar pattern, which is the alternating sequences of riffles and pools separated by diagonal fronts migrating downstream. Unfortunately, only the early stages of the free bar development can be simulated. As the simulation progresses in time, spurious bed oscillations destroy the bars, as shown in Figure F-7. By t = 1169 s, oscillation similar to the ones observed in Figure F-5 start to appear, but in a two-dimensional pattern. The oscillations continue to grow and make the model break down. 207 Depth FREE BARS, t = 908 s Qin = 0.0300 Qout = 0.0300 10.045 0.045 Time 19m 29.000s Figure F-7. Contour plots of computed free alternate bars at two different time steps: t = 908 s and t= 1189 s. By observing the results shown in Figures F-5 and F-7, it seems the model is able to simulate the initial stages of the bed evolution but, at certain point in time, strong spurious bed oscillations caused by the lack of upwinding take over the solution. This is also true for the FD SOC scheme shown in Figure F-4. In short, the model is able to simulate forced bars driven by curvature effects as demonstrated in Chapters 3 and 4, but it cannot simulate free migrating bars because of the G F E M employed. 208 Sediment delta into Check Dam Arminini and Larcher (2001) performed laboratory experiments of sediment deposition upstream of open slit check dams. These types of dams are built in steep mountain streams to contain debris and sediment flows. Open dams have windows or vertical slits that allow a portion of water and sediment to pass through the dam, delaying the water-sediment hydrograph. Since flows in mountain streams are usually supercritical, a hydraulic jump may form upstream of the dam (Fig. F-8) causing sediment deposition. The experimental runs were performed in a 10 m long and 0.5 m wide channel with a steep slope of 5%, covered with uniform sediment with median size Dso = 3.0 mm transported mainly as bedload. Sediment moved up to the jump, where it deposited in the subcritical reach. The depositional pattern propagated in both directions, upstream in the supercritical reach and 209 downstream in the subcritical reach. The sediment front and the hydraulic jump advanced downstream until it practically reached the location of the dam. Busnelli et al. (2001) applied a one-dimensional Finite Volume morphological model with transcritical flow capabilities to reproduce some of the experiments of Arminini and Larcher (2001). For a steady-state case the flow rate was 6.68 L/s, the normal depth 0.0178 m. Since their numerical model solves for subcritical boundary conditions (known upstream inflow and downstream water surface elevation), a fictitious pool with subcritical flow conditions was introduced upstream, from where a S2 profile develops towards the supercritical normal depth. The downstream boundary condition was defined assuming critical flow conditions for a vertical slit with an opening of 0.10 m. The sediment transport equation selected was the Meyer-Peter Muller equation. The upstream part of the channel was assumed as fixed bed. This test case, where transcritical flow occurs, is a fairly challenging problem for 2D models. Most morphological models with transcritical flow capabilities are limited to ID only (Spinewine and Zech 2002, Papanicolaou in press, Cui and Wilcox in press, Cui et al. in press a,b). River2D-MOR was applied to compute first the initial steady-state flow using a similar approach to the fictitious pool of Busnelli et al. (2001). As shown in Figure F-8, the slope of the first 2 m of the channel was assumed flat. Roughness height was assumed to be 0.872 mm and the downstream water depth as 0.077 m. The Froude number in the supercritical reach was 1.8 and the hydraulic jump was located at 8.8 m from the inlet, in reasonable agreement with the experimental values. The results of the morphological simulation are shown in Figure F-9. Unfortunately, the results are very poor and do not agree with the experimental data. The advancing sediment front coinciding with the location of the hydraulic jump was qualitatively correct, but its extension was significantly under-predicted. In the experiment, after 19 minutes, the sediment deposit covered the last 6 m of the flume, with a maximum height of almost 0.1 m. In the numerical model, the sediment deposit hardly extends over 1 m, even after 60 minutes of simulation (Fig. F-9). 210 The bed profiles exhibit notorious spurious oscillations and, even worse, a large erosion hole is predicted in the upstream end of the sediment pulse, which was not observed in the experiments (Fig. F-9). The scour hole produces a corresponding dip in the water surface profile and a second hydraulic jump, which might be in part responsible for preventing the normal growth of the sediment deposit. Figure F-9. Computed bed and water surface profiles for check dam simulation. The under-prediction of sediment deposition upstream of the dam and the scour hole have also been observed in qualitative tests under full subcritical flow conditions (not shown). Therefore, it is unlikely that the poor performance of the model is due to the occurrence of transcritical flow conditions; it is more likely a consequence of intense convection, as suggested by the migrating nature of the sediment front. The oscillations may be a consequence of the convection dominated bed front, as suggested by Bradford and Katapodes (2000). 211 S U M M A R Y The sediment continuity (Exner) equation, which governs the morphodynamics of alluvial channels, can be interpreted as the combination of two physical transport mechanisms: diffusion (dispersion) and convection (translation). In certain cases, one of those mechanisms may dominate over the other. For highly diffusive problems, the so-called "centered" numerical schemes, such as central Finite Differences or Galerkin Finite Elements Methods (GFEM), are appropriate to solve the Exner equation. The applications shown in Chapter 2 (bed aggradation/degradation, knickpoint scour) and Chapters 3 and 4 (forced point and alternate bars driven by curvature) were successfully simulated by a G F E M because in those applications convection was not dominant. However, for alluvial channels where the translation or migration of bedforms is important (e.g. free alternate bars, advancing sediment fronts, sediment waves); in other words, where convection is the dominant mechanism, the G F E M breaks down and becomes unstable. This is because convection requires the so-called "upwind" numerical schemes that exploit the upstream or upwind information to compute the fate of the moving front. One possibility to correct this problem could be replacing the traditional Galerkin scheme with a type of Petrov-Galerkin scheme. For example, the Streamline Upwind Petrov-Galerkin scheme (SUPG) currently used by the hydrodynamic model of River2D is an attractive possibility. Further research is needed to study the most appropriate scheme. The applications explored in this Appendix could be used as future test cases for that purpose. 212 APPENDIX G QUALITATIVE SIMULATIONS OF DYKE B R E A C H AND DAM BREAK This Appendix shows some qualitative exploratory tests that were made to assess the potential of the model as a tool to conduct future research related to the rapid morphological changes caused by dyke breach or dam break phenomena. These are fairly challenging problem because of their rapid and transient nature, normally involving transcritical flow and sharp moving fronts. Most hydrodynamic models, even popular one-dimensional models such as HEC-RAS, are intended for gradually varied flow and cannot deal with this type of hydraulic problem. Given the exceptional transcritical flow capabilities of River2D-Morphology, it has the potential to become a useful tool for future research in this field that could eventually lead to a practical engineering program for applications to dam breaks over movable beds or similar applications. Two idealized test cases in straight channels are explored as detailed below: dyke breach caused for flow overtopping and sudden dam break. 213 D Y K E B R E A C H This idealized test case is a modification of the knickpoint migration experiment simulated in Chapter 2. Instead of having the channel split into two succesive reaches at different elevations separated by the over-steepened reach (see Fig. 10), the channel reach upstream from the knickpoint was lowered, such that the over-steepened 10%-slope reach became the downstream face of a dyke, as shown in Figure G - l . Only a short 0.10 m reach upstream from the knickpoint was left to represent the top of the dyke. The upstream face of the dyke was also made with a 10% slope reach. A l l other conditions remain the same as in the original simulation (see Chapter 2 for details), except for the inflow discharge. Two runs were made, one with the discharge Q = 0.488 L/s and the other with no inflow Q = 0 L/s. It is assumed that at time t = 0 water is already flowing over the dyke, the same assumption made for the knickpoint test case. Figure G - l . Sketch of initial conditions for idealized dyke breach simulation. 214 Constant inflow Q = 0.588 L/s The computed bed and water surface profiles at three different simulation times are shown in Figure G-2. At t = 0, downstream from the knickpoint the water surface profile is identical to the previous knickpoint simulation (Fig. 9). However, upstream the water surface approaches an almost horizontal backwater profile ( M l curve), as typically found upstream of spillways or weirs. Figure G-2 shows that most of the changes happened during the first minute of the simulation, when most of the height of the dyke was eroded, which caused a release of a large volume of the water impounded upstream. The additional flow released by the initial breach of the dyke flooded the downstream reach, causing a transient increase in downstream water levels that partially drown the hydraulic jump and pushed it downstream, as shown in the curve label "W.S 1 min" in Figure G-2. The transient increase in the flow rate during the first minute can also be seen in Figure G-3. The peak discharge quickly doubles the initial value. After about 4 minutes, the discharge is almost back to its original value. The discharge hydrographs at two different locations, upstream and downstream from dyke, are shown in Figure G-2 to illustrate the phase lag between them. After intense and rapid changes during the first minutes, the sediment is slowly transported downstream. By t = 160 min, water levels have receded considerably upstream, the dyke has become a large bump in the bottom of the channel. Notice in Figure G-2 that the location of the hydraulic jump and the sediment front coincide. This is because the high sediment transport capacity along the supercritical reach sharply drops as the flow become subcritical past the hydraulic jump (Busnelli et al. 2001, Arminini and Larcher 2001, Belial et al. 2003) In practice, the downstream increase of flows is one of the most important effects of breaches in dams and dykes. The increased flow rates are usually accompanied by severe flooding and damage downstream. 215 Figure G-2. Computed bed and water surface (W.S.) profiles for dyke breach simulation. 1.4 i 0.0 4 - - , - , , , 0 60 120 180 240 Time (s) Figure G-3. Unsteady flow conditions caused by the breach both upstream (x = 10 m) and downstream (x = 13 m) from dyke. 216 No upstream inflow Q = 0 L/s The computed bed and water surface profiles at three different simulation times are shown in Figure G-4. At / = 0, the initial water surface profile is identical to the previous test case (Figure G-2), afterwards the changes in water and bed levels are also similar (Figure G-2) but slightly different. The long term bed changes in this case are less pronounced because there is no inflow of water to maintain sediment transport. After 60 s bed changes almost stop as water levels and velocity drop towards zero. The flood hydrograph also shows a decreasing tendency as shown in Figure G-5. 217 10 10.5 11 11.5 12 12.5 13 Distance (m) Figure G-4. Computed bed and water surface (W.S.) profiles for dyke breach simulation with no inflow (Q = 0). 1.2 n 0 60 120 180 240 Time (s) Figure G-5. Unsteady flow conditions caused by the breach both upstream (x = 10 m) and downstream (x = 13 m) from dyke, with no upstream inflow (Q = 0). 218 D A M B R E A K A sketch of the idealized dam break test case is shown in Figure G-6. In a straight channel 10 m long and 1 m wide, it is assumed that a vertical dam is located at the middle (x = 5 m), spanning over the entire width. A horizontal sediment layer of 1 m is assumed to exist upstream from the dam. The initial water levels are 1 m and 3m and upstream from the dam respectively (loosely inspired by Leal et al. (2001) and Leal et al. (2003) experiments). The upstream inflow was Q = 1 m3/s. The mesh was made of 1408 nodes and 2566 elements; nodes spacing was about 0.10 m. Sediment properties were assumed similar to the previous tests. In the model the vertical walls of water and sediment at the dam were simulated as inclined between the nodes of a single element. At t = 0 the dam is assumed to break suddenly. DAM Q = 1 m3/s ^ > \ 2 m / * iir"iiiii iiiiiiiiii^Si:^"iHsSsjsZj?sjs&j}-sjsjs* ^ «^ *0*\>*}»'5»,*V.»VV»V,«V,»V,V^ •*-^\\^\y\^\\\\\i*\y\\\^ / v 1 m i Figure G-6. Initial configuration of the idealized dam break problem. The computed bed profiles every second, for the first 8 seconds of the simulation are shown in Figure G-7. The sediment layer upstream from the dam is rapidly eroded and deposited as an 219 advancing front downstream. The speed of the front seems to increase with time. By 8 seconds the front has reached the downstream end of the flume (where it got reflected). The computed bed profiles upstream from the dam (erosion) are smooth. However, the depositional front advancing downstream exhibits some oscillations as shown in Figure G-7) and more clearly in Figure G-8. These oscillations are similar to those observed in convective problems (Appendix F), such as the deposition upstream from a check dam (Fig. F-9). They are probably of similar origin. Figure G-9 shows both the bed and water surface profiles at 5 different times and Figure G-10 shows the flood hydrograph downstream from the dam at x = 6 and x = 9 m. The maximum water level downstream from the dam was 1.8 m (Fig. G-9), while the peak discharge exceeded 4 m3/s. The flood wave moved very fast, traveling 4 m (JC = 9 m) in about 0.8 s (Fig. G-10). As observed before, the sediment front traveled downstream in conjunction with the hydraulic jump (Fig. G-9), which is in agreement with the analysis of Ferreira at al. (2001). 0 2 4 6 8 10 12 Distance (m) Figure G-7. Evolution of bed profiles every second for the dam break test case. 220 Figure G-8. Bed topography at t = 0 and t = 8 s for the dam break test case. 5.0 -, 0 1 2 3 4 5 6 7 8 Time (s) Figure G-10. Downstream flood hydrographs caused by dam break. S U M M A R Y Preliminary tests on the rapid morphological changes caused by dyke overtopping or sudden dam break led to reasonable results, at least at the qualitative level. The highly transient flood water and sediment hydrographs were simulated by the model, although the accuracy of the prediction is unknown. Some bed oscillations were observed in advancing sediment fronts, probably because of the convection problems explained in Appendix F. Nevertheless, the results suggested that the model developed in this research is a promising tool for future research regarding dam break or similar rapid varied flows over alluvial beds. 222
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Two-dimensional finite element river morphology model
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Two-dimensional finite element river morphology model Vásquez, José Alfredo 2006
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Two-dimensional finite element river morphology model |
Creator |
Vásquez, José Alfredo |
Publisher | University of British Columbia |
Date Issued | 2006 |
Description | Alluvial streams are those that flow over a movable bed that they have formed by the transport of sediment, either suspended in the water current or dragged along the bottom. Bedload represents the fraction of that sediment that is transported along the channel bottom. Since bedload transport is associated with the removal (erosion) or deposit of sediment over the bed, it is the main transport mechanism responsible for the morphological changes in alluvial channels. Channel changes made by sediment deposition increase the risk of river flooding, block navigation canals, reduce the life span of reservoirs and modify the course of streams. While erosion can undermine the foundation of built structures (e.g. bridges, pipelines, bank revetments), scour riparian habitats or induce drawdown in the water table levels, just to mention a few effects. For those reasons, bedload transport in alluvial streams is a topic of high interest for both river engineers and fluvial morphologists. In this thesis, a two-dimensional numerical model based on the Finite Element method was developed for the morphodynamic simulation of scour and deposition in straight and curved alluvial channels. The model uses unstructured meshes formed by elements of triangular shape. The model was developed by coupling a solver for the bedload sediment continuity equation to existing fixed-bed hydrodynamic models to dynamically update the channel bed elevation as sediment is transport by the flow. The model was successfully applied to reproduce experimental observation of bed degradation and aggradation along straight flumes; as well as scour and deposition along curved channels. Data from a meandering river was also applied to test the model for full-scale applications. For the tests performed, the model proved stable and accurate. There are two versions of the model. A conventional vertically averaged (VA) model and a novel vertically averaged and moment (VAM) model. The VAM is a quasi-3D flow model that provides more detail of the vertical structure of the flow than a VA model for curved channels. The model has some unique capabilities for modeling rapidly varied flows and transcritical flow (supercritical flow and hydraulic jumps) over alluvial beds. This is probably the first two-dimensional river morphology model based on unstructured triangular Finite Elements that has combined capabilities for bedload transport, secondary flow effect in bends, transverse bed slope and transcritical flow. The results of this thesis suggest that the model has the potential for future morphodynamic applications to dam-breaks, dam removal and steep mountain streams. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0063259 |
URI | http://hdl.handle.net/2429/18195 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2006-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2006-130415.pdf [ 17.41MB ]
- Metadata
- JSON: 831-1.0063259.json
- JSON-LD: 831-1.0063259-ld.json
- RDF/XML (Pretty): 831-1.0063259-rdf.xml
- RDF/JSON: 831-1.0063259-rdf.json
- Turtle: 831-1.0063259-turtle.txt
- N-Triples: 831-1.0063259-rdf-ntriples.txt
- Original Record: 831-1.0063259-source.json
- Full Text
- 831-1.0063259-fulltext.txt
- Citation
- 831-1.0063259.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0063259/manifest