A TWO DIMENSIONAL HYDRODYNAMIC RIVER MORPHOLOGY AND GRAVEL TRANSPORT MODEL by Stephen Kwan PhD., University of Liverpool, 1997 MCS., Regent College, 2008 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (VANCOUVER) July 2009 © Stephen Kwan, 2009 Abstract A depth-averaged two-dimensional hydrodynamic-morphological and gravel transport model is presented. Based on River2D by Prof Peter Steffler et al. (2002) of the University of Alberta, River2D-Morphology (R2DM) is capable of simulating flow hydraulics, sediment transport for uni-size and mixed size sediment, and morphological changes of a river over time. Changes in bed elevation are calculated by solving Exner’s equation for sediment mass conservation with the mixed size sediment transport function of Wilcock and Crowe (2003). R2DM is capable of exploring the dynamics of grain size distributions, including fractions of sand in sediment deposits and on the bed surface of rivers and streams. Results have been verified with experimental data for aggradation and downstream fining (Paolo, 1992; Seal, 1992; Toro-Escobar, 2000); for degradation (Ashida and Michiue, 1971); and qualitatively with two dimensional flow (Tsujimoto, 1987). Application of R2DM to actual rivers has shown that it is ready for general release as a tool to assist engineers with river restoration projects, the design of hydraulic structures and with fish habitat modeling. ii Table of Contents Abstract ...................................................................................................................... ii List of Figures ............................................................................................................ vii List of Equations ....................................................................................................... xii List of Symbols ......................................................................................................... xiv Acknowledgements ................................................................................................ xvii 1.0 Introduction ..........................................................................................................1 1.1 Thesis Objective and Outline ..................................................................................... 2 2.0 River2D .................................................................................................................3 2.1 River2D Hydrodynamic Model .................................................................................. 4 2.2 Limitations................................................................................................................. 6 2.3 River2D Morphology Version 1 ................................................................................ 6 2.3.1 R2DM’s Secondary Flow Correction Algorithm ................................................. 9 2.3.2 Limitations of River2D Morphology (Version 1) .............................................. 11 3.0 River2D Morphology Version 2 ............................................................................ 16 3.1 Simulation of Moving Bed-forms ............................................................................ 16 3.2 Mixed Size Sediment Transport............................................................................... 20 3.3 Threshold of Movement in Mixed Sized Sediments ................................................ 22 3.3.1 Gradation Independence ................................................................................. 22 3.3.2 Equal Mobility .................................................................................................. 23 3.3.3 Experimental Observations.............................................................................. 23 3.4 The Gravel Transport Algorithm ............................................................................. 30 iii 3.5 Calculation of Roughness Height, Ks ....................................................................... 33 3.6 Calculation of Radius of Curvature ......................................................................... 33 3.7 The Graphical User Interface .................................................................................. 39 4.0 Results: Mixed Size Sediment Transport Model ................................................... 42 4.1 Verification of the Transport Model ....................................................................... 42 4.1.1 Simulation Procedure ...................................................................................... 42 4.1.2 Results .............................................................................................................. 43 4.2 Aggradation in a Straight Flume ............................................................................. 45 4.2.1 Simulation Procedure ...................................................................................... 47 4.2.2 Results of Narrow Channel Simulations, Runs 1-3 .......................................... 48 4.2.3 Results of 1D Wide Channel Simulations, Runs 4 – 6 ...................................... 54 4.2.4 Results of 2D Wide Channel Simulations, Runs 4 – 6 ...................................... 57 4.2.5 Discussion of Results ........................................................................................ 58 4.3 Degradation in a Straight Flume ............................................................................ 59 4.3.1 Simulation Procedure ...................................................................................... 59 4.3.2 Results .............................................................................................................. 59 4.3.3 Discussion......................................................................................................... 62 4.4 Dispersion of a Sediment Pulse ............................................................................... 63 4.4.1 Simulation Procedure ...................................................................................... 64 4.4.2 Results and Discussion for Run 2 ..................................................................... 64 4.4.3 Results and Discussion for Run 3 ..................................................................... 69 4.4.4 Summary .......................................................................................................... 73 5.0 Two Dimensional Simulations Using R2DM .......................................................... 74 5.1 Scour and Deposition in a Curved Flume ................................................................ 74 5.2.1 Simulation Procedure ...................................................................................... 74 5.1.2 Results .............................................................................................................. 75 iv 5.1.3 Discussion......................................................................................................... 77 5.2 Scour and Deposition in a Variable Width Flume ................................................... 78 5.2.1 Simulation Procedure ...................................................................................... 79 5.2.2 Results and Discussion ..................................................................................... 79 5.3 Application of R2DM to a Natural River ................................................................. 81 5.3.1 Simulation Procedure ...................................................................................... 81 5.3.2 Results and Discussion ..................................................................................... 82 5.4 Summary ................................................................................................................. 83 6.0 Conclusions and Further Research ....................................................................... 84 6.1 Conclusion ............................................................................................................... 84 6.2 Thesis Contribution ................................................................................................. 84 6.3 Further Research ..................................................................................................... 85 6.3.1 Two Dimensional Validation Tests ................................................................... 85 6.3.2 Stability Problems ............................................................................................ 85 6.3.3 Non-Erodible Areas .......................................................................................... 86 6.3.4 Large Scale Simulations of Rivers..................................................................... 86 6.3.5 Modeling of Suspended Load .......................................................................... 87 7.0 References .......................................................................................................... 88 Appendix: Running R2DM ......................................................................................... 92 Step 1: Run a Steady State Hydrodynamic Simulation ................................................. 92 Step 2: Specify the Boundary Conditions for a Morphological Simulation ................... 92 Step 3: Select the Wilcock and Crowe Sediment Transport Equation ........................... 93 Step 4: Run Morphology ............................................................................................... 94 Step 5: Viewing Results with R2DM .............................................................................. 94 Vector Data ............................................................................................................... 94 v Contour Plots ............................................................................................................ 95 Output Files ............................................................................................................... 96 vi List of Figures Figure 1. Flow chart of the computational procedure of R2DM ....................................... 7 Figure 2. Observed and simulated bed profiles at 40 minutes due to sediment overloading four times the equilibrium level (after Vasquez et al., 2005). .................................. 8 Figure 3. Observed and simulated bed profiles at 10 hours for bed degradation experiment due to sediment shut off for two correction coefficients Cr=1.0 and Cr=0.4 (after Vasquez et al., 2005). .......................................................................... 8 Figure 4. Diagram illustrating how secondary flow is formed (after Vasquez, 2005). ...... 9 Figure 5. Typical cross section at a natural river bend (after Henderson, 1966). ............. 10 Figure 6. Deviation of bed shear stress around bends. ..................................................... 10 Figure 7. Illustration of the stability problems associated with migrating sand wave in R2DM version 1.0. .................................................................................................. 12 Figure 8. Sediment flux indicated by circle is a mathematical error caused by the secondary correction algorithm. .............................................................................. 14 Figure 9. Sediment flux vectors if secondary flow correction is not used. ....................... 14 Figure 10. Bed shear deviation angle for a section of the Fraser River using the velocity gradient based curvature algorithm (Equation 10). ................................................. 15 Figure 11. Example of sediment flux into and out of an element in the mesh. ................ 17 Figure 12. Simulation of migrating bed-form using up-winding. In this simulation an upwinding coefficient of 0.7 is used............................................................................ 19 Figure 13. Effect of varying the up-winding coefficient on the simulation of a moving bed-form. In this case an up-winding coefficient of 0.3 is sufficient to damp out the oscillations. .............................................................................................................. 19 Figure 14. The particle-weight effect: large particles are harder to move because they are heavier (after Southard, 2006). ................................................................................ 21 Figure 15. The hiding effect: larger particles are fully exposed to the flow whereas smaller particles are sheltered by the larger particles (after Southard, 2006). ........ 21 Figure 16. The rollability effect: larger particles can easily roller over smaller ones whereas small particles cannot easily roll over larger ones (after Southard, 2006). 21 Figure 17. Graph of τci/τc50 vs Di/D50 for gradation independence and equal mobility .... 24 vii Figure 18. Variation of reference shear stress as a function of sand content. ................... 25 Figure 19. The effect of sand on the transport rate of sediment. ...................................... 26 Figure 20. Similarity collapse of all fractional transport observations. Source: Wilcock and Crowe (2003), used with permission. ............................................................... 27 Figure 21. Flow chart showing calculation of bed load transport rate............................. 29 Figure 22. Diagram to illustrate sediment flux entering and leaving element with area, AE ................................................................................................................................. 30 Figure 23. Conceptual model of a gravel bedded river during aggradation. The bed load mixes with the surface layer to form a new gravel size distribution. The surface layer thickness is assumed to be a constant thickness Ls (after Cui, 2007). ............ 32 Figure 24. Conceptual model of a gravel bed river during degradation. The surface layer mixes with the sub-surface layer to form a new grain size distribution. The surface layer thickness is assumed to be a constant thickness Ls (after Cui, 2007). ............ 32 Figure 25. Cumulative discharge in a “wavy” flume with total discharge = 12 l/s. ......... 33 Figure 26. Illustration of stream tube (after Steffler et al. , 2002). ................................... 35 Figure 27. Discharge across an elemental stream tube (after Steffler et al., 2002). ......... 35 Figure 28. Diagram to show three points that lie on the same streamline. ....................... 37 Figure 29. Bed shear deviation angle in a typical meandering river (flow is from right to left). When the channel bends towards the right the bed deviation angle is positive. When the channel bends to the left the bed deviation angle is negative. ................ 38 Figure 30. Bed deviation angle in a section of the Fraser River calculated using the streamline based curvature algorithm. ..................................................................... 39 Figure 31. “Sediment Options” dialogue box used to set up a morphological simulation in R2DM ...................................................................................................................... 40 Figure 32. “Set Grain Size Distribution” dialogue box used for specifying additional information for a mixed size sediment transport simulation. .................................. 41 Figure 33. Experimental set up used by Wilcock et al. to measure transport rates for different sand/gravel mixtures at various constant water discharge rates. .............. 42 Figure 34. Section of the one dimensional mesh used to simulate the Wilcock et al. (2001) experiments.............................................................................................................. 43 viii Figure 35. Plot of dimensionless transport function W* vs τ/τr for Wilcock et al. (2001) experiments and R2DM simulations. J06, J14, J21 and J27 represent mixtures containing 6, 14, 21 and 27% sand content. ............................................................ 44 Figure 36. Grain size distribution of the feeds used for the SAFL experiments (after ToroEscobar, 2000). ψ =1 corresponds to the sand/ gravel boundary (D=2 mm) .......... 46 Figure 37. Experiment setup of the SAFL experiments (after Cui, 2007) ........................ 47 Figure 38. 2D mesh used for SAFL wide channel runs 4 and 5. ...................................... 48 Figure 39. Simulated bed profiles shown with experimental results for run 1 shown with final water surface elevation. ................................................................................... 49 Figure 40. Simulated bed profiles shown with experimental results for run 2 shown with final water surface elevation. ................................................................................... 49 Figure 41. Simulated bed profiles shown with experimental results for run 3 shown with final water surface elevation. ................................................................................... 50 Figure 42. Simulated and observed grain size, ψ ( = log2D) for run 1 (sand included). . 51 Figure 43. Simulated and observed gravel grain size (ψ) for run 2. ................................ 51 Figure 44. Simulated and observed gravel grain (ψ) sizes for run 3 (sand excluded). .... 52 Figure 45. Simulated and observed sand fraction for run 1. Experimental data obtained from Cui, 2007......................................................................................................... 53 Figure 46. Simulated and observed sand fraction for run 2. Experimental data obtained from Cui, 2007......................................................................................................... 53 Figure 47. Simulated and observed sand fraction for run 3. Experimental data obtained from Cui, 2007......................................................................................................... 54 Figure 48. Simulated bed profiles shown with experimental results for run 4 (1-D simulation). .............................................................................................................. 55 Figure 49. Simulated bed profiles for run 5 (1-D simulation). ........................................ 55 Figure 50. Characteristic gravel grain sizes for run 4 (1-D simulation). ......................... 56 Figure 51. Characteristic gravel grain sizes (ψ) for run 5 (1-D simulation). ................... 57 Figure 52. Perspective views of final bed surface of (a) SAFL run 4 (after Toro-Escobar et al., 2000) and simulated results using R2DM. .................................................... 58 Figure 53. Simulated and measured scour depths for Ashida and Michiue (1971) experiment. .............................................................................................................. 60 ix Figure 54. Simulated change in sand fraction 1m and 10 m from inlet. .......................... 61 Figure 55. Simulated change in D90 1m and 10 m from inlet. ......................................... 61 Figure 56. Simulated and measured surface layer distribution at 10 m from inlet at 600 minutes (after Wu, 2001). ....................................................................................... 62 Figure 57. Grain size distribution of sediment in feed and pulse (feed uses same distribution as run 2). After Cui et al., 2003. .......................................................... 63 Figure 58. Observed long profile of average bed elevation for run 2 Source: Cui et al., 2003 (used with permission). .................................................................................. 65 Figure 59. Simulated long profile of bed elevation for run 2. ......................................... 65 Figure 60. Observed sediment transport rate for run 2. Source: Cui et al., 2003 (used with permission). ............................................................................................................. 66 Figure 61. Simulated sediment rate for run 2 using R2DM. ............................................. 67 Figure 62. Mean surface grain size for run 2. Source: Cui et al., 2003 (used with permission). ............................................................................................................. 68 Figure 63. Simulated mean surface grain size for run 2. .................................................. 68 Figure 64. Observed long profile of average bed elevation for run 3. Source: Cui et al., 2003 (used with permission). .................................................................................. 69 Figure 65. Simulated long profile of bed elevation for run 3. .......................................... 70 Figure 66. Observed sediment transport rate for run 3. Source: Cui et al., 2003 (used with permission). ............................................................................................................. 71 Figure 67. Simulated transport rate for run 3 using R2DM. ............................................. 71 Figure 68. Observed mean surface grain size for run 3. ................................................... 72 Figure 69. Simulated mean surface grain size for run 3. .................................................. 73 Figure 70. Mesh used for the “Laboratory of Fluid Mechanics” flume............................ 75 Figure 71. Scour and deposition in a curved flume for sand bed and gravel bed. ............ 76 Figure 72. Bed change for gravel bed simulation. ............................................................ 76 Figure 73. Velocity vectors of the flow and bed-load transport for gravel bed simulation. ................................................................................................................................. 77 Figure 74. Transverse bed profiles in sinusoidal varying channel of Tsujimoto (1987). After Vasquez et al., 2005. ...................................................................................... 78 Figure 75. Mesh used for simulation of sinusoidal varying channel. ............................... 79 x Figure 76. Transverse bed profiles in channel with sinusoidal varying width. ............... 80 Figure 77. Bed elevation in a sinusoidal varying flume. UW=0.9, ks=2.5, kT=0.1. ......... 80 Figure 78. The Fraser River near Chilliwack showing bed elevation and outline of water’s edge.......................................................................................................................... 81 Figure 79. Bed changes in metres after 24 hours. ............................................................. 82 Figure 80. Bed Shear deviation angle and sediment flux vectors. .................................... 83 Figure 81. The “Sediment Options” dialogue box ............................................................ 93 Figure 82. The “Wilcock and Crowe Settings: Set Grain Size Distribution” dialogue box. ................................................................................................................................. 94 Figure 83. The “Vector Plot” dialogue box ...................................................................... 95 Figure 84. The “Colour/Contour” dialogue box ............................................................... 95 Figure 85. Where to select CSV files in the “Run Morphology” dialogue box................ 96 xi List of Equations Equation 1. Conservation of mass ...................................................................................... 4 Equation 2. Conservation of momentum in x direction ...................................................... 5 Equation 3. Conservation of momentum in y direction ...................................................... 5 Equation 4. Discharge intensity in x direction .................................................................... 5 Equation 5. Discharge intensity in y direction .................................................................... 5 Equation 6. Bed load transport equation ............................................................................. 6 Equation 7. Direction of bed shear stress.......................................................................... 11 Equation 8. Bed shear deviation angle. ............................................................................. 11 Equation 9. Secondary flow correction coefficient........................................................... 11 Equation 10. Radius of curvature...................................................................................... 13 Equation 11. The wave equation ....................................................................................... 16 Equation 12. First order up-wind equation ....................................................................... 16 Equation 13. Up-winding scheme used in R2DM ............................................................ 17 Equation 14. Flux through side 12 .................................................................................... 18 Equation 15. Flux through side 23 .................................................................................... 18 Equation 16. Flux through side 31 .................................................................................... 18 Equation 17. Change in bed elevation .............................................................................. 18 Equation 18. Dimensionless shear stress .......................................................................... 22 Equation 19. Condition for gradation independence. ....................................................... 22 Equation 20. Relationship between normalized critical shear stress and normalized grain size (after Southard, 2006). ..................................................................................... 22 Equation 21. Condition for equal mobility (after Southard, 2006). .................................. 23 Equation 22. Wilcock and Crowe (2003) hiding function. ............................................... 24 Equation 23. Exponent for Wilcock and Crowe (2003) hiding function. ......................... 24 Equation 24. Shear stress .................................................................................................. 25 Equation 25. The Dimensionless Wilcock and Crowe (2003) transport function ............ 26 Equation 26. Dimensionless transport rate for fraction i .................................................. 26 Equation 27. Volume of fraction i in element with area, AE ............................................ 31 Equation 28. Volume of fraction i .................................................................................... 31 xii Equation 29. New surface layer fraction........................................................................... 31 Equation 30. R2DM bed roughness formula .................................................................... 33 Equation 31. Equation of a streamline. ............................................................................. 34 Equation 32. The stream function ..................................................................................... 34 Equation 33. Differential of the stream function. ............................................................. 34 Equation 34. Differential of the stream function in terms of flow.................................... 34 Equation 35. Discharge in elemental stream tube. ............................................................ 36 Equation 36. Change in stream function equals discharge through it............................... 36 Equation 37. Flow between 2 streamlines. ....................................................................... 36 Equation 38. Equation of a circle through points A (x1,x2), B (x2,y2) and C (x3, y3) in matrix format (Source: The Math Forum @ Drexel University). ........................... 38 Equation 39. Engelund-Hansen sediment transport equation. .......................................... 74 Equation 40. Width of flume ............................................................................................ 78 Equation 41. Suspended load sediment equation. ............................................................. 87 xiii List of Symbols A = coefficient in secondary flow correction ac = centripetal acceleration C90 = multiply D90 with this factor to get ks C = Chezy friction coefficient (m1/2/s) c = wave celerity D50 = size of 50th percentile grain size Dsm = median surface grain size Di = grain size in fraction i dz = change in bed elevation Fi = proportion of fraction i in the surface size distribution Fb = proportion of fraction i in bed load size distribution Fs = proportion of sand in surface Fss = proportion of fraction i in subsurface size distribution fbi = bed load fraction g = gravitational acceleration h = water depth ks = effective roughness height in River2D kT = transport rate factor in R2DM L = length of flume Ls = surface layer thickness qb = transport rate per unit width qbi = transport rate of size fraction i per unit width qbT = total bed load transfer rate qs = specific volumetric sediment flux qsIN = upstream sediment supply rate (m2/s) qw = specific water discharge qx = sediment flux in x direction qy = sediment flux in y direction rc = radius of curvature Sd = distance between nodes for calculation of radius of curvature xiv Sox = bed slope in x direction Soy = bed slopes in y direction Sfx = friction slope in x direction Sfy = friction slope in y direction s = specific gravity of sediment SAFL = St. Anthony Falls Laboratory t = time in s T = Final time u* = shear velocity u = velocity component in x direction v = velocity component in y direction Wi* = dimensionless transport rate of size fraction i Wr* = reference value of dimensionless transport rate (=0.002) zb = bed elevation ρ = water density γ = specific weight = ρ∗g κ = von Karman’s constant λ = porosity φ = τ/τri = -log2D ψ = log2D σ = stream function τ = shear stress τci = critical shear stress of size fraction i τr = reference shear stress τri = reference shear stress of size fraction i τrm = reference shear stress of mean size of bed surface τi* = dimensionless Shields stress for size fraction i τrm* = reference dimensionless Shields stress for mean size of bed surface τ*c50 = dimensionless critical Shields stress for mean grain size τxx = component of turbulent stress tensor in xx xv τxy = component of turbulent stress tensor in xy τyx = component of turbulent stress tensor in yx τyy = component of turbulent stress tensor in yy ξd = tailgate elevation in SAFL experiment xvi Acknowledgements I offer my enduring gratitude to Dr Robert Millar who gave me the opportunity to study here at UBC, to Dr Jose “Pepe” Vasquez for making this project possible, and to my fellow students in the hydrotechnical group, especially Arturo Haro, Marcel Luthi, Alana Smiarowski who beta tested the program. I owe particular thanks to Antony Pranata and Dennis Hodge for helping me understand the C++ programming language. I am also thankful to Dr. Gary Parker, Dr Joanna Curran, and Dr Yantao Cui for willingly giving me their data and experimental results. Finally, special thanks are owed to my wife Nelly and our children, Abigail and Ariana for their love and support. xvii 1.0 Introduction Numerical models provide the basis for understanding the hydrodynamic and geomorphic conditions in river ecosystems and are a valuable tool towards solving complex issues in river and environmental engineering. Using numerical models provides the ability to simulate possible scenarios under altered hydraulic or watershed conditions, allowing scientists and engineers to provide solutions to existing problems and to anticipate future issues before they happen. Most of the sediment transport models used in river engineering are one dimensional, especially those used for long-term simulation of a long river reach (e.g Hoey and Ferguson, 1994; Ferguson and Church, 2009). One dimensional models generally require the least amount of field data for calibration and testing. The numerical solutions are more stable and require less computational time than two dimensional simulations but are generally not suitable, however, for simulating the effects of curvature, channel width and multiple reaches or other two dimensional local phenomena. A recent example of a one dimensional model is the TUGS (The Unified Grain Sand) model developed by Cui (2007). This uses the surface-based bed load equation of Wilcock and Crowe (2003) and links the grain size distributions in the bed load, surface layer and subsurface with the gravel transfer function of Hoey and Ferguson (1994) and Toro-Escobar et al. (1996), a hypothetical sand transfer function, and hypothetical functions for sand entrainment/infiltration from/into the subsurface. It has been validated using St. Anthony Falls Laboratory (SAFL) narrow channel flume experiments and the flushing flow experiments of Wu and Chou (2003). An advantage of the TUGS model, is its ability to model the evolution of sand in gravel deposits (few models can do this). This is important to understand because the habitat of aquatic species can be seriously affected by the transport of fine sediments (e.g. in salmonid bearing rivers, adult salmonids select locations with favourable hydraulic conditions and appropriate grain size distributions to lay their eggs). A limitation of the TUGS model, though, is that it is principally a research tool and not generally available for others to use. Parker on the other hand, has made his one dimensional sediment transport models free available on the internet.1 On his website there are a variety of models that can be used to simulate aggradation and degradation on both gravel and sand bed rivers under constant flow or transient flow conditions. However, a limitation of these models is that only channels with an initial constant bed slope can be simulated and the output options are limited because the program is written as an Excel macro. The abovementioned problems are just a few reasons why we need a numerical model that is 1) capable to simulate one and two dimensional problems; 2) capable of modelling uni-size and mixed size sediment transport (especially gravel); 3) able to model the evolution of gravel and sand; 4) user friendly and 5) available as freeware. At this time only a few numerical sediment transport models have all the features. River2DMorphology version 2 is one of them. 1.1 Thesis Objective and Outline This thesis will present the technical merits River2D-Morphology (R2DM). A brief overview River2D and R2DM version 1 is presented in Section 2. Section 3 presents an introduction to the transport of mixed size sediments and how it is computed in R2DM. In Section 4 the transport model is verified using the flow measurements of Wilcock et al. (2003), the St. Anthony Falls experiments (Paola et al, 1992; Seal et al, 1997; Toro-Escobar et al, 2000), the degradation experiments of Ashida and Michiue (1971) and the sediment pulse experiments of Cui and Parker (2003). Qualitative results of mixed size sediment transport in curved flumes and natural rivers are given in Section 5. Finally, a summary and the contributions of this thesis are presented in Section 6 along with some suggestions for further work. 1 Dr Gary Parker’s code is available at http://vtchl.uiuc.edu/people/parkerg/excel_files.htm 2 2.0 River2D Numerical hydrodynamic models have become popular in recent years due to the advent of inexpensive, powerful personal computers and are now routinely used in research and engineering practice. There are many commercially available hydrodynamic models based on finite difference, finite volume, and finite element methods (e.g. Mike21, CCHE2D, Fluent PHOENICS, Star-CCM+ etc) and each numerical method has advantages and disadvantages. For instance, it has been widely argued that finite volume methods offer the best stability and efficiency while finite element methods offer the best geometric flexibility. Thus, for complex river bed topographies, finite element methods may be the most suitable. (For a more complete review see Islam, 2009). River2D (Steffler et al., 2002) is a two-dimensional, depth averaged finite element model developed for use on natural streams and rivers. It solves the basic mass conservation equation and two components of momentum conservation. Outputs from the model are two velocity components (in x and y directions) and a depth at each node. Velocity distributions in the vertical are assumed to be uniform and pressure distributions are assumed to be hydrostatic. Three-dimensional effects, such as secondary flows in curved channels, are not calculated. Additional Modules In additional to solving the hydrodynamics, River2D has an ice module and fish habitat module. The ice module, models the flow of under a floating ice cover with known geometry. This is important because ice can affect the flow hydraulics by increasing the magnitude of shear stress on the flow (Steffler et al., 2002). The fish habitat module is based on the Weighted Usable Area (WUA) concept used in the PHABSIM family of fish habitat models (Steffler et al., 2002). WUA is calculated according conditions flow conditions (velocity, depth and channel substrate) that fish species prefer. 3 Data Requirements of River2D The input data required by River2D include channel bed topography, roughness and transverse eddy viscosity distributions, boundary conditions, initial flow conditions and a mesh that can capture the flow variations. Bed Topography Obtaining an accurate representation of bed topography is the most critical, difficult, and time consuming aspect of two dimensional hydrodynamic modeling (Steffler and Blackburn, 2002) and so simple cross-section surveys are generally not adequate. The developers of River2D recommend that a minimum of one week of field data collection is required per study site followed by post processing of the field data with a quality digital terrain model before being used as input for the 2D model. Boundary Conditions Boundary conditions take the form of a specified total discharge at the inflow and fixed water surface elevations or rating curves at the outflow. It is important to locate flow boundaries some distance from areas of interest is important to minimize the effect of boundary condition uncertainties. Initial conditions are also important because they can significantly reduce the total run time and in some cases make the difference between a stable run and an unstable one. 2.1 River2D Hydrodynamic Model Like most 2D models, River2D calculates the hydrodynamics based on the 2D vertically averaged St. Venant equations: ∂h ∂q x ∂q y + + =0 ∂t ∂x ∂y Equation 1. Conservation of mass 4 ∂q x ∂ (uqx ) ∂ (vq x ) g ∂h 2 1 ∂ (hτ xx ) ∂ (hτ xy ) + + + = gh(Sox − S fx ) + + ρ ∂x ∂x ∂x ∂y 2 ∂x ∂y Equation 2. Conservation of momentum in x direction ∂q y ∂x + ∂ (uq y ) ∂x + ∂ (vq y ) ∂y + g ∂h 2 1 ∂ (hτ yx ) ∂ (hτ yy ) = gh(Soy − S fy ) + + ρ ∂x 2 ∂y ∂y Equation 3. Conservation of momentum in y direction Where h is the water depth, u and v are the vertically averaged velocities in x and y respectively, and qx and qy are the discharge intensities related to the velocity components through: q x = uh Equation 4. Discharge intensity in x direction q y = vh Equation 5. Discharge intensity in y direction Sox, Soy are the bed slopes; Sfx, Sfy are the friction slopes; τxx, τxy, τyx, τyy are the components of turbulent stress tensor; g is the gravitational acceleration and ρ is the water density. The assumptions in the equations are: 1. The pressure distribution in hydrostatic thus limiting the accuracy in areas of steep slopes and rapid changes in bed slopes. 2. the horizontal velocities are constant with depth and so information on secondary flows and circulation are not available 3. Coriolis and wind forces are assumed negligible. 5 2.2 Limitations A limitation of River2D is that it can only compute the hydrodynamics for a fixed bed. This limits application to situations where sediment transport is not significant or not of interest. In natural rivers, both the water surface and bed elevations can change considerably according to flow conditions or the presence of obstacles in the river. 2.3 River2D Morphology Version 1 R2DM was developed in 2005 by Jose Vasquez at the University of British Columbia so that morphodynamic changes of a river bed could be simulated (Vasquez, 2005). It exists as a separate module that links to the River2D hydrodynamics by solving the bed load transport continuity equation: (1 − λ ) ∂zb ∂qsx ∂qsy + + =0 ∂t ∂x ∂y Equation 6. Bed load transport equation where qsx and qsy are the components of volumetric rate of bedload transport per unit length in x and y. λ is the porosity of the bed material, t is time and zb is the bed elevation. The bedload transport, qs, can be computed using either of the four empirical equations available: Engelund-Hansen, Meyer-Peter-Müller, Van Rijn, and an empirical formula (Kassem and Chaudrey, 1998). Each equation is generally applicable for uni-size sand or gravel but not gravel-sand mixtures. To run a morphology model, it is first necessary to obtain a converged steady state solution for the hydrodynamic model using River2D. The saved output of this (the CDG file) is then used as the initial conditions for a transient morphological simulation. The additional input data required are the sediment transport equation, sediment feed rate, sediment parameters (D50 and porosity), outflow bed elevation boundary conditions and secondary flow correction parameters. At every time step during the simulation, the hydrodynamic model computes water depth, h and the velocity components, u and v. This information is then used to calculate the bed load sediment transport for mixed size sediment, qs using Equation (6). The method is shown in the flow chart in Figure 1. 6 Figure 1. Flow chart of the computational procedure of R2DM 7 Version 1 has been successfully tested in cases of bed elevation changes in straight flumes (Vasquez et al. 2005, 2006), scour and deposition in laboratory bends (Vasquez et al. 2005), meandering rivers (Vasquez et al. 2005) and knickpoint migration (Vasquez et al. 2005) under transcritical flow conditions. The results of some of these simulations are reproduced here for convenience (Figures 2-3). Figure 2. Observed and simulated bed profiles at 40 minutes due to sediment overloading four times the equilibrium level (after Vasquez et al., 2005). Figure 3. Observed and simulated bed profiles at 10 hours for bed degradation experiment due to sediment shut off for two correction coefficients Cr=1.0 and Cr=0.4 (after Vasquez et al., 2005). 8 2.3.1 R2DM’s Secondary Flow Correction Algorithm The difference between flows in straight channels and those in curved channels is the presence of a centripetal acceleration in the latter. Assuming a simple 1D streamline, the centripetal acceleration, ac = u2/rc (where rc is the radius of curvature of a streamline, and u is the velocity) increases from zero at the river bed to a maximum close to the surface of the water. At a cross section in the bend of a channel the surface layers of water move outwards due to the centripetal force (Figure 4). In order to satisfy mass conservation, water in the lower layers close to the bed must move inwards. When this secondary flow circulation combines with the primary flow it generates a 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. In natural rivers this flow pattern causes scour sediment from the outside of the curve and deposit it on the inside, forming a cross section like the one shown in Figure 5. (Henderson, 1966). Figure 4. Diagram illustrating how secondary flow is formed (after Vasquez, 2005). 9 Figure 5. Typical cross section at a natural river bend (after Henderson, 1966). Figure 6. Deviation of bed shear stress around bends. 10 The secondary flow causes the direction δ of the bed shear stress to deviate from the direction of the mean depth averaged flow velocity by a deviation angle δs (Figure 6). Since a 2D model cannot simulate this 3D effect, a semi-empirical secondary flow correction is introduced that assumes a local fully developed curved spiral flow (Struiksma et al. 1985): v u δ = arctan − δ s Equation 7. Direction of bed shear stress. where δ is the direction of the bed shear stress vector, arctan(v/u) is the direction of the depth averaged velocity, and δs is given by: h r c δ s = arctan A Equation 8. Bed shear deviation angle. A is parameter of order 10 and can be estimated as a function of von Karman coefficient, κ and Chezy’s roughness coefficient, C, and is defined by: A= g 2 1 − κ 2 κC Equation 9. Secondary flow correction coefficient. 2.3.2 Limitations of River2D Morphology (Version 1) While version 1 has been verified with flume experiments, it does have some limitations that restrict its use. The four main ones are: 1) its inability to simulate moving bed-forms; 2) its inability to simulate multiple size fractions and development of a 11 surface armour layer; and 3) its inability to correctly calculate secondary flow correction for natural rivers and 4) the absence of a user-friendly graphical user interface. Unable to Simulate Moving Bed-forms Version 1 could only simulate bed elevation changes in channels that were not dominated by advection, such as moving bed-forms (Vasquez, 2005). This was because the numerical discretization method (a conventional Galerkin Finite Element Method) used by R2DM is not capable of solving advection dominated problems (Vasquez, 2005). This means that mathematical instabilities may arise during the simulation of large bed forms that migrate downstream (e.g. alternate free bars, sediment pulses and pro-grading deltas). An example of such a simulation is shown in Figure 7 which shows the time evolution of a hypothetical sand wave and the development of mathematical instability. Figure 7. Illustration of the stability problems associated with migrating sand wave in R2DM version 1.0. 12 Simulations Limited to Uni-size Sediment Version 1 could only simulate sand-bed rivers because the traditional sediment transport equations of Meyer-Peter Muller, Engelund-Hansen, Van Rijn and the empirical formula all assume that the sediment is effectively of a single size. This presents a problem because there are many applications, particularly in gravel bed rivers, where multiple size fraction capability is necessary. Furthermore, most gravel rivers are characterized by a coarse surface armour layer which plays a significant role on regulating transport rates, limiting degradation and scour depths (Parker et al., 1982). Unstable Method to Calculate Radius of Curvature River 2D uses an unstructured mesh and so the algorithm developed by Vasquez (2005) to calculate curvature is not well suited to complex 2D geometries or natural river bed. As an example, Figures 8 and 9 show the sediment flux vectors for simulations on a natural river with and without secondary flow correction respectively. In Figure 8, it can be seen that with secondary flow correction, the vectors switch abruptly 90 degrees in different regions. This is because of the way the algorithm calculates the radius of curvature rc. Version 1 of R2DM computes rc using velocity gradients according to the formula (Vasquez, 2005): 3 (u 2 + v 2 ) 2 rc = ∂u ∂u ∂v ∂v − uv − v 2 + u 2 + uv ∂x ∂y ∂x ∂y Equation 10. Radius of curvature For initially flat beds, this procedure gives good results for simple geometries with uniform bed topography, however, in natural river beds, local velocity spikes could develop that could lead to unrealistic values of rc (Vasquez, 2009). This is evident in the example shown in Figure 8 because the bed shear deviation angle is seen to vary from -45o to 45o indicating that the calculated values of rc are either too small or indeterminate. This then creates unrealistic high values of δs (~+/-89o) in Equation 8 as rc approaches zero which the secondary flow correction algorithm limits to +/-45o (Figure 10). 13 Figure 8. Sediment flux indicated by circle is a mathematical error caused by the secondary correction algorithm. Figure 9. Sediment flux vectors if secondary flow correction is not used. 14 Figure 10. Bed shear deviation angle for a section of the Fraser River using the velocity gradient based curvature algorithm (Equation 10). No Graphical User Interface Finally version 1 was never actually released as an official version because all the parameters (e.g porosity, sediment feed rate, boundary conditions etc.) were hard coded into the program. So whenever a certain parameter was changed, the program had to be re-compiled and run from within the C++ environment (this required users to have both knowledge of C++ and a Microsoft C++ compiler). 15 3.0 River2D Morphology Version 2 The main objective of this research work was to address these four limitations of R2DM (v1). This section will discuss how these limitations were addressed. 3.1 Simulation of Moving Bed-forms The simulation of moving bed-forms (advective transport) was made possible by the introduction of up-winding into the numerical scheme. Up-winding is a form of numerical discretization methods for solving hyperbolic partial differential equations and utilize an adaptive or solution-sensitive finite difference stencil to numerically simulate more properly the direction of propagation of information in a flow field. First order schemes are the simplest methods and introduce severe numerical diffusion into the solution where large gradients exist (Patankar, 1980). For example, consider the following hyperbolic partial differential wave equation: ∂z ∂z +c =0 ∂t ∂t Equation 11. The wave equation Where c = wave celerity. z in+1 = z in + c ∆t n z i − z in−1 ∆x [ ] Equation 12. First order up-wind equation Equation 12 is considered an up-wind solution because it uses information from the upwind node (n-1) and not the downstream node (n+1). Second order schemes have improved spatial accuracy over first order schemes but are less diffusive (Patankar, 1980). Third order schemes are less diffusive compared to the second-order and introduce slight dispersive errors in the region where the gradient is high. 16 A simple first order up-winding method has been incorporated into R2DM that calculates the sediment flux in each element by taking into account the flux flowing in from neighbouring elements: Q = (1 − UW ) * Qdownstream + UW * Qupstream Equation 13. Up-winding scheme used in R2DM where 0 ≤ UW ≤ 1 is the up-winding weighting factor. UW=1 means full up-winding, and UW=0 means no up-winding is used (conventional Galerkin Finite Element Method). Figure 11. Example of sediment flux into and out of an element in the mesh. Figure 11 illustrates an example in which an element with nodes 1, 2, and 3 has sediment entering through sides 1-2 and 1-3, while leaving by side 2-3. Sediment fluxes through each element side will be computed as: 17 Q + Q3 Q12 = (1 − UW ) * 1 + UW * Q A 2 Equation 14. Flux through side 12 Q + Q3 Q23 = (1 − UW ) * 2 + UW * Q1 2 Equation 15. Flux through side 23 Q + Q3 Q13 = (1 − UW ) * 1 + UW * Q13 2 Equation 16. Flux through side 31 The average bed change ∆zb at the element during a time interval ∆t can therefore be calculated using: ∆zb = ∆t Q + Q13 − Q23 * 12 (1 − λ ) A Equation 17. Change in bed elevation where A= area of element, λ = porosity The effect of using up-winding is illustrated in Figures 12 and 13. Figure 12 shows the results of the same hypothetical sand wave as the one shown in Figure 7 but with up-winding. Figure 13 shows how the results for various values of UW for this simulation. 18 Figure 12. Simulation of migrating bed-form using up-winding. In this simulation an up-winding coefficient of 0.7 is used. Figure 13. Effect of varying the up-winding coefficient on the simulation of a moving bed-form. In this case an up-winding coefficient of 0.3 is sufficient to damp out the oscillations. 19 3.2 Mixed Size Sediment Transport R2DM version 2 is capable of simulating mixed size sediment transport using the transport algorithm of Wilcock and Crowe (2003). Mixed size sediment transport differs from uni-size sediment transport because the transport of sediment size fractions (a specified range of sizes within the size distribution) is usually dependent on other sediment size fractions. So, when dealing with mixed size sediment transport it is more appropriate to look at the fractional transport rates. The fractional transport rate, qbi, is the transport rate in mass per unit time, per unit cross-stream width of the flow for each fraction in the mixture of the transported sediment (q denotes the unit transport rate, subscript i denotes the ith fraction in the mixture and b denotes the bed load). The fractional transport rates can vary between two extremes: gradation independence and equal mobility (Parker et al., 1982). Gradation independence is when the transport rate of each size fraction is independent of the presence of all the other size fractions. So if every grain in the mixture acted as though it were surrounded by grains of similar size to itself, the critical Shields number (see Equation 18) would be the same for all grains in the surface, and the boundary shear stress required to move grain size Di would increase linearly with size Di (i.e. a grain that is twice the size of another grain would require twice the boundary shear stress to move). Equal mobility is when the ratio of the fractional transport rate of a given size fraction to the proportion of the given size fraction in the bed sediment is the same for all the size fractions. The condition of equal mobility may seem confusing because common sense would indicate that coarser, heavier particles would be more difficult for the flow to transport than the finer particles. This condition, known as the particle-weight effect, is however offset by 2 other conditions: 1) the hiding-sheltering effect which occurs when smaller particles are sheltered from the flow of the force by the larger particles which are more exposed to the flow (Einstein, 1950); and 2) rollability (Southard, 2006). This occurs because larger particles can roll over easily over the bed of smaller particles, but smaller particles cannot roll easily over the bed of larger particles. These two effects are an essential factor in mixed size sediment transport (see Figures 14-16). 20 Figure 14. The particle-weight effect: large particles are harder to move because they are heavier (after Southard, 2006). Figure 15. The hiding effect: larger particles are fully exposed to the flow whereas smaller particles are sheltered by the larger particles (after Southard, 2006). Figure 16. The rollability effect: larger particles can easily roller over smaller ones whereas small particles cannot easily roll over larger ones (after Southard, 2006). 21 3.3 Threshold of Movement in Mixed Sized Sediments The threshold for movement of a specific size fraction in a mixed-size sediment is different to that of a uni-size or very well sorted sediment. This is because the movement threshold in mixed size sediments is dependent on a combination of the hiding-sheltering effect and the rollability effect. To illustrate how the threshold of movement varies in mixed size sediments we will examine the two extremes, gradation independence and equal mobility and then look at some actual experimental observations. 3.3.1 Gradation Independence For gradation independence the threshold for each fraction is the same as if it were a uni-size sediment, so τ*ci, the Shields dimensionless stress: τ ci γDi ( s − 1) = constant for all fractions Equation 18. Dimensionless shear stress where τci= critical shear stress of fraction i , γ= specific weight, s=specific gravity, Di = size of fraction i Thus a grain of a given size D within a mixture has exactly the same mobility as it would have if the bed were composed entirely of size D. Algebraically this can be expressed as follows: τ*ci = constant for all fractions i Equation 19. Condition for gradation independence. Setting τ*ci = τ*c50 the threshold value for the 50th percentile of the sediment mixture gives: τ ci γDi ( s − 1) = τ c 50 γD50 ( s − 1) ⇒ D τ ci = i τ c 50 D50 Equation 20. Relationship between normalized critical shear stress and normalized grain size (after Southard, 2006). 22 3.3.2 Equal Mobility With equal mobility, each size fraction has the same movement threshold as each other, and so the particles move at the same threshold value of τ0: τci = constant for all fractions i Setting τci = τc50 the threshold value for the 50th percentile of the sediment mixture gives: τ ci =1 τ c 50 Equation 21. Condition for equal mobility (after Southard, 2006). 3.3.3 Experimental Observations To illustrate what actually happens in reality we can look at the experimental study on the transport of mixed size sediments by Wilcock and Crowe (2003). Their study showed that sediment mixtures behave somewhat in between gradation independence and equal mobility and can be clearly seen in Figures 17 which shows the normalized critical shear stress vs normalize grain size. Gravel smaller than D50 tend to equal mobility while those larger than D50 tend more towards gradation independence (Figure 17). 23 Figure 17. Graph of τci/τc50 vs Di/D50 for gradation independence and equal mobility Hiding Function This curve represents the “hiding function” since it shows that finer grains are hard to move (because they are hidden by the larger grains) and larger grains are easier to move (because they tend to protrude more into the flow). It can be represented by the power function: τ ci Di = τ c 50 D50 b Equation 22. Wilcock and Crowe (2003) hiding function. b= 0.67 1 + exp(1.5 − Di / D50 ) Equation 23. Exponent for Wilcock and Crowe (2003) hiding function. where Di is grain size in fraction i, D50 is the median surface grain size. 24 Variation of Shields number with sand content The results of this study also demonstrate that sand plays an important role in determining the gravel transport rate. By plotting the dimensionless Shields number versus the percentage of sand on the bed surface they found a relationship that could be expressed by the following equation: τ ci* = 0.021 + 0.015 exp(−20 Fs ) Equation 24. Shear stress Shown graphically in Figure 18, the trend shows that the critical shear stress, τ*ci decreases with increasing sand content but beyond 20% there is no more significant decrease. This suggests that the sand increases the transport rate of the sediment by either making the grain particles protrude more into the flow or by increasing the rollability effect (Figure 19). Figure 18. Variation of reference shear stress as a function of sand content. 25 Figure 19. The effect of sand on the transport rate of sediment. Transport Function When the transport rates of each experiment were plotted as a function of normalized shear stress,τ/τri (see Figure 20) a distinct trend was observed which could be expressed as: 0.002φi7.5 φ < 1.35 4.5 * Wi = 0.894 φ ≥ 1.35 14 * 1 − φ 0.5 i Equation 25. The Dimensionless Wilcock and Crowe (2003) transport function where φ = τ , Fi = proportion of i on the bed surface, u* = shear velocity, qbi = τ ri volumetric transport rate per unit width of size i, s = ratio of sediment to water density, g = gravity, and Wi* is the dimensionless transport rate defined by: Wi* = ( s − 1) gqbi Fiu*3 Equation 26. Dimensionless transport rate for fraction i 26 Figure 20. Similarity collapse of all fractional transport observations. Source: Wilcock and Crowe (2003), used with permission. 27 It should be noted that Equation (25) is simply a best fit through all the data points and for τ/τri < 5 the variation in transport rate can be as high as two or three orders of magnitude. Thus individual measured data points may need to be scaled by a transport factor kT = 0.1 to 10 times the estimate provided by Equation (25). In R2DM the procedure for calculating the mixed size sediment transport is shown in Figure 21. 28 Figure 21. Flow chart showing calculation of bed load transport rate. 29 After the total transport rate, qbT is calculated, it is then used in Equation (6) to compute the bed elevation ∆zb in a given time step dt. The bed elevation is updated as zbnew = zbold + ∆zb ; and then the surface and subsurface grain distributions are updated according to the flow of each grain size fraction into and out of an element using the gravel transport algorithm 3.4 The Gravel Transport Algorithm The gravel transport algorithm is used to update the grain size distribution of the armour (surface) layer. It assumes that the active surface layer remains at a constant thickness, Ls, set by the user (set approximately equal to 2 x D90). At each time step the grain size distribution for the surface layer distribution is recalculated according to the volume of sediment entering or leaving the element. Figure 22 illustrates how sediment enters and leaves a typical element. Figure 22. Diagram to illustrate sediment flux entering and leaving element with area, AE During aggradation, the net volume of sediment entering an element is positive (see Figure 23) and the new fraction of grain size i in the surface layer of the element can be determined: 30 Vi = (Q12i + Q23i + Q31i )dt + ( Ls − dz ) * (1 − λ ) Fi AE Equation 27. Volume of fraction i in element with area, AE where Vi = volume of fraction i in the surface layer (Ls), Fi = surface layer fraction i, and Ls = surface layer thickness (~2 x D90), Q12i = volume of sediment in fraction i entering or leaving through side 12 of the element per unit time. During degradation (∆z < 0) sediment leaves the element and the surface layer mixes with the substrate to maintain a constant Ls (see Figure 24). The new volume of fraction i is therefore: Vi = (1 − λ ) AE [( Ls − dz ) * Fi + dz * Fsi ] Equation 28. Volume of fraction i where Fsi is the fraction of grain size i in the substrate. Thus, the new surface layer fraction at time step, t+∆t is: Fsnew = Vi Vi = VTotal Ls A(1 − λ ) Equation 29. New surface layer fraction This formulation accounts for the changes in surface layer composition during aggradation or degradation and satisfies mass continuity. 31 Figure 23. Conceptual model of a gravel bedded river during aggradation. The bed load mixes with the surface layer to form a new gravel size distribution. The surface layer thickness is assumed to be a constant thickness Ls (after Cui, 2007). Figure 24. Conceptual model of a gravel bed river during degradation. The surface layer mixes with the sub-surface layer to form a new grain size distribution. The surface layer thickness is assumed to be a constant thickness Ls (after Cui, 2007). 32 3.5 Calculation of Roughness Height, Ks In River2D (the hydrodyanamic model), once the roughness height, ks, is defined in the computational domain, it remains unchanged throughout the simulation since the bed morphology remains constant. For simulations where the bed and/or surface grain distribution change over time, this method of defining ks is therefore not suitable. When the mixed Wilcock and Crowe sediment transport equation is selected in R2DM, the roughness height is recalculated after the surface distribution is updated according to the following Manning-Strickler formulation: k s = C 90 D90 Equation 30. R2DM bed roughness formula where D90 is the size of the surface material such that 90% is finer and C90 is a dimensionless constant for a given simulation (typically ~2-3.5) (Parker, 2008). 3.6 Calculation of Radius of Curvature Version 2 uses the streamlines generated by the Cumulative Discharge function of River2D to compute the radius of curvature at each node in the computational domain. The Cumulative Discharge function works by computing the streamtubes from one side of the channel where the discharge = 0 to the other side where the discharge = total discharge (Figure 25). Figure 25. Cumulative discharge in a “wavy” flume with total discharge = 12 l/s. 33 Mathematically the cumulative discharge is calculated as follows: For two dimensional flow the equation of a stream line is q x dy − q y dx = 0 Equation 31. Equation of a streamline. Defining a scalar function such that qx = ∂σ ∂y and qy = − ∂σ ∂x Equation 32. The stream function where σ(x,y) is a scalar function called the stream function. Considering the total differential of σ: dσ = ∂σ ∂σ dx + dy ∂x ∂y Equation 33. Differential of the stream function. and substituting in equation 32: dσ = q x dy + q y dx Equation 34. Differential of the stream function in terms of flow. Therefore Equation (31) = Equation (34) when dσ = 0. This also illustrates that σ is constant along any streamline. 34 Figure 26. Illustration of stream tube (after Steffler et al. , 2002). Figure 27. Discharge across an elemental stream tube (after Steffler et al., 2002). 35 In Figure 26, we have two streamlines, σ1 and σ2, separated by a distance, ds in a twodimensional flow field. The streamlines confine the flow into a stream tube since there can be no flow across a streamline. The discharge in the elemental stream tube across ds (Figure 27) can be calculated as: dQ = q x dy + q y dx Equation 35. Discharge in elemental stream tube. Therefore dQ = dσ Equation 36. Change in stream function equals discharge through it. Hence the change in σ across this elemental stream tube equals the discharge through it. If Equation (36) is integrated between σ1 and σ2, the difference in the stream function between any two streamlines is equal to the flow between the streamlines: 2 Q12 = ∫ dσ = σ 2 − σ 1 1 Equation 37. Flow between 2 streamlines. Due to the relationship between stream function and discharge, if σ1 = 0 at one side of a channel, then the value of σ2 at the other side of the channel will be equal to the total discharge across the channel. 36 Figure 28. Diagram to show three points that lie on the same streamline. The algorithm to compute the curvature from the cumulative discharge works in the following way: 1. Determine the cumulative discharge in a node. 2. Search through entire mesh for two points that lie on the same streamline separated by a distance Sd = 0.25 – 0.5 x rc (the minimum radius of curvature), one upstream of the node and one downstream (Figure 28). For nodes close to the inflow or outflow no points exist upstream or downstream respectively so the algorithm chooses an arbitrary high value to ensure that the calculated radius is high. 3. Determine the equation of the circle that passes through these three coordinates by solving the following matrix: 37 Equation 38. Equation of a circle through points A (x1,x2), B (x2,y2) and C (x3, y3) in matrix format (Source: The Math Forum @ Drexel University).2 4. Determine centre of the circle and the radius. 5. Determine direction of bed deviation calculated using Equation (8). Figure 29 shows the bed deviation angle in a typical meandering river. Figure 29. Bed shear deviation angle in a typical meandering river (flow is from right to left). When the channel bends towards the right the bed deviation angle is positive. When the channel bends to the left the bed deviation angle is negative. This method of calculating the radius of curvature from the streamlines is found to be highly stable and works for all geometries. This is illustrated in Figure 30 which shows the bed deviation angle in the same section of Fraser River as in Figure 10 but using the improved streamline based curvature algorithm. 2 Actual formula available online at http://mathforum.org/library/drmath/view/55166.html. 38 Figure 30. Bed deviation angle in a section of the Fraser River calculated using the streamline based curvature algorithm. 3.7 The Graphical User Interface The major improvement to version 1 is the addition of the “Sediment Options” dialogue box (Figure 31) to run morphological simulations. This allows the user to specify parameters such as the sediment equation, sediment properties (D50, porosity), boundary conditions, and parameters associated with two dimensional flow.3 3 For further information on these boundary conditions please refer to R2DM user manual by Kwan (2009) 39 Figure 31. “Sediment Options” dialogue box used to set up a morphological simulation in R2DM The addition of the dialogue box makes it easy to run morphological simulations without the need of editing the code and recompiling the software. It also makes R2DM readily available for other user to do river morphology simulations for uni-size or mixed size sediments If the Wilcock and Crowe selected then the user must enter additional data by clicking the “W/C settings” button. This opens the “Set Grain Size Distribution” dialogue box shown in Figure 32. 40 Figure 32. “Set Grain Size Distribution” dialogue box used for specifying additional information for a mixed size sediment transport simulation. To run a mixed size sediment transport simulation it is necessary to specify how the grain sizes are distributed in the sediment feed, the initial surface layer of the river bed, and the substrate since these distributions will affect the sediment transport and the morphological changes in the river bed. Other factors that will affect the transport rate are the transport rate factor, kT , C90 factor, and the depth of the surface layer, Ls. 41 4.0 Results: Mixed Size Sediment Transport Model 4.1 Verification of the Transport Model The computational method to calculate mixed size sediment transport was verified using the results of Wilcock et al. (2001) in their study of the transport of mixed sand and gravel. In their study, a wide range of transport rates for different sand/gravel mixtures (containing 6, 14, 21, 27 and 34% sand) in a laboratory flume were measured to determine the effect of sand on the transport of gravel. Forty eight experimental runs were conducted with water discharge rates spanning a threefold range in the experimental set up shown in Figure 33. The flume was 7.9 m long, 0.6 m wide and the gravel mixture was screened flat on the bed. Figure 33. Experimental set up used by Wilcock et al. to measure transport rates for different sand/gravel mixtures at various constant water discharge rates. 4.1.1 Simulation Procedure Using the same experimental and hydraulic conditions, these experiments were modeled using R2DM. Since the experiments were one dimensional in nature, a simple 1D mesh could be used to simulate them. The mesh used 188 triangular elements and 190 nodes with the average distance between nodes set to 0.6 m (Figure 34) and the flume was lengthened to 50 m to minimize boundary condition effects. Boundary conditions and the roughness ks were adjusted until the flow depth and shear stress matched those of Wilcock et al. The corresponding total transport rate was 42 recorded for each experiment and W* and the dimensionless transport rate was plotted as a function of bed shear stress. The reference shear stress at W* = 0.002 was found by interpolating or extrapolating the data. Figure 34. Section of the one dimensional mesh used to simulate the Wilcock et al. (2001) experiments. 4.1.2 Results The solution was not sensitive to up-winding (UW=0 and UW=0.9 gave the similar results), kT was unadjusted (= 1) and ks was varied between 0.1-0.3. All the data was collapsed onto a single curve by plotting W* vs τ/τr shown in Figure 35 with experimental results and the Wilcock and Crowe function. It can be seen that there is excellent agreement between the transport rates simulated by R2DM and experimental results. This indicates that there are no computational errors in the algorithm for the calculation of sediment transport rate for mixed size sediment. 43 Figure 35. Plot of dimensionless transport function W* vs τ/τr for Wilcock et al. (2001) experiments and R2DM simulations. J06, J14, J21 and J27 represent mixtures containing 6, 14, 21 and 27% sand content. 44 4.2 Aggradation in a Straight Flume The most comprehensive set of experiments to date on gravel transport have been performed at the St. Anthony Falls Laboratory (SAFL) in Minnesota by Parker and his co-workers (Paola et al., 1992; Seal et al., 1997; Toro-Escobar et al., 2000). The performance of the R2DM was examined by simulating the three SAFL downstream fining narrow flume runs and three wide channel flume runs (Paola et al., 1992; Seal et al., 1997; Toro-Escobar et al., 1996). These six experiments were performed at for the purpose of reproducing downstream fining of a gravel-sand mixture in response to bed aggradation. The channel used for the experiments had a depth of 1.83 m, a width of 2.74 m, length 60 m and with an initial concrete-bottom. For runs 1,2,3 and 6, the flume was narrowed to 0.305 m using a partition. The flume is ponded at its downstream reach by setting a constant water surface elevation at the downstream end, which drives channel aggradation and downstream fining. The water flow rate, Qw, sediment feed mixture, Qs, and the tailgate elevation ξd were set at constant values for each experiment; grain size of the sediment had a specific weight of 2.65. A summary of the relevant parameters for the all the runs are given in Table 1 and the sediment feed distributions are shown in Figure 27. During the course of the experiments a depositional wedge (Figure 37), having a mildly concave profile and distinct avalanche front, was formed which progressed downstream with time. Channel bed elevations were constantly monitored and the sediment was sampled at the end of each run using both the Klingeman method and a standard Wolman pebble count. 45 Table 1. Conditions for experimental runs in SAFL flume. Source: Toro-Escobar et al. (2000). Parameters Units Width m 0.3 0.3 0.3 2.7 2.7 0.3 initial slope % 0.2 1 1 1 1 1 Qw l/s 49 49 49 440 110 12.2 Hd m 0.4 0.45 0.5 0.45 0.15 0.15 Qs kg/min 11.3 5.65 2.83 25.2 35 8.5 qw m /s 2 0.163 0.163 0.163 0.163 0.041 0.041 qs 2 2.37E-04 1.18E-04 5.93E-05 5.87E-05 8.15E-05 1.78E-04 689 1379 2753 2776 500 228 33 33 33 33 55 55 m /s Run 1 Run 2 Run 3 Run 4 Run 5 Run 6 ratio of water to sediment sand content % Figure 36. Grain size distribution of the feeds used for the SAFL experiments (after Toro-Escobar, 2000). ψ =1 corresponds to the sand/ gravel boundary (D=2 mm) 46 Figure 37. Experiment setup of the SAFL experiments (after Cui, 2007) The final deposits of the SAFL experiments were wedge shaped and had a maximum thickness of approximately 0.7 m and a length of 38 m. The wedge was characterized by a steeply dipping front and sand filled the interstices of the gravel deposit. The wedge shaped deposit sloped gently downstream with a concave profile until an abrupt increase in slope at the front. The width-to-depth ratios of the flows for runs 1 – 3 were unrealistically low compared to natural rivers to ensure that no appreciable transverse bed topography would develop. Runs 4 – 6 used a wider channel of 2.7 m so that the width-to-depth ratio was more typical of natural rivers (about 18) to allow for the formation of transverse bed topography and sediment patches. 4.2.1 Simulation Procedure The mixed sediment transport feature of R2DM was developed to model river beds and so it was necessary to specify a grain size distribution for the bed even though the initial bed of the flume was concrete. To model the concrete bed, a surface distribution which consisted mostly of grains 16 mm and above was therefore used to simulate the concrete bed. The narrow flume was discretized using 318 triangular elements and 320 nodes in a simple 1D mesh similar to the one in Figure 34. The average distance between nodes was set to 0.35m. The results could be fine tuned by adjusting: the depth of surface layer 47 (Ls in Equation 21), C90, and the transport rate factor, kT. Good agreement with observed values for runs 1-3, could be obtained using the following values: Ls =0.05; C90=4; kT = 0.1-0.2, UW=0.9. The wide channel runs were discretized using 462 triangular elements and 282 nodes (Figure 38). The average distance between the nodes was set to 0.8 m. The lateral variations in the wide channel runs were introduced by altering the bed by either 1) increasing the elevation of random nodes by 0.5 mm or 2) by apply secondary flow correction. It was also found that slight lateral variations in the flow arose without making these changes because of the irregularity of the mesh. Good agreement with observed values for runs 4 and 5, could be obtained using the following values: Ls =0.05; C90=2.5; kT = 0.1-0.2 and UW=0.9. Figure 38. 2D mesh used for SAFL wide channel runs 4 and 5. 4.2.2 Results of Narrow Channel Simulations, Runs 1-3 Bed Elevation Model results of runs 1 – 3 are shown in Figures 39 – 41 for bed elevations. 48 Figure 39. Simulated bed profiles shown with experimental results for run 1 shown with final water surface elevation. Figure 40. Simulated bed profiles shown with experimental results for run 2 shown with final water surface elevation. 49 Figure 41. Simulated bed profiles shown with experimental results for run 3 shown with final water surface elevation. The bed elevations were fine tuned to fit the experimental results by adjusting the kT (transport factor) and C90 (which subsequently alters the bed roughness, ks) with the surface layer kept at a constant thickness of 0.05 m (this value is used by Cui, 2007). It was also found that different initial surface layer distributions (with sand content varying from 0 to 90%) had little effect on the morphological changes of the bed. The initial substrate distribution was set to the same values as the feed and had no effect on the results because no degradation occurred. For runs 1-3, nk =4, and kT=0.2 was found to give the best match to experimental results while. For runs 4 and 5, good matches were found with C90 =2.5, and kT=0.1. Grain Size The characteristic grain sizes along the length of the flume are shown in Figures 42 – 44. The results for run 1 were obtained from Paola et al.(1992) and give ψ without sand. runs 2 and 3 were taken from Seal et al. (1997) and give ψ values with sand included. 50 Figure 42. Simulated and observed grain size, ψ ( = log2D) for run 1 (sand included). Figure 43. Simulated and observed gravel grain size (ψ ψ) for run 2. 51 Figure 44. Simulated and observed gravel grain (ψ ψ) sizes for run 3 (sand excluded). It can be seen that R2DM can reproduce grain size distributions. Observed and simulated gravel characteristic sizes decrease in size downstream. Sand Fraction Figures 45 – 47 show the simulated and observed results for sand fractions. 52 Figure 45. Simulated and observed sand fraction for run 1. Experimental data obtained from Cui, 2007. Figure 46. Simulated and observed sand fraction for run 2. Experimental data obtained from Cui, 2007. 53 Figure 47. Simulated and observed sand fraction for run 3. Experimental data obtained from Cui, 2007. The simulated sand fractions show that they are approximately between 5-10% of observed values in both time and space with the exception of two high sand fractions for the time interval between 6 and 14 hours in the experiment. The observed high sand fraction values are samples from a small amount of sand deposit downstream of the main depositional front (Cui, 2007). 4.2.3 Results of 1D Wide Channel Simulations, Runs 4 – 6 Bed Elevation R2DM simulated results of runs 4 – 5 are shown in Figures 48 and 49 for bed elevations. R2DM could not reproduce run 6 because the ratio of water discharge to sediment feed rate was too low causing mathematically instabilities to arise after about 3 hours of simulated time. 54 Figure 48. Simulated bed profiles shown with experimental results for run 4 (1-D simulation). Figure 49. Simulated bed profiles for run 5 (1-D simulation). 55 Grain Size Figures 50 and 51 show the simulated and observed results for characteristic gravel grain sizes, ψ . The grain sizes for the 1-D wide channel runs are under-predicted by approximately a factor of 1.5. This may be due to high sediment feed rate but without further experimental data it is difficult to verify this hypothesis. Figure 50. Characteristic gravel grain sizes for run 4 (1-D simulation). 56 Figure 51. Characteristic gravel grain sizes (ψ ψ) for run 5 (1-D simulation). 4.2.4 Results of 2D Wide Channel Simulations, Runs 4 – 6 Plan views of the final bed surface for runs 4 and 5 are shown in Figure 52 alongside with the experimental results of Toro-Escobar et al. (2000). R2DM was not able to simulate run 5 in 2D because large lateral flows caused the solution to become unstable after 3 hours of simulation time. 57 Figure 52. Perspective views of final bed surface of (a) SAFL run 4 (after Toro-Escobar et al., 2000) and simulated results using R2DM. 4.2.5 Discussion of Results It has been shown that R2DM is capable of simulating the propagation of aggradation wedges, predicting both the deposition amounts and the advance rate of the wedge. Results of the wide channel simulations indicate that R2DM can reproduce bed elevations and downstream fining. The slight variation between the one dimensional and two dimensional simulations of run 4 indicate that one dimensional simulations are accurate provided that there is not much lateral variation in the flow. While the wide channel simulations demonstrate that R2DM can simulate two dimensional sediment transport, these experiments do not adequately test R2DM’s two dimensional capabilities because of the low lateral variations in the flow. 58 The fact that run 6 could not be simulated is not seen as a major problem because the sediment feed rate in this experiment is a lot higher than typical rivers. For example, in the Fraser River, the ratio of water discharge to sediment input is typically much greater than 100,000.4 4.3 Degradation in a Straight Flume Bed degradation can be observed in a flume by shutting off the upstream sediment supply. Ashida and Michiue (1971) performed such experiments in a flume to investigate bed degradation and armouring downstream of a dam. The experimental flume was 20 m long, 0.8 m wide, a slope of 1% and was filled with gravel with a median size of 1.5mm. Water was fed at the inlet at 0.0314 cubic metres per second and the outflow was kept constant at 0.06 m. 4.3.1 Simulation Procedure The flume was discretized using a one dimensional mesh of 52 triangular elements with 50 nodes. The average distance between nodes was set to 0.8 m. Good agreement with experimental results was found with UW=0, C90 =3 and kT = 5, Ls=0.07. 4.3.2 Results The experiments showed that the total scour process could be divided into two phases: 1) from 0-100 minutes, the bed was intensively scoured; 2) after 100 minutes, the scouring rate was reduced and an armouring layer was formed gradually. This has been reproduced and can be seen in Figure 53 which shows the measured and simulated scour depths at 7 m, 10 m and 13 m from the inflow (the actual results were obtained from Wu, 2001). Figures 54-55 show that over time the smaller grains leave the flume and the bed slowly coarsens up. The bed reaches equilibrium when all the small grains are taken by the flow leaving only the coarse grains behind. This can also be seen in Figure 56 which shows the simulated and measured surface layer distribution at 10 m from the inlet. 4 Using highest sediment transport rates reported by Ferguson and Church (2009) and lowest discharge on record from Water Survey of Canada at Mission. 59 Figure 53. Simulated and measured scour depths for Ashida and Michiue (1971) experiment. 60 Figure 54. Simulated change in sand fraction 1m and 10 m from inlet. Figure 55. Simulated change in D90 1m and 10 m from inlet. 61 Figure 56. Simulated and measured surface layer distribution at 10 m from inlet at 600 minutes (after Wu, 2001). 4.3.3 Discussion There is good agreement between simulated and measured results for the bed changes. The surface layer distribution at 10 m from the inlet shows good agreement for grain sizes less than 5mm but simulated grain sizes in the upper range are approximately 1.5 times larger than measured. However, with only one profile to validate the results it is difficult to comment on the accuracy of the degradation algorithm of R2DM. 62 4.4 Dispersion of a Sediment Pulse Sediment often enters mountain rivers in discrete pulses associated with landslides and debris flows. R2DM was used to simulate a sediment pulse to see whether it could capture the effects of dispersion. To observe how they evolve, Cui et al. (2003) simulated sediment pulses using a flume 45 m long, 0.5 m wide with a water discharge of 9 l/s, sediment feed rate of 45 g/min and slope of 0.0108. The following experimental procedure was adopted: 1) the bottom of the flume was covered with sediment; 2) the flume was allowed to reach a “pre-pulse” equilibrium; 3) the flow is temporarily halted so that a sediment pulse (3.5 cm high and 60 cm long for runs 2 and 3) could be formed at the upstream end of the flume; 4) the flow was started again and the flume was allowed to re-equilibriate. In this study, R2DM was used to simulate runs 2 and 3 which used a sediment pulse comprised of gravel. The grains size distribution of the pulse for run 2 was identical to that of the feed sediment and run 3 used a sediment pulse with grain size distribution that was coarser than the feed (Figure 57). Figure 57. Grain size distribution of sediment in feed and pulse (feed uses same distribution as run 2). After Cui et al., 2003. 63 4.4.1 Simulation Procedure The flume was discretized using a one dimensional mesh of 300 triangular elements with 302 nodes. The average distance between nodes was set to 0.3 m. It was not possible to simulate the same conditions as in the experiments because R2DM uses the same grain size distribution for the initial bed material throughout the computational domain. The simulation therefore used the same sediment mixture for the pulse and bed material upstream and downstream of it. Agreement with experimental results were obtained for run 2 using UW=0.9, C90 =5, kT =0.1, Ls=0.01 and for run 3 using UW=0.9, C90 =5, kT =0.1, Ls=0.01. 4.4.2 Results and Discussion for Run 2 Bed Elevation The experimental observation in Figure 58 shows that the mode of pulse deformation is dispersive with a small degree of translation. The simulated results (Figure 59) are in reasonable agreement with observed results. Better results of bed elevation could be obtained by increasing kT but this increased sediment rates two orders of magnitude more than measured values. 64 Figure 58. Observed long profile of average bed elevation for run 2 Source: Cui et al., 2003 (used with permission). Figure 59. Simulated long profile of bed elevation for run 2. 65 Sediment Transport Rate Figure 60 shows the observed sediment transport rate in g/min. Ciu et al. (2003) report that the sediment samplers were uncalibrated and so the measured values should be “interpreted as order-of magnitude estimates” (Cui et al., 2003). Observed sediment transport rates went from 15 times the sediment feed rate at 5 minutes (45 g/min) to over 40 times at 27 minutes. After 8 hours the transport rates returned to values close to the sediment feed level. Simulated results (Figure 61) were notably higher by an order of magnitude than the observed results. However, at the upstream part of the pulse, the model predicts very low values of sediment transport rates comparable to what is observed in the experiments. Similarly, at 15 m the model predicts the peak in sediment transport rate observed in the experiments. Figure 60. Observed sediment transport rate for run 2. Source: Cui et al., 2003 (used with permission). 66 Figure 61. Simulated sediment rate for run 2 using R2DM. Grain Size Figures 62 and 63 show the arithmetic mean size of the surface layer for observed and simulated results respectively. A direct comparison of the results cannot be made the simulation used the same sediment material for the pulse and initial bed (this is a limitation of R2DM v2). The experimental results show that at 36 min the surface layer becomes significantly finer than the initial condition (identified as pre-pulse on the Figure 64) and then returns to the initial condition after 5 hours 12 mins. The simulated results also show a similar trend: at 36 min the surface layer becomes finer and then gradually coarsens up. 67 Figure 62. Mean surface grain size for run 2. Source: Cui et al., 2003 (used with permission). Figure 63. Simulated mean surface grain size for run 2. 68 4.4.3 Results and Discussion for Run 3 Bed Elevation Figure 64 shows the long profiles of bed elevation for run 3. After 7 hours the pulse has not completely been dissipated. It is evident that the coarser pulse material used in run 3 has caused the dissipation rate to slow down considerably. This is also evident in the simulated results shown in Figure 65. Figure 64. Observed long profile of average bed elevation for run 3. Source: Cui et al., 2003 (used with permission). 69 Figure 65. Simulated long profile of bed elevation for run 3. Sediment Transport Rate Figures 66 and 67 show the sediment transport rates of observed and simulated results respectively. Similar to run 2, the simulated results were notably higher by an order of magnitude than the observed results. 70 Figure 66. Observed sediment transport rate for run 3. Source: Cui et al., 2003 (used with permission). Figure 67. Simulated transport rate for run 3 using R2DM. 71 Grain Size Figures 68 and 69 show the arithmetic mean size of the surface layer for observed and simulated results respectively. Observed results show that there is not much significant variation in the mean surface grain size along the length of the flume and also over time. Simulated results show that the surface is becoming coarse over time with slight downstream fining. Figure 68. Observed mean surface grain size for run 3. 72 Figure 69. Simulated mean surface grain size for run 3. 4.4.4 Summary R2DM has been successfully used to capture the effects of dispersion of a sediment pulse. The results are not ideal but this was expected because of the following limitations of the model: 1. a different grain size distribution and for the sediment pulse and surrounding bed material could not be modeled; and 2. the thickness of the sediment pulse could not be specified; The model over predicted transport rates by approximately one order of magnitude. This could be reduced by decreasing kT but it would also have significantly reduced the changes to bed elevation. 73 5.0 Two Dimensional Simulations Using R2DM The advantage of R2DM over other models is its ability to model two dimensional sediment flux in curved and variable width channels. Three examples are given in this section but only qualitative analyses were performed since no experimental data was found for gravel transport in curved or shaped flumes. 5.1 Scour and Deposition in a Curved Flume There are very few published results of experiments that observe the effects of gravel transport in curved flumes and so it is not possible at this time to verify the 2D aspects of mixed size sediment transport with R2DM. Therefore a qualitative analysis of a simulation of a gravel bed in a curved flume will be compared to a simulation that uses a sandy bed for which experimental data exists. 5.2.1 Simulation Procedure The simulation uses the Laboratory of Fluid Mechanics (LFM) flume (Struiksma, 1985) to reproduce the results of Vasquez et al. (2005b). The flume was discretized using a mesh of 3831 elements with 2150 nodes (Figure 70). The bed elevation of both the inflow and outflow sections was kept constant and to minimize boundary effects the length of the straight sections were increased 10 m. The simulation was modeled using the Engelund-Hansen (Engelund and Hansen, 1967) sediment transport equation with D50=0.78 mm, β=0.4, and UW=0. qs = 0.05C 2 g g ( s − 1) D503 τ *1.5 Equation 39. Engelund-Hansen sediment transport equation. where g=gravitational force, τ*=Shield’s number, C=Chezy friction factor. For the gravel bed simulation, the parameters used were: UW =0, kT=1, Ls = 0.05, and C90=2.5 with the “SAFL” sand poor grain size distribution shown in Figure 36. 74 Figure 70. Mesh used for the “Laboratory of Fluid Mechanics” flume. 5.1.2 Results The results of the sand bed and gravel bed simulations are shown in Figure 71. Since it appears that there is not much change in bed elevation for the gravel bed simulation, bed change for the gravel bed simulation is shown in Figure 72. The results of both simulations are qualitatively similar but the sand bed simulation shows a greater bed change than the gravel bed simulation because the fine sediment (D50=0.78 mm) is more mobile than gravel. 75 Figure 71. Scour and deposition in a curved flume for sand bed and gravel bed. Figure 72. Bed change for gravel bed simulation. 76 5.1.3 Discussion The results show that the model is correctly simulating the sediment transport. The curvature generates a secondary flow in the transverse direction from the outer bank, where scour occurs to the inner bank where a point bar forms. R2DM uses an algorithm to account for this secondary flow by making the sediment transport direction deviate from the depth-averaged flow velocity direction. This is shown in Figure 73 which shows the velocity vectors of the flow and bed-load transport. Figure 73. Velocity vectors of the flow and bed-load transport for gravel bed simulation. 77 5.2 Scour and Deposition in a Variable Width Flume Experiments by Tsujimoto (1987) in a sinusoidally varying channel show that, at equilibrium, bed topography varies quite significantly along the channel. Figure 74 shows the transverse bed profiles at the widest and narrowest sections of a channel that varies according to: B ( x) = B0 − B1 sin(ωx) Equation 40. Width of flume where ω = 2π / L where B0 is the mean channel width of the sinusoidal channel, B1 is the amplitude of the sinusoidal oscillation, and ω is the wave number, and L is the sinusoidal channel wavelength. Figure 74. Transverse bed profiles in sinusoidal varying channel of Tsujimoto (1987). After Vasquez et al., 2005. 78 5.2.1 Simulation Procedure For this simulation, the following parameters were used: Q = 0.056 m3/s, slope = 0.001, Ls = 0.05, ks = 2.5, kT = 0.1, the “SAFL” initial grain size distribution (Figure 36) and a simulation period of 1000 s using a time step of 1 s. The computational mesh, the channel was assumed to be 2 m long containing three complete wavelengths and two 0.40 m long sections with constant 0.23 m widths upstream and downstream. The node spacing was approximately 5 cm in the upstream and downstream sections and 1 cm in the variable width region. The total number of nodes was 2957 with 5571 elements (Figure 75). No secondary flow correction was used because the channel is straight with no bends to cause the generation of secondary flows. Figure 75. Mesh used for simulation of sinusoidal varying channel. 5.2.2 Results and Discussion The computed results at t=1000 s using at time step of 1 s are shown in Figures 67and 68. In contrast to the results of Vasquez (2005b), no smoothing of the bed topography was necessary. By comparing Figures 76 and 77 it is clear that the numerical model can 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. 79 Figure 76. Transverse bed profiles in channel with sinusoidal varying width. Figure 77. Bed elevation in a sinusoidal varying flume. UW=0.9, ks=2.5, kT=0.1. 80 5.3 Application of R2DM to a Natural River R2DM can be used to model the sediment transport in the Fraser River during a high flow event of 6000 m3/s. Figure 78 shows the original bed elevation of a section of the Fraser River near Chilliwack. Figure 78. The Fraser River near Chilliwack showing bed elevation and outline of water’s edge. 5.3.1 Simulation Procedure The bed file was created using data imported from a GIS software package with a resolution of 20 m. The mesh consisted of 12340 elements with 6400 nodes. The simulation used the following parameters: UW =0, kT=0.1, Ls = 0.05, C90=2.5 with the “SAFL” sand poor initial surface layer distribution (Figure 36), and a simulation period of 24 hours using a time step of 60 s. Secondary flow correction was used but was not perhaps necessary because the curvature is very low (rc~ 500 – 1000 m) for the majority of the reach. 81 5.3.2 Results and Discussion The bed changes after a period of 24 hours are shown in Figure 79. The results show that the bed change is minimal in the majority of the simulated reach but there are a few areas where the bed degrades by up to 0.154 m and one area where the bed aggrades by 0.27 m. In Figure 80, the sediment flux vectors are shown for upper section of the reach where the channel bends to the left. Since the radius of curvature is high (~ 500 m) the secondary flow effects are minimal causing the sediment vectors to deviate only slightly (< -8o) from the main flow. These results are reasonable but, in the absence of detailed experimental results, it is not possible to comment on the accuracy of this simulation. Figure 79. Bed changes in metres after 24 hours. 82 Figure 80. Bed Shear deviation angle and sediment flux vectors. 5.4 Summary R2DM has been successfully used to capture the effects of scour and deposition in two dimensional flows. More specifically: 1. In curved flumes the secondary flow correction algorithm correctly simulates the transport of sediment from the outer bank to the inner bank. 2. The effects of contraction and expansion in a channel can be simulated. 3. The curvature algorithm can accurately calculate the radius of curvature of bends on natural rivers. 83 6.0 Conclusions and Further Research 6.1 Conclusion In this thesis, the river morphology model developed by Vasquez (2005) was upgraded to model mixed size sediment transport. This model has been successfully applied to straight alluvial channels, curved alluvial channels, and to natural rivers. In Section 4, the one dimensional simulations showed R2DM to be stable and reliable and able to capture the effects of 1) downstream fining; 2) bed aggradation and degradation; 3) and dispersion and translation of mixed size sediment. In Section 5, the two dimensional simulations in curved flumes demonstrate R2DM’s ability to simulate scour and deposition of gravel beds in bends and channels with variable width. The model successfully replicated point bars and scour pools along the bends. The results of the analysis of the Fraser River suggest that R2DM can be used to simulate aggradation and degradation in gravel bedded rivers provided that there is adequate data to calibrate the model. 6.2 Thesis Contribution As a result of this study, a two dimensional finite element depth averaged model for river morphology and gravel transport was developed with the following features: 1. the ability to model uni-size and mixed size sediment; 2. the ability to model the evolution of gravel and sand; 3. the ability to model sediment sorting; 4. the ability to model moving bed-forms using up-winding; 5. a new method to calculate the radius of curvature in unstructured meshes; and 6. a user friendly graphical user interface. However, R2DM does have limitations that need to be addressed through future work. These include: 1. stability problems associated with wet and dry areas; 84 2. stability problems associated with sharp bends in side walls; 3. it cannot explicitly model bank erosion; 4. all the bed is considered erodible since there is no provision yet for non-erodible layers (such as rocky outcrops and concreted areas); 5. all the bed has to be initialized to the same grain size distribution; and 6. the model is only valid for bed-load (suspended load is ignored). The aforementioned limitations indicate the early stage of development of this model and it is expected that these problems will be resolved with further research. 6.3 Further Research 6.3.1 Two Dimensional Validation Tests This work has shown that the two dimensional simulations of mixed size sediment in curved and shaped flumes show good qualitative results. However, this is not sufficient to sufficiently validate the model and so further tests must be performed in order to confirm R2DM’s ability to model two dimensional mixed size sediment transport. Unfortunately, the limited amount of available experimental work on gravel transport in shaped flumes prevents us from doing this. Further research work would therefore involve performing experiments in variable width flumes, curved flumes and meandering flumes. This work will also enable us to further refine the gravel transport module so that it can better capture the effects of vertical sorting and downstream fining. 6.3.2 Stability Problems In all the models, except the Fraser River simulation, water covered the entire computational domain and no dry areas existed. Under these conditions the model was stable and produced good results. However, for low flow conditions some areas of the computational domain can become dry which can subsequently cause computational instability and create unusual bed topography. Similarly, it was found that sharp bends in side walls can also cause instabilities to occur. Both these problems were inherited from 85 version 1 and are preventing application of the model to its full extent. While smoothing algorithms can be coded into the model to prevent such instabilities, it does not solve the root cause of the problem and so further research is needed to find a solution. 6.3.3 Non-Erodible Areas In R2DM versions 1 and 2 all the bed is considered to be erodible and there is no provision to model non-erodible areas such as rock outcrops or where the natural bed has been replaced by concrete. A feature that allows for non-erodible areas will be created in the next version of R2DM. 6.3.4 Large Scale Simulations of Rivers Once the stability issues have been resolved, it will be possible to use R2DM to model bed morphology and gravel transport in large scale rivers to determine how the bed will respond to different scenarios such as high flow events or when a large portion of gravel is removed. This type of analysis would help engineers in the design of flood protection dykes and would also be particularly useful in the study of fish habitat since the development of bars, island and secondary channels due to gravel deposition can affect the quality of aquatic habitat. This is because: 1) the reworked stable gravels provide spawning substrate for chum and pink salmon; 2) back channels and slackwater zones provide rearing habitat for several salmonid species, sturgeon and other fishes; and 3) the multiple channels provide extensive bankline where hiding zones and drop-in food sources occur (Church, 1999). Therefore, any man made changes to the channel or gravel removal may disrupt this cycle of habitat renewal since the amount of “new” gravel entering the reach each year is small. R2DM can be used to estimate safe levels of gravel removal to maintain flood safety whilst preserving aquatic habitat. 86 6.3.5 Modeling of Suspended Load R2DM is currently limited to bed-load dominated channels. To further improve the model to include suspended load would require solving the convection diffusion equation for suspended sediment (Equation 41) with a mechanism to exchange sediment between the bed load and suspended load. ∂c ∂c ∂c ∂ ∂c ∂ ∂c + u + v − Ex − E y = S ∂t ∂x ∂y ∂x ∂x ∂y ∂x Equation 41. Suspended load sediment equation. where c=concentration, u and v are the velocity components in x and y, Ex and Ey = diffusion co-efficients, S=source term. To model the suspended load is not a simple task since the equation is far more complex than the bed-load transport continuity equation. Equation (41) may have to be discretized and solved using a different method to the bed-load equation and so this may also mean increased computational time. 87 7.0 References Ashida, K. and Michiue, M. (1971). An investigation of river bed degradation downstream of a dam, Proc. 14th Congress of the IAHR, pp. 245-255. Church, M. (1999), Sedimentation and Flood Hazard in the Gravel Reach of Fraser River: Progress Report, University of British Columbia, available on line at http://www.geog.ubc.ca/fraserriver/reports/report1999.pdf. Cui, Y., Parker, G., (2003). Sediment pulses in mountain rivers: 1. Experiments. Water Resources, Water Resources Research, 39(9): 1239-1251. Cui, Y., Parker, G., Pizzuto, J., Lisle, T., E., (2003). Sediment pulses in mountain rivers: 2. Comparison between experiments and numerical predictions, Water Resources, Water Resources Research, 39(9): 1240. Cui, Yantao (2007). The Unified Gravel-Sand (TUGS) Model: Simulating Sediment Transport and Gravel/Sand Size Distributions in Gravel-Bedded Rivers, Water Resources Research, 43. Einstein, H. A., 1950, The Bed-load Function for Sediment Transportation in Open Channel Flows, Technical Bulletin 1026, U.S. Dept. of the Army, Soil Conservation Service. Engelund, F., and Hansen, E. (1967). “A Monograph on Sediment Transport in Alluvial Streams,” Technical Forlag, Copenhagen, Denmark Ferguson, R. and Church, M. (2009). A critical perspective on 1D modeling of river processes: gravel load and aggradation in lower Fraser River, Water Resources Research (in review). 88 Hey, R. D., Gravel-bed Rivers: Form and Processes, Gravel-Bed Rivers, eds Hey, R. D., Bathurst, J. C., Thorne, C. R., Wiley, 1982. Islam, A. K. M. S. (2008). PhD Thesis: Gravel transport and morphological modeling for the lower fraser river, BC, University of British Columbia. 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. Kwan, S. (2009). R2DM User Manual. Neill, C. R. and Hey, R. D., Gravel-bed Rivers: Engineering Problems, Gravel-Bed Rivers, eds Hey, R. D., Bathurst, J. C., Thorne, C. R., Wiley, 1982. Paola, C., Parker, G., Sinha, S. K., Southard, J. B., Wilcock, P. R. (1992). Downstream fining by selective deposition in a laboratory flume, Science, 258: 1757-1760. Parker, G., Klingeman, P.C., McLean, D.G., Bedload and size distribution in paved gravel-bed streams. American Society of Civil Engineers, Proceedings, Journal of the Hydraulics Division 108 HY4 (1982), pp. 544–571 Parker, G., Transport of Gravel and Sediment Mixtures in Sedimentation Engineering: Process, Measurements, Modeling and Practice, ed. Marcelo H. Garcia, 2008. Patankar, Suhar, Numerical Heat Transfer and Fluid Flow, Taylor Francis, 1980. Rozovskii, J. L. (1957). Flow of water in bends of open channel. Academic Science of the Ukraine SSR, Kiev. 89 Seal, R., C. Paola, G. Parker, and B. Mullenbach (1995), Laboratory experiments on downstream fining of gravel, narrow channel runs 1 through 3: supplemental methods and data, External Memorandum M-239, St. Anthony Falls Laboratory, Univ. of Minnesota. Seal, R., Paola, C., Parker, G., Southard, J. B., Wilcock, P. R. (1997). Experiments on Downstream Fining of Gravel: Part 1 Narrow Channel Runs, Journal of Hydraulic Engineering, 126(3): 198-208. Soni, J. P., Garde, R. J., and Ranga Raju, K. G. (1980). Aggradation in streams due to overloading, Journal of Hydraulic Engineering, ASCE, 106: 117-132. Southard, John (2006). An Introduction to Fluid Motions, Sediment Transport, and Current-generated Sedimentary Structures. MIT Open Course Ware, Massachussets Institute of Technology, available on line at http://ocw.mit.edu/OcwWeb/Earth-Atmospheric--and-Planetary-Sciences/12-090Fall-2006/CourseHome/index.htm. Steffler, P. M. & Blackburn, J. (2002). River2D: Two-dimensional depth averaged model of river hydrodynamics and fish habitat. Introduction to depth averaged modeling and user’s manual. Edmonton, University of Alberta. Struiksma, N. (1985). Prediction of 2-D bed topography in rivers, Journal of Hydraulic Engineering, 111(8): 1169-1182. Tsujimoto, T. (1987). Non-uniform bed load transport and equilibrium bed profile. Proceedings of the XXII IAHR Congress, Lausanne, Switzerland, Fluvial Hydraulics, Vol 1, 177-182. Toro-Escobar, C., Paola, C., Parker, G., Wilcock, P. R., Southard, J. B., (2000). Experiments on Downstream Fining of Gravel. II: Wide and Sandy Runs, Journal of Hydraulic Engineering, 90 Van Rijn, L.C. (1984). Sediment Transport, Part I: Bed load transport, Journal of Hydraulic Engineering, ASCE, no 10. Vasquez, J. A. (2005). PhD Thesis: Two Dimensional Finite Element River Morphology Model. University of British Columbia. Vasquez, J. A., Millar, R. G. & Steffler, P. M. (2005a). River2D Morphology, Part I: Straight Alluvial Channels, 17th Canadian Hydrotechnical Conference, Edmonton, 17-19 August 2005. Vasquez, J. A., Millar, R. G. & Steffler, P. M. (2005b). River2D Morphology, Part II: Curved Alluvial Channels, 17th Canadian Hydrotechnical Conference, Edmonton, 17-19 August 2005. Vasquez, J. A. (2009), Personal communication. Wilcock, P. R. and Crowe, J. C. (2003). Surface-based transport model for mixed-sized sediment, Journal of Hydraulic Engineering, 129(2), 120-128. Wilcock, P. R., Kenworthy, S. T., Crowe, J. C. (2001). Experimental study of the transport of mixed sand and gravel, Water Resources Research, 37(12): 3349-3358. Wong, M. and Parker, G. (2006). Re-Analysis and Correction Of Bedload Relation Of Meyer-Peter and Müller using their own Database. Journal of Hydraulic Engineering, 132(11): 1159-1168. Wu, W. (2001) CCHE2D sediment transport model, Technical Report No. NCCHE-TR2001-3, National Center for Computational Hydroscience and Engineering, The University of Mississippi. 91 Appendix: Running R2DM This section gives details of how to run a mixed size sediment transport simulation using R2DM. Step 1: Run a Steady State Hydrodynamic Simulation 1. Calibrate the hydrodynamic model for a given flow scenario by adjusting the roughness values (bed values for open water, bed) until model results agree with observed data. 2. Run a model using the steady solver until it reaches steady state. 3. Save the CDG file. Step 2: Specify the Boundary Conditions for a Morphological Simulation To run a morphological model the boundary conditions must be specified. These are: 1. the properties of the sediment (e.g porosity, angle of repose) 2. the sediment transport equation, 3. boundary conditions and 4. factors which affect the mathematical stability (e.g. up-winding factors, secondary flow etc). These properties can be defined under Options->Sediment Options: 92 Figure 81. The “Sediment Options” dialogue box Step 3: Select the Wilcock and Crowe Sediment Transport Equation After selecting Wilcock and Crowe, press the W/C settings button to set additional initial conditions and parameters. These are: 1. Surface layer distribution of computational domain 2. Substrate distribution 3. Grain size distribution of sediment feed 4. Initial surface layer depth 5. Transport function weighting factor (0-100) 6. ks factor (1-10) These can be set by selecting W/C settings button on the “Sediment Options” dialogue box. 93 Figure 82. The “Wilcock and Crowe Settings: Set Grain Size Distribution” dialogue box. Step 4: Run Morphology 1. Select the “Run Morphology” option (it replaces the Run Transient option) under the “Hydrodynamics” menu item. 2. Select the type and frequency of transient output for the simulation (if any is desired). Step 5: Viewing Results with R2DM Vector Data The information related to bed morphology can be displayed using the same output options as River2D. In the “Vector Plot” dialogue box the “Sediment Transport” option replaces the “Discharge Intensity”: 94 Figure 83. The “Vector Plot” dialogue box Contour Plots In the “Contour Plot” dialogue box, the following options have been added: • Curvature • qsx (x-component of sediment flux) • qsy (y-component of sediment flux) • qs (sediment flux) • Bed Change Figure 84. The “Colour/Contour” dialogue box If the Wilcock and Crowe sediment transport equation is used, then the following parameters can be displayed: 95 • D50 (Substrate) • D50 (Surface) • D50 (Bed load) • D90 (Substrate) • D90 (Surface) • D90 (Bed load) • Sand Fraction If the D90 is selected the percent finer distribution will be displayed when the mouse is moved over the flow field. Output Files CSV files containing the bed morphology, ks, grain sizes (D10, D30, D50, D70, D90), gravel grain sizes (GD5, GD16, GD30, GD50, GD84, GD90), depth and sand fraction at various time intervals can be obtained by selecting the output options button on the “Run Morphology” dialogue box (see Figure 59). Figure 85. Where to select CSV files in the “Run Morphology” dialogue box. 96
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A two dimensional hydrodynamic river morphology and...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A two dimensional hydrodynamic river morphology and gravel transport model Kwan, Stephen 2009-12-31
pdf
Page Metadata
Item Metadata
Title | A two dimensional hydrodynamic river morphology and gravel transport model |
Creator |
Kwan, Stephen |
Publisher | University of British Columbia |
Date | 2009 |
Date Issued | 2009-07-31T20:00:24Z |
Description | A depth-averaged two-dimensional hydrodynamic-morphological and gravel transport model is presented. Based on River2D by Prof Peter Steffler et al. (2002) of the University of Alberta, River2D-Morphology (R2DM) is capable of simulating flow hydraulics, sediment transport for uni-size and mixed size sediment, and morphological changes of a river over time. Changes in bed elevation are calculated by solving Exner’s equation for sediment mass conservation with the mixed size sediment transport function of Wilcock and Crowe (2003). R2DM is capable of exploring the dynamics of grain size distributions, including fractions of sand in sediment deposits and on the bed surface of rivers and streams. Results have been verified with experimental data for aggradation and downstream fining (Paolo, 1992; Seal, 1992; Toro-Escobar, 2000); for degradation (Ashida and Michiue, 1971); and qualitatively with two dimensional flow (Tsujimoto, 1987). Application of R2DM to actual rivers has shown that it is ready for general release as a tool to assist engineers with river restoration projects, the design of hydraulic structures and with fish habitat modeling. |
Extent | 3908928 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Collection |
Electronic Theses and Dissertations (ETDs) 2008+ |
Date Available | 2009-07-31 |
Provider | Vancouver : University of British Columbia Library |
DOI | 10.14288/1.0063132 |
URI | http://hdl.handle.net/2429/11562 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2009-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- ubc_2009_fall_kwan_stephen.pdf [ 3.73MB ]
- Metadata
- JSON: 1.0063132.json
- JSON-LD: 1.0063132+ld.json
- RDF/XML (Pretty): 1.0063132.xml
- RDF/JSON: 1.0063132+rdf.json
- Turtle: 1.0063132+rdf-turtle.txt
- N-Triples: 1.0063132+rdf-ntriples.txt
- Original Record: 1.0063132 +original-record.json
- Full Text
- 1.0063132.txt
- Citation
- 1.0063132.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
Canada | 16 | 10 |
United States | 13 | 62 |
China | 10 | 0 |
Peru | 8 | 4 |
Mexico | 6 | 3 |
Brazil | 5 | 4 |
Germany | 5 | 38 |
India | 4 | 3 |
France | 4 | 3 |
Ethiopia | 2 | 0 |
Ecuador | 2 | 2 |
Finland | 1 | 0 |
Spain | 1 | 1 |
City | Views | Downloads |
---|---|---|
Unknown | 34 | 63 |
Ashburn | 7 | 0 |
Coyoacan | 6 | 3 |
Shenzhen | 4 | 0 |
Québec | 4 | 0 |
Beijing | 4 | 0 |
Vancouver | 2 | 0 |
Cuenca | 2 | 2 |
San Jose | 2 | 0 |
Montreal | 2 | 0 |
Wellington | 1 | 1 |
Coogee | 1 | 0 |
Olympia | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
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}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0063132/manifest