UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Application of a two-dimensional hydrodynamic model to the Fraser River Gravel Reach Yusuf, Faizal 2001

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2001-0559.pdf [ 27.94MB ]
Metadata
JSON: 831-1.0063721.json
JSON-LD: 831-1.0063721-ld.json
RDF/XML (Pretty): 831-1.0063721-rdf.xml
RDF/JSON: 831-1.0063721-rdf.json
Turtle: 831-1.0063721-turtle.txt
N-Triples: 831-1.0063721-rdf-ntriples.txt
Original Record: 831-1.0063721-source.json
Full Text
831-1.0063721-fulltext.txt
Citation
831-1.0063721.ris

Full Text

Application of a two-dimensional hydrodynamic model to the Fraser River Gravel Reach by Faizal Yusuf B .A .Sc , University of Waterloo, 1999 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE D E G R E E OF MASTER OF APPLIED SCIENCE in THE FACULTY OF G R A D U A T E STUDIES DEPARTMENT OF CIVIL ENGINEERING We accept this thesis as conforming to the required standard The University of British Columbia August, 2001 © Faizal Yusuf, 2001 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Cc*i / j5<4j/V> CCrth The University of British Columbia Vancouver, Canada DE-6 (2/88) Application of a 2D Hvdrodynamic Model to the Fraser River Gravel Reach ABSTRACT This thesis investigates the potential analytical uses of two-dimensional (2D) models for large braided river systems such as the Fraser River Gravel Reach, which extends from Laidlaw to Mission, B.C., and examines the limitations and advantages of a 2D modeling approach. The 2D depth-averaged model, River2D, was applied to two sections of the gravel reach: a 4.5 km section of the Fraser River at the Agassiz-Rosedale Bridge and an 8.5 km section near the Harrison River confluence and Minto Island. The secondary objective of this work was to characterize the hydraulics in the modeled reaches. River2D was used to examine various gravel extraction scenarios, investigate bank erosion issues, estimate superelevation of the water surface around bends and to determine local depths and velocities for use in habitat delineation and mapping studies. The two sections of the river studied in this thesis have been the focus of several engineering studies by local consultants over the past few years. Some of the recommendations made in those studies to alleviate local flooding risks and erosion concerns were investigated with River2D. The key steps in a typical 2D modeling study were examined in detail beginning with the development of a digital elevation model (DEM) from topographic and bathymetric survey data, which had been previously collected, followed by the development of a River2D model file and model calibration. The most important step in 2D modeling was found to be obtaining an accurate representation of the bed topography. Channel roughness, which is represented by the roughness length, ks, in River2D, was estimated with the surface dso of the bed material. In this study, River2D has been shown to be technically sound in terms of its hydrodynamic formulation. Further development of this program is recommended to include a morphodynamic module specifically designed to address issues in the gravel reach. ii Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach T A B L E OF C O N T E N T S Abstract » Table of Contents iii List of Tables vi List of Figures vii List of Symbols xi Acknowledgements xiii 1. Introduction 1 1.1 Overview 1 1.2 Objectives 5 1.3 Scope 5 2. Literature Review 7 2.1 Hydrodynamic Models 7 2.1.1 General 7 2.1.2 Governing Equations 7 2.1.3 Numerical Solution 15 2.1.4 Hydrodynamic Model Examples and Applications 19 2.1.5 River2D 23 2.2 Flow Resistance 27 2.3 Turbulence Models 31 2.4 Acoustic Doppler Current Profilers (ADCP) 33 3. Fraser River Gravel Reach 36 3.1 Setting 36 3.2 Hydrology 38 3.3 Current Issues 39 4. Survey Data 42 iii Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 4.1 Floodplain Topography 42 4.2 River Bathymetry 43 4.3 A D C P Survey 43 5. River2D Sensitivity Analysis 45 5.1 Roughness Length 45 5.2 Eddy Viscosity Coefficient 51 5.3 Computational Node Spacing 51 6. Model Development and Calibration: Agassiz-Rosedale Reach 57 6.1 Agassiz-Rosedale Reach Overview 57 6.2 DEM Generation 60 6.2.1 General 60 6.2.2 Agassiz-Rosedale Reach DEM 62 6.3 River2D Model Development 67 6.4 Calibration 70 6.4.1 1999 Freshet 71 6.4.2 A D C P Data :....72 6.4.3 Aerial Photograph 77 6.5 Summary 79 7. Applications: Agassiz-Rosedale Reach 80 7.1 Superelevation at Agassiz-Rosedale Bridge 80 7.2 Mean Annual Flood 81 7.3 Estimation of Grain Size Spatial Distribution 87 7.4 Design Flood 89 7.5 Gravel Removal Scenarios 93 7.5.1 Localized Removal 94 7.5.2 Large-Scale Removal: Powerline Island 99 7.6 Habitat information at Big Bar 104 8. Model Development and Calibration: Harrison River Confluence/Minto Island Reach 110 iv Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 8.1 Harrison River Confluence/Minto Island Reach Overview 110 8.2 DEM Generation 113 8.3 River2D Model Development 114 8.4 Calibration 115 8.4.1 1999 Freshet 116 8.4.2 A D C P Data 118 8.5 Summary 124 9. Applications: Harrison River Confluence/Minto Island Reach 125 9.1 Harrison River Confluence Energy Loss 125 9.2 Design Flood 127 9.3 Harrison Bar Pilot Channel 130 9.4 Dredging of Harrison Bar 137 9.5 Island 22 Bank Erosion 142 10. Discussion and Conclusions 151 11. Recommendations 159 12. References 162 v Application of a 2D Hvdrodynamic Model to the Fraser River Gravel Reach LIST OF T A B L E S Table 3.1: Hydraulic characteristics at the mean annual flood 38 Table 5.1: ks values used in the sensitivity analysis shown with corresponding Manning's n 46 Table 5.2: Node spacing sensitivity analysis 53 Table 6.1: 1999 freshet calibration results in the Agassiz-Rosedale reach 72 Table 6.2: Flow splits (%) in the Agassiz-Rosedale reach 76 Table 6.3: Comparison of water levels measured during the A D C P survey in the Agassiz-Rosedale reach 76 Table 7.1: Comparison of 1D hydraulic parameters at the Agassiz-Rosedale Bridge for the mean annual flood 82 Table 8.1: 1999 freshet calibration results in the Harrison River confluence/Minto Island reach 118 Table 8.2: Flow splits (%) in the Harrison River confluence/Minto Island reach 122 Table 8.3: Comparison of water levels measured during the A D C P survey in the Harrison River confluence/Minto Island reach 123 vi Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach LIST OF F IGURES Figure 1.1: Fraser River basin and gravel reach study area 4 Figure 2.1: Definition sketch for 2D depth averaged free surface flow 9 Figure 3.1: Fraser River Gravel Reach 37 Figure 5.1: Average velocity versus ks 46 Figure 5.2: Average depth versus ks 47 Figure 5.3: Velocity change from ks = 0.04 m to ks = 2.04 m 49 Figure 5.4: Depth change from ks = 0.04 m to ks = 2.04 m 49 Figure 5.5: Mean average velocity error versus node spacing 54 Figure 5.6: Mean average depth error versus node spacing 55 Figure 6.1: Agassiz-Rosedale reach 58 Figure 6.2: Fraser R.: Agassiz-Rosedale Bridge looking upstream at Q = 6700 m 3/s : 59 Figure 6.3: Fraser R.: Agassiz-Rosedale Bridge looking upstream at Q = 650 m 3/s 59 Figure 6.4: Agassiz-Rosedale reach contour plot 65 Figure 6.5: Surface plot of the upper half of Agassiz-Rosedale reach (looking downstream) 66 Figure 6.6: Surface plot of the lower half of Agassiz-Rosedale reach (looking upstream) 67 Figure 6.7: A D C P transect locations in the Agassiz-Rosedale reach 73 Figure 6.8: Simulated depth versus A D C P depth for the Agassiz-Rosedale reach.. 74 Figure 6.9: Simulated velocity versus A D C P velocity for the Agassiz-Rosedale reach 74 Figure 6.10: Depth contours superimposed over aerial photo of the Agassiz-Rosedale reach (March 7, 2001; Q=500 m3/s) 78 Figure 6.11: Velocity contours superimposed over aerial photo of the Agassiz-Rosedale reach (March 7, 2001; Q=500 m3/s) 78 Figure 7.1: Bed elevation contours for the mean annual flood at Agassiz 83 Figure 7.2: Depth contours for the mean annual flood at Agassiz 83 Figure 7.3: Water surface elevation contours for the mean annual flood at Agassiz 84 vii Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Figure 7.4: Velocity contours for the mean annual flood at Agassiz 84 Figure 7.5: Discharge per unit width contours for the mean annual flood at Agassiz 85 Figure 7.6: Froude number contours for the mean annual flood at Agassiz 85 Figure 7.7: Shear stress contours for the mean annual flood at Agassiz 86 Figure 7.8: Discharge per unit width vectors for the mean annual flood at Agassiz . 86 Figure 7.9: Design flood profile for the Agassiz-Rosedale reach 91 Figure 7.10: Local gravel removal on Powerline Island 95 Figure 7.11: Localized gravel removal on Powerline Island: Pre-removal depth contours 96 Figure 7.12: Localized gravel removal on Powerline Island: Post-removal depth contours 96 Figure 7.13: Localized gravel removal on Powerline Island: Pre-removal velocity contours 97 Figure 7.14: Localized gravel removal on Powerline Island: Post-removal velocity contours 97 Figure 7.15: Localized gravel removal on Powerline Island: Pre-removal discharge vectors 98 Figure 7.16: Localized gravel removal on Powerline Island: Post-removal discharge vectors 98 Figure 7.17: Location of Ferry Island and Island "B" 99 Figure 7.18: Proposed gravel removal area to reduce downstream bank erosion .100 Figure 7.19: Powerline Island gravel removal: Contour plot of the change in velocity due to extraction 102 Figure 7.20: Powerline Island gravel removal: Pre-removal discharge vectors 103 Figure 7.21: Powerline Island gravel removal: Post-removal discharge vectors.... 103 Figure 7.22: Location of Big Bar 105 Figure 7.23: Big Bar habitat information: Depth contours for 700 m 3/s 106 Figure 7.24: Big Bar habitat information: Velocity contours for 700 m 3/s 106 Figure 7.25;- Big Bar habitat information: Depth contours for 1,500 m 3/s 107 Figure 7.26: Big Bar habitat information: Velocity contours for 1,500 m 3/s 107 viii Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Figure 7.27: Big Bar habitat information: Depth contours for 3,000 m 3/s 108 Figure 7.28: Big Bar habitat information: Velocity contours for 3,000 m 3/s 108 Figure 7.29: Big Bar habitat information: Depth contours for 6,000 m 3/s 109 Figure 7.30: Big Bar habitat information: Velocity contours for 6,000 m 3/s 109 Figure 8.1: Harrison River confluence/Minto Island reach 111 Figure 8.2: Photo taken at the downstream end of Harrison River, looking towards the Fraser River 112 Figure 8.3: Harrison River confluence/Minto Island reach contour plot 114 Figure 8.4: Harrison River confluence/Minto Island reach roughness regions 117 Figure 8.5: A D C P transect locations in the Harrison River confluence/Minto Island reach 119 Figure 8.6: Simulated depth versus A D C P depth for the Harrison River confluence/Minto Island reach 120 Figure 8.7: Simulated velocity versus A D C P velocity for the Harrison River confluence/Minto Island reach 121 Figure 9.1: Water surface elevation contours in the Harrison River confluence area for the design flood 126 Figure 9.2: Design flood profile for the Harrison River confluence/Minto Island reach 128 Figure 9.3: Location of Harrison Bar pilot channel 131 Figure 9.4: Design flood depth contours for the Harrison Bar pilot channel: Pre- and post-excavation 133 Figure 9.5: Design flood velocity contours for the Harrison Bar pilot channel: Pre-and post-excavation 134 Figure 9.6: Design flood profile along the main channel for the Harrison Bar pilot channel 135 Figure 9.7: Location of Harrison Bar dredge area 138 Figure 9.8: Design flood depth contours for the Harrison Bar dredge: Pre- and post-excavation 139 Figure 9.9: Design flood velocity contours for the Harrison Bar dredge: Pre- and post-excavation 140 ix Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Figure 9.10: Design flood profile along the main channel for the Harrison Bar dredge 141 Figure 9.11: Island 22 bank erosion 143 Figure 9.12: 1999 freshet discharge vectors for the Queen's Bar gravel deposit proposed excavation 146 Figure 9.13: 1999 freshet velocity contours for the Queen's Bar gravel removal: Pre-and post-excavation 147 Figure 9.14: 1999 freshet velocity contours for the Queen's Bar gravel removal: Change in velocity due to extraction 148 Figure 9.15: 1999 freshet depth contours for the Queen's Bar gravel removal: Change in depth due to extraction 148 x Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach ^  LIST OF S Y M B O L S x, y coordinates in the horizontal plane [m] z elevation with respect to a defined datum [m] H depth of flow [m] U velocity [m/s] u, v depth averaged velocity components [m/s] w vertical velocity component [m/s] Q discharge [m3/s] q discharge per unit width [m2/s] qx, qy discharge per unit width components [m2/s] R hydraulic radius [m] g acceleration due to gravity [m/s2] p density of water [kg/m3] S o bed slope [m/m] Sox, S o y bed slope components [m/m] Sfx, Sfy friction slope components [m/m] to mean bed shear stress [N/m2] Tox, T0y bed shear stress components [N/m2] TXX, Txy, Tyx, Tyy transverse shear stress components [N/m2] T* dimensionless bed shear stress TC* critical dimensionless bed shear stress U* shear velocity [m/s] C Chezy roughness coefficient [m1 / 2/s] C s dimensionless Chezy roughness coefficient ks effective roughness length [m] n Manning's roughness coefficient [s/m1 / 3] f Darcy-Weisbach friction factor [dimensionless ] Vxx\ Vxy', vyy', vyy' depth averaged turbulent exchange coefficient components [m2/s] xi Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach v' eddy viscosity coefficient [m2/s] s user-defined River2D coefficient that controls v' [dimensionless] c/50 median grain size diameter [m] K constant m constant T transmissivity [m2/s] S storativity [dimensionless]; and specific gravity of sediment [dimensionless] z b ground surface elevation in groundwater flow equation [m] k kinetic energy of the turbulent motion per unit mass [m2/s2] / length scale [m] N total number of survey points AH superelevation [m] r radius of curvature [m] W width of river [m] a velocity coefficient [dimensionless] c wave celerity [m/s] xii Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach A C K N O W L E D G E M E N T S I would like to thank my research supervisor, Dr. Robert Millar, for his valuable insight and guidance throughout this study. His patience and helpfulness are greatly appreciated. I am very grateful to Colin Rennie for providing the A D C P data including the velocity profile information which he calculated and also for his contributions to this work through discussions. Darren Ham was extremely helpful in providing the bathymetric and topographic survey data in the requested formats while Ron Henry at the Ministry of Environment, Lands and Parks supplied several reports, maps and useful information. A special thanks to Dr. Noboru Yonemitsu for his assistance in tracking down River2D modeling papers. The Natural Sciences and Engineering Research Council of Canada provided funding which helped make this work possible. xiii Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 1. INTRODUCTION 1.1 Overview The simulation of flow in streams and rivers provides useful information for various purposes. For example, water surface profiles can be generated for a range of discharges which provides the basic design criteria for flood protection measures such as dyking systems. Similarly, the adequacy of existing flood control structures can also be assessed. Open channel flow simulation can also aid in the evaluation of fish habitat and river diversion works and the analysis of local flow hydraulics including flow around a bend or around an obstruction such as a bridge pier. The investigation of sediment transport and water quality issues can often be greatly assisted by computer simulations as well. Numerous computer models are available to simulate the flow in streams and rivers. These models can differ greatly in their complexity and numerical formulation therefore a particular model's applicability to a certain modeling scenario requires adequate knowledge of the model's inherent assumptions, data requirements, output and limitations. An important distinction in open channel flow is made between steady and unsteady flow. Flow is steady if the depth and velocity of flow do not change in magnitude and direction with respect to time. In natural rivers, steady flow conditions are rare; however, in some cases the temporal variations of the flow may be small and over short time intervals, the flow may be assumed to be steady. If a steady flow assumption can be made for a channel or reach, iterative analytical procedures exist to calculate water surface profiles. Although the calculations can be done by hand, for irregular channel cross-sections and when water surface elevations at more than a few locations are required, which is often the case, a water surface profile can be generated with a software program such as H E C - R A S . Most river engineering problems deal with unsteady or transient flows, where the flow depth and/or velocity change with time. For example, a steady flow model could not be used for a braided reach with multiple flow paths and a model capable 1 Application of a 2D Hvdrodvnamic Model to the Fraser River Gravel Reach of simulating unsteady flow conditions would be required. Unsteady flow is described by a set of differential equations that consist of an equation of continuity and equation(s) of conservation of momentum. Since these equations cannot be solved analytically, hydrodynamic models which employ numerical solution techniques are used. Hydrodynamic models differ in the numerical method used (e.g. finite difference, finite element and finite volume), the solution approach of the numerical method (implicit or explicit), the computational grid employed (rectilinear, curvilinear, or mesh) and the number of spatial dimensions. The distinction between a fixed and mobile bed is also of importance since many morphological and sediment transport issues require the simulation of mobile bed conditions. One-dimensional (1D) hydrodynamic models are frequently used in river engineering although technological advances in the areas of personal computers and computational software have enabled two-dimensional (2D) models to become more common. Two-dimensional hydrodynamic models allow local velocity and depth distributions to be calculated which is a significant step from the average velocity and water surface elevations at each cross-section typically solved for in 1D models. While the potential analytical uses of 2D models far outweigh those of 1D programs, so do the computational challenges which accompany them. The difficulty in modeling both supercritical and subcritical flows and dealing with wetting and drying areas are two of the potential problems that must be addressed. Several different types of models have been applied to the lower Fraser River for various purposes. This section of the Fraser River represents a diverse set of issues that require careful management to maintain the river's ecological quality while keeping the interests of various stakeholders in mind. The ecological significance of the Fraser River Gravel Reach from Laidlaw to Mission (Figure 1.1) is evident by at least 28 species of fish in addition to several types of mammals, amphibians and birds that are found there. This reach is an important migration 2 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach route for several species of salmon and steelhead that spawn in up-river tributaries. In addition, several species offish use the complex habitat offered by the gravel reach for spawning and rearing. The gravel reach also serves as a route for navigation, a popular location for recreation and fishing and the traditional territory of several First Nations communities. There are significant erosion and flooding risks to both aboriginal and non-aboriginal communities along the river as well. Historically, gravel has been extracted from this part of the Fraser River for commercial use in construction projects. The management and decision-making processes regarding these issues has been aided by the use of hydrodynamic models. The focus of the various modeling studies on the lower Fraser River has been as diverse as the issues which the river presents. ONE-D, which was developed by Environment Canada, is employed operationally on the lower Fraser River for the simulation of flow in tidal-affected regions where flow reversals temporarily occur in several reaches. The unsteady fixed-bed model, MIKE 11 has also been applied to the lower Fraser River as a general predictive tool for assessing hydrodynamic conditions for vessel navigation and safety (Scott et al., 1999). The design flood profiles in the gravel reach were recently updated using MIKE 11 as well (LIMA, 2000). A morphodynamic module in MIKE 11 was used to investigate the effects of a proposed excavation in the vicinity of the Harrison River mouth to reduce the flood risk in that area of the gravel reach which is characterized by complex hydraulics (Water Management Consultants, 2001). A 2D model, TELEMAC-2D, has been applied in the tidally influenced lower 40 km of the Fraser River, below the gravel reach, to provide information on current patterns in certain key reaches (Scott et al., 1999). 3 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach The focus of this study is the application of a 2D hydrodynamic model to the Fraser River Gravel Reach. Prior to this work, the gravel reach has not been studied with such a computational tool and therefore the potential analytical uses and limitations of this modeling approach have not been fully investigated. This study will hopefully provide the initial steps to guide local engineering practice in adopting 2D modeling techniques in the analysis of river hydrodynamics and fish habitat. The 2D depth-averaged model, River2D, which was developed at the University of Alberta, was selected for use in this study. River2D is public domain software and was therefore readily available for this research project. Furthermore, the model has been proven to be technically sound through its application to several modeling scenarios, as described in Section 2.1.5. 1.2 Objectives The primary objectives of this study are to investigate the potential analytical uses of 2D models for large braided river systems such as the Fraser River Gravel Reach and to examine the limitations and advantages of a 2D modeling approach. The secondary objective of this work is to characterize the hydraulics in areas of interest. Specifically, River2D will be used to examine various gravel extraction scenarios, investigate bank erosion issues, estimate superelevation of the water surface around bends and to determine local depths and velocities for use in habitat delineation and mapping studies. 1.3 Scope Two sections of the gravel reach were selected for the application of River2D. The first was a 4.5 km section of the Fraser River at the Agassiz-Rosedale Bridge. The upstream end of this reach consists of two distinct channels which convey approximately equal amounts of flow converging into one main channel that flows under the Agassiz-Rosedale Bridge. Mid-channel bars with vegetated islands divide the river, which is about 500 m wide at the bridge, creating side channels. 5 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach The second section of river that was modeled was an 8.5 km reach near the Harrison River confluence and Minto Island. Significant features in this reach include very complex flow hydraulics at the mouth of the Harrison River, where the main channel of the Fraser River is abruptly redirected through a 90° bend as it flows into the outcropping bedrock of Harrison Knob which results in strong secondary currents and large energy losses. Minto Island divides this reach into two distinct channels. This reach presents many more computational obstacles than the Agassiz-Rosedale reach and has been the focus of several engineering studies by local consultants over the past few years. Some of the recommendations made in those studies to alleviate local flooding risks and erosion concerns were investigated with the 2D model. It should be made clear that this thesis is not intended to provide a thorough 2D hydraulic assessment of the Fraser River Gravel Reach. The primary focus of this study is on the investigation of 2D models as analytical tools to assess local hydraulics. In some cases, the results of the 2D simulations and analysis are not conclusive with respect to the Fraser River flow hydraulics and should be considered preliminary. In those instances further work would be required to investigate and potentially validate the results. A literature review is included in this thesis to provide relevant background information on hydrodynamic models, open channel flow and the use of acoustic Doppler current profilers to measure velocity distributions that can aid in 2D model calibration. The Fraser River Gravel Reach is described in terms of its main morphological, hydrological and hydraulic features. A description of the survey data that was used to create the digital elevation models is also included. The results of a sensitivity analysis on River2D parameters that are controlled by the user are presented and discussed. The 2D model development and calibration for both the Agassiz-Rosedale and Harrison River confluence/Minto Island reaches are described in detail and several applications of River2D in both reaches are presented. 6 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 2. LITERATURE REVIEW 2.1 Hydrodynamic Models 2.1.1 General This section on hydrodynamic models is intended to provide general background information on the theoretical considerations and assumptions made in these models as well as to describe some of the various types of models used and their applications. In Section 2.1.2, the governing equations of continuity and conservation of momentum are presented. The assumptions made in the development of these equations and their applicability to simulating flow in different situations are discussed. A general discussion of the different numerical solution techniques used in these models is included in Section 2.1.3. Since a detailed review of the numerical solution schemes (e.g. finite difference and finite element) is not attempted, relevant literature is cited. Section 2.1.4 provides information on several different hydrodynamic models that are currently used in practice and describes some of their applications. The 2D model used in this study, River2D, is described in Section 2.1.5. 2.1.2 Governing Equations The governing equations for open channel flow consist of an equation of continuity or conservation of mass and one or more equations of conservation of momentum. The number of .spatial dimensions considered determines the number of equations of the conservation of momentum. One of the main assumptions made in the derivation of these equations is that the vertical acceleration of the water particles are negligible in comparison to the total acceleration so the streamlines are essentially parallel. This allows hydrostatic pressure distributions to be assumed. Also, the effects of boundary friction can be accounted for with resistive forces described by empirical relationships such as the Manning or Darcy-Weisbach equations. Bed slopes in the principal flow direction are assumed to be small, such that sine s tanG (slope < 10%), where 6 is the angle the bed makes with the 7 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach horizontal plane in the principal flow direction (Ghanem et al., 1995b). The density of water is also assumed to be homogeneous. Since this study focuses on 2D hydrodynamic models, the governing equations for two-dimensional free surface flow will be presented. These equations can be obtained by integrating the three-dimensional Reynolds equations over the depth of flow (e.g. Weiyan, 1992; Stefflerand Jin, 1993). Alternatively, a control volume approach can be used to derive these equations, which are collectively known as the St. Venant equations, as described in Steffler (2000). Coriolis forces, which are body forces due to the rotation of the earth, are important in atmospheric, ocean and lake dynamics but insignificant in most rivers (Millar and Barua, 1999). Therefore, these forces are not included in the form of the St. Venant equations presented and shear stresses induced by the wind on the free surface are assumed to be negligible also. In the control volume approach the basic principles of conservation of mass and momentum are applied to a prismatic vertical water column (Figure 2.1) that is bounded by the bed from the bottom and the free surface at the top. This approach yields one equation of continuity or the conservation of mass and two equations for the conservation of momentum in the x and y directions, shown on the following page. 8 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach H Figure 2.1: Definition sketch for 2D depth averaged free surface flow (after Ghanem et al., 1995a) Conservation of mass: dt dx. dy [2.1] Conservation of x-direction momentum: d g ' d ^ .+ (uqx) + -(vqx) + ±\ —{H*) = gH(S0X -Sfx) + - ^T(HT„) ot ox dy 2\dx ) p\dx 1 d + • P dy • J [2.2] Conservation of y-direction momentum: diy d 1 Kuh y , dt dx dy ^(H2)) = gH(Soy - Sfi) + -(^r(HTyx) dy J P\dx + • P d dy [2.3] 9 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach In the above equations, x and y represent the cartesian coordinates in the horizontal plane while t refers to time. The depth of flow (m) is denoted by H while the components of discharge per unit width (m2/s) in the x and y directions are qx and <7y. The dependent variables H, qx and ay are functions of the independent variables x, y and t. The depth averaged velocity components (m/s) are u and v and are related to qx and qy through u = qx/H and v = q/H. The acceleration due to gravity (m/s2) is denoted by g and the density of water (kg/m3) is p. Sox and Soy represent the bed slope components in equations [2.1] to [2.3]. A bed resistance model can be used to determine the friction slope components, Sfx and Sfy, by relating the bed shear stress components, TOX and r o y , to the magnitude and direction of the depth-averaged velocity (Steffler, 2000). A two-dimensional form of the Chezy equation is used to solve for Sfx and Sfy, shown below. The assumption made is that frictional resistance formulae for steady 1D flow are applicable to unsteady 2D flow. V 2 2 U + V Sfx = = [2.4] * PgH Cs2gH Toy V-yJu2 + V2 Sfy=-^-= " ' [2.5] fy PgH Cs2gH The roughness parameter in equations [2.4] and [2.5] is the dimensionless Chezy coefficient, C s and can be related to the effective bed roughness height (m), ks, and the depth of flow through the following relation: Cs =5.751og 11— [2.6] 10 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach The two-dimensional form of the St. Venant equations also include four depth averaged transverse shear stresses (N/m2) that are caused by turbulent flow interactions. These are represented by rxx, zxy, ryx and ryy and can be determined through a simplified transverse shear model that relates the turbulence stresses to mean flow properties. An example of a transverse shear stress equation is shown below: c%. dv dy dx [2.7] In equation [2.7], v x y ' is a depth averaged turbulent exchange coefficient (m2/s) and its magnitude is a function of the structure of the turbulence. If an isotropic turbulence assumption is made, i.e. v = vxx'= vxy'= vyx'= vyy', then v'can be defined as the eddy viscosity coefficient. If the dominant turbulence generation mechanism is assumed to be bed shear, then v will depend on the depth of flow and bed shear stress which yields equation [2.8] which includes a coefficient, e. For rivers, reasonable values for s range from 0.2 to 1.0 (Steffler, 2000). H^u2+v2 v'=s— 2.8 C . The bed resistance and transverse shear models described in equations [2.4] to [2.8] are similar to the formulation in River2D. In order to obtain local solutions for H, qx and qy, values for the effective roughness length, /cs, would be required at each computational point. The specific data requirements to obtain a solution to the St. Venant equations presented here are discussed in Section 2.1.5. In many other hydrodynamic models, parameters such as Manning's n or the Darcy-Weisbach friction factor, f, are used to describe energy losses due to boundary roughness. More discussion on open channel flow roughness parameters is provided in Section 2.2 and a brief review of turbulent stresses in rivers is included in Section 2.3. 11 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach In the 1D form of the St. Venant equations the flow is approximated with uniform velocity at each cross-section. Further simplifications are that the free surface is taken to be a horizontal line across each section so the centrifugal effect due to channel curvature is ignored (Water Modelling Section, 1988). Since cross-stream velocities are ignored, only the continuity and x-direction conservation of momentum equations are required. In a three-dimensional (3D) model, the velocity varies in all three directions, therefore, a vertical velocity component, w, is introduced. To obtain a solution in this case, a fourth equation is required: the conservation of momentum in the z-direction. The continuity equation conserves mass from section to section and controls storage routing in the St. Venant equations. The conservation of momentum equations consist of a series of terms for the local and convective acceleration and pressure, gravity and friction forces, which are shown below. Equation [2.2] which is the x-direction conservation of momentum equation can be re-written ignoring the transverse shear stress terms as follows: dix d d . . ~t + -z(uvx) + -z-(vqI) ot ox qy + ox - gH(Sa sfx) = o [2.9] Local Convective Pressure Gravity Friction acceleration acceleration force force force term terms term term term [2.10] [2.11] [2.12] 12 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach When both acceleration terms and all three force terms are included in the conservation of momentum equation, then it is referred to as the dynamic wave equation [2.10] which simulates unsteady, non-uniform flow. If the acceleration terms are omitted, the resulting form is the diffusion wave equation [2.11]. If further simplifications are made and the pressure force term is ignored, the remaining two terms represent the kinematic wave equation [2.12], which simulates uniform flow conditions. Most river engineering problems deal with unsteady or transient flows, where the flow depth and/or velocity change with time. Flow is considered uniform if the depth of flow does not change with distance along a channel, however uniform flow only occurs in long, prismatic channels that have constant cross-sections and slopes. Since these conditions do not occur in most natural channels, the flow is usually non-uniform, resulting in an acceleration or deceleration of the flow. Given that unsteady, non-uniform flow conditions are the most likely to occur in natural channels, it would imply that the full dynamic wave equation would be required to simulate a flood wave. In practice, the simplified models, i.e. diffusion and kinematic wave, are frequently used due to analytical and numerical advantages. It should be noted that the local acceleration term containing 3ns often retained in the kinematic and diffusion wave models to account for unsteady flow conditions. The applicability of these equations to certain modeling situations requires an examination of the acceleration and force terms shown in equation [2.9]. Henderson (1966) presented a scale analysis where he investigated the magnitudes of different slope terms based on typical values taken from an actual river in a steep alluvial region. The values he presented were based on a very fast rising flood where the acceleration terms would be relatively large. The slope terms he examined were S 0 , — , — — and and their magnitudes in m/km were approximately 4.92, 0.095, dx g dx g dt 0.035 and 0.0095, respectively. These four slope terms correspond to the gravity 13 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach force, pressure force, convective acceleration and local acceleration terms in equation [2.9]. This scale analysis indicates that the three latter slope terms are very small compared to the bed slope, therefore the gravity force term would be the most significant in steep rivers (S 0 > «0.002). In this case, the kinematic wave equation, which accounts for only pure translation and no attenuation of the flood wave, would probably suffice in approximating flow conditions. For bed slopes flatter than about 0.002, the pressure force terms becomes significant and may approach the same order of magnitude as the gravity force term. Under these conditions, the diffusion wave equation would be required to more accurately model the flow. The difference between the diffusion wave and dynamic wave approximation is the exclusion of the convective acceleration term, which is sensitive to velocity changes in the river. Therefore, where velocity changes would be significant, the full dynamic wave approximation should be used. Henderson (1963) concluded that for steep slopes the kinematic wave equation adequately describes the flood wave, for gentle slopes the diffusion wave equation is an appropriate approximation and for intermediate slopes the full dynamic wave equation best describes the flood wave. This criterion is difficult to apply inpractice since gentle, intermediate and steep slopes were not well defined. Choi et al. (1993) also investigated the applicable ranges of the kinematic and diffusion wave equations based on channel bed slopes, a dimensionless parameter describing depth increase at the upstream boundary and Froude numbers. They found that the relative magnitude of the pressure term in mild slopes was large, suggesting that the kinematic wave equation may not be applicable for mild channel bed slopes. In response to increasing Froude numbers, the applicable range of the kinematic wave equation increases while the applicable range of the diffusion wave equation decreases. The applicable range of both models increase with increasing channel bed slope. In general, the applicable range of diffusion wave models was found to be wider than that of kinematic wave models for open channel flows. 14 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach The above discussion on the various forms of the conservation of momentum equation suggests that while uniform flow conditions are rare, the kinematic wave equation can be applied to simulate the flow in relatively steep rivers (S 0 > -0.002). The choice between the diffusion and dynamic wave equation is not as straightforward and would likely require some level of engineering judgement. 2.1.3 Numerical Solution The St. Venant equations presented in the previous section are an example of differential equations that cannot be solved analytically. Therefore, to obtain a solution numerical methods are used to convert the continuous differential equations into finite algebraic expressions that can be solved numerically. Before applying a numerical scheme to solve the St. Venant equations for a particular reach or channel, the study area has to be discretized. The modeled area is composed of an infinite number of points that would require an infinite number of equations; therefore, a mesh or grid is established over the study area consisting of a finite number of points in space and time. Each point in the mesh or grid represents a location where the algebraic approximations of the St. Venant equations are applied and solved to obtain the depth and velocity components. The computational grid can be either rectilinear, curvilinear, or a mesh composed of polygons. A rectilinear computational grid defines the x and y co-ordinates relative to an established geographical frame of reference such as UTM (Millar and Barua, 1999). This type of grid is not suitable for simulating flow in curved river channels since the grid dimensions or aspect ratios cannot be changed. Curvilinear grids, which are better suited to natural river channels, allow the x and y co-ordinates to be transformed into stream-wise and cross-stream directions. Complex flow patterns can be more accurately modeled using a mesh consisting of elements of varying sizes and shapes (Steffler, 2000). In areas that are particularly important or rapidly varying areas, including boundaries, smaller elements can be defined. The vertices 15 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach of the elements, which are usually triangles, are the nodes or computational points in the model. In 3D applications, the mesh usually consists of bricks or tetrahedrons rather than triangular elements (Viessman and Lewis, 1996). The most common numerical methods used in hydrodynamic modeling are finite difference and finite element, although finite volume methods are also used. The finite difference method utilizes a Taylor series expansion to convert the St. Venant equations to algebraic finite difference equations. Descriptions of the use of the finite difference method in solving the St. Venant equations can be found in Wood (1993) and Graf (1998). Either rectilinear or curvilinear grids can be used with a finite difference scheme, whereas finite element methods require a mesh to be defined over the modeled area and therefore offer more geometrical flexibility. Finite element methods use numerical integration to solve integral equations that have been converted from the governing differential equations. A solution can be obtained through a weighted residual technique that involves using a trial function to obtain an approximate solution (Steffler, 2000). The trial function is specified and has a number of adjustable degrees of freedom. Introducing trial functions for the unknown variables (e.g. H, u and v for a 2D model) into the equations of continuity and conservation of momentum will result in a residual, as the equations are not exactly equal.to zero. Since the goal is to make the residuals as small as possible, they are multiplied by a weighting function and integrated over the element area with the result.set to zero. A separate weighting function, which is adjusted until the residuals are minimized, is used for each degree of freedom of the trial function. A detailed review of different numerical method formulations is beyond the scope of this thesis. Further information on the formulation and application of the finite element method in computational hydraulics can be found in Davis (1975), Walters and Casulli (1998) and Vreugdenhil (1989). Other numerical methods that are used in hydrodynamic models are finite volume (Demuren, 1993) and the method of characteristics (French, 1985; Wood, 1993). 16 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Discretization schemes such as finite difference and finite element result in a set of nonlinear algebraic equations which are solved using either explicit or implicit methods. In the explicit method, the unknown values are solved sequentially along a time line from one distance point to the next (Chow, 1998). Since each value is calculated independently, matrix solutions are not required, leading to faster computer solutions. Implicit methods are more complex mathematically but generally more stable for larger time steps (Millar and Barua, 1999). Matrix solutions are required to simultaneously solve all the unknown values on a certain time line as all variables at a new time are considered to depend on each other as well as the values at the previous time step. The main.problem with explicit methods is that the time step, At, is restricted in order to keep the solution stable. The Courant condition {At < Ax/c, where Ax is the distance between computational nodes and c is the wave celerity equal to (gH)05) must be satisfied in order to achieve numerical stability, which refers to the propagation and magnification of round-off errors and other errors associated with the discretization approximations. Although explicit formulations require satisfaction of the Courant condition, it alone does not necessarily guarantee stability. Barnett (1976) identified four categories of instability in unsteady open channel flow computations, of which only Courant and non-homogeneous instability usually need to be considered. Non-homogeneous instability arises from inadequate treatment of terms that are not homogeneous with the first order differential terms in the St. Venant equations. Unlike explicit methods, the satisfaction of the Courant condition is not a necessary requirement to achieve numerical stability for implicit formulations. For a finite element mesh, the Courant condition is related to the time taken by a shallow water wave to travel from one node to the next. If the mesh spacing is reduced in any area, the time step for the entire simulation must be decreased (Steffler, 2000). 17 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Numerical solution schemes are prone to other possible computational problems besides numerical instability, such as numerical dispersion (Maghsoudi and Simons, 1993), which refers to the artificial spreading of the flood wave due to round-off errors. An example of this would be in the case of a kinematic wave which simulates pure translation. No attenuation of the flood peak would be expected and any dispersion that occurs would be numerical dispersion. The convergence of numerical methods, that is its sequence of solutions along a timeline should asymptotically approach a fixed value, is another consideration (Thompson, 1993; Wood, 1993). A numerical model's consistency refersto whether the sequence of solutions converges to the solution of the differential equations which govern the physical phenomenon being modeled. The solution of the St. Venant equations for natural rivers introduces additional computational challenges. Many hydrodynamic models are not capable of simulating either supercritical or transcritical flows which describe the existence of both subcritical and supercritical flow simultaneously within the computational domain (Meselhe et al., 1993). In rivers, local areas of supercritical flow may develop over gravel bars when they become submerged at high flows or at other local high bed elevations. Severe width contractions can also induce localized supercritical flow conditions. The wetting and drying of bars, islands and floodplain areas is also an important consideration in hydrodynamic models (Ghanem et al., 1995a). Rapid drying in certain areas can potentially result in the model crashing if the drying occurs within a time period less than the simulated time step. To address this problem, some models turn cells or elements on and off and insert no-flow boundaries between them based on a minimum depth criteria or they simply change the fluid properties at very low depths so that a very thin layer of fluid is always present (Steffler, 2000). An alternative approach is to utilize groundwater flow equations in the model so that a common free surface is calculated, both above the ground in the channel and below the ground on dry surfaces. 18 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Generally, as the level of discretization of the modeled area increases, i.e. the computational grid or mesh gets finer, the numerical method approximation of the governing differential equations gets better (Viessman and Lewis, 1996). However, accuracy gained through a finer discretization comes at the expense of the increased computational resources required to obtain a solution in a reasonable amount of time. For all hydrodynamic models, boundary conditions must be specified in order to obtain a unique solution to the governing equations. The flow, water surface elevation or a stage-discharge relationship can usually be entered as a set of boundary conditions. For example, typical boundary conditions for hydrodynamic models may be a discharge specified at the upstream boundary and a fixed water surface elevation.specified at the downstream boundary. Often flow boundaries are located some distance upstream and downstream from areas of interest to account for lateral variation in the discharge or water level. In addition, the initial conditions at each computational node throughout the modeled region need to be defined at the beginning of the first time step. In some models an initial estimate of the upstream water level will allow the model to compute initial depths at all computationalpoints withinthe modeled domain if a fixed water level downstream has already been specified. 2.1.4 Hydrodynamic Model Examples and Applications Although steady flow conditions are rare in natural rivers, in some cases the temporal variations of the flow may be small thereby allowing such an assumption to be made over short time intervals. Generally, the steady flow assumption tends to be conservative, especially in flood plain mapping studies (Environment Canada, 2000). Backwater models such as HEC-2 and H E C - R A S are available to compute water surface profiles based on steady, gradually varied flow principles for channels with varying and irregular cross-sections. These models were developed by the Hydrologic Engineering Center of the US Army Corps of Engineers and they 19 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach iteratively solve the 1D steady-state energy equation and evaluate energy losses due to friction based on the Manning equation (US Army Corps of Engineers, 1991). HEC-2 can simulate either subcritical or supercritical flow and compute water surface profiles through structures such as bridges, culverts and weirs. HEC-RAS, which includes a graphical user interface, was developed subsequent to HEC-2. Several computational improvements were made in H E C - R A S including utilizing the momentum equation in situations where the water surface profile becomes rapidly varied, such as hydraulic jumps, flow around bridge piers and at stream junctions. (US Army Corps of Engineers, 1997). Since the majority of river engineering problems involve transient flow conditions, hydrodynamic models based on the governing equations and solution techniques described previously are often employed. The one-dimensional hydrodynamic model, ONE-D, was developed by Environment Canada in 1970 and has been used to model flow in several Canadian rivers (Environment Canada, 2000). The model, which includes an optional water quality routine, uses an implicit finite-difference scheme to solve the 1D St. Venant equations for single or multi-channel networks and estuaries. ONE-D has been applied successfully to the Fraser River (Morse et al., 1991) and to the Mackenzie River Delta (Fassnacht, 1997). The ONE-D model is employed operationally on the lower Fraser River for the simulation of flow in tidal-affected regions where flow reversals temporarily occur in several reaches (Environment Canada, 2000). ONE-D-SED is an extended version of the original fixed-bed model and includes a mobile bed morphodynamic module that can simulate sediment transport processes in river networks (Morse et al., 1991). Another commonly used 1D hydrodynamic model is MIKE 11, which was developed by the Danish Hydraulic Institute (DHI). MIKE 11 uses a six-point implicit finite difference scheme to obtain a solution for the 1D St. Venant equations. In addition to the hydrodynamic module, MIKE 11 includes components that can simulate rainfall-runoff, advection-dispersion, water quality and cohesive and non-20 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach cohesive sediment transport. Applications of MIKE 11 include flood forecasting and reservoir operation, simulation of flood control measures, operation of irrigation and surface drainage systems, design of channel systems and tidal and storm surge studies in rivers and estuaries (DHI Water and Environment, 2001). The DHI have also developed the 2D models, MIKE 21 and MIKE 21C. The main difference between these two models is that MIKE 21 uses a rectilinear computational grid, while MIKE 21C employs a curvilinear grid. A morphodynamic module is also included in the model, which allows for the simulation of advection-dispersion, sediment transport and bank erosion. MIKE 11 was used in conjunction with MIKE21C to simulate flow and sediment transport conditions in the Brahmaputra River. Design flood levels and the sediment budget were estimated with MIKE 11 while local hydraulics and river morphology were investigated with MIKE21C (Olesen, 1992; Engrob and von Lany, 1994). The design flood profile for the Fraser River Gravel Reach was recently updated using MIKE 11 as well (UMA, 2000). This model was also used to investigate mitigation scenarios involving dyke raising and gravel extraction that would alleviate the flood risk. Subsequently, Water Management Consultants (2001) utilized the morphodynamic module in MIKE 11 to simulate the effect of a proposed excavation in the Harrison Bar area to reduce the upstream design flood profile. The T E L E M A C models use semi-implicit finite element techniques to solve either the 2D or 3D St. Venant equations. These models include a mobile bed morphodynamic module which is applicable to gravel-bed rivers (Millar and Barua, 1999). T E L E M A C 2D and 3D have been used in a wide variety of modeling studies around the world, For example, T E L E M A C 2D has been used to determine the expected flow pattern near a captive jetty at Pipavav, India and to assess the impact of associated dredging on flow in the surrounding areas for several different scenarios. T E L E M A C 3D was used to determine the impact of new bridge piers on 21 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach the flow through surrounding channels and the Pearl Estuary in Hong Kong (HR Wallingford, 1998). Another example of a 3D hydrodynamic model is Delft3D, which employs an implicit finite difference scheme to obtain solutions for either 2D or 3D problems. H3D, which evolved from an oceanographic model of the Department of Fisheries and Oceans, uses a semi-implicit finite difference scheme with a rectilinear computational grid. While 3D hydrodynamic models are not as commonly used as 1D or even 2D models, they do have significant potential uses in the study of natural rivers. For example, a 3D model would allow the turbulent flow structure in a meandering channel to be examined so the secondary currents around the bend could be estimated (Demuren, 1993). Scott et al. (1999) applied MIKE11 and T E L E M A C 2D to the Fraser River as part of a forecasting system for water levels and river currents. The forecasts were intended to provide hydrodynamic conditions for vessel navigation and safety on the Fraser River. MIKE11 was used to simulate approximately 250 km of channels in the Fraser River delta and lower part of the river to obtain water level estimates. MIKE11 also provided the boundary conditions for T E L E M A C 2D, which was used to simulate the hydrodynamics in the lower 40 km of the river in order to provide information on current patterns in certain reaches of interest. There are several other examples of applications of 2D hydrodynamic models to specific modeling situations. A 2D depth averaged finite element model, RMA-2V, was used to predict flow conditions for the design of a proposed fish barrier on the Sacramento River (Shrestha et al., 1993). Simulation results were used to estimate riprap size on the barrier face to ensure stability. In another study conducted by Berger et al. (1993), the hydraulic performance of high velocity channels subject to supercritical flow conditions was evaluated using the 2D model, HIVEL2D. The goal of this work was to eventually develop a 2D numerical model that could simulate flow through man-made flood control channels so the potential location of oblique 22 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach standing waves and hydraulic jumps could be predicted and the superelevation of the water surface in channel bends could be determined. 2.1.5 River2D River2D is an implicit finite element two-dimensional, depth averaged hydrodynamic model based on a conservative Petrov-Galerkin upwinding formulation (Ghanem et al., 1995a; Hicks and Steffler, 1992). This model, which was developed at the University of Alberta, runs on a Windows platform and is public domain software. The form of the St. Venant equations and the bed resistance and transverse shear models presented in Section 2.1.2 describe the mathematical representation of natural river processes simulated by River2D. The model includes the ability to handle subcritical, supercritical and transcritical flows. Ghanem et al. (1995a) successfully tested River2D's ability to simulate supercritical flow and to predict the locations and heights of standing waves resulting from a channel constriction. In River2D, the effective roughness length, ks, is specified at every computational node and is the primary calibration parameter. A secondary calibration parameter is the user defined coefficient, e, shown in Equation [2.8] which controls the eddy viscosity. This coefficient is applied globally and is set to a default value of 0.5. More discussion on roughness parameters and turbulence models for open channel flow is provided in Sections 2.2 and 2.3, respectively. Furthermore, a sensitivity analysis on these River2D calibration parameters is presented in Section 5.0. For 2D models, calibration can be achieved by comparing simulated water surface elevations and velocities to measured values. For large river systems, such as the Fraser River, this velocity information can be obtained by using acoustic Doppler current profilers, which are described in Section 2.4. The most critical step in the modeling process using a 2D hydrodynamic model such as River2D is to obtain an accurate representation of the bed 23 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach topography. Survey data is used to create a digital elevation model (DEM) which provides the bed topography input data in the form of x, y and z coordinates for River2D. In addition to the bed topography data, the input file also requires ks to be defined at every node. A bed topography module is included in River2D to graphically edit and refine the input file. The computational discretization, i.e. the finite element mesh, can be created from the bed topography file using the mesh generation program, R2D_Mesh (Steffler, 2000; Ghanem etal. , 1995a). The computational boundaries, including inflow and outflow sections, are graphically defined at this stage. The modeled domain is then filled with nodes at a spacing specified by the user before the mesh generation utility is used to perform a triangulation. Only triangular elements are included in the mesh as these elements are considered the simplest possible in two dimensions and result in the minimum execution time for a given number of nodes (Steffler, 2000). In regions of high interest or where flow variations are large, for example near a sloping bank, more closely spaced nodes can be used to minimize errors. Once a satisfactory mesh has been created and boundary conditions specified, the River2D simulations can proceed to solve for the depths and velocities at all the computational nodes within the modeled domain. The boundary conditions required by River2D are typically a discharge at upstream sections and a fixed water surface elevation at downstream sections. The option of a depth-discharge relationship based on Manning's equation is also available as the downstream boundary condition. This relationship is shown below and is applied to each point across the section. q = KHm [2.13] In equation [2.13], q is the discharge per unit width (m2/s) and K a n d m are constants. Based on the Manning's equation, K = ^—^- and m = 5/3. This option for n the downstream boundary condition is slightly better behaved numerically and can 24 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach be used to get close to the final solution before changing to the known fixed water surface elevation (Steffler, 2000). The difficulty of dealing with wet and dry areas in computational hydraulics was briefly described in Section 2.1.3. In locations where the depth of flow becomes very shallow or where there is no water over a part of the modeled area, River2D utilizes groundwater flow equations. In these areas the continuity equation is replaced by equation [2.14] so a continuous free surface which extends both above and below the ground is calculated. In the above equation, T is the transmissivity (m2/s), which is set by default to 0.1, and S is the storativity of the artificial aquifer, which is taken as unity and describes the volume of water absorbed or expelled per unit change in head and area. The ground surface elevation (m) is represented by zb. The actual groundwater discharge is negligible since the transmissivity is set to a low value although it can be changed by the user. Ghanem et al. (1995a) describes several tests performed with River2D to examine and verify its ability to adequately simulate flow over wet and dry areas in various scenarios. A fish habitat component is also included in River2D and is based on a weighted usable area (WUA) concept (Bovee, 1982). The fish habitat model uses the computed depths and velocities from the hydrodynamic simulations for its calculations. The general approach of this component is derived from the Physical Habitat Simulation System or PHABSIM, which requires precise values of depth and velocity to produce relationships between streamflow and usable habitat area for different life stages of various fish species (Ghanem et al., 1995a). 25 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Examples of the application of River2D include the Kananaskis River, which is a braided mountain river ( S 0 « 0.01) in Southwestern Alberta (Christison et al., 1999). The study site was a 2 km long by 0.5 km wide anabranching reach with simulated discharges ranging from 2 m 3/s to 50 m 3/s. Wooded islands and gravel bars divide this portion of the gravel-bed Kananaskis River. Typical channel dimensions were as high as 10 m in width and 1 m in depth. The channel topography was represented by almost 9,000 survey points which were organized into 900 feature lines that defined the geometric characteristics of the river (e.g. thalweg, top and bottom of banks, ridges and pools). Computer resource limitations required that the study reach be sub-divided into fourteen linked sub-reaches as over 44,000 nodes were used to represent the channel topography. Water levels and velocities were surveyed in one of the sub-reaches to provide the model's calibration data. Although the calibration data was limited, it did indicate that the model was performing reasonably well. The Elbow River near Calgary was used as the study site in a comparison between River2D and a 1D model in simulating flow through a small habitat stream (Waddle et al., 2000). A low (4 m3/s) and a high (21 m3/s) discharge were simulated for the study site, which is 315 m long, 35-50 m wide and has a mean slope, S 0 « 0.0037. The channel bed consists of cobbles and boulders. For the 1D simulations, a backwater model (WSP) was coupled with the IFG4 program to distribute the velocities across the channel based on a set of point values of Manning's n calculated at a calibrated flow. The results for velocity and water surface elevation were comparable between the 1D model, River2D and measured values. However, River2D was found to more accurately capture complex flow patterns where significant transverse flow was present. Lacey (2001) investigated the use of River2D for the assessment and design of instream channel restoration works. He applied River2D to a side channel of the Chilliwack River to assess the morphological and hydraulic effects of instream large woody debris and boulder structures. The modeled section of the side channel was 26 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach approximately 320 m long by less than approximately 20 m wide. Computed depth and velocity measurements were found to compare favourably to measured values. 2.2 Flow Resistance Resistance to flow in open channels is due to several factors including the nature and scale of the boundary material and distortions of the channel alignment. Additional flow resistance is encountered in locations where the flow patterns change significantly which result in increased energy dissipation. Carling (1996) described four types of flow resistance: grain, form, internal distortion and spill resistance. Grain resistance is generated by individual sand and/or gravel particles on the channel bed while form resistance is due to grain protrusion and the presence of bedforms and bars. Internal distortion resistance is associated with banks and changes in channel alignment. An example of spill resistance, which is typically associated with rapid flow conditions, is the flow over boulders in a mountain stream. A broader division of the total resistance can be made into simply grain and form resistance, where the grain resistance is defined as above and form resistance takes into account all factors contributing to flow separation and subsequent eddy losses. There are several different roughness parameters which can be used to describe the flow resistance in a reach. Two such parameters are Manning's n and Chezy C. The commonly used Manning equation [2.15] is an empirical one-dimensional formula. U = -R2,3S01'2 [2.15] n In the above equation, U is the average velocity (m/s) in the channel, R is the hydraulic radius (m) and n is the Manning's roughness coefficient (s/m 1 / 3). Values for Manning's n for a range of channel types are tabulated in Chow (1959). The relationship between Manning's n and the Chezy coefficient (m 1 / 2/s) is shown below 27 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach and the dimensionless form, C s , which appears in equations [2.4] to [2.6] and [2.8], is found by dividing C by y[g . R U 6 C = — [2.16] n A value of the resistance can also be estimated from an examination of velocities in an open channel since the velocity profile is dependent upon the roughness length, ks. Turbulent flow over a rough boundary in an open channel can be described by the depth-integrated law of the wall, shown below. - ^ = 5.751og| v ks j [2.17] All the variables in equation [2.17] have been defined previously except U*, which is the shear velocity (m/s) and can be expressed as a function of the mean bed shear stress (N/m2), T0, as follows: U* = ^ = JgHS~0 [2-18] Another measure of roughness is provided by the Darcy-Weisbach friction factor, f, which can be related to the mean velocity and the Chezy coefficient as shown in equation [2.19]. (7 C — = -$== I- [2.19] The work of Nikuradse (1933), who studied the behaviour of f and ks in circular pipe flow, led to the following relation by Williamson (1951): 28 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 1/3 / = 0.113 [2.20] The applicability of equation [2.20] to open channel flow situations is justified since the effects of cross-sectional shape are not significant (see Henderson, 1966). Therefore, substituting equations [2.16] and [2.20] into equation [2.19] yields a relationship between the Manning's n coefficient and the roughness length, ks of the following form: n = 0.038kj'6 [2.21] Equation [2.21] is very similar to Strickler's formula (1923) which relates Manning's n to the median diameter (m) of the bed material, d50: n = 0.042d50V6 [2.22] Equations [2.21] and [2.22] provide a measure of the grain resistance since they relate Manning's n directly with sediment size. Millar (1999) investigated the division of grain and form resistance for gravel-bed rivers. By calculating the grain resistance for 176 gravel reaches, he defined a relation [2.23] between ks and dso which formed a lower bound for the overall flow resistance. ks = 5.9d50 [2.23] Other investigators, such as Bray (1982) have suggested similar relationships between ks and a characteristic grain size diameter: ks = 6.Sd50 [2.24] 29 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Nikuradse (1933) originally defined ks for uniform boundary sediment to be approximately equal to the sediment diameter, d, which is supported by the similarity between equations [2.21] and [2.22]. For non-uniform boundary sediment ks would be expected to be approximately equal to the mean or median grain diameter, however equations [2.23] and [2.24] suggest that ks is actually a multiple of a characteristic grain size. The effect of the multiplier, which is empirically determined, is to produce the actual energy loss due to all forms of resistance when applying equations such as [2.17] to rough open channel flow. Clifford et al. (1992) attempted to provide a physical explanation of the roughness length multiplier for gravel-bed rivers. They reasoned that the multiplier takes into account small-scale form resistance associated with irregular particle groupings such as cluster bedforms. The treatment of the roughness parameter is one of the fundamental differences between 1D and 2D models. For example, in 1D hydrodynamic models Manning's n is often used to account for all forms of resistance, including both grain and form roughness. One-dimensional models assume that many of the two-dimensional effects in natural rivers are accounted for by the roughness term. A consequence of this assumption is that in locations where energy dissipation is unusually high, the 1D roughness parameter, e.g. Manning's n, may have to be increased to unrealistic values to simulate the increased energy losses. Conversely, in 2D models the resistance term only.represents direct bed shear and should not be used to account for form resistance such as losses associated with expansion and contraction in planform geometry (Waddle et al., 2000). The roughness parameter should be used to represent grain resistance only and therefore observations of bed material and bedform size are usually sufficient to establish reasonable roughness estimates for the application of 2D models (Steffler, 2000). Based on the above discussion, when simulating flow conditions in a gravel-bed river using a 2D hydrodynamic model ks can generally be taken to be approximately equal to c/50, although the presence of bedforms on the channel bottom would require the roughness length to be increased. Furthermore, equations 30 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach [2.21] and [2.22] can be used to provide an estimate of the 2D Manning's n, which is essentially a measure of the grain resistance. Prestegaard (1983) studied the various forms of resistance in twelve straight gravel-bed streams. She concluded that bar resistance contributes between 50 to 75% of the total resistance in these situations, while the remainder is due to grain resistance. The significance of bar resistance in unvegetated gravel-bed rivers was also studied by Hey (1988). He stated that bar form resistance was due to the accelerations and decelerations in the flow between pools and riffles so he analyzed the hydraulics between pool-riffle sequences to quantify the contributions of grain and bar form resistance. In the application of a 2D hydrodynamic model, the effects of bar form resistance would already be accounted for by the model's physical formulation and would not have to be considered when estimating ks values. 2.3 Turbulence Models This section is not intended to provide a rigorous explanation of the turbulent flow structure in open channels, which is well beyond the scope of this thesis. Instead, some additional background information on turbulence models will be presented. A much more thorough explanation of the relationships between mean flow properties and the shear stresses caused by turbulence is provided in the A S C E Task Committee on Turbulence Models (1988) and Rodi (1984). If an attempt would be made to accurately model the actual turbulent structure in a river, a three-dimensional computational grid would be required with spacings smaller than the smallest turbulent motion and a time step less than that associated with the fastest eddies. The computational requirements to obtain such a solution are extreme and therefore one of two alternate approaches can be taken to solve this problem. The first method is to use statistical theories on the structure of turbulence based on observations and empirical correlations while the second approach employs a semi-empirical analysis of the effects of turbulent motions on 31 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach mean flow properties (Lane, 1998). The second method is the one commonly used in hydrodynamic models. The transverse shear stresses, which are commonly referred to as Reynolds shear stresses, described in Section 2.1.2 represent the lateral momentum transfer in a river that is due to turbulence. In the 2D St. Venant equations, these transverse shear stresses represent four additional unknown terms. Since they can never be zero in turbulent flow and there have not been additional equations developed to solve for these quantities, turbulence models are used which relate the transverse shear stresses to flow parameters which are either known or can be readily solved for. These turbulence models are typically referred to as turbulence closure schemes since they allow the governing equations for shallow open-channel flow to become closed when considering the effects of turbulence. Most turbulence closure schemes are based on a Boussinesq (1877) approximation which allows an estimation of the eddy viscosity coefficient to be made. Zero, one and two-equation turbulence models are available for this purpose (Lane, 1998). Zero-equation models specify a constant eddy viscosity or determine the eddy viscosity based on a mixing length hypothesis. One-equation models attempt to account for the convective and diffusive transport of the turbulent velocity scale by incorporating the kinetic energy of the turbulent motion per unit mass, k. The velocity scale for large-scale turbulent motion is represented by 4k . More complex two-equation models also consider the effect of transport processes on the length scale, /, in addition to the velocity scale. The relatively straightforward turbulence model in River2D is an example of a zero-equation model where the horizontal turbulence is simulated by including a constant eddy viscosity coefficient which allows the transverse shear stresses to be related to the mean flow velocity gradients. Other hydrodynamic models employ more sophisticated turbulence closure schemes such as T E L E M A C 2D, which provides the choice between the ks, Smagorinsky or Elder turbulence models. A 32 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach description of a two-equation turbulence model is described in Lane and Richards (1988) who tested and applied a 2D hydrodynamic model to a multi-thread reach of a proglacial stream. The model used in that study calculated the eddy viscosity from k and a mixing length, /. Since the turbulence stress terms usually have a stabilizing effect on the numerical solution process, they are sometimes exploited in some models to stabilize schemes that would be otherwise unstable through the addition of excessive and unrealistic values (Ghanem et al., 1995a). Generally, the overall solution of the St. Venant equations for open channel flow is not sensitive to the selection of the user defined coefficient, e in equation [2.8] since the effect of the transverse shear stresses are small (Steffler, 2000). The limits of this coefficient for rivers range between 0.2 and 1.0 in order to keep the eddy viscosity within a reasonable range. Fisher et al. (1979) proposed that for flow in open channels, the value of the eddy viscosity, v', which represents isotropic turbulence, can be found by: v'=(0.14±0.07)(7*/f [2.25] As previously mentioned, the eddy viscosity is a constant that is applied globally in River2D. The sensitivity of the user defined coefficient, and consequently the eddy viscosity, is investigated in Section 5.0. 2.4 Acoustic Doppler Current Profilers (ADCP) Calibration of 2D hydrodynamic models can be done through the comparison of both modeled water levels with measured gauge readings and modeled velocities with observed values. For smaller streams, conventional current meters can be used to obtain velocity measurements throughout the study reach, however for larger rivers such as the Fraser River this approach is not feasible. In these 33 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach situations, acoustic Doppler current profilers (ADCP's) can be used to measure velocity along several transects in the study reach. A D C P ' s use the Doppler shift, which is the change in the observed sound pitch that results from relative motion, of an acoustic signal to measure 3D water velocity (RD Instruments, 1996). In general terms, if a source of sound is moving relative to a receiver, the frequency of the sound at the receiver is shifted from the transmit frequency. These frequencies can be related to the relative velocity between the source and the receiver by accounting for the speed of sound in water. A D C P ' s use the Doppler shift by emitting short acoustic pulses at a fixed frequency through the water column and listening to echoes returning as a result of some of the acoustic energy being scattered back to the instrument from small particles suspended in the water. A key assumption is that these small particles move at the same velocity as the water. Only a small fraction of the sound is reflected back to the A D C P , and this amount is Doppler shifted. RD Instruments (1996) provide a description of the principles of operation and potential sources of error associated with ADCP 's . Four acoustic beams are used by the A D C P to compute three velocity components. These instruments are able to measure velocity profiles by dividing the water column into uniform segments or bins. Computed velocities are averaged over the entire depth range of each bin and are transformed to an earth-referenced coordinate system using an internal fluxgate compass. The A D C P can be mounted on the side of a moving vessel to allow velocity profiles to be computed along desired transects in a river. While short transmit pulses are used to obtain vertical velocity information, longer pulses are used for the A D C P ' s bottom tracking capability. Bottom tracking is used to correct water column velocities for boat motion and it can also provide depth information that can be used to compare with model results. Information obtained 34 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach from bottom tracking can be combined with the velocity profile data to obtain estimates of discharge at each surveyed transect. Some of the computed velocity profile information should be rejected since the echo from a hard surface such as the river bottom is much stronger than the echo from suspended particles in the water. For A D C P ' s with a 30° beam transducer angle, the bottom 15% of the data is likely biased towards zero due to the channel bottom and should be ignored. Another source of error that must be considered is due to the random error of the velocity measurements. For example, velocity errors based only on a single ping can range from a few mm/s to as high as 0.5 m/s (RD Instruments, 1996). Averaging of several pings reduces the standard deviation of the velocity error by the square root of the number of pings. An example of the use of A D C P ' s to measure water velocity is described in Chu et al. (1993) who used a BroadBand A D C P to survey tidal currents in two inlets on the south shore of Long Island, New York. In another study, Rennie (2001) used an A D C P ' s bottom tracking capability to develop a technique to measure the apparent velocity of bedload. The study site for this work was the Fraser River Gravel Reach. 35 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 3. F R A S E R RIVER G R A V E L R E A C H 3.1 Setting The Fraser River is the largest in British Columbia with a length of 1370 km from its headwaters in Mount Robson Provincial Park to its mouth at the Strait of Georgia. The river has a watershed of approximately 250,000 km 2 , including most of south-central British Columbia (see Figure 1.1), although a drainage diversion in the upper Nechako basin reduced the effective drainage area to 233,000 km 2 in 1952. The drainage basin is bounded by the Coast Mountains to the west and the Rocky and Columbia mountains to the east. After flowing northwestward from its source, the Fraser River turns abruptly in a southerly direction in the vicinity of Prince George and continues through the Coast Mountains, where it has created a deep gorge. At Yale, the river exits from the canyon, where the gradients are relatively steep, and flows 190 km through the lower Fraser Valley to the sea. From Yale to Laidlaw the Fraser River flows in a single channel confined by bedrock walls, landslide debris and Pleistocene terraces. Downstream of Laidlaw is a 55 km reach that extends to Sumas Mountain and is characterized by numerous mid-channel islands and low-order braiding. This section of the river is referred to as the gravel reach (Figure 3.1). Several gravel bars are present in this reach including lateral bars that are attached to the shore and to vegetated islands, point bars at bends and mid-channel bars in areas of flow expansion (McLean et al., 1999). The main channel meanders between the right and left banks within the active braided zone. The distribution of subsurface channel material is bimodal, containing both coarse gravel and medium sand fractions. The median diameter of the gravels is typically between 25 and 30 mm. With the exception of the sand sizes, which are absent, the surface material that is exposed on the bars has a similar composition as the subsurface material. 36 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Figure 3.1: Fraser River Gravel Reach The river changes from a braided gravel reach to a single-thread sand-bed channel at Sumas Mountain, which is near Mission. The river flows in this single channel for 50 km to the head of the alluvial delta at New Westminster. The Fraser River is tidally influenced below Sumas and the tidal range at Mission varies from a few centimeters during the spring freshet to over 1 m during the highest winter tides. Hydraulic characteristics for the mean annual flood at three gauging stations in the lower Fraser River are provided in Table 2.1. The three Water Survey of Canada gauging stations listed are located at Hope (08MF005), Agassiz (08MF035) and Mission (08MH024). Hope is located about halfway between Yale and Laidlaw while Agassiz is located approximately 17 km downstream of Laidlaw. 37 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Table 3.1: Hydraulic characteristics at the mean annual flood (after Church and McLean, 1994) Station Q Depth Width Velocity Froude Slope Surface (m3/s) (m) (m) (m/s) # d5o (mm) Hope 8,770 10.1 268 3.2 0.32 0.60 x 10"3 100 Agassiz 8,760 6.6 512 2.6 0.32 0.48 x 10"3 42 Mission 9,790 12.6 540 1.5 0.13 0.05 x 10"3 0.38 3.2 Hydrology The dominant hydrological event is spring snowmelt since much of the Fraser River basin has elevations above 1,000 m. The freshet typically begins in April with the peak flows occurring in late May, June and early July. The long-term mean flow at Hope is 2,830 m 3/s and the mean annual floods at the three gauging stations in the gravel reach are shown in Table 3.1. The discharge at Hope and Agassiz is basically identical which indicates no local inflow between those stations. The increased discharge at Mission is primarily due to inflows from the Harrison and Chilliwack Rivers (see Figure 1.1). These two rivers add about 4.6% to the Fraser River's drainage area while increasing the mean flow by about 18%. The largest flood on the Fraser River occurred in 1894 and had a peak discharge of approximately 17,000 m 3/s at Hope. At the time, the lower Fraser Valley was sparsely populated and thus damage was not severe. However in June 1948, another very large flood (Q = 15,200 m 3/s at Hope) resulted in $20 million of damage since the adjacent lands had become a highly developed agricultural area with increasing commercial, industrial and suburban residential developments. These large floods prompted dyke construction along much of the lower Fraser River. The design flow for this flood protection work was based on the 1894 flood of 17,000 m 3/s at Hope which has an estimated return period of 200 years. There are no flow records of earlier or larger flood for the Fraser River than the 1894 event. 38 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 3.3 Current Issues Between Laidlaw and Sumas Mountain, the bed of the Fraser River is rising as the river slows due to a decrease in gradient which results in gravel deposition. In this reach, the movement and deposition of gravel from upstream sources is increasing the flooding risk in some areas as the capacity of the existing dykes to contain floodwaters may have been reduced. In addition, the movement of gravel within this reach results in the river channel being shifted which has resulted in increased bank erosion in some areas that leads to a loss of land and threatens the dyking system. To manage these risk, possible strategies are to raise or set-back the dykes, armour the banks or attempt to re-align the channel by reducing sharp bends which creates a backwater effect and raises the upstream flood profile. Another option is to remove gravel from areas of aggradation in an attempt to increase the river's conveyance. There is an increasing demand for gravel resources for construction applications in the Lower Mainland and Fraser Valley, however gravel mining in the Fraser River Gravel Reach could have adverse affects on the river's morphology and ecosystem. Although it would seem that gravel removal would help to alleviate the flooding risk and also to satisfy the demand from potential commercial extractors, the careful management of the river's resources is key to maintaining its ecological character. Furthermore, the effects of gravel extraction on the river cannot be predicted precisely. For example, it is difficult to determine how the removal of gravel in a particular area would impact on the hydraulics at a downstream location or even how much gravel would have to be removed to mitigate the flood hazard in areas of concern. The hydraulic modeling study of the Fraser River Gravel Reach that was recently completed by UMA Engineering Ltd. (2000), focused on updating the design flood profile. In that study, which identified several areas along the river where the existing dykes are providing inadequate protection, various mitigation options were 39 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach investigated based on combinations of dyke raising and gravel extraction. While the 1D model used for that work provides preliminary estimates on overall gravel volumes that would have to be extracted to alleviate the flooding risk on a reach-wide scale, the impacts of such mining operations on the river's morphology and ecology, which would likely be substantial, cannot be inferred. The complex ecology of the Fraser River Gravel Reach is a result of the diverse habitat that have been formed and influenced by the river's natural processes such as gravel transport and deposition and the annual spring freshet. The river provides habitat for several species of mammals, amphibians and birds in addition to at least 28 species of freshwater fish (Fraser Basin Council, 2001). Many of these fish species support in-river recreational and aboriginal fisheries and ocean and estuary commercial fisheries. The gravel reach is also an important migration route for several species of salmon and steelhead and is also used by several fish species for spawning and rearing. Two habitat classification systems have been developed or are in the process of being developed to provide useful habitat information that can aid in the management of the river's resources. The Department of Fisheries and Oceans developed the first classification system while the Department of Geography at UBC is currently working on a classification system based on year-round field sampling (Church et al., 2000) and relating gathered information to previously conducted channel morphology and sediment budget studies (McLean et al., 1999 and McLean and Church, 1999). First Nations communities have several interests in the river as well, ranging from their rights to certain traditional and archaeological sites to economic interests that include fishing and the river's gravel resources. The use of this reach as a recreational area and a navigation route also require consideration in the river's management. A 2D hydrodynamic model, such as River2D, has the potential to aid in the decision-making process as different alternatives to alleviate local flooding and/or 40 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach erosion risks can be investigated and possible downstream impacts of any in-river work can be assessed. As described in Section 1.2, this thesis investigates the potential analytical uses of 2D hydrodynamic models through the application of River2D to selected sections of the Fraser River Gravel Reach. The survey data used in the modeling process is described in the following section. 41 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 4. S U R V E Y DATA 4.1 Floodpla in Topography A floodplain topography survey of the entire dyked corridor of the lower Fraser River was conducted on March 16-19, 1999 by Terra Surveys Ltd. using a helicopter equipped with a laser unit (UMA Engineering Ltd., 2000). All above-water areas were surveyed with this technique including floodplains, islands and most of the gravel bars, which were exposed since the river was at a relatively low stage during the survey. The laser beam could not penetrate water therefore a different survey method using a depth sounder was used to obtain bathymetric information. This survey is described in Section 4.2. The laser unit measured the distance to ground surfaces and it's position in space was tracked by a differential G P S . This data was combined to obtain the topographic information. Any angular deviations of the laser beam were noted and accounted for and video equipment was used to simultaneously record the ground image being targeted by the laser. Points created due to laser reflections from vegetation were deleted so only the actual ground surface was represented. The data was collected in a series of transects which were approximately 200 m apart in the gravel reach. Most of the surveyed transects were perpendicular to the corridor of the dyked river but not necessarily perpendicular to the direction of flow in the individual channels. UMA Engineering Ltd. (2000) reviewed the quality of this data by comparing it to previously surveyed dyke crest elevations and to the design profile of the C P Railway line across Seabird Island. They concluded that the topographic data was generally reliable and although the average error is likely quite low, individual points may be as much as 0.25 m in error. 42 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 4.2 River Bathymetry Publ ic W o r k s and Government Serv i ces C a n a d a ( P W G S C ) conducted a hydrographic survey of the gravel reach from a boat using a depth sounder and differential G P S from August 15-25, 1999. The bathymetric data w a s col lected in t ransects which were 200 m apart, similar to the topographic data, however the locations of the transects from these surveys are offset by up to 100 m in some areas. The main channel and all other navigable s ide channe ls were surveyed, resulting in 515 measured cross-sect ions. The survey points ac ross each transect were very c losely spaced , approximately 0.5 m apart on average. The data indicates that the channel bed is generally smooth in the gravel reach. Severa l c ross sect ions do not extend to the full wetted portion of the cross-sect ion s ince depth soundings were not obtained near the banks at some locations. The data points were tied to the N A D 8 3 datum and latitude and longitude coordinates were used although they were later converted to U T M . The data w a s processed and checked for quality by P W G S C although compar isons could not be made with previous channel surveys as the annual spring freshet continually modif ies bed elevat ions. 4.3 A D C P Survey A s part of the hydrographic survey, Publ ic W o r k s C a n a d a used an A D C P to collect velocity direction and magnitude information in the Fraser River Grave l R e a c h during July 30 to August 11, 1999. In addition to depth averaged velocity vectors, the A D C P data enabled the d ischarge distributions between the main channel and major s ide channels to be determined. The velocity information is based on a single ping at each vertical profile s ince no time averaging w a s done during the survey. G P S w a s used to track the boat 's posit ion. Every survey line was dupl icated and s o m e were run three or four t imes to obtain more than one d ischarge est imate for each transect. A lso , the boat w a s 43 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach unable to sample near the bank which would likely lead to a bias towards lower discharge estimates. A differential G P S was used to estimate water surface elevations along each of the A D C P transects, which can be compared to model results. 44 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 5. RIVER2D SENSITIVITY ANALYSIS 5.1 Roughness Length The primary calibration parameter in 2D models is the roughness parameter, which in the case of River2D is the roughness length, ks. A sensitivity analysis was performed on this parameter to determine its effect on the simulated depths and velocities. The study reach for the sensitivity analysis was directly downstream of the Agassiz-Rosedale Bridge and was approximately 2 km long by an average of approximately 800 m wide. All model parameters and conditions were kept constant while the ks value was varied for each run. A node spacing of 25 m was used and the boundary conditions consisted of a discharge of 3,000 m 3/s at the upstream section and a fixed water surface elevation of 13.5 m at the downstream section. The depth and velocity at 100 pre-determined points distributed throughout the reach were extracted from each run to allow comparisons to be made. These 100 points were scattered throughout both the main and side channel and also included deeper areas in the thalweg as well as shallow areas near the banks. Since the surface d50 at Agassiz is 42 mm (Table 3.1) and ks can generally be taken to be equal to dso (Section 2.2), the base case for this analysis was assumed to be ks - 0.04 m. No data on the spatial distribution of the dso grain sizes was available so uniform values of ks were applied to each node in the modeled region. Table 5.1 shows the different ks values which were used to assess the sensitivity of this parameter. The percentage change in ks relative to the base case of 0.04 m is shown in the first column while the third column shows the equivalent Manning's n value found by using equation [2.21]. 45 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Table 5.1: ks va lues used in the sensitivity analysis shown with corresponding Manning 's n % Change ks (m) n -90 0.004 0.015 -50 0.02 0.020 0 0.04 * 0.022 +100 0.08 0.025 +500 0.24 0.030 +1,000 0.44 0.033 +2,000 0.84 0.037 +5,000 2.04 0.043 * - base c a s e Figures 5.1 and 5.2 show the average velocit ies and depths, respectively, of the 100 points plotted as a function of the ks value. 1.40 j ^ 1-35 <r tf) 1.10 -i , 1 , 1 0.0 0.5 1.0 1.5 2.0 2.5 ks (m) Figure 5.1: Ave rage velocity versus ks 46 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 4.2 Figure 5.2: Ave rage depth versus ks The mean error (ME) of the velocit ies in each simulation, relative to the base case of ks = 0.04 m, was calculated as shown in equation [5.1]. The velocity at each point is denoted by U while / var ies from 1 to N, where N represents the number of points (in this c a s e 100) where the depth and velocity are extracted. Therefore, UikS# represents the velocity at the /* point for the ks value of the current simulation and Uiks=o.o4 represents the velocity at the / h point for a ks value of 0.04 m. ME = -Ui 100% /b=0.04 N [5.1] The mean depth error of each simulation was calculated similarly to the velocity errors. The mean velocity error var ies from almost +2% for ks = 0.004 m to 47 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach - 1 0 % for ks = 2.04 m while the mean depth error ranges from approximately - 2 % for ks= 0.004 m to almost +14% for ks= 2.04 m. To assess the spatial sensitivity of ks, the velocity and depth differences between the base case {ks = 0.04 m) and the most extreme case (/cs=2.04 m) were calculated and used to create the contour plots shown in Figures 5.3 and 5.4, respectively. The Agassiz-Rosedale Bridge is located at approximately x = 589,000 and the flow direction is from east to west. The blanked out portion in the center of the contour plots represent the vegetated island on Big Bar which divides the river at this section into the main channel to the north and a side channel to the south. Also shown on the figures are the locations of the 100 data points used for the comparisons. The green areas in Figures 5.3 and 5.4 indicate an increase in either the depth or velocity, while yellow regions represent decreases. An increase in ks should result in a decrease in velocities throughout the modeled region which is seen in Figure 5.3 except along part of the north bank and in most of the side channel. The average velocities are reduced by up to approximately 0.3 m/s in some locations in the main channel. Figure 5.4 shows an increase in depths due to the increase in ks, as would be expected. The depth change gradually varies from about +1 m at the upstream end to about -0.25 m at the downstream section. The depth decrease in the vicinity of the outflow section can likely be attributed to boundary condition effects as a fixed water surface elevation is used as the downstream boundary condition. While it is difficult to explain the velocity increases in the side channel, the contour plots do satisfy the basic conservation of mass principles. Figure 5.4 shows that there is a larger depth increase in the main channel than in adjacent areas of the side channel, however this is balanced by a larger velocity decrease in the main channel than in the side channel. The net result should be an identical discharge through the reach between the ks = 0.04 m and ks = 2.04 m simulations. 48 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 5452200-1—1 587000 587400 587800 588200 588600 X(m) Figure 5.3: Velocity change from ks = 0.04 m to ks = 2.04 m 5452200 587000 587400 587800 588200 588600 X(m) Figure 5.4: Depth change from ks= 0.04 m to ks = 2.04 m 49 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach The results of this sensitivity analysis indicate that while /c scan potentially have a very large range of values for gravel-bed rivers, the effects on depth and velocity, which are significant, are probably less than would be expected. A 5,000% increase in ks, which is perhaps not an unreasonable value in some localized areas since there could possibly be very large boulders in parts of the thalweg, results in an average change in velocities and depths of only 10% and 14%, respectively. Furthermore, a 5,000% increase in the base value of ks corresponds to a Manning's n of 0.043, which is quite possible for large gravel-bed rivers (Chow, 1959). In terms of the spatial sensitivity of ks, the somewhat curious anomalies in Figures 5.3 and 5.4 aside, large increases in the boundary roughness can be expected to decrease velocities by an amount approximately proportional to the depth at a given location. For example, in deeper areas, such as in the thalweg, the velocity decrease is greater than near the banks or in side channels. On the other hand, the depth changes appear to be more uniform as Figure 5.4 indicates a gradual variation in the depth difference in only the streamwise direction. Ghanem et al. (1995b) also found that the simulation results from River2D are not overly sensitive to the exact ks value. They modeled a 50 m wide reach at a discharge of 14.6 m 3/s and determined that a 100% increase in ks leads to approximately an 8% increase in calculated depths. In this sensitivity analysis, a 100% increase in ks resulted in an average depth increase of only 1%. This difference is due to the relative depths between the two studies. Ghanem et al. (1995b) modeled a reach with depths in the order of 1 m where the influence of boundary roughness would be much greater than in the Fraser River downstream of the Agassiz-Rosedale Bridge where the average depths are approximately 4 m for a discharge of 3000 m 3/s, but considerably deeper in the thalweg. This analysis and discussion would suggest that the boundary roughness is not of primary importance in 2D hydrodynamic models. An accurate representation of the bed topography would likely be a more important factor in obtaining quality simulation results than the precise definition of local ks values. This further supports 50 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach the assertion made in Section 2.2 that equating ksto dsowould suffice in most 2D modeling applications. 5.2 Eddy Viscosity Coefficient The same study reach and boundary conditions used to assess the sensitivity of ks was also used to investigate the effects of the eddy viscosity coefficient, a. A node spacing of 25 m and ks = 0.04 m were used for all the runs. The default value for e is set to 0.5 in River2D and the reasonable range of values given for this coefficient is between 0.2 and 1.0 (Steffler, 2000). Therefore, values of 0.2, 0.5, 0.8 and 1.0 were used to determine the model's sensitivity to this parameter. The effects of the eddy viscosity coefficient were found to be insignificant since the average velocity between the different simulations, which were calculated from the 100 points as previously described, only varied in the fourth decimal place and the average depth only varied in the third decimal place. Furthermore, the mean velocity and depth errors between ^=1.0 and s = 0.5 were only 0.13% and 0.03%, respectively. The spatial sensitivity of this parameter was not examined as its overall influence is of little consequence. Based on this analysis, the default value of e = 0.5 was used for all subsequent River2D simulations. 5.3 Computational Node Spacing Ideally the hydrodynamic simulations should produce a solution which is independent of the computational node spacing. Lane and Richards (1998) investigated the effects of different grid densities in the application of the 2D model S T R E M R to a proglacial braided reach in Switzerland and found a point spacing of 0.05 m would yield grid-independent solutions for that particular modeling application. The average density of bed elevations obtained through surveys was 0.5 m, therefore to avoid modeling an artificially constructed surface resulting from the effects of bed sampling, a grid spacing of 0.10 m was used. In critical areas, such as zones of flow recirculation, a reduced spacing was used. 51 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Previous studies involving River2D have used coarser grid spacings. For example, Waddle et al. (2000) created a mesh with an average node spacing of approximately 4 m with a finer spacing used near channel boundaries and coarser spacing on in-stream islands. They used 3,302 nodes to model a reach that was 315 m long by 35-50 m wide (area « 13,400 m2). The distributed topographic survey data was collected at an average overall spacing of approximately 6 m. Lacey (2001) created a DEM with a 0.5 m spacing, which was similar to the distances between surveyed points, for a side channel of the Chilliwack River. The study reach was 320 m long by approximately 20 m wide on average (area « 6,400 m2), and the finite element mesh had a spacing of 1 m. Due to computational restrictions, the reach was divided into upper and lower sections which were modeled separately. For the Fraser River Gravel Reach, a study section for 2D modeling purposes may be 5 km long by 500 m wide, or 2,500,000 m 2. In many locations this is likely a conservative estimate since the main channel, side channels and vegetated islands occupy a much greater distance in certain areas and also a portion of the floodplains would have to be included in the finite element mesh for modeling higher flows. Since any potential study area in the Fraser River Gravel Reach would be orders of magnitude larger than the previous study reaches modeled with River2D, computational restrictions make it impossible to have node spacings as fine as was used in prior studies. A practical limit on the number of nodes for a River2D finite element mesh is in the order of 10,000 (Steffler, 2000), although this limit is dependent on the computer's available memory and processing capabilities. For the 5 km by 500 m reach, an average node spacing of 16 m would result in 9,766 nodes. The nodes could be non-uniformly distributed throughout the modeled region so more closely spaced nodes are used in areas of high interest or flow variation. A sensitivity analysis on River2D node spacing was conducted using the same 2 km study reach (area « 1,876,000 m2) downstream of the Agassiz-Rosedale 52 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Bridge that was used to investigate the effects of ks and e. Similar boundary conditions were also used. The mesh generation program that was developed for use with River2D was used to create several meshes with varying node spacing. Uniform node spacing was used for each mesh to allow for easier comparisons and each node was assigned a ks value of 0.04 m. The results of the different simulations, including the time in minutes required to reach steady-state conditions for each run, are shown in Table 5.2. The simulation times, which are relatively low given the number of nodes, were achieved using a moderately fast P C with a Windows 2000 operating system, AMD Athlon 700 MHz processor and 256 MB RAM. Initially only 128 MB RAM was installed on the computer which ran under a Windows 98 operating system at the time. Doubling the memory allowed a greater number of nodes to be defined in any given mesh and also decreased simulation times. Changing to a Windows 2000 operating system vastly improved the model's stability as the program crashed far less frequently than before. Table 5.2: Node spacing sensitivity analysis Node spacing (m) # of nodes Simulation time (minutes) 10 >22,000 * n/a 12 15,951 170 25 3,849 8 30 2,735 4 35 2,043 3 40 1,606 2 50 1,067 1 100 325 0.5 * - mesh generation program crashed during triangulation of nodes The finest node spacing that could be used for this reach was 12 m as the mesh generation program would crash while attempting a triangulation for finer 53 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach uniform node spacings. The 12 m spacing simulation was used as the base case since the finest possible discretization should provide the most accurate geometric representation of the bed. The mean velocity and depth errors of the 100 points used for comparison for each simulation relative to the base case is plotted as a function of the node spacing in Figures 5.5 and 5.6, respectively. Figure 5.5: Mean average velocity error versus node spacing 54 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 2.0 -0.5 ~\ , , , , 1 0 20 40 60 80 100 Node spacing (m) Figure 5.6: Mean average depth error versus node spacing Figure 5.5 indicates that the mean average velocity changes by less than 1% when using a uniform node spacing of up to approximately 35 m. Coarser node spacings result in a decrease in model accuracy, as the mean velocity error is almost - 8 % for a uniform node spacing of 100 m. Node spacing has less of an effect on depth than velocity since the mean depth error is less than 2% for even the coarsest mesh considered in this analysis. These results would suggest that a general rule of thumb for applying River2D to a large river system such as the Fraser River Gravel Reach would be to use uniform node spacings less than or equal to about 35 m, although the finest possible mesh should be used given the computational limitations. While maintaining a node spacing less than or equal to 35 m may yield decent results on a reach-wide scale in most 2D modeling applications for the Fraser River Gravel Reach, such a relatively coarse mesh could lead to localized regions where the errors are higher than a few percent. For example, in areas where 55 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach shallow flow passes over an isolated mesh node, excessively high local Froude numbers may occur. These anomalies are an artifact of the calculation and require refinement of the mesh in those areas to avoid high velocities (Waddle et al., 2000). The solution method in River2D is based on a weighted average process where the nodal value at any particular point reflects the surrounding conditions. Where there are large variations in the flow, such as near a sloping bank or in the vicinity of an isolated bar in a stream, a finer node spacing is often necessary to minimize discretization errors (Steffler, 2000). Although a finer mesh will usually mean increased accuracy of model results, computational resources will often restrict how close nodes can be placed to one another, especially when modeling larger rivers. In general, a greater emphasis should be placed upon the quality of the survey data and how it is used to represent the channel bathymetry, which is typically done through the creation of a quality digital elevation model (DEM). The development of the DEM and River2D model for the Agassiz-Rosedale reach is described in the next section. 56 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 6. MODEL D E V E L O P M E N T AND CALIBRATION: AGASSIZ -ROSEDALE REACH 6.1 Agassiz-Rosedale Reach Overview The floodplain topography and channel bathymetry survey information, described in Section 4.0, were combined into a single data set containing transects approximately 200 m apart covering the dyked corridor of the gravel reach. This data was used to create the DEM's for selected reaches in the lower Fraser River. The first study reach, shown in Figure 6.1, selected for 2D modeling was a 4.5 km section centered at the Agassiz-Rosedale Bridge. Photos of the Agassiz-Rosedale Bridge looking upstream the Fraser River at flows of 6,700 and 650 m 3/s are shown in Figures 6.2 and 6.3, respectively. This part of the Fraser River was selected as the first application reach for River2D for several reasons. The first is that it presents fewer complications in terms of 2D modeling compared to more complex downstream reaches. In addition, there are several hydraulic and morphological features that make this reach attractive for 2D hydrodynamic modeling investigations such as the bend in the river at the Agassiz-Rosedale Bridge which creates a superelevation of the water surface. There are current gravel extraction and bank erosion issues in this reach that can be investigated with a 2D model as well. Some of the hydraulic characteristics at the bridge section are listed in Table 3.1, which indicates the width of the river is approximately 500 m at the bridge, although at higher flows large portions of the floodplain are occupied by floodwaters. On the north floodplain, the Kent A and Kent B dykes are set back much further from the north bank than the Chilliwack dyke is from the south bank. In certain locations, the Kent A dyke is as far as 1,200 m from the river's edge. The upstream and downstream boundaries of the modeled reach are shown in Figure 6.1, which indicates that there are two distinct flow branches at the upstream end. The aerial photograph in Figure 6.1 was taken in 1996 at a relatively low flow so many of the 57 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach gravel bars are exposed but the subsequent annual freshets have modified the bars to some degree since then. There are also two vegetated islands in the Agassiz-Rosedale reach which create side channels that become inundated at higher flows. Three water level gauges are also identified in Figure 6.1 as these were used during calibration to compare simulated water levels to measured values. 58 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Figure 6.2: Fraser R.: Agassiz-Rosedale Bridge looking upstream at Q = 6,700 m 3/s Figure 6.3: Fraser R.: Agassiz-Rosedale Bridge looking upstream at Q = 650 m 3/s 59 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 6.2 DEM Generation 6.2.1 General A digital elevation model, or DEM, is a three-dimensional surface generated from topographic and/or bathymetric survey data using digital terrain modeling software. The goal is to obtain an accurate and continuous representation of the channel bed and floodplain surface given a finite number of survey points, which in this case were collected in a series of transects approximately 200 m apart. In 2D hydrodynamic modeling applications, x, y and z coordinate data can be extracted at user-specified intervals from the DEM for use in the design of a finite element mesh. The creation of a digital elevation model involves the application of some type of interpolation algorithm which "fills in" the gaps between surveyed points by estimating grid point elevations. Mathematical and image methods are two general approaches that are commonly used to represent a surface. The mathematical methods involve fitting a continuous three-dimensional function to survey points that allows the elevation at any other point to be determined by evaluating the function at that point. Basically the elevations, z, are expressed as a function of horizontal coordinates, x and y. Alternatively, image methods explicitly give the elevation at some set of points with no functional dependence on horizontal coordinates. There are numerous examples of mathematical interpolation methods, two of which are Lagrange and cubic spline interpolation. The Lagrange method of interpolation will find a polynomial of degree N through A/+1 survey points which allows a curve to be plotted through the sampled points. A drawback to this method is that the interpolated curve can begin to oscillate between the sample points and deviate significantly from a linear interpolation. Cubic spline interpolation fits a piecewise curve through the sampled points and ensures that the first and second derivatives of each function are equal at the intersection of two pieces of the fitted curves. This helps to ensure that the pieces are smooth and continuous with no 60 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach breaks or sudden changes throughout the curve resulting in a closer approximation to a linear interpolation. The image method can employ a regular grid or elevation matrix, which has the difficulty of matching the grid to the terrain complexity, or can utilize an irregular grid based on a triangulated irregular network (TIN). In a TIN, the surface is represented as a sheet of edge-connected triangular elements based on Delauney triangulation of irregularly-spaced control points. Some advantages of a TIN include the ability to accommodate breaklines such as ridges along triangle sides and also the capability to vary the size and shape of the triangles according to terrain features. Other interpolation techniques can be derived from the field of geostatistics. Geostatistics is a collection of statistical methods that were traditionally used in geosciences and describe spatial autocorrelation among sample data and use it in various types of spatial models. Kriging is a method of interpolation used in geostatistics that utilizes a variogram to express spatial variation while minimizing the error of predicted values based on the spatial distribution of estimated values. A variogram is a three-dimensional plot of the semivariances as a function of distance from a point and the semivariance is the measure of degree of spatial dependence between survey points. The smaller the distance between points, the smaller the semivariance and larger distances lead to a larger semivariance. Essentially, variograms measure how quickly things change on average and are strongly dependent on direction. Two main kriging methods are point, or ordinary kriging, and block kriging. The difference is that in point kriging, the value of a point is estimated from a set of nearby sample values using weighting factors while in block kriging the value of a block is estimated instead. Since the variogram changes with distance, the weights depend on the known sample distribution and for ordinary kriging the sum of the weights is equal to one. A detailed review of geostatisitcal methods, and specifically 61 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach kriging, is beyond the scope of this thesis. The reader is referred to Isaaks and Srivastava (1989) and Journel (1989) for further information. 6.2.2 Agass iz -Roseda le Reach D E M The digital terrain modeling software, Surfer 7.0 by Golden Software, Inc. was used to generate all DEM's for the Fraser River Gravel Reach in this study. Surfer offers several different gridding techniques to convert irregularly spaced xyz survey data to a regularly spaced, rectangular array of z values. Some of the methods include kriging, radial basis function, nearest neighbour, triangulation with linear interpolation and polynomial regression. Each of these gridding methods has several specific user-controlled options to further modify and tailor the interpolation process to the particular application. Given the numerous number of gridding possibilities available, kriging and the radial basis function were focused upon to generate a DEM because these two methods generally tend to produce the best quality grid data, including cases when larger data sets (i.e. > 1,000 observations) are being analyzed (Golden Software, Inc., 1999). When applying either kriging or the radial basis function gridding techniques to the survey data it was determined through experimentation that anisotropy, which refers to the preference of data points to have a higher degree of continuity in one particular direction, was significant. Generally, anisotropy can be described using an analogy with an ellipse which helps to represent the anisotropy ratio and angle. The anisotropy ratio is the relative weighting applied in the interpolation process and is simply the maximum range divided by the minimum range while the anisotropy angle defines the direction of the major axis in degrees. The preferred orientation of the survey points was found to be in the streamwise direction. When anisotropy was not specified (i.e. ratio = 1.0), several "holes" would appear in the channel bed. These holes could be clearly viewed when contour plots of the gridded data were created in Surfer. The occurrence of these holes, which is probably mainly due to the relatively large spacing of 200 m between surveyed transects, was greatly 62 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach reduced by incorporating mild anisotropy ratios of 2.0 and setting the anisotropy angle to coincide with the channel alignment in the streamwise direction. Higher anisotropy ratios were investigated however they did not offer significant improvement over milder ratios. Additional problems in the DEM were noticeable on the banks of the Fraser River where an undulating pattern would develop rather than linear smooth features. The inclusion of anisotropy ratios partially helped to improve this phenomenon although the definition of breaklines along channel banks further improved the results. Breaklines are used to define and preserve linear features in a channel such as banks, thalwegs or other breaks in the slope. When the gridding algorithm encounters a breakline, it calculates the z value of the nearest point along the breakline and uses that value in conjunction with nearby survey points to compute grid node values. If a grid point lies on a breakline, then the value of the breakline takes precedence over that point. For the Agassiz-Rosedale reach, breaklines were defined for the north and south banks as well as the thalwegs by plotting individual cross-sections in Excel using the bathymetric survey data and identifying the x, y and z coordinates at the breaks in slope which correspond to those linear channel features. The dykes on the north and south banks were also included as breaklines. There are several different.types of gridding alternatives within the radial basis function group of interpolation methods. Two types that were examined more closely than others were the multiquadratic and thin plate spline functions. The details of these functions will not be discussed in this thesis, however the reader is referred to Carlson and Foley (1991) for a more complete overview. Gridding using the radial basis functions partially resolved the problems associated with holes in the channel bed and undulating banks by trial and error adjustment of gridding options and by incorporating anisotropy and breaklines, however the resulting DEM was still not considered acceptable for use in 2D modeling. 63 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Point, or ordinary kriging, was found to provide an adequate representation of the channel bathymetry using a few very simple options. As described previously, breaklines were used and the anisotropy ratio was set to 2.0 and the anisotropy angle was set to coincide with the river's angle in the streamwise direction. For the Agassiz-Rosedale reach this required the survey data to be gridded in two different sections, as changes in the channel alignment require adjustment of the anisotropy angle. A very simple linear variogram model, which includes the anisotropy details, was used with a slope equal to 1.0. Based on the results of the computational node spacing sensitivity analysis presented in Section 5.3 and the 200 m spacing between surveyed transects, a grid node spacing of 25 m was considered adequate for the Agassiz-Rosedale reach. The gridding process requires the definition of a search ellipse, which defines the local group of survey points to consider when interpolating each grid node. The option of using all the data for each interpolation is also available in Surfer but was not used. Since the transects are approximately 200 m apart, the major axis of the search ellipse was set to 250 m and the minor axis made equal to 125 m. This search ellipse would ensure that the survey points along the transects directly upstream and downstream of each grid node would be included in each interpolation. The relative dimensions of the search ellipse (major axis + minor axis = 2.0) were purposely made to coincide with the anisotropy ratio. Similarly, the angle of the search ellipse was also set equal to the anisotropy angle. The gridding process in Surfer produces rectangular arrays of data which includes areas outside the study reach. Since additional points unnecessarily increase computation time, the extraneous portions of the grid files were blanked so only the study reach and portions of the floodplain contained within the dykes were included. A contour plot of the grid data is shown in Figure 6.4. An interesting feature of the bathymetry is a large scour hole directly upstream of the bridge in the thalweg. 64 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 5453000 5452500 5452000 >• 5451500 5451000 5450500 5450000 Agassiz-Rosedale Bridge LEGEND - Dykes — Banklines/lslands Contour interval = 1 m 587000 587500 588000 588500 589000 589500 590000 590500 591000 591500 X J m ) Figure 6.4: Agassiz-Rosedale reach contour plot A three-dimensional surface plot of the grid data is shown in Figures 6.5 and 6.6. Surfer does not allow surface plots from different grid files to be combined therefore the upstream and downstream halves of the reach are shown separately for illustrative purposes. The dykes and constant elevation lines of 5, 10 and 15 m are shown in color and the light blue arrows indicate the flow directions into and out of each section. The piers for the Agassiz-Rosedale Bridge, which are approximately 10 m long by 4.4 m wide, were not included in the DEM since an individual pier is smaller than the grid size used (25 m by 25 m) to represent the channel surface. Furthermore, River2D would likely not be able to compute the hydraulic effects of the bridge piers in sufficient detail given the 25 m resolution. 65 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 66 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Figure 6.6: Surface plot of the lower half of Agassiz-Rosedale reach (looking upstream) 6.3 River2D Model Development The x, y and z coordinates were extracted from the DEM at the same 25 m intervals as the grid data. Grid data that was blanked in Surfer was deleted in Excel (the x and y coordinates remain for blanked regions), leaving 16,036 data points representing the study area. The data was organized in Excel using five columns with one node per each row, according to the formatting requirements of a River2D .bed file. The first column represented the node number, which ranged from 1 to 16,036, while the second, third and fourth columns represented the x, y and z data, respectively. The final column represented the ks value, which was set to 0.04 m, for each node. The ks values could not be varied according to channel, vegetated 67 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach island and floodplain areas at this stage since Surfer does not have the capability to link additional attributes, other than a z value, to individual grid nodes. River2D_Bed, which is the bed topography file editor developed for use with River2D, allows the roughness to be defined by region using a graphical approach; however, this is a crude feature and was not considered necessary at this point since the creation of a "fully calibrated" 2D model is not the primary goal of this study. It was assumed that flow over floodplains would only be significant at discharges much greater than the mean annual flood and therefore the exact ks specification in non-channel areas would not be critical for the simulation of discharges considerably lower than the design flood. The .bed file was opened in R2D_Mesh, which is the mesh generation program for River2D. The first step in the mesh program is to draw the external computational boundary around the region to be modeled. Next, the inflow and outflow sections are specified and a discharge at the upstream section and fixed water level at the downstream section are assigned. River2D does not allow non-uniform flow or water level distributions to be specified across inflow or outflow sections so it is important to locate these boundaries some distance away from the actual areas of interest to account for any inherent errors in the lateral distribution of boundary conditions. For the Agassiz-Rosedale reach there are two distinct inflow channels. Through experimentation (i.e. by creating and running a model which included several kilometers of river upstream of the actual inflow section) it was determined that the flow split between the two channels is approximately 0.53*Q in the north arm and 0.47*Q in the south arm. This flow split seemed to be consistent over a wide range of flows and was therefore used for all simulations in this reach. Since the two inflow sections are adjacent to one another, River2D adds up the discharge in each arm and then redistributes the total flow across the entire inflow section which makes the initial flow assignments in the north and south arms trivial. To avoid this problem, a very short no-flow section was specified between the two inflow sections which forces River2D to treat each inflow section separately and not sum and redistribute all flows into the reach. 68 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach The next step in the mesh program is to specify the boundary node spacing, which is followed by filling the interior area of the computational boundary with nodes. For the Agassiz-Rosedale reach, uniform node spacings of either 25 m or 30 m were used. For lower discharges (i.e. Q > mean annual flood), a 25 m node spacing was used while the simulation of discharges larger than the mean annual flood required a 30 m spacing. This distinction arose since the maximum node limit of approximately 10,000 has to be maintained (see Section 5.3) while including a greater extent of floodplain area in the modeled domain when simulating higher flows. Once the modeled region has been filled with nodes, a triangulation is performed which creates a finite element mesh consisting of triangular elements. The quality of the mesh can be tracked by a quality index, QI, which should be between 0.2 and 0.5 for an acceptable mesh (Steffler, 2000). After triangulation, the QI is typically quite low, however values of greater than 0.3 can be attained after smoothing the mesh a few times. Smoothing makes the triangles more regular in shape which creates a more gradual transition between different sized triangles. At this stage, the mesh and all boundary conditions can be saved in a .cdg file, which is the file format required by River2D. The "Save As River2D Input File" command automatically creates the .cdg file with the necessary formatting. Prior to the creation of this file, the program prompts the user for an initial inflow water surface elevation estimate. The definition of breaklines is possible and recommended in the finite element mesh to preserve the linear features of the channel topography. However, if breaklines are defined in the mesh, the time required to perform a triangulation increases substantially. For example, the triangulation of 10,000 nodes may take in the order of 10 minutes using an Athlon 700 MHz computer with 256 MB RAM when no breaklines are defined. However when only two sets of breaklines (e.g. north and south dykes) are included for the same problem, the triangulation process is unsuccessful as the computer crashes after an excessively long period of time. Obviously including breaklines for all linear channel features such as the banks and 69 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach thalwegs would further compound this problem. Therefore, no breaklines were included in the finite element mesh in this study. The .cdg file contains a line which controls whether a dynamic or diffusion wave solution will be used. The program defaults to the dynamic wave approximation, however the user can manually change this to a diffusion wave solution using a text editor. Millar and Barua (1999) showed that there is no significant attenuation or lag in the timing of the flood peak between Hope and Mission. This would suggest steady-state conditions through the reach and therefore, a kinematic wave approximation could be used to model the flow. However, the braided reach is relatively flat and comprised of multiple flow paths which would require an unsteady simulation and therefore either a diffusion or dynamic wave approximation would be necessary. In this study, the difference between the dynamic and diffusion wave approximations as applied to the Fraser River Gravel Reach was not investigated in detail. Section 2.1.2 provides background and discussion on this topic, however future studies should involve the investigation of the sensitivity of River2D results to the choice between the two equations for the gravel reach. The simulations were run in River2D until steady state conditions were achieved. The convergence of the solution towards steady state was tracked by observing the "Net Outflow" and "Solution Change". The "Net Outflow" should approach zero at steady state while the "Solution Change" should become sufficiently small (« 0.00001). The model results were extracted and post-processing and further analysis was done using Excel and Surfer. 6.4 Calibration In order to investigate potential applications of the 2D model in the Fraser River Gravel Reach, the model had to be calibrated to some degree to ensure that the simulation results were reliable. The process and data requirements of a 2D 70 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach hydrodynamic model calibration and verification were examined given the available data. Two types of calibration data were considered. The first was measured water levels from the 1999 freshet and the second was A D C P data collected during August of 1999. Since the goal was not to obtain a fully calibrated model but rather to investigate 2D modeling processes and possible applications, less attention was given to fine tuning model parameters, especially on a local scale. 6.4.1 1999 Freshet The peak flow during the 1999 freshet occurred on June 23, 1999 and was measured to be 11,100 m 3/s at Hope. Measured water levels at the three gauges (Kent 5, 22 and Chwk 2) shown in Figure 6.1 were available to compare to model results. Since there is basically no attenuation in the flood peak between Hope and Mission (Millar and Barua, 1999), a discharge of 11,100 m 3/s was used as the upstream boundary condition. The flow was split between the north and south inflow arms as described in the previous section. The downstream water level for this discharge was obtained from the MIKE11 profiles generated in the hydraulic modeling study conducted by UMA Engineering Ltd. (2000a). The goal was to match the simulated water surface elevations to the measured values at the three gauge locations as closely as possible. A range of ks values were used from 0.02 m up to 0.28 m, however a global value of ks = 0.04 m was found to give the best results. Table 6.1 displays the results of this part of the calibration. Shown are the observed and simulated water levels as well as the difference between them. In addition, the MIKE 11 results (UMA, 2000) are also shown for comparative purposes. 71 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Table 6.1: 1999 freshet calibration results in the Agassiz-Rosedale reach Gauge # MIKE 11 Observed River2D Difference (m) (m) (m) (m) 22 17.22 17.00 17.02 +0.02 Kent 5 17.20 16.86 16.86 0.0 Chwk 2 16.75 16.78 16.88 +0.10 The River2D results are very good showing a maximum error of 0.10 m for the three gauges. Although these results are encouraging, additional measured water levels would be useful to further assess the model's performance since it is difficult to consider this a thorough calibration based on the comparison of only three water surface elevations. 6.4.2 A D C P Data The A D C P survey through the Agassiz-Rosedale reach was conducted on August 9-10, 1999. The discharge at Hope during this time was 6100 m 3/s and therefore this flow was used as the upstream boundary condition. The downstream water level was not known for the simulated discharge; therefore the depth-discharge outflow condition (Equation 2.13) was used to obtain an estimate of the water surface elevation at the outflow section, which was approximately 14.8 m. Applying this outflow condition requires the Manning's n value for the reach to be known. Manning's n was calculated to be 0.030 for this reach based on the 1D hydraulic geometry at Agassiz, shown in Table 3.1. All nodes were assigned a ks value of 0.04 m for this simulation. The simulation results were compared with the A D C P data and no adjustment of model parameters were made to obtain a better match. Therefore, in strict terms it is probably incorrect to refer to this part of the study as "calibration". Figure 6.7 shows 72 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach the location of the A D C P transects as red lines in the Agassiz-Rosedale reach that were used to compare to model results. Figure 6.7: A D C P transect locations in the Agassiz-Rosedale reach Since the A D C P data was based on single pings, which means no averaging was done in the field, spatial averaging was performed in Excel. There was one measurement made approximately every meter so every ten values of velocities and depths along each transect were averaged to reduce errors associated with individual point measurements. The simulated depths and velocities were extracted at the x, y coordinates of the averaged A D C P data. Figures 6.8 and 6.9 show scatter plots of the simulated versus measured depths and velocities, respectively. Also shown on the scatter plots are the 1:1 lines. 73 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach ADCP Depth (m) Figure 6 . 8 : Simulated depth versus A D C P depth for the Agassiz-Rosedale reach 0.0 0.5 1.0 1.5 2.0 2.5 3.0 ADCP Velocity (m/s) Figure 6 . 9 : Simulated velocity versus A D C P velocity for the Agassiz-Rosedale reach 74 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach The mean error (ME) for the velocity and depth comparisons were calculated and are shown in Figures 6.8 and 6.9. This parameter, which indicates the dispersion of the data about the 1:1 line, was determined as follows: z mod; UADCPi ME = u 100% ADCPi J [6.1] In the above equation the modeled and A D C P velocities are represented by Umod and UADCP, respectively and N is the total number of points included in the calculation. The depth error was calculated similarly as the velocity errors. The mean depth error is very low at -2.6% while more variability is seen in the velocity comparisons (ME = 8.8%). All the averaged A D C P depth data was used for the comparisons however two sections of the averaged A D C P velocity data were filtered. The A D C P velocities were much lower than model results in transect 19B and the portion of transect 17 that crosses the gravel bar shown in Figure 6.7. The depths were very shallow in these two areas and the survey notes even indicated that due to shallow depths along transect 19B, the discharge estimate, which is based on measured velocities, may not be valid. The remaining data shows a decent correlation between measured and simulated velocities. The A D C P data was also used to compare the flow splits between the main and side channel in transects 18A/B and 19A/B, shown in Table 6.2. The MIKE 11 results are also shown for comparative purposes. The relative proportion of the flow distribution between adjacent channels, which is probably more relevant, is shown rather than absolute values. 75 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Table 6.2: Flow splits (%) in the Agassiz-Rosedale reach Transect A D C P MIKE 11 River2D 18A 83.5 85.8 83.4 18B 16.5 14.2 16.6 19A 96.8 91.4 93.9 19B 3.2 8.6 6.1 The results from both MIKE 11 and River2D are very comparable with the A D C P values. The River2D results are basically identical to the A D C P estimates for transects 18A/B and are within 3% for transects 19A/B. The water levels measured along each line during the A D C P survey using a differential G P S and a high speed boat are compared to the model results in Table 6.3. The measured water levels were given as constant values along each line while the River2D results show some lateral variation in the water surface elevations. Therefore, the simulated water levels across each transect were averaged so they could be directly compared with the measured values. Table 6.3: Comparison of water levels measured during the A D C P survey in the Agassiz-Rosedale reach Transect Measured Simulated Difference (m) (m) (m) 17 14.91 15.19 +0.28 18A 15.11 15.37 +0.26 18B 15.11 15.28 +0.17 19A 15.6 15.67 +0.07 19B 15.6 15.75 +0.15 76 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach The River2D simulated water levels are consistently higher than the measured values. Since the accuracy of the measured water levels is affected by superelevation, it is difficult to assess the reliability of the measured values. In addition, the survey notes indicated that the water level at line 19B was actually measured 500 m away. Therefore, this.comparison of water levels does not offer conclusive information on the accuracy of the Agassiz-Rosedale River2D model. 6.4.3 Aerial Photograph An alternate method of calibrating 2D models was attempted which involved superimposing model results over air photos of the study area. This method would allow the lateral or two-dimensional extent of the water surface to be compared to the wetted area shown on the photos. Of course the simulation would have to be run for the discharge that occurred when the photos were taken. Black and white aerial photographs of the gravel reach were taken on March 7, 2001 when the discharge at Hope was approximately 500 m 3/s. River2D was used to obtain depth and velocity contours for the Agassiz-Rosedale reach for this discharge. The contour plots were superimposed onto the aerial photographs and are shown in Figures 6.10 and 6.11. Due to differences in scale, the contour plots had to be manipulated in order to get a satisfactory match of the water's edge with the air photos. Given the level of manipulation required to obtain the match shown in Figures 6.10 and 6.11, it is difficult to refer to this process as calibration. The concept is still potentially useful in 2D model calibration and warrants further investigation in future studies. 77 The contour plots indicate some interesting features in this reach. Figure 6.10 shows the deep scour hole just upstream of the Agassiz-Rosedale Bridge where the depths are in excess of 12 m for the low discharge of 500 m 3/s. The velocity contours indicate that velocities in the thalweg upstream of the bridge are much less than 1 m/s. The higher velocities towards the downstream part of the reach are likely an anomaly due to boundary condition effects at the outflow section. The contour plots show that the side channel to the south of the main channel in the downstream half of the reach is dry at 500 m 3/s although the air photo clearly indicates that there is water present. This discrepancy could be due to the relatively 78 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach short length of river modeled in this simulation. It is possible that sediment buildup at the mouth of the side channel has restricted flows at lower discharges and that most of the water in the side channel is actually backwatered from downstream of the island/gravel bar (Big Bar) which divides the river at this location. The modeled reach in River2D does not include the downstream portion of Big Bar therefore no water can be backwatered into the side channel in the simulation. 6.5 Summary In this section of the report an overview of the Agassiz-Rosedale reach was provided. General background information on the development of DEM's was discussed and the details regarding the Agassiz-Rosedale DEM generated in this study were presented. The development of the River2D model for this reach, including assumptions and problems which had to be overcome, were described as well. The comparison of model results to measured water levels and A D C P data indicates that an adequate level of calibration has been achieved for the purposes of this study and therefore further modeling scenarios can be investigated. Several applications of River2D in the Agassiz-Rosedale reach are presented, in the following section. 7 9 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 7. APPLICATIONS: AGASSIZ -ROSEDALE REACH 7.1 Superelevation at Agassiz-Rosedale Bridge The superelevation of the water surface at the Agassiz-Rosedale Bridge is evident by the measured levels for the 1999 freshet peak flow at gauge 22 and Kent 5, shown in Table 6.1. The measured values indicate that the superelevation at the bridge section is 0.14 m however MIKE 11 simulated a lateral variation in water level of only 0.02 m (UMA Engineering Ltd., 2000). This low value is to be expected as 1D models assume a constant water surface across each section. Equation [7.1], shown below, can be used to estimate the superelevation (AH) at the outer bend at a given location, based on the radius of curvature (r), depth-averaged velocity (U), width of river (W) and acceleration due to gravity (g). A velocity coefficient (a) is also incorporated into equation [7.1] and can be taken as approximately 1.05. all2 AH = ^—W [7.1] gr At the Agassiz-Rosedale Bridge, the width and radius of curvature of the Fraser River are approximately 500 m and 3,000 m, respectively. The average velocity in the main channel was approximately 3 m/s during the 1999 freshet, based on the River2D simulations. Substituting these values into equation [7.1] provides an estimate of 0.16 m for the superelevation at the Agassiz-Rosedale Bridge, which is very comparable to the difference between observed water levels at gauges 22 and Kent 5. Based on the modeled water levels at the two gauge locations at either end of the bridge, River2D simulates a superelevation of 0.16 m for the 1999 freshet. This simulated water level difference is identical to the value calculated with equation [7.1] and within 2 cm of the observed difference. This example illustrates 80 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach the usefulness of 2D models at estimating the superelevation of water surfaces around river bends. 7.2 Mean Annual F lood The mean annual flood at Agassiz, which is 8,760 m 3/s, was selected as a discharge to be used to provide an example of the hydraulic parameters that can be computed with 2D models and how they can be presented graphically. A mesh was created with a total of 8,103 nodes, which were spaced at 30 m intervals. All nodes were assigned a ks value of 0.04 m. The upstream constant discharge boundary condition was 8,760 m 3/s while the downstream water level (15.7 m) was estimated with the depth-discharge boundary condition described in Section 2.1.5. The simulation required 137 iterations to reach steady state conditions. The following parameters at each node were extracted from the River2D simulation and saved to a .csv file format which can be opened in Excel for further analysis: > bed elevation, z > x discharge intensity, qx > bed roughness, ks > y discharge intensity, qy > depth, H > Froude number, Fr Using the equations presented in Section 2.0, several additional parameters were calculated at each node in an Excel spreadsheet, including the u and v velocity components, overall velocity magnitude (LO, friction slopes in the x and y directions (Su and Sfy), mean bed shear stress (r0) and the discharge per unit width (q). r 0 and q were determined from taking the square root of the sum of the squares of their x and y components. The water surface elevation was simply calculated as z + H at each node, although River2D also provides these computed values. 81 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Surfer was used to create grid files and contour plots of the model results, which allowed the computed values across the bridge section to be digitized. The digitized values were used to calculate the average velocity, depth and Froude number. A comparison was made between the 1D hydraulic parameters at the Agassiz-Rosedale Bridge given by Church and McLean (1994) and the River2D results. Table 7.1 shows the cross-section averaged parameters computed from the River2D results and the previously published values. Table 7.1: Comparison of 1D hydraulic parameters at the Agassiz-Rosedale Bridge for the mean annual flood Parameter Church and McLean (1994) values Calculated from River2D results % Difference Average velocity (m/s) 2.6 2.52 -3.1% Average depth (m) 6.6 6.80 +3.0% Average Froude # 0.32 0.31 -3.6% The calculated values are within 3.0% to 3.6% of the published figures which further indicates that the model is adequately simulating the hydraulics through this reach. One of the main advantages of a 2D hydrodynamic model is the ability to determine the spatial variation of hydraulic parameters and also to visualize the direction of flow at any point. To provide an example of the possible plots that can be created using 2D model results, Surfer was used to create contour maps of several parameters for the mean annual flood at Agassiz. Figures 7.1 to 7.7 show contour plots of the bed elevation, depth, water surface elevation, velocity, discharge per unit width, Froude number and the mean bed shear stress, respectively while a plot of the discharge vectors is shown in Figure 7.8. 82 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 587000 587500 588000 588500 589000 589500 590000 590500 591000 XJm) Figure 7.1: Bed elevation contours for the mean annual flood at Agassiz 83 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 84 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 587000 587500 588000 588500 589000 589500 590000 590500 591000 X M Figure 7.5: Discharge per unit width contours for the mean annual flood at Agassiz 587000 587500 588000 588500 589000 589500 590000 590500 591000 XJm) Figure 7.6: Froude number contours for the mean annual flood at Agassiz 85 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 5452500 5452000 S 5451500 >-5451000 5450500 Discharge Vectors • * > 0.1 m A 2 / s 60 m A 2 / s 587000 587500 588000 588500 589000 589500 590000 590500 591000 Xjm) Figure 7.8: Discharge per unit width vectors for the mean annual flood at Agassiz 86 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 7.3 Estimation of Grain Size Spatial Distribution Currently no information is available on the spatial distribution of the d5o grain' sizes and so a uniform value of ks is used throughout the modeled region. Therefore, an analysis was considered that would allow the spatial variation of the d50 of the bed material to be estimated. This analysis was based on an assumption that at some particular flow, all the grains on the wetted surface of the channel bed are at the threshold of entrainment. McLean et al. (1999) have suggested that this threshold condition occurs near flows of 5,000 m 3/s at Agassiz, at which significant gravel transport commences. Unfortunately there are some uncertainties inherent in these assumptions, which prevented any dso distribution contour plots from being produced; however, the methodology of this analysis will still be presented as the potential for future work is possible. The methodology that was proposed was to simulate a discharge of 5000 m 3/s through the Agassiz-Rosedale reach and use the computed velocities and depths to estimate the d50 at each node by equating the dimensionless shear stress, r*, to the critical dimensionless shear stress, TC*. The threshold of particle mobilization is represented by the critical shear stress. Equations [7.2] to [7.4], shown below, were used to derive an expression relating r c*to the depth averaged velocity (U), depth (H) and dso at each computational node. All variables shown in the following equations have been previously defined with the exception of S, which is the specific gravity of the sediment. [7.2] f 12.2// 47 = 2.031og [7.3] V 50 [7.4] 87 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Substituting d 5 0 for ks and rearranging allows equation [7.3] to be expressed as follows: f = log 0.243 r\22H} V ^ 5 0 J [7.5] Equating r*to z-c*and combining equations [7.2], [7.4] and [7.5] yields equation [7.6]. The acceleration due to gravity, g, and the specific gravity of the sediment, S, have been taken as 9.81 m/s 2 and 2.65, respectively. 0 .00187£ / 2 log V ^50 J [7.6] d 50 Successfully applying equation [7.6] requires that the dimensionless critical shear stress, also known as the Shields parameter, be specified at each node. The Shields parameter is often taken as 0.06 for uniform sediment and 0.03 for a natural poorly sorted channel bed with no movement. McLean et al. (1999) found that the Shields parameter can reach values as high as 0.07 to 0.09 during flood conditions at Agassiz, which are almost twice as large as the accepted values. These relatively high values can be explained by the effects of particle imbrication and protrusion (Church etal. , 1998). For the simulated discharge of 5000 m 3/s, a r c*of 0.048, which is a compromise between uniform sediment and a poorly sorted bed, was initially proposed to be applied to all nodes. However, the Shields parameter will likely vary considerably throughout the reach and applying one value globally will produce misleading results. In addition, the original assumption that every particle on the wetted surface is at the threshold of entrainment at a particular flow is highly unlikely 88 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach since protruding grains may become entrained at discharges lower than 5,000 m 3/s while other particles may remain stationary up to considerably larger flows, perhaps due to clustering and shape effects. Another potential difficulty with this analysis is that the computed depths and velocities from the model, which are to be used to calculate individual dso sizes, are based on a uniform ks value, in this case equal to 0.04 m. This problem could be overcome by employing an iterative procedure where a uniform ks value could be used in River2D to provide initial depths and velocities to calculate d50 sizes with equation [7.6]. These calculated grain sizes could then be used to revise the ks values in the River2D model and update the depths and velocities in order to re-apply equation [7.6] to calculate new dso sizes. This process would be repeated until there was little change in the computed grain sizes. This analysis could potentially provide valuable information on the spatial distribution of grains sizes, which is not feasible to obtain through conventional sampling methods due to the river's size. However, to successfully apply the methodology discussed would require information regarding the spatial variability of the Shields parameter. At the very least, the variation of this parameter between gravel bars, riffles and in the thalweg in any proposed reach would be necessary. The assumption that all particles are at the threshold of entrainment at a given flow, i.e. 5000 m 3/s at Agassiz, would require additional investigation also. 7.4 Design Flood One of the main issues in the Fraser River Gravel Reach deals with flood control, which prompted the design flood profile to be recently updated using MIKE11 (UMA Engineering Ltd., 2000). The MIKE11 design flood profiles have been accepted as the criteria for flood protection along the gravel reach. For comparative purposes, River2D was used in this study to compute a water surface profile for the design flood for the Agassiz-Rosedale reach. The design flow of 89 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 17,000 m 3/s at Hope was applied as the upstream boundary condition and the downstream water level of 17.91 m was obtained from the MIKE11 results. A mesh spacing of 30 m was used with an initial value of ks = 0.04 m assigned to all nodes, including those on the floodplain. The water levels along the thalweg were used to generate the flood profile. The resulting design flood profile was considerably lower than the MIKE11 results. At the downstream end, the profiles matched exactly which was expected since MIKE11 provided the downstream boundary condition for River2D, however the difference between the two simulations increased considerably in the upstream direction. For example, at the bridge the difference was almost 0.3 m. While the global ks values of 0.04 m have yielded good results thus far, it is apparent that for such a large discharge, which would result in significant flow over floodplains and islands, a much higher roughness value should be assigned to non-channel areas. Since Surfer does not have the capability of linking ks, or any other parameter besides a single z (bed elevation in this case), to grid nodes, the "Roughness by Region" command in River2D's bed configuration program was used to spatially vary the roughness. This command is very limited as the boundaries of complex polygons tend to become distorted. Therefore, in order to preserve the desired shape of the polygons, i.e. the selected modeled area, polygons had to be limited to five or six sides. This makes it extremely difficult to accurately define the boundaries of islands or floodplain areas with this command. Ideally, a GIS can be used to define different areas within the modeled region and assign ks values locally. The bed configuration program could then be bypassed as x, y, z and ks data could be extracted directly from the GIS. Since no other option was available in this study, the non-channel areas were crudely defined with the "Roughness by Region" command. The north and south floodplains were represented by several four-sided polygons and the boundaries of the two largest islands in the reach were approximated with simple polygons as well. The ks values 90 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach for the channel areas were kept at 0.04 m, however the floodplain and island roughnesses were increased substantially. Through trial and error, a ks value of 2.0 m for all non-channel areas was found to give a fairly decent match between the River2D and MIKE11 design flood profiles, shown in Figure 7.9. Also shown in this figure is the initial River2D generated profile using a /Vsof 0.04 m for non-channel areas. Although the Agassiz-Rosedale model extends almost 2 km upstream of the bridge, the profiles for the upper 1.5 km of the reach are not shown. Despite increasing the floodplain and island roughness, the River2D flood profile still deviates significantly from the MIKE11 profile in this upper region. This difference is possibly due to upstream boundary condition effects in the 2D model as any assumed errors in the lateral distribution of the discharge at the inflow section require some distance to dissipate. Therefore, the design flood profiles in Figure 7.9 are shown for approximately 2.5 km along the thalweg in the Agassiz-Rosedale reach. The chainages correspond to the distances shown in Figure 6.1. 19.0 17.8 J r - - 4 - r - - 4 ~ I 129+000 129+500 130+000 130+500 131+000 131+500 Chainage MIKE11 ks=0.04 m for non-channel areas - — ^ — ks=2.0 m for non-channel areas Figure 7 . 9 : Design flood profile for the Agassiz-Rosedale reach 91 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach While it is possible to obtain fairly accurate water surface profiles over short distances with River2D, the use of a 1D model for this purpose is recommended. It would not be feasible to use a 2D model to generate the design flood profile for the entire gravel reach as it would have to be done in several different simulations and the resulting profiles would have to be spliced together. Conversely, a 1D model can generate a profile for the entire gravel reach in a single simulation. In addition, accurate boundary conditions would be required for each 2D simulation. Without detailed water level information along the river, this information would not be available. If a calibrated 1D model already exists, as for the gravel reach, water level boundary conditions can be obtained from 1D simulation results for use in more localized 2D modeling efforts. This simulation was the first one attempted using relatively high ks values (>0.5 m) combined with very large discharges. It was found that the inclusion of these higher roughness values led to some nonlinear instabilities in the model. As previously described, the convergence of any simulation towards steady state conditions can be tracked with the "Net Outflow" and "Solution Change" values in River2D. When regions of high ks were included in the modeled area, the "Solution Change" underwent frequent fluctuations rather than steadily approach a sufficiently small value. This caused the computational time step to undergo frequent alternating periods of reduction and gradual increase. Eventually the "Net Outflow" did reach a small enough discharge to be acceptable and the simulation was stopped at that point. Examining the .log file after the simulation had been stopped indicated that the instabilities were very localized and restricted to only a few nodes in the modeled area with higher ks values. Overall, the simulation had reached steady state, as the solution was not changing significantly at the vast majority of the nodes. At the few problem nodes, negative water surface elevations or extremely high velocities were not uncommon. In addition, the area surrounding these nodes were usually a combination of adjacent wet and dry regions. 92 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Since these problem areas were very localized, the overall simulation results were still accepted and the results at the nodes causing the instabilities were ignored in further analysis. The exact cause of these instabilities is uncertain, however it is likely that the problem is also related to the coarse node spacing (30 m) used in the finite element mesh. If finer node spacings are specified in the vicinity of the few problem nodes that appear in any particular simulation, then new nonlinear instabilities develop at other floodplain or island nodes that have high ks values. Given the practical limitation on the total number of nodes, it is not possible to use finer node spacings in all areas with high ks values so the localized instabilities cannot be avoided. The simulation results can be accepted as having reached steady state conditions once the net outflow converges to a discharge close to zero. 7.5 Gravel Removal Scenarios Historically there has been several gravel removal operations in the lower Fraser River, however the British Columbia Assets and Land Corporation have regulated these extractions since 1974. Also, the Department of Fisheries and Oceans (DFO) has required permits for any gravel extraction since 1980. McLean and Church (1999) have estimated that the average gravel load transported past the Agassiz-Rosedale Bridge is approximately 200,000 tonnes per year, which corresponds to a volume of about 108,000 m 3 annually, all of which is deposited within the gravel reach. By comparison, an average of 120,000 m 3 of gravel is mined annually from this reach. About 80% of the total gravel extracted from the river is done at two locations in the Minto side channel near Chilliwack. While gravel is a valuable commercial resource some studies have suggested that removal could also help alleviate local flooding risks and possibly reduce bank erosion problems. UMA Engineering Ltd. (2000) used MIKE11 to simulate two gravel extraction scenarios that would increase the river's conveyance and reduce the design flood profile in areas of concern. The total volumes of gravel extracted in these two scenarios were approximately 17,000,000 and 28,000,000 m 3 Such large-scale 93 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach removals would likely have significant effects on the river's morphology and ecology. There are many questions regarding how much gravel can be removed before there are significant morphological and ecological effects in the lower Fraser River. Some of those questions can be partially answered with the aid of a 2D morphological model. The fixed bed model, River2D, can be used to assess the impact on local hydraulics in the river due to varying volumes of gravel removal. In this section, two gravel removal scenarios are considered. The first is a localized gravel removal. (« 130,000 m3) on Powerline Island, upstream of the Agassiz-Rosedale Bridge. The second scenario looks at large-scale gravel scalping of most of Powerline Island, as recommended by Northwest Hydraulic Consultants Ltd. (1999), to reduce bank erosion downstream of the bridge. 7.5.1 Local ized Removal A relatively small rectangular area, slightly larger than 200 m by 100 m, on Powerline Island (Figure 7.10), about 1 km upstream of the Agassiz-Rosedale Bridge was focused upon to determine the effects on river hydraulics of a localized gravel removal. This removal site is located at the edge of a gravel bar, just north of the main channel. Figure 7.10 shows the upstream and downstream limits of the reach used in this simulation. The downstream boundary was taken at the bridge as no effects on river hydraulics were expected to occur downstream of that point. Two River2D model files were created for this section of river, with the only difference between them being that the bed elevations in one of these files were reduced in the proposed removal area to represent the gravel extraction. Surfer was used to modify the z elevations in the proposed removal area in the Agassiz-Rosedale reach grid files. The depths of gravel removal were kept fairly constant at about 5 m and the total volume of gravel removed was 132,000 m 3. A mesh was created with nodes spaced at 25 m intervals, which were all assigned a ks of 0.04 m. The mean annual flood of 8,760 m 3/s was selected as the upstream boundary 94 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach condition and the downstream water level was set to 16.3 m, which was estimated from previous simulations. Both the pre- and post-removal scenarios were simulated with River2D and the results were extracted for comparison. Figure 7.10: Local gravel removal on Powerline Island Pre- and post-removal depth contours are shown in Figures 7.11 and 7.12, respectively while pre- and post-removal velocity contours are shown in Figures 7.13 and 7.14, respectively. Discharge per unit width vectors are shown for the pre- and post-removal conditions in Figures 7.15 and 7.16. The contour and vector plots indicate that this localized gravel removal does not significantly affect the hydraulics in the main channel adjacent to the extraction site or at any other location in the reach. The pre- and post removal contours and vector plots show that the depths, velocities and direction of flow are virtually identical before and after gravel removal everywhere in the reach other than the removal site itself. In the excavated area the depths are about 5 m deeper, which was expected, and the flow forms a whirlpool as water circulates around the perimeter of the pit. The velocity in the center of the excavation is lower than during pre-removal conditions. This analysis suggests that localized gravel removals in the lower Fraser River have negligible impact on the 95 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach flow hydraulics outside of the excavated area. Therefore, such gravel removals are likely not a feasible method to alleviate flooding and erosion risks in this river. 5452000 5451800 5451600 E 5451400 >-54512001 5451000 5450800 589000 589200 589400 589600 589800 590000 590200 590400 590600 590800 591 XJm) 000 16 m 14 m 12 m 10 m 8 m •6 m 4 m 2 m 0 m -2 m Figure 7.11: Localized gravel removal on Powerline Island: Pre-removal depth contours 5452000 5451800 5451600 5450800 R e m o v a l area —J 589000 589200 589400 589600 589800 590000 590200 590400 590600 590800 591000 X(rn) 16m 14m 12 m 10 m 8 m 6 m 4 m : 2 m 0 m •2 m Figure 7.12: Localized gravel removal on Powerline Island: Post-removal depth contours 96 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 97 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 5452000 5451800 5451600 5451000 5450800 Discharge Vectors > 1 m A 2/s 56.5 m A 2/s * is tr" *" «* it & iS * / iS S is *" d *• nT j r * *r* *" i f * • * * 7 , * V * ' * / V S s * 589000 589200 589400 589600 589800 590000 590200 590400 590600 590800 591000 X (na) Figure 7.15: Localized gravel removal on Powerline Island: Pre-removal discharge vectors 5452000 5451800 5451600 Discharge Vectors 5451000 5450800 1 m A 2/s 56 m A 2/s Remova l area 589000 589200 589400 589600 589800 590000 590200 590400 590600 590800 591000 X ( m ) Figure 7.16: Localized gravel removal on Powerline Island: Post-removal discharge vectors 98 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 7.5.2 Large-Scale Removal: Powerline Island Northwest Hydraulic Consul tants Ltd. (1999) prepared a report for the C h e a m Indian Band regarding an erosion problem on the left bank of the Fraser River downstream of the Agass i z -Roseda le Bridge. The river is eroding portions of the left bank at both Ferry Island and Island "B" , shown in Figure 7.17. Cont inued erosion in these areas may threaten the dykes behind the two islands. Figure 7.17: Location of Ferry Island and Island " B " (after Northwest Hydraulic Consul tants Ltd., 1999) Northwest Hydraulic Consul tants Ltd. (1999) suggested that the erosion on the left bank is primarily a result of the growth of the point bar on the downstream end of Powerl ine Island, which causes a progressive shift of the thalweg towards the left bank. A lso , the large gravel stockpile on Powerl ine Island accelerates the growth of the point bar. Al though this is a natural process in a braided river, the shifting of the thalweg increases the angle of attack on the left bank which leads to the erosion problem. 99 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach To reduce the bank erosion, Northwest Hydraulic Consultants Ltd. (1999) proposed that the gravel bars on Powerline Island be scalped so the effective flow area during high flows is increased. Theoretically this would lead to a reduction in velocities and possibly prevent further shifting of the thalweg towards the left bank, which would help to minimize erosion. The proposed removal area, which would extend vertically 2 to 3 m below the low water level, is shown in Figure 7.18 Figure 7.18: Proposed gravel removal area to reduce downstream bank erosion (after Northwest Hydraulic Consultants Ltd., 1999) River2D was used to assess this gravel removal scenario to determine whether it would reduce erosion on the left bank downstream of the bridge. Since the annual peak freshet flow causes most of the erosion damage, the 1999 peak flow of 11,100 m 3/s was used in this analysis. The low water level in the main channel adjacent to Powerline Island was estimated to be approximately 11.75 m by simulating a discharge of 500 m 3/s through the Agassiz-Rosedale reach. As with the previous gravel removal scenario, two River2D files were created, with the only difference being that one had the bed elevations lowered in the proposed gravel removal area. Surfer was used to digitize the limits of the gravel bars and three heavily vegetated areas within the gravel-covered region on Powerline Island. This allowed the coordinates of the gravel bar areas on the island to be isolated since only the bed elevations in those areas were to be lowered. The 100 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach area of gravel removal was 265,000 m 2 and all bed elevations in this region were set to 9 m, which corresponded to an average depth below the low water level of 2.75 m. The average cut at the grid nodes was 4.42 m and the total volume of gravel removed was 1,170,000 m 3. A mesh was created for both the pre- and post-removal cases with nodes spaced 30 m apart and assigned ks values of 0.04 m. The reach limits and boundary conditions used in this analysis were similar to those used in the 1999 freshet calibration simulation, i.e. Q = 11,100 m 3/s at the inflow section and water surface elevation = 16.3 m at the outflow section. The simulation results were extracted and Excel and Surfer were used to perform further calculations, prepare plots and compare results. A contour plot showing the difference in velocities between the pre- and post-removal cases is shown in Figure 7.19. The floodplain portions of the reach have been blanked so only the channel areas are shown. The removal boundary is shown in red and the three interior red lines within the larger boundary represent vegetated areas where no gravel was extracted. Increases in velocities are represented by shades of blue while velocity decreases are shown in green. Velocities are increased in the excavated zone by up to 1 m/s in certain areas. The velocities are increased by up to 0.5 m/s in the northern part of the main channel and part of the side channel near the upstream end of Powerline Island as a portion of the flow diverges into the excavated region due to the increased conveyance capacity. The gravel removal also results in a velocity increase up to 0.5 m/s in the streamlines downstream of the excavation and over portions of Big Bar Island, towards the downstream end of the reach. There is a small side channel, known as Camp Slough, between Ferry Island and Island "B". There has been some erosion at the entrance of this slough and the contour plot indicates that gravel removal on Powerline Island would increase velocities by up to 0.5 m/s in the slough itself and possibly even higher at its entrance. 101 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 5450500 | 587000 587500 588000 588500 589000 589500 590000 590500 Figure 7.19: Powerline Island gravel removal: Contour plot of the change in velocity due to extraction In the remainder of the reach, the gravel removal on Powerline Island results in a noticeable decrease in channel velocities. In the thalweg adjacent to the excavation, velocities are decreased by up to 1 m/s. Velocity decreases up to 0.5 m/s are seen in the remainder of the channel areas adjacent to and downstream of Powerline Island. Discharge per unit width vectors are shown in the channel areas for the pre-and post -removal cases in Figures 7.20 and 7.21, respectively. The vector plots indicate that gravel removal on Powerline Island decreases the maximum unit discharge in the thalweg by only 2.5 m 2/s, which corresponds to a 4.2% decrease. Other than in the excavated area, there is not a significant change in flow direction anywhere in the reach. With a fixed-bed 2D model such as River2D, it is not possible to determine whether this gravel extraction scenario would prevent future shifting of the thalweg towards the left bank. A 2D morphodynamic model would be required to address that question. 102 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 5452500 5452000 E. 5451500 > 5451000 5450500 Discharge Vectors 1 m A 2 /s 59 m A 2/s Left bank Camp Slough 587000 587500 588000 588500 589000 589500 590000 590500 591000 X J m ) Figure 7.20: Powerline Island gravel removal: Pre-removal discharge vectors 5452500 5452000 'E 5451500 >-5451000 5450500 Discharge Vectors 1 m A 2 /s 56.5 m A 2/s Removal boundary shown in black Left bank Camp Slough 587000 587500 588000 588500 589000 589500 590000 590500 591000 X J m ) Figure 7.21: Powerline Island gravel removal: Post-removal discharge vectors 103 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach This analysis does indicate that scalping of the gravel surface on Powerline Island would reduce velocities in the adjacent channel, which would help to reduce the erosion along the left bank of the river in this reach. However, it is unlikely that this reduction in channel velocities alone would alleviate the erosion risk entirely. Furthermore, the modeling results suggest that the proposed gravel removal would increase velocities both at the mouth of and in Camp Slough, which could accelerate existing local erosion of bank material. No significant change in flow direction in the main channel as a result of the excavation could be detected from the vector plots. 7.6 Habitat information at Big Bar The Department of Geography at UBC is currently developing a habitat classification system based on year-round sampling of the abundance, location and habitat offish and invertebrates in the gravel reach (Church et al., 2000). This classification system is comprised of a hierarchy of three levels of classification units: sub-reaches, pool-bar-riffle and bar-edge habitat. The first classification unit divides the river from Hope to Mission into five sub-reaches, which is useful for strategic planning for fisheries management. The intermediate level of classification identifies the pool-bar-riffle sequences in each sub-reach, which can help provide guidance to field studies and the operational management of fisheries and habitat protection in the river (Fraser Basin Council, 2001). The finest classification level, which deals with bar-edge habitat, is based on morphology, substrate, flows and depth. To aid in the bar-edge habitat classification, 2D hydrodynamic models can be used to obtain local depths and velocities around the gravel bars for a range of flows. Although the finest classification level has not yet been applied to any bars in the Agassiz-Rosedale reach, in this study River2D was used to provide an example of the type of data that can be obtained from a 2D model to aid in bar-edge habitat classification. A range of flows was simulated around Big Bar, which is located 104 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach towards the downstream end of the modeled reach (Figure 7.22). The 25 m DEM data was used to create a finite element mesh with nodes approximately 20 m apart. The upstream boundary, which is shown in Figure 7.22, was defined at the bridge section. Figure 7.22: Location of Big Bar The discharges simulated were 700 m 3/s, 1,500 m 3/s, 3,000 m 3/s and 6,000 m 3/s, as requested by the UBC Department of Geography. Downstream water levels for each flow were estimated using the depth-discharge relationship boundary condition. The simulated depths and velocities were extracted from River2D and Surfer was used to interpolate the data to a 10 m grid around Big Bar. The finer resolution is necessary as a 25 m spacing may be too coarse to define some bar-edge habitat features. Contour plots of depths and velocities were prepared in Surfer and are shown for the simulated flows in Figures 7.23 to 7.30. The green lines represent channel boundaries and the brown lines represent the edge of the gravel bar. Depth and velocity contours are only shown for areas on and surrounding Big Bar. 105 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 5452000 5451750 5451500 £ > 5451250 5451000 5450750 14.0 m 3.0 m 2.0 m 1.0 m J 0.0 m 587000 587250 587500 587750 588000 588250 588500 588750 X J m ) Figure 7.23: Big Bar habitat information: Depth contours for 700 m 3 /s 5452000 5451750 5451500 E, 5451250 5451000 5450750 -12.00 m/s 1.50 m/s 1.00 m/s 10.50 m/s •0.00 m/s 587000 587250 587500 587750 588000 588250 588500 588750 XJm) Figure 7.24: Big Bar habitat information: Velocity contours for 700 m 3 /s 106 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 5452000 5451750 5451500 E 5451250 5451000 5450750 5.0 m 4.0 m 3.0 m 2.0 m 1.0 m 0.0 m 587000 587250 587500 587750 588000 588250 588500 588750 X ( m ) Figure 7.25: Big Bar habitat information: Depth contours for 1,500 m /s 5452000 5451750 5451500 E > 5451250 5451000 5450750 2.00 m/s 1.50 m/s 1.00 m/s 0.50 m/s •0 .00 m/s 587000 587250 587500 587750 588000 588250 588500 588750 XJm) Figure 7.26: Big Bar habitat information: Velocity contours for 1,500 m 3/s 107 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 5452000 5451750 5451500 > 5451250 5451000 5450750 7.0 m 6.0 m 5.0 m 4.0 m 3.0 m 2.0 m .0 m 0.0 m 587000 587250 587500 587750 588000 588250 588500 588750 XJm) Figure 7.27: Big Bar habitat information: Depth contours for 3,000 m 3/s 5452000 5451750 5451500 E. >-5451250 5451000 5450750 2.50 m/s 2.00 m/s 1.50 m/s 1.00 m/s 10.50 m/s ' 0 . 00 m/s 587000 587250 587500 587750 588000 588250 588500 588750 XJm) Figure 7.28: Big Bar habitat information: Velocity contours for 3,000 m 3/s 108 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 5452000 5451750 5451500 E, >-5451250 5451000 5450750 8.0 m 7.0 m 6.0 m 5.0 m - - 4 0 m J-M 3.0 m 2.0 m 1.0 m 0.0 m 587000 587250 587500 587750 588000 588250 588500 588750 XJm) Figure 7.29: Big Bar habitat information: Depth contours for 6,000 m 3/s 5452000 5451750 5451500 E >-5451250 5451000 5450750 i3.00 m/s 2.50 m/s 12.00 m/s i 1.50 m/s 1.00 m/s 10.50 m/s •0.00 m/s 587000 587250 587500 587750 588000 588250 588500 588750 XJm) Figure 7.30: Big Bar habitat information: Velocity contours for 6,000 m 3/s 109 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 8. MODEL DEVELOPMENT AND CALIBRATION: HARRISON RIVER CONFLUENCE/MINTO ISLAND REACH 8.1 Harrison River Confluence/Minto Island Reach Overview The second reach that was modeled using River2D was an 8.5 km long section that stretches from upstream of the confluence of the Harrison and Fraser Rivers to downstream of Minto Island (Figure 8.1). The confluence of the two rivers is characterized by complex flow hydraulics. The main channel of the Fraser River flows into the outcropping bedrock at Harrison Knob and gets abruptly redirected in a southerly direction through a 90° bend, which creates strong secondary currents and significant energy losses. The backwater effect not only raises upstream water levels that increase the flooding risk, but also reduces outflows from the Harrison River. The Harrison River drains through a large lake and contributes approximately 7.6% of the peak discharge (design flood) of the Fraser River at Hope, but its drainage basin accounts for only 3.6% of the tributary area above the Fraser River at Hope. Another interesting feature of the Harrison River is its clarity compared to that of the Fraser River. Figure 8.2 shows the abrupt change from the clear waters of the Harrison River to the sediment-laden Fraser River. In the center of this study reach is Minto Island, which is shown in Figure 8.1 and divides the river into two distinct channels. The division of flow into these two channels affects the rate of downstream progression of gravel bars and also bank erosion at critical downstream locations. Dykes on the north and south banks are also identified in Figure 8.1. The Kent D and Nicomen Island dykes provide flood protection for areas on the north side of the river while the Chilliwack, Young Road and Wing dykes protect the land on the south side. The distance from the dyke on one side of the river to the dyke on the opposite bank is almost 2.5 km in some parts of the reach. Four water level gauges are shown in Figure 8.1, which provide measured water levels that can be compared to simulation results. 110 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Figure 8.1: Harrison River confluence/Minto Island reach 111 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Figure 8.2: Photo taken at the downstream end of Harrison River, looking towards the Fraser River This reach presents many more computational obstacles than the Agassiz-Rosedale reach due to a combination of the complex flow hydraulics at the confluence of the Harrison and Fraser Rivers and also due to its size. The practical limitation on the total number of computational nodes in River2D makes obtaining meaningful 2D simulation results for this reach a difficult task. To successfully model the flow through this section of the Fraser River, the reach was divided into two sections, separated by the dashed blue line shown in Figure 8.1. Two inflow sections, for the Fraser and Harrison Rivers, were specified in the upper half of the reach, which has been the focus of several engineering studies by local consultants over the past few years. Some of the recommendations made in those studies to alleviate local flooding risks and erosion concerns will be investigated with River2D. 112 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 8.2 D E M Generation The DEM for this reach was generated using techniques similar to those established for the Agassiz-Rosedale reach. Point kriging was used in Surfer 7.0 to convert the bathymetric and topographic survey data into regularly spaced grid data at 25 m intervals. The same linear variogram model and search ellipse described in Section 6.2.2 was used in this gridding application. The anisotropy ratio was set to 2.0 throughout the study area and the anisotropy angle was varied in relation to the orientation of the channels. Significant changes in the river's alignment required a corresponding adjustment in the anisotropy angle to ensure holes in the channel bottom or other anomalies in the DEM would not develop. Four separate grid files were created corresponding to the following sections of the reach: the upper 3.5 km of the Fraser River, Harrison River, main channel just downstream of the sharp bend and Minto Island including adjacent channels and downstream areas. A contour plot created from all four grid files is shown in Figure 8.3. Breaklines were defined for all dykes and along the edge of Harrison Hill, as shown in Figure 8.3; however they were not defined for the thalwegs or any other banks. This decision was made since the quality of the DEM seemed adequate for preliminary modeling investigations and the task of defining breaklines for this reach would be an extremely time consuming task. In future studies, the proper definition of breaklines for all major breaks in slope is recommended. The contour plot indicates that there are two deep scour holes in the main channel just downstream of the abrupt bend. There is another significant scour hole in the main channel downstream of Minto Island where the thalweg crosses over towards the south bank. 113 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 8.3 River2D Model Development As with the DEM generation, a similar procedure as for the Agassiz-Rosedale reach was used to develop the River2D model for the Harrison River confluence/Minto Island reach. The River2D .bed file was created from the DEM grid node x, y and z data. This data, spaced at 25 m intervals, was extracted from Surfer and opened in Excel where node numbers and a ks value of 0.04 m were added for all points. The resulting .bed file, which had over 38,000 points and took an excessively long period of time to load into the mesh generation program (> 30 minutes), was split into two separate files representing the upper and lower part of the reach as described in Section 8.1. Unlike the Agassiz-Rosedale reach, higher ks values were defined for islands and floodplains in this reach even for the simulation of lower flows (i.e. < mean annual flood). The bed topography file editor was used to define several four and 114 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach f ive-sided polygons over islands and floodplains. The ks va lues in these regions were determined through calibration, descr ibed in Sect ion 8.4.1. The bed files were opened in the mesh generation program where the computational boundary, inflow and outflow sect ions and assoc ia ted boundary condit ions were defined. The upper half of the reach had two inflow sect ions, corresponding to the Fraser and Harrison River inflows, and one fixed water level outflow sect ion at the intermediate reach boundary, whose location is approximately shown in Figure 8.1. The upstream boundary condit ions for the lower half of the reach were obtained from the upper half simulations. Simulated d ischarges in the main channel and Minto channel were defined as upstream boundary condit ions for the lower half model and a single fixed water level was used as the downstream boundary condition. Through trial and error it was determined that the finest uniform node spac ing that could be practically defined in the mesh was 35 m for both the upper and lower sect ions of the reach. This spac ing, which results in a maximum number of nodes of less than 11,000 in both upper and lower halves, is relatively coarse although it should still yield decent simulation results based on the findings of the sensitivity analysis conducted in Sect ion 5.3. It was not possible to define breakl ines in the mesh due to computational restrictions. The .cdg files were created as the final step in the mesh generation program and were subsequent ly opened in River2D. A s before, upon completion of the simulations, model results were extracted for post-processing and analysis using a combination of Exce l and Surfer. 8.4 Calibration Similar to the Agass i z -Roseda le reach, two sets of calibration data were used for this reach: measured water levels from the 1999 freshet and A D C P data col lected during August of 1999. This calibration process is only intended to ensure 115 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach the model is producing relatively accurate results so that various 2D modeling applications can be investigated with some degree of confidence. 8.4.1 1999 Freshet The peak flow during the 1999 freshet on the Fraser River, which occurred on June 23, 1999 and was measured to be 11,100 m 3/s at Hope, was used as a calibration event for this study reach. Measured water levels at the four gauges (Kent 8, Chwk 7, 40 and Chwk 10) shown in Figure 8.1 were available to compare to model results. The upstream boundary conditions for the upper half of this reach were discharges of 11,100 m 3/s and 500 m 3/s for the Fraser and Harrison River inflows, respectively. The Harrison River discharge of 500 m 3/s was estimated by UMA Engineering Ltd. (2000) to correspond to the timing of the 1999 peak flow in the Fraser River. A fixed downstream water level of 11.0 m, which was specified at the intermediate reach boundary, was obtained from measurements from the Chwk 7 and 40 gauges. Obviously the water levels measured at these two gauges for the peak flow could no longer be used to compare with simulation results. The upstream boundary conditions for the lower half of the reach, which includes Minto Island, the adjacent channels and channel and floodplain areas downstream of Minto Island, were obtained from the simulated discharges in the upper reach. The inflows for this part of the reach were 7,700 m 3/s and 3,900 m 3/s for the main channel and Minto channel, respectively. The downstream water level was fixed at 9.45 m, which was obtained from the MIKE11 water surface profiles (UMA Engineering Ltd., 2000a). The goal was to match the simulated water surface elevations to the measured values at the four gauge locations as closely as possible. Since the water levels from gauges Chwk 7 and 40 were used to provide boundary conditions, comparisons could only be made at the Kent 8 and Chwk 10 gauges. Several combinations of ks values were attempted for island and floodplain areas until the 116 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach water levels at the two gauges matched closely. The channel k s was maintained at 0.04 m since the upstream end of this reach is only 8 km downstream of the Agassiz-Rosedale bridge where the d50 of the surface sediment is 0.042 m. It is unlikely that the dso of the grains on the channel bed would change significantly in that relatively short distance. Through trial and error, the best match in water levels was obtained when the vegetated portions of Minto Island and two floodplain areas were assigned ks values of 2.0 m. This was the same roughness value used for the non-channel areas in the Agassiz-Rosedale reach. The two floodplain areas that were given higher roughnesses were located on the south bank near the inflow section and on the north bank near the outflow section. The roughness regions, which could only be approximately defined using the bed topography file editor, are shown in Figure 8.4. Figure 8.4: Harrison River confluence/Minto Island reach roughness regions 117 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Table 8.1 displays the observed and simulated water levels as well as the difference between them. In addition, the MIKE 11 results from the UMA Hydraulic Modeling Study (2000a) are also shown for comparative purposes. Although water levels can only be compared at two gauges, the simulation results are very good with differences of only -0.03 m and -0.01 m. Table 8.1: 1 9 9 9 freshet calibration results in the Harrison River confluence/Minto Island reach Gauge # MIKE 11 Observed River2D Difference (m) (m) (m) (m) Kent 8 11.84 11.88 11.85 -0.03 Chwk 7 10.94 11.02 *** -40 10.96 10.97 *** -Chwk 10 9.85 9.75 9.74 -0.01 8.4.2 ADCP Data The A D C P survey through the Harrison River confluence/Minto Island reach was conducted on August 4-5, 1999. The transects in the upper half of the reach were surveyed on August 5, 1999 when the measured flow for the Fraser River at Hope was about 6,800 m 3/s. This discharge was used as the upstream boundary condition for the Fraser River inflow and the Harrison River inflow boundary condition was set to 1,360 m 3/s, which was the average discharge estimate from the three A D C P transects on the Harrison River. The downstream water level for the upper half of the reach was set at 9.30 m, which was obtained from the water level estimates made during the A D C P survey. The inflows for the lower half of the reach were 5,700 m 3/s and 2,460 m 3/s, which were assigned to the main channel and Minto channel, respectively. This flow split was determined from the simulation results from the upstream portion of the reach. The downstream water level was set 118 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach to 7.90 m, which was also obtained from the water level estimates made during the A D C P survey. The ks values that were assigned in the 1999 freshet calibration were maintained for the A D C P simulations. Figure 8.5 shows the location of the A D C P transects as red lines in the Harrison River Confluence/Minto Island reach that were used to compare to model results. Figure 8.5: A D C P transect locations in the Harrison River confluence/Minto Island reach 119 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach The River2D results were compared with the A D C P data and no adjustment of model parameters were made to obtain a better match. As with the Agassiz-Rosedale reach A D C P data, the A D C P measurements for this study reach were based on single pings, which means no averaging was done in the field, therefore spatial averaging was performed in Excel. There was one measurement made approximately every meter, so every ten values of velocities and depths along each transect were averaged to reduce errors associated with individual point measurements. The simulated depths and velocities were extracted at the x, y coordinates of the averaged A D C P data. Figures 8.6 and 8.7 show scatter plots of the simulated versus measured depths and velocities, respectively. The results from the upper and lower halves of the reach have been combined and the 1:1 lines are also shown on the scatter plots. ADCP depth (m) Figure 8 .6 : Simulated depth versus A D C P depth for the Harrison River confluence/Minto Island reach 120 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 3.5 co 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 A D C P velocity (m/s) Figure 8.7: Simulated velocity versus A D C P velocity for the Harrison River confluence/Minto Island reach The mean depth error is very low at 1.7% while more variability is seen in the velocity comparisons (ME = -4.3%), which is a similar trend as the Agassiz-Rosedale reach comparisons. All the averaged A D C P depth data was used for the comparisons however certain sections of the averaged A D C P velocity data were filtered. The A D C P velocities were significantly lower than model results in relatively shallow depth areas (< « 3 m), which usually coincided with gravel bars. Therefore, the data in these areas were omitted from the velocity comparisons. The A D C P data was also used to compare the flow splits between the main channel and Minto channel in transects 8B2/8A3, 8B1/8A1 and 7A/7B, shown in Table 8.2. The MIKE 11 results are also shown for comparative purposes. The relative proportion of the flow distribution between adjacent channels, which is probably more relevant, is shown rather than absolute values. 121 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Table 8.2: Flow splits (%) in the Harrison River confluence/Minto Island reach Transect A D C P MIKE 11 River2D 8B2 67.2 61.1 69.1 8A3 32.8 38.9 30.9 8B1 66.0 59.2 68.4 8A1 34.0 40.8 31.6 7A 92.2 83.7 84.1 7B 7.8 16.3 15.9 The results from both MIKE 11 and River2D are comparable with the A D C P values. The largest discrepancy in Table 8.2 occurs in transect 7A and 7B, where both River2D and MIKE11 simulate approximately 8% more flow along line 7B, which is a side channel near the downstream end of this reach. The two computer simulation results are likely more accurate in this case since the A D C P survey notes indicate that there was a gap in that sampled cross-section due to shallow water, which would have reduced the discharge estimate along that line. The water levels measured along each line during the A D C P survey using a differential G P S and a high speed boat are compared to the model results in Table 8.3. The measured water levels were given as constant values along each line while the River2D results show some lateral variation in the water surface elevations. The simulated water levels across each transect were averaged so they could be directly compared with the measured values. 122 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Table 8.3: Comparison of water levels measured during the A D C P survey in the Harrison River confluence/Minto Island reach Transect Measured Simulated Difference (m) (m) (m) 9 10.2* 10.62 +0.42 8A3 9.54 9.50 -0.04 8B2 9.46 9.51 +0.05 HARRIS3 9.9** 10.37 +0.47 HARRIS2 9.9*** 11.0 +1.1 8A1 8.47 8.61 +0.14 8B1 8.72 8.73 +0.01 7A 8.1 8.07 -0.03 7B 8.30*** 8.12 -0.18 Notes: * - line measured between two closely spaced gravel bars, poor line? ** - elevations measured 300 m away *** - elevations measured 900 m away **** - elevations measured 800 m away The River2D simulated water levels are within 0.14 m of the measured values, except in the following cases. Along transects HARRIS2, HARRIS3 and 7B the measured water levels were taken several hundred meters away, as noted in Table 8.3, and consequentially the difference between these and the simulated levels are very high, especially for the two lines in the Harrison River. The measurement along line 9 was taken between two closely spaced gravel bars and the A D C P survey notes indicated that the line may be poor due to the boat wandering. This may help to explain the difference of 0.42 m between the measured and simulated water level along this transect. 123 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 8.5 Summary In this section an overview of the Harrison River confluence/Minto Island reach was provided. A description of the DEM generation and River2D model development for this reach was briefly discussed since detailed information on these topics was previously provided for the Agassiz-Rosedale reach. The comparison of model results to measured water levels from the 1999 freshet and A D C P data indicates that an adequate level of calibration has been achieved for the purposes of this study and therefore further modeling scenarios can be investigated. Several applications of River2D in the Harrison River confluence/Minto Island reach are presented in the following section. 124 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 9. APPLICATIONS: HARRISON RIVER CONFLUENCE/MINTO ISLAND REACH 9.1 Harrison River Confluence Energy Loss Some of the features of the complex flow hydraulics that occur in the confluenceof the Fraser and Harrison Rivers have been previously described, such as the abrupt bend in the main flow as the Fraser River strikes the outcropping bedrock at Harrison Knob. Strong secondary currents in this area characterize the three-dimensional flow patterns. The River2D model developed for this reach cannot simulate the eddies that form in this region because the node spacing of the finite element mesh is 35 m, which is too coarse to accurately simulate flow hydraulics at the eddy scale. To properly model the complex flow patterns in the confluence of the two rivers, a 3D hydrodynamic model would be required to resolve the secondary flow patterns. Although a 2D model cannot be used to analyze the complex flow hydraulics in this part of the river, it can be used to investigate the large amount of energy that is dissipated. Backwater conditions encourage sedimentation upstream of the bend and the secondary currents and associated turbulence that occur in the bend result in high energy losses, which leads to a steep drop in the water surface. The MIKE11 water surface profiles (UMA Engineering Ltd., 2000a) indicated a drop in water level from approximately 13.4 m to 12.6 m in this section of the river for the design flood simulation. The distance this energy expenditure occurs over is about 750 m, which gives a very steep hydraulic grade line slope of approximately 0.0011. The design flood was simulated using River2D for this reach (see Section 9.2) and water surface elevation contours were plotted, shown in Figure 9.1. 125 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 545400 5453000 > 5452000 5451000 5450000 576000 577000 578000 579000 X(m) L E G E N D 13.6 m 13.4 m 13.2 m 13.0 m 12.8 m 12.6 m '12.4 m Flow direction High energy loss area Banklines/lsland boundary Contour interval = 1 m Figure 9 . 1 : Water surface elevation contours in the Harrison River confluence area for the design flood The 2D model, which has the advantage of providing information on the spatial variation in hydraulic parameters, shows similar results as MIKE11. The water surface drops from about 13.4 m to 12.5 m in a distance of approximately 1,000 m, which gives a hydraulic grade line slope of 0.0009. Unlike the MIKE11 profile, which shows a linear drop in energy, River2D indicates that the energy is dissipated in two distinct stages. This is seen clearly in the design flood profile (Figure 9.2), presented in the following section. In the first stage, the water level drops 0.55 m in about 350 m (HGL slope = 0.0016). This zone of high energy expenditure is followed by a 250 m length of river where the water surface profile is 126 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach relatively flat. Examination of the bed topography indicates that there is a very deep scour hole in this region with local elevations approaching - 7 m, likely caused by the erosive power of the high kinetic energy in the river just upstream. A second significant drop in water surface elevation occurs downstream of this scour hole where the water level falls approximately 0.35 m in 400 m (HGL slope = 0.00088). Both MIKE11 and River2D provide similar values for the overall energy loss through the 1 km long reach of the Fraser River immediately downstream of the Harrison River mouth, although the roughness values in the 1D model had to be increased significantly in order to do so. Since MIKE11 has no other means of simulating such a large energy loss through an abrupt bend, the Manning's n values had to be raised to a maximum of 0.10 in parts of the reach in order to match measured water levels both at the mouth of the Harrison River and at downstream locations (UMA Engineering Ltd., 2000a). On the other hand, the ks values assigned during the calibration process were not changed in River2D, which was able to simulate the high energy expenditure without any manipulation of model parameters. This result illustrates how the roughness parameter in 2D models only represents losses due to direct bed shear and all other forms of energy dissipation are accounted for by the model itself. This application provides an example of how River2D can be successfully applied to predict unusually high energy losses in gravel-bed rivers. 9.2 Design F lood As with the Agassiz-Rosedale reach, River2D was used to compute the water surface profile for the design flood through the Harrison River confluence/Minto Island reach. The simulation was performed in two stages, with many of the boundary conditions provided by the MIKE11 results. For the upper half of the reach, the two inflows that were specified were 17,000 m 3/s and 1,300 m 3/s for the Fraser and Harrison Rivers, respectively. The estimation of the Harrison River inflow corresponding to the timing of the Fraser River design flood is described in a 127 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach report on the hydrologic and hydraulic investigations of the two rivers by UMA Engineering Ltd. with Ward & Associates and Quick (2000a). The downstream water levels, obtained from the MIKE11 generated profiles, were 12.4 m for the main channel and 12.6 m for Minto channel. The inflows for the lower half of the reach, determined from the upper half simulations, were 12,000 m 3/s and 6,300 m 3/s for the main channel and Minto channel, respectively. The MIKE11 results provided the downstream water level of 11.3 m. The node spacing in the finite element mesh was 35 m and the roughness values indicated in Figure 8.4 were not adjusted. The water levels along the thalweg of the main channel were used to generate the flood profile, shown in Figure 9.2. The upstream and downstream profiles are represented by separate lines and the MIKE11 profile is also shown for comparative purposes. The chainages correspond to the distances shown in Figure 8.1 c o "5 > M m u m t Z2 cn •_ 0) 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 11.4 11.2 111 2 ^ " +000 112+000 113+000 114+000 115+000 116+000 117+000 118+000 Chainage (m) MIKE11 •River2D u/s •River2D d/s Figure 9.2: Design flood profile for the Harrison River confluence/Minto Island reach 128 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach At the downstream end of both the upper and lower halves of the reach, the match between River2D and MIKE11 is very good. This is largely due to the MIKE11 water levels being used as the downstream boundary conditions in the River2D simulations. As discussed in the previous section, the two models give very comparable results through the high energy loss zone just downstream of the Harrison River mouth. Although the 2D model for this reach extends about 2.5 km upstream of the Harrison River mouth, the profiles for the upper 2 km of the reach are not shown. The River2D profile, which is very flat upstream of the confluence, deviates significantly in this region from the MIKE11 profile, which continues to increase in the upstream direction. In both the lower part of this study reach and in the Agassiz-Rosedale reach, the River2D profiles are lower by about 0.1 m or more in some areas in the upper 1.5 to 2 km. Those differences may be due to assumed errors in the lateral distribution in the inflow discharge in the 2D model, which require some distance to dissipate before the q distribution more closely approximates the actual flow patterns in the river. However, it is unlikely that this possibility can fully explain the large difference between the two profiles at the upper end of this reach. This difference was not able to be resolved in this study. This application illustrates some of the limitations of using 2D models to generate water surface profiles over long distances for large rivers and further supports the assertion that 1D model are better suited for this purpose. The simulated water levels by River2D upstream of the confluence should be investigated in future work. A 2D model could be developed for the Harrison confluence reach to include 2 km of river upstream of the existing inflow section which would help to examine boundary condition effects on water levels through the area. This model would have to have a coarser mesh node spacing unless some of the existing reach area is excluded. The Minto channel outflow section could be moved further upstream which would remove a large part of the side channel and possibly also a portion of Minto Island from the model. In addition, the MIKE11 profiles could be used to specify a fixed water level as the upstream boundary 129 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach condition, which would cause the discharge to become a variable that the model computes. 9.3 Harrison Bar Pilot Channel UMA Engineering Ltd. (2000) developed gravel extraction scenarios to assess possible flood hazard mitigation strategies that could be used in the gravel reach. One of these scenarios included the excavation of a 600 m wide channel across Harrison Bar, which is located near the mouth of the Harrison River. This channel excavation was part of a maximum gravel excavation scenario where gravel areas in the reach would be excavated down to either the low water level or 1 m or 2 m below the low water level. The relief, or pilot channel, was to be excavated to a depth of 2 m below the low water level and would require the clearing of some floodplain forests in the area. In the MIKE11 study, the pilot channel not only helped to reduce water levels in the main channel, since much of the flow was conveyed through the newly formed channel, but also significantly reduced the steep drop in the design flood profile downstream of the Harrison River mouth. In this study, River2D was used to assess the effects of excavating such a pilot channel across Harrison Bar in more detail. The location of the modeled pilot channel is shown in Figure 9.3. The length of the excavated channel is approximately 2.5 km and is 430 m wide on average, smaller than the 600 m width suggested by UMA Engineering Ltd. The total excavated area would be 112.6 ha, of which 53.7 ha, or almost 48%, would be forested regions. 130 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Figure 9.3: Location of Harrison Bar pilot channel For this application, simulations were only performed for the upper half of the reach. The first step was to determine low water levels so the excavation depths could be established. The model was run for discharges of 500 m 3/s for the Fraser River and 100 m 3/s for the Harrison River. From examination of long-term flow records of the two rivers, these flows were among the lowest that occurred during February and March in the driest years. Precise definition of the long-term minimum flows was not necessary as only approximate low water levels were required for the purposes of this analysis. From the River2D simulation, the water surface elevations in the main channel at the upstream and downstream ends of the proposed pilot channel were 7.4 m and 3.9 m, respectively. It is interesting to note that this low discharge simulation reached steady state conditions with no instability whatsoever. In all previous simulations where regions of high ks had been included, a portion of the total discharge would flow over the non-channel areas resulting in local 131 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach instabilities, although on a reach-wide scale steady state conditions would be eventually achieved. The areas of local instability were in the vicinity of adjacent wet and dry regions. In this case, all the flow was restricted to channel areas and consequently no stability issues arose which suggests that the instability may be caused by a combination of high ks regions and wet/dry areas. The DEM bed elevations in the excavation area were modified using Surfer. At the upstream end of the pilot channel, the bed elevations were set to 5.4 m while at the downstream end of the channel the elevations were reduced to 1.9 m. The elevations between these points were adjusted to create a smooth bed slope of about 0.0014, which is about an order of magnitude steeper than the bed slope of the main channel in this reach. The average and maximum cuts in the excavated area were 5.18 m and 8.44 m, respectively. The total volume of gravel removed to form the pilot channel was 5,835,000 m 3. The x, y, z data from the modified DEM were used to create a new finite element mesh, with nodes spaced at 35 m as before. All the nodes in the pilot channel were assigned a ks of 0.04 m. The design flood was simulated for this scenario using the boundary conditions described in the previous section. The model results were extracted and contour plots were created showing the effects of the pilot channel excavation. Contour plots of depth and velocity for both pre- and post-excavation conditions are shown for the design flood in Figures 9.4 and 9.5, respectively. Figure 9.6 shows the effects of the pilot channel excavation on the design flood profile along the main channel in comparison to the pre-excavation condition and also to the MIKE11 profile. 132 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach a.) depth contours without pilot channel: 576000 576500 577000 577500 578000 578500 579000 579500 I X ( m ) b.) depth contours with pilot channel: 576000 576500 577000 577500 578000 578500 579000 579500 I X j rn) \ Figure 9.4: Design flood depth contours for the Harrison Bar pilot channel: Pre- and post-excavation 133 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach a.) velocity contours without pilot channel: 576000 576500 577000 577500 578000 578500 579000 579500 Xjm) b.) velocity contours with pilot channel: , , J 1 5Lio , , , 576000 576500 577000 577500 578000 578500 579000 579500 Xjrn) Figure 9.5: Design flood velocity contours for the Harrison Bar pilot channel: Pre-and post-excavation 134 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach c o CD > Q) <D o (f) (D -#-» CD 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 115+500 116+000 116+500 117+000 Chainage (m) R 2 D : without pilot ( :hannel R 2 D : with pilot cha nnel 117+500 118+000 •MIKE11 •River2D: without pilot channel •River2D: with pilot channel Figure 9.6: Design flood profile along the main channel for the Harrison Bar pilot channel The contour plots show that the depths in the pilot channel are greater than 8 m at the upstream end and greater than 10 m at the downstream end, while the pilot channel velocities range between 2 and 2.5 m/s. While it is difficult to determine whether the effects on depth are significant based solely on the contour plots, the effects on velocity are very noticeable both in the main channel and also in Minto channel. On Harrison Bar, at the mouth of the Harrison River, the pilot channel causes a reduction in velocity of about 0.5 m/s, and a similar reduction in velocity occurs in the steep reach downstream of the Harrison River mouth. An even greater decrease in velocities occurs in Minto channel where a reduction of approximately 1.5 m/s occurs following the pilot channel excavation. Another effect due to the pilot channel is a reduced backwater effect on the Harrison River, which is evident by the post-excavation velocity contours that indicate higher velocities at the 135 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach downstream end of the Harrison River. This is due to reduced flows in the Fraser River main channel at the confluence. Figure 9.6 indicates that the pilot channel reduces the large energy loss just downstream of the Harrison River mouth by 0.4 m. Water levels upstream of the confluence are also reduced by almost 0.4 m. In this part of the river, the effects of the pilot channel excavation seem beneficial from both flood hazard mitigation and bank erosion perspectives. The reduction in velocities in the main channel and Minto channel would help to reduce some erosion concerns while the lowering of the design flood profile upstream of the confluence would benefit flood control. It would also appear that the pilot channel would reduce the possibility of a future avulsion in this reach that would result in Minto channel becoming the main channel. This potential avulsion would considerably increase erosion risks along the south bank. However, the ecological and morphological effects of such a major excavation cannot be determined. A 2D hydrodynamic model capable of mobile bed simulations would be useful to assess the morphological impacts downstream of the pilot channel. UMA Engineering Ltd. (2000) suggested that if a much smaller channel than modeled in their hydraulic modeling study, or in this study for that matter, were excavated across Harrison Bar, than the energy potential due to the difference in upstream and downstream water levels would be large enough to enlarge the channel on its own. This high-energy channel could potentially grow to become the main channel of the river and remove the steep drop in the water surface profile. The main drawback to the high-energy channel concept is that bed material that would be eroded from the natural enlargement of the channel would be deposited at downstream locations, which would likely lead to additional instability and channel shifting during the annual freshets. The high-energy channel concept cannot be investigated with the fixed- bed River2D model, as sediment transport capabilities would be required. 136 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 9.4 Dredging of Harrison Bar The Harrison Bar pilot channel concept was investigated by Water Management Consultants (2001) in a preliminary study for the City of Chilliwack. The study focused on the hydraulic, environmental and morphological issues which would arise with a bar excavation in the area. The goal of improving existing dyke protection against the design flood by reducing upstream water levels and stabilizing the river such that the possibility of Minto Channel becoming the main channel decreases initiated this study. Water Management Consultants (2001), the City of Chilliwack and BC Environment decided that the pilot, or high-energy, channel across Harrison Bar would lead to undesirable downstream morphological impacts. Therefore, the proposed solution was to excavate, or dredge, the downstream or northern tip of Harrison Bar, which would act as a sediment trap and likely lessen downstream impacts. The proposed Harrison Bar dredge would decrease the bend's radius of curvature and increase the conveyance through the region due to an increase in cross-sectional area, which should lower water levels. The rationale or potential benefit of reducing the radius of curvature would hopefully be to decrease the strength of secondary currents, which is one the main causes of the high energy loss through the bend. A reduction in energy loss would reduce upstream water levels. The morphological module in MIKE11 was used by Water Management Consultants to simulate the effects of the dredging of Harrison Bar. The location of this proposed excavation is shown in Figure 9.7. The areal extent and depths of excavation simulated in the MIKE11 morphological model were used to modify the DEM in the study reach in order to investigate this gravel removal scenario with River2D. The original DEM developed for the Harrison River confluence reach was modified in Surfer by reducing bed elevations in the proposed excavation area. The bed elevations were reduced to 5.5 m and 4 m at the upstream and downstream ends of the excavation, respectively. According to Water Management Consultants, these elevations are 137 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach below low winter water levels in those locations. The average and maximum cuts were 3.01 m and 6.03 m, respectively. In total, 950,000 m 3 of gravel was removed over an area of 31.4 ha. Figure 9.7: Location of Harrison Bar dredge area The x, y, z data from the modified DEM were used to create a new finite element mesh, with nodes spaced at 35 m as before. The design flood was simulated for this scenario using the similar boundary conditions as described in Section 9.2. The model results were extracted and contour plots were created showing the effects of the pilot channel excavation. Contour plots of depth and velocity for both pre- and post-excavation conditions are shown for the design flood in Figures 9.8 and 9.9, respectively. Figure 9.10 shows the effects of dredging of Harrison Bar on the design flood profile along the main channel in comparison to the pre-excavation condition, the pilot channel profile and also the MIKE11 profile. 138 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach a.) depth contours without Harrison Bar dredge: 5454000 5453500 5453000 5452500 >- 5452000 5451500 5451000 5450500 16m 14m 12 m 10m 576000 576500 577000 577500 578000 578500 579000 X(m) 579500 b.) depth contours with Harrison Bar dredge: 5454000 5453500 5453000 5452500 >• 5452000 5451500 5451000 5450500-576000 576500 577000 577500 578000 578500 579000 579500 Xjm] Figure 9.8: Design flood depth contours for the Harrison Bar dredge: Pre- and post-excavation 139 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach a.) velocity contours without Harrison Bar dredge: 576000 576500 577000 577500 578000 578500 579000 579500 XJm] b.) velocity contours with Harrison Bar dredge: 576000 576500 577000 577500 578000 578500 579000 579500 XJm] Figure 9.9: Design flood velocity contours for the Harrison Bar dredge: Pre- and post-excavation 140 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach o +-< ro > a> <x> o •E CO i _ ro 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 — • , • m R2D * if * •m I — R2D: pilot chann 115+500 116+000 116+500 117+000 Chainage (m) 117+500 118+000 - - MIKE11 •River2D •River2D: Pilot Channel •River2D: Dredged Bar Figure 9.10: Design flood profile along the main channel for the Harrison Bar dredge The contour plots indicate that dredging of Harrison Bar results in a slight increase in velocities upstream and downstream of the excavated area. The velocities within the dredged area are reduced, although not significantly, due to the increased flow area. In addition, the depth of flow is higher in the excavated area after dredging, as would be expected. Impacts of the dredging on Minto Channel are minimal since only a modest decrease in velocities is observed from the contour plots. Figure 9.10 indicates that the proposed dredging would reduce water levels for the design flood in the Harrison Bar area by only 0.12 m, which is similar to the value simulated by the MIKE 11 model. The energy loss through the bend is also reduced by approximately 0.12 m, although the drop still occurs in two distinct stages, as described in Section 9.1. 141 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Overall, dredging Harrison Bar does not have a significant impact in this part of the reach. The modest decreases in upstream water levels and in Minto Channel velocities that would be achieved through the proposed dredging would provide little benefit towards minimizing flood hazards and bank erosion. The hydraulic control section for this reach is clearly the steep reach, where the large energy loss occurs, downstream of Harrison Bar. The reduction of flows in that part of the river is key to reducing upstream water levels which would aid flood mitigation efforts, as was illustrated by the pilot channel 2D investigation presented in the previous section. Excavating the bar upstream of the bend does not adequately address the high energy loss that takes place downstream. The initial 1D morphological modeling study of the dredging of Harrison Bar was performed because it was thought that this excavation offered a compromise between reducing flood risks in the reach while also minimizing downstream morphological impacts. Although the latter may be true, the limited flood control benefits that this excavation strategy would provide probably do not justify its future implementation. 9.5 Island 22 Bank Eros ion Following the 1999 freshet, a bank erosion problem along an 800 m length of unprotected bank at Island 22, which is located on the south bank of the river downstream of Minto Island (Figure 9.11) developed and has since rapidly increased. This erosion threatens a municipal campground and more importantly the Wing Dyke, which provides flood protection for much of the centre of the City of Chilliwack. The main channel of the Fraser River, which flows on the north side of Minto Island, crosses over to the opposite side of the river just downstream of the island and impinges on the left bank, resulting in the erosion. Two consulting firms, Northwest Hydraulic Consultants Ltd. (2001) and UMA Engineering Ltd. (2001) were retained by the City of Chilliwack to assess this bank erosion problem. A series of cross-sections were surveyed downstream of Minto Island in August 1999, January 2001 and March 2001. The successive cross-142 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach sections show the thalweg migrating laterally towards the left bank and the channel becoming deeper. Between August 1999 and March 2001 the bank line had moved back 20 m to 30 m due to the erosion, however the future rate of migration of the thalweg will more closely match the rate of bankline recession since the portion of the left bank below water has been scoured considerably (UMA Engineering Ltd., 2001). 1996 air photo: 2000 air photo: Figure 9.11: Island 22 bank erosion 143 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach UMA Engineering Ltd. (2001) predicts that the bankline will recede a further 30 m to 60 m due to the 2001 freshet, while Northwest Hydraulic Consultants (2001) estimate an erosion of 70 m, although that figure is based on an average 2001 freshet flow of 9000 m 3/s, which is considerably higher than that experienced this year. At the current rate of erosion, the Wing Dyke may be threatened and could fail in 2003 or 2004, or possibly even next year (2002) if a large freshet were to occur. The main cause of the Island 22 bank erosion has been attributed to a rapidly growing gravel deposit on Queen's Bar, which is on the opposite side of the main channel (UMA Engineering Ltd., 2001 and Northwest Hydraulics Consultants, 2001). Figure 9.11 shows portions of two air photos, taken in 1996 and 2000. The 1996 air photo shows the downstream half of the study reach and indicates the location of the Queen's Bar deposit, which began accumulating gravel in 1999. The portion of the 2000 air photo shown is of only a small section of the river but the newly formed Queen's Bar gravel deposit is evident. The gravel deposit has continued to grow since last year and has reduced the channel's cross-sectional area, which increases velocities. It has been widely suggested that the deposit also redirects the flow of the main channel so the river impinges more directly on the unprotected portion of the left bank. Both consultant reports suggest that as a result of the gravel deposit growth, the thalweg migrates further towards the opposite side of the river and scours the left bank. . UMA Engineering Ltd. (2001) described two possible contributing factors to the rapid growth of the Queen's Bar gravel deposit. The first is the propagation of a shift in the meander pattern in the main channel below the Harrison River that began in the 1960's and 1970's. The meander pattern has naturally progressed downstream such that it now actively erodes part of a gravel bar on the right bank of the main channel, opposite the downstream end of Minto Island, and likely deposits this material on Queen's Bar, creating the rapidly growing gravel lobe. Another possible factor are reduced flows in Minto Channel due to increased sediment 144 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach deposition at its entrance. This has increased the momentum of the main channel flow while decreasing the Minto Channel flow momentum, which results in the flow being gradually rotated in a counter clockwise direction at their confluence downstream of Minto Island. Therefore, the increasing angle of attack of the main channel flow on the left bank and gravel deposition on Queen's Bar in the region of relatively lower velocities can likely be partially attributed to the Minto Channel reduced flows. River2D can be a valuable tool in further investigations of this complex bank erosion problem. For example, Northwest Hydraulics Consultants (2001) proposed an excavation of the gravel bar deposit on Queen's Bar to reduce Island 22 bank erosion. River2D could be used to determine the effects of this proposed excavation, including an analysis of shear stresses along the left bank. Based on a bar survey conducted in January 2000, Northwest Hydraulic Consultants proposed the excavation of approximately 4.25 ha on Queen's Bar down to 2 m below the water level at the time of excavation. The DEM generated in this study, which is based on August 1999 survey data, was modified to represent the Queen's Bar deposit excavation down to a level of just under 3 m, which was similar to the preliminary assumption made by Northwest Hydraulics. The average and maximum cuts in the DEM were 2.37 m and 3.03 m, respectively and the total volume of gravel removed was approximately 100,000 m 3. By comparison, sections through the proposed excavation prepared by Northwest Hydraulics showed a maximum cut of almost 5 m. This indicates that between August 1999 and January 2000 about 2 m of gravel had deposited in this area. If this potential solution were implemented today the total volume of gravel removal would obviously be much greater than 100,000 m 3, especially considering that the Queen's Bar deposit continued to grow at a fairly rapid rate following the January 2000 bar survey. Since no survey information collected following the extensive reach-wide August 1999 survey was obtained for this study the overall effects of the proposed Queen's Bar excavation cannot be determined at this time. Instead, River2D was 145 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach used to compare the effects of the gravel removal with respect to the August 1999 bed conditions. The pre-excavation velocities and depths were obtained from the calibration simulations performed for the lower half of the Harrison River confluence/Minto Island reach using the 1999 freshet peak flow while the post-excavation conditions were represented by the modified DEM described in the preceding paragraph. The modified DEM was used to create a mesh that was used to simulate the 1999 freshet through the lower half of the reach using inflows of 7,700 m 3/s and 3,900 m 3/s for the main channel and Minto channel, respectively and the downstream water level was set to 9.45 m. These were identical to the boundary conditions used for the 1999 freshet calibration simulations (Section 8.4.1). A q vector plot for the entire lower half of the Harrison River confluence/ Minto Island reach for the 1999 freshet peak flow is shown in Figure 9.12 for the post-excavation conditions. Contour plots of the pre- and post-excavation velocities are shown in Figure 9.13 while the change in velocities in the vicinity of Queen's Bar due to the proposed gravel removal is shown in Figure 9.14. A similar contour plot indicating depth changes due to the excavation is presented in Figure 9.15. 5452000 5451500J 5451000 5450500 5450000 5449500 5449000 5448500 Discharge Vectors -* 3* 1 m/sA2 43 m/sA2 s v / * r -i I 7 * * /*; 7 I * A Excavation area 571500 572000 572500 573000 573500 574000 574500 575000 575500 576000 576500 577000 577500 XJm) Figure 9.12: 1999 freshet discharge vectors for the Queen's Bar gravel deposit proposed excavation 146 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach a.) velocity contours prior to the excavation of the Queen's Bar gravel deposit (August 1999 survey data): 571500 572000 572500 573000 573500 574000 574500 575000 575500 576000 576500 577000 577500 X(m) b.) velocity contours following the excavation of the Queen's Bar gravel deposit: X(m) Figure 9.13: 1999 freshet velocity contours for the Queen's Bar gravel removal: Pre-and post-excavation 147 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 545050a 5450000^  5449500 5449000 5448500 572500 573000 573500 574000 574500 575000 X(m) 0.2 m/s 0.1 m/s 0.0 m/s -0.1 m/s -0.2 m/s -0.3 m/s -0.4 m/s -0.5 m/s Figure 9.14: 1999 freshet velocity contours for the Queen's Bar gravel removal: Change in velocity due to extraction 5450500 5450000^  S 5449500^  >-5449000 5448500 Triangle represents removal boundary 2.0 m 0.2 m 0.0 m 572500 573000 573500 574000 574500 575000 " - 0 . 2 m X(m) Figure 9.15: 1999 freshet depth contours for the Queen's Bar gravel removal: Change in depth due to extraction 148 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach The velocity contours in Figure 9.13 indicate that channel velocities are between 2 and 2.5 m/s along the left bank downstream of Minto Island for the 1999 freshet and the discharge vectors indicate that the thalweg is situated extremely close to the north side of Island 22. From examination of the pre- and post-excavation channel velocities in Figure 9.13, the effects of the excavation appear to be very minimal. Figure 9.14 confirms this fact as the velocity differences are typically between ±0.1 to 0.2 m/s. The slightly higher velocity differences at the downstream end of the reach should be ignored due to boundary condition uncertainties. The largest velocity difference of about -0.5 m/s occurs within the excavated region and there is a slight increase in velocity just upstream and downstream of the removal area. No other significant trends in velocity changes are evident and the minimal differences elsewhere likely fall within the level of accuracy of the computational model given the relatively coarse 35 m node spacing. The depth differences are also quite low as the Queen's Bar gravel removal only results in depth changes of less than ± 0.2 m in the majority of the channel area shown in Figure 9.15. The only exceptions are the excavated area where the depths are increased by more than 2 m, as would be expected, and at the downstream boundary. Although these 2D simulation results do not indicate the overall effects of the Queen's Bar gravel removal since the latest survey data was not used, the procedure for conducting such an analysis using River2D and the types of figures which would be useful in the investigation of this bank erosion problem haven been presented and discussed. In future studies, the recent March 2001 survey data could be used to create a DEM and River2D model to allow the two-dimensional hydraulics related to the erosion problem to be analyzed in greater detail and more thoroughly assess the effectiveness of potential solutions such as the proposed Queen's Bar gravel deposit excavation. Other relevant analyses that could be performed with River2D would be to assess the effect of reduced flows in Minto Channel. A series of simulations could 149 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach be performed with varying proportions of the total flow divided between the main channel and Minto Channel and q vectors could be plotted showing the direction of flow in the channel adjacent to Island 22. Shear stresses along the left bank could also be calculated for the varying flow splits. A 2D morphological model of the reach would be even more useful as it would enable analysis of the effects of reduced Minto Channel flows on the growth of the Queen's Bar gravel deposit over time. A mobile bed model would also allow the downstream propagation of the meander shift in the main channel to be studied more closely. Such modeling studies would provide valuable insight and lead to a more thorough understanding of the Island 22 bank erosion problem. Based on the various gravel removal scenarios conducted in this study it would seem unlikely that the proposed removal would significantly reduce bank erosion at Island 22, let alone stop it altogether. The localized gravel removal investigated in the Agassiz-Rosedale reach (Section 7.1.5) considered an excavation of slightly more than 130,000 m 3 over an approximate area of 2 ha. The hydraulic effects of this gravel removal were limited to the excavation area and basically no noticeable change.in depths or velocities were evident in the adjacent main channel. In the case of the Queen's Bar gravel deposit, the current excavation volume could very well be in excess of 200,000 m 3 covering an area about twice as large as the localized removal scenario described. It is improbable that such a gravel extraction would appreciably reduce channel velocities in the thalweg for a river as large as the Fraser. Furthermore, gravel removal in this area would likely be a short-term solution of questionable effectiveness since the excavated area would still be prone to deposition in subsequent years. 150 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 10. DISCUSSION AND CONCLUSIONS The recommended modeling approach for the Fraser River Gravel Reach suggested by Millar and Barua (1999) in the hydrodynamic model selection study involved the use of two different models to satisfy the various modeling objectives that were defined. For the purpose of updating design flood profiles and identifying deficiencies in the dyking system, the development of a 1D unsteady network model for the entire reach was considered the most feasible solution. Depth-averaged 2D hydrodynamic models were recommended for local applications to investigate issues related to flow hydraulics, including gravel extraction scenarios, local scour and riverbank erosion and habitat assessment. The first part of the recommended modeling approach was satisfied by UMA Engineering Ltd. (2000 and 2000a) who recently completed a hydraulic modeling study of the gravel reach in which they updated its design flood profile using the 1D hydrodynamic model, MIKE11. This study begins to address the second part of the suggesting modeling strategy as the depth-averaged 2D model, River2D was applied to two sections of the gravel reach to investigate the potential analytical uses of 2D models for large braided river systems; to examine the limitations and advantages of a 2D modeling approach; and to characterize the hydraulics in reaches of interest. The two sections of the gravel reach that were studied with River2D were a 4.5 km section of the Fraser River at the Agassiz-Rosedale Bridge and an 8.5 km reach near the Harrison River confluence and Minto Island. The key steps in a typical 2D modeling study were investigated in detail beginning with the development of a digital elevation model (DEM) from topographic and bathymetric survey data, which had been previously collected. The most important step in any 2D modeling study is obtaining an accurate representation of the bed topography. Point kriging in Surfer 7.0, incorporating a simple linear variogram model with relatively mild anisotropy ratios, was found to be a reliable gridding technique to convert survey data collected in transects, which is not the 151 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach most suitable form for 2D modeling applications, into rectangular arrays of point elevations. The quality of the DEM is enhanced by defining breaklines for most of the significant breaks in slopes, such as the tops of banks and thalwegs. Breakline definition for channel features can be a time consuming task since in this study individual cross-sections had to be plotted in Excel in order for the x, y, z points corresponding to breaks in slope to be determined. Surfer is limited in this respect since it does not offer the ability to efficiently perform such tasks from a three-dimensional surface. The use of a GIS, or perhaps another surface modeling program such as Softdesk, would aid the DEM generation process in terms of both time and likely quality as well. Only the north bank along Harrison Hill was defined as a breakline within the channel areas in the Harrison River confluence/Minto Island reach due to the extreme time requirements. However, all the dykes in both reaches studied were defined as breaklines, thus preserving them as continuous, linear features. The size of the Fraser River, which is considerably larger than any previous channel modeled with River2D, required that the node spacing be set to 25 m to 35 m for most of the 2D modeling applications considered such that the practical limit of total number of nodes was not exceeded. This limitation was found to be in the order of 10,000 nodes for a given simulation, although it is largely dependent on computational resources. The resolution of the DEM grid was 25 m for both modeled reaches, which was considered appropriate based on the 200 m spacing between surveyed transects. In most cases the 25 m spacing was maintained in the finite element mesh, which is basically the computational node spacing; however, this distance had to be increased to 35 m for the Harrison River confluence/Minto Island reach due to its large size. A sensitivity analysis of the computational node spacing indicated that there would not be significant changes on a reach-wide scale for uniform spacings of up to 35 m, as the mean average velocity and mean average depth were found to vary by less than 1% and 1.5%, respectively. Although it is possible, and recommended, to define breaklines in the finite element mesh, due to computational restrictions this was not possible in this 152 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach study. Specifically, the mesh generation would crash during triangulation of nodes when breaklines were included. It is uncertain whether the use of the most "advanced" PC available on the market today would resolve this issue. A relationship was derived between Manning's n and the roughness length, ks, that is almost identical to Strickler's equation, which relates Manning's n to the median grain size diameter, dso- This suggests that ks, which is the primary calibration parameter in River2D, can be estimated directly with the surface dso of the bed material. This hypothesis was tested and verified in the Agassiz-Rosedale reach where the surface dso is 0.042 m. The calibration process involved comparing measured water levels at a few gauges for the 1999 peak freshet flow to simulated levels. A very close match between observed and simulated water levels was obtained when ks was set to 0.04 m at all nodes within the channel areas. The roughness length was not varied spatially within the channel areas since no further information regarding the surface dso was available. Non-channel areas such as floodplains and vegetated islands were assigned ks values of 2.0 m, which was required in order to match the design flood profile simulated by River2D with the MIKE 11 profile computed by UMA Engineering Ltd. Similar ks values for channel and non-channel areas were used in the Harrison River confluence/Minto Island reach. Assigning ks values locally using the bed topography file editor in River2D is crude and not very accurate as only simple polygons can be defined to represent areas where the roughness differs from the global values previously set in the bed topography files. Furthermore, the program only provides a visual display of the rectangular array of nodes, which doesn't differentiate between channel and non-channel areas. The ability to overlay channel, island and gravel bar boundaries, as in a GIS, on top of the array of nodes would be very useful. Currently, non-channel areas can only be defined with the aid of displayed x, y coordinates corresponding to the cursor location. Since Surfer 7.0 does not have the capability to link more than one attribute other than a z coordinate (elevation) to each grid node in the DEM, the 153 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach crude River2D ks definition was used in this study. The approach used was to assign all nodes a ks value of 0.04 m in Excel and then to graphically define polygon regions that approximately covered the non-channel areas where higher roughnesses were considered appropriate. The other River2D parameter that can be adjusted by the user is a coefficient that controls the eddy viscosity. Through a sensitivity analysis, this parameter was found to have negligible effects on model results and therefore it was concluded that the default value of 0.5 can be used in all subsequent simulations. Although the eddy viscosity does not noticeably affect simulation results, it is required in River2D to achieve closure of the St. Venant equations since no other equations are available to solve for the four transverse shear stresses at each computational node. Comparisons made between measured and simulated water levels for the 1999 freshet and the ADCP-derived and simulated flow splits in adjacent channels indicated that River2D was adequately simulating flow hydraulics through the modeled reaches. The relative proportions of flow in adjacent channels were deemed more important than the actual flow magnitudes. A comparison of A D C P and simulated velocities showed a poorer correlation in both reaches than similar comparisons of depths obtained from the A D C P ' s bottom tracking capability and River2D. Although some spatial averaging on the A D C P data was performed in Excel, the field measurements were based on only a single ping at each vertical, which can result in errors ranging from a few mm/s to as high as 0.5 m/s (RD Instruments, 1996). Time averaging is essential to reduce these errors to obtain more accurate velocities to use in 2D model calibration. Depths obtained from the A D C P bottom tracking proved to be useful as they compared very favourably with model results. The water levels measured with the differential G P S and high speed boat were in relatively close agreement with simulated levels in the Harrison River confluence/Minto Island reach, however they were lower than model results by up to 0.28 m in the Agassiz-Rosedale reach. The cause of this discrepancy is uncertain, although it could possibly be related to equipment setup issues in the boat during the 154 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach survey. Additional gauged water levels in both reaches would have been very useful to provide a more thorough calibration of the River2D results. Several applications of River2D were attempted in both reaches that illustrate the potential benefits of utilizing 2D models. A series of contour plots for the mean annual flood in the Agassiz-Rosedale reach provide an example of how River2D results can be used to gain valuable information regarding the spatial variation of various hydraulic throughout modeled regions. From a sediment transport or morphological perspective, the ability to calculate shear stresses from model results is very useful. River2D was shown to accurately simulate the superelevation at the Agassiz-Rosedale Bridge section, which was not properly modeled by MIKE11. The analysis of several different gravel extraction scenarios indicated that the hydraulic impact of localized removals of the order of 100,000 m 3 is limited to the excavated area and no significant change in depth, velocity or flow direction should be expected outside this region. The model was also used to determine local depths and velocities for a range of flows around Big Bar in the Agassiz-Rosedale reach. This enabled contour plots to be created for use in habitat delineation and mapping studies, such as the one currently being conducted by the Department of Geography at the University of British Columbia. One area where 1D models are superior to 2D models is in the generation of design flood profiles. MIKE11 was able to provide a design flood profile for the entire 65 km gravel reach without having to splice results from several different simulations whereas River2D was used to produce comparable profiles for only a few kilometers in any given simulation. River2D produced a similar profile for the lower 2-3 km of each modeled section although the upstream portions were consistently lower than MIKE11 water levels. This deviation could possibly be due to assumed errors in the upstream boundary condition as a uniform discharge is assigned across the inflow sections and some distance may be required for the flow distribution to achieve the actual lateral variation in the flow pattern. The MIKE11 profiles were used to provide the downstream boundary conditions for use in the 2D 155 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach simulations. This is another benefit of the recommended modeling strategy of utilizing both 1D and 2D models to study the Fraser River Gravel Reach. Having a calibrated 1D unsteady flow model such as MIKE11 allows the boundary conditions for local applications of River2D to be relatively easily determined. Without the MIKE11 results, successful application of River2D would have been difficult since measured water levels for the range of flows considered in the modeled reaches are limited, and in many cases non-existent. Due to the coarse node spacing required, River2D was not able to provide detailed information on the complex flow patterns that occur at the mouth of the Harrison River. A 3D model would likely be required to properly model the strong secondary currents in this part of the river. However, the large energy loss that takes place at the bend downstream of the confluence of the Harrison and Fraser Rivers was simulated without having to artificially increase the roughness parameter. This is significant since Manning's n had to be increased to values up to 0.10 in the MIKE11 model to match measured water levels at the mouth of the Harrison River. The treatment of the roughness parameter is one of the fundamental differences between 1D and 2D hydrodynamic models. In 1D models, where no other means of dissipating energy is available, the roughness parameter is used to account for all forms of energy losses whereas in 2D models, the roughness parameter only represents.losses due to direct bed shear (grain roughness) and all other energy abstractions (form roughness) are accounted for by the two-dimensional formulation of the model itself. This helps to explain why ks can be simply equated to the median grain size diameter, d s o , in 2D modeling applications. The large-scale scalping of gravel bars on Powerline Island in the Agassiz-Rosedale reach that was recommended by Northwest Hydraulics Consultants to alleviate erosion along the left bank downstream of the bridge would reduce velocities by up to 0.5 m/s or more in the main channel. This may provide some reduction in bank erosion although the preliminary analyses performed in this study indicate no change in the direction of flow in the main channel. Additional 156 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach investigations, possibly with the use of a 2D morphodynamic model, would be required to determine whether the future migration of the thalweg towards the left bank would be avoided by the proposed gravel excavation. The River2D simulations and analysis showed that the pilot channel concept across Harrison Bar is superior to the bar dredging proposal in reducing flood risks in the gravel reach upstream of the Harrison River confluence. The Harrison Bar dredging scenario would result in very limited benefits from both flood risk mitigation and bank erosion reduction perspectives while a pilot channel could potentially reduce upstream water levels during the design flood by approximately 0.4 m. Considering the hydraulics in this reach, the control section is the steep reach downstream of the Harrison River mouth where the large energy drop occurs. The pilot channel re-directs much of the main channel flow away from this section while the bar dredging merely results in flows being accelerated even more in the steep reach. However, a pilot channel excavation would likely lead to significant morphological and ecological impacts both in the immediate vicinity and at downstream locations. Some of the potential morphological impacts could be investigated with a 2D model capable of mobile bed simulations. Bank erosion along Island 22, downstream of Minto Island, is a complex problem that would benefit from a detailed study of the morphological processes in the reach. The most up to date survey data could be used with River2D to assess the effect of the Queen's Bar gravel deposit excavation on channel velocities. However, to gain a more thorough understanding of the development of the bank erosion and to determine whether proposed solutions may in fact alleviate erosion concerns and not compromise dyke safety, would be greatly assisted by 2D morphological modeling. The wide variety of modeling applications that were investigated with River2D allowed the benefits of a 2D modeling approach to be examined while also identifying some of the limitations that must be considered. The node spacing that 157 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach must be used to successfully model a given reach in the lower Fraser River is perhaps coarser than would be initially desired; however this study has shown that decent 2D model results can still be achieved at this scale using the methods described in this thesis. On a reach-wide scale, the model was able to achieve steady-state solutions without numerical stability concerns, although when flow over non-channel areas with high ks values (floodplains and vegetated islands) occurred, some local instabilities would arise. Even though these instabilities were very localized at only a few nodes throughout the modeled region and did not prevent the solution at the remaining nodes from achieving steady-state conditions, the investigation of the causes of this instability is recommended. One of the main drawbacks associated with River2D is the lack of quality pre-and post-processing software. For example, the bed topography file editor is of limited use due to its functional deficiencies previously described. Furthermore, incorporating GIS-type functionality such as the ability to overlay channel and island boundaries would enhance both the bed topography and mesh generation programs. In the latter case, this would make it easier for the user to vary the node spacing spatially by assigning fewer nodes in non-channel areas. The lack of post-processing capabilities in River2D can be overcome by utilizing both Excel and Surfer 7.0 to analyze and process model results, as was done in this study. It was stated in the opening section.of this thesis that this work is not intended to provide a thorough hydraulic assessment of the gravel reach but rather to focus on the use of 2D models to study large river systems such as the Fraser. While this work has hopefully initiated the process of increased 2D modeling in local river engineering studies, the further development of a model such as River2D should also be considered. Several references throughout this thesis have been made to 2D morphological modeling and how it could potentially benefit the study of the Fraser River Gravel Reach. Since River2D has been shown to be technically sound in terms of hydrodynamic simulations, the development of a morphodynamic module specifically designed to address issues in the gravel reach would be very beneficial. 158 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 11. RECOMMENDATIONS Based on the information presented and discussed in this thesis, the following recommendations are made to assist future 2D hydrodynamic modeling studies of the Fraser River Gravel Reach. 1. A GIS or other surface modeling software that allows for the graphical definition of breaklines for all linear channel features should be used to generate DEM's. At the very least, breaklines should be defined for all tops and bottom of banks (where applicable), thalwegs and dykes. 2. A GIS should be incorporated into future 2D modeling studies so that the roughness parameter can be easily varied spatially according to the type of surface (e.g. channel, gravel bar, vegetated island, floodplain, etc.). GIS interfacing would allow DEM information to be directly input into River2D and model results could be graphically displayed on the digital terrain surface. 3. Additional water level gauges should be considered in reaches of interest to provide more calibration data. 4. A D C P data should be collected with sufficient time averaging to minimize horizontal velocity errors. 5. If possible, use the existing calibrated MIKE11 model of the gravel reach to obtain water levels to be used as boundary conditions for 2D modeling applications. 6. Investigate the use of aerial photos for use in 2D model calibration by comparing the lateral extent of the water surface at a given flow. 159 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 7. Investigate the sensitivity of 2D model results to the choice of the dynamic and diffusion wave solutions. 8. Investigate the local instabilities that arise in River2D simulations when flow over high ks regions occurs. Determine whether this instability is also related to the coarse node spacing required in 2D models of the gravel reach. 9. Investigate the cause of the flat water surface profile upstream of the Harrison River mouth for the design flood simulated by River2D. 10. Conduct field studies to establish the typical values of the Shields parameter for gravel bars, riffles and main and side channels in the gravel reach. This would allow River2D to be used to obtain information on the spatial distribution of grain sizes throughout modeled reaches. To apply the iterative procedure described in this thesis the threshold of entrainment at a particular flow criterion should be investigated further. 11. River2D should be used to determine the spatial variation of hydraulic parameters, including depth, velocity, Froude number, shear stress and unit discharge, in reaches of interest. 12. River2D should be used to estimate the superelevation of the water surface around channel bends. 13. River2D should be used to assess the hydraulic effects of proposed gravel excavations. The model can also be used to determine the volume of gravel that must be removed from the river, and from what particular location, to achieve desired changes in depth and/or velocity. 160 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach River2D should be used in habitat classification studies to provide local depths and velocities around bar-edges. 14. 15. A morphodynamic module for River2D should be developed specifically designed to address the issues in the Fraser River Gravel Reach. Among numerous other potential uses, this module could provide valuable information regarding the possible impacts of in-stream work such as the pilot channel concept and also help to better understand the processes that led to the development of the bank erosion problem at Island 22. In addition to these recommendations, many of the 2D modeling techniques discussed in this report should be adopted as they have been proven to lead to satisfactory simulation results. 161 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach 12. R E F E R E N C E S A . S . C . E . Task Committee on Turbulence Models in Hydraulic Computations, 1988. Turbulence modelling of surface water flow and transport. A.S.C.E. Journal of Hydraulic Engineering, 114(9): 1052-1073. Barnett, A .G . , 1976. Numerical stability in unsteady open channel flow computations. Proceedings of the International Symposium on Unsteady Flows in Open Channels, University of Newcastle-Upon-Tyne, England, 12-15 April, 1976. Berger, R.C. and Stocksill, R.L., 1993. A 2-D numerical model for high velocity channels. Proceedings of the A.S.C. E. Hydraulic Engineering 1993 Conference, 1085-1090, San Francisco, CA, 25-30 July, 1993. Boussinesq, J . , 1877. Essai sur la theorie des eaux courantes. Memoires presentes pars divers savants a I'Academie des Sciences, Paris, France. Bray, D.I., 1982. Flow resistance in gravel-bed rivers. In Gravel-bed Rivers, 109-133, Hey, R.D., Bathurst, J .C. and Thorne, C.R., (Eds.), John Wiley and Sons, Chichester, England. Bovee, K.D., 1982. A guide to stream habitat analysis using the Instream Flow Incremental Methodology. U.S.D.I., Fish and Wildlife Service, Office of Biological Services, Instream Flow Information Paper No. 23, 248 pp. Carling, P.A., 1996. In-stream hydraulics and sediment transport. In River flows and channel forms: 160-184, Petts, G.E. and Calow, P., (Eds.), Blackwell Science, Oxford, England. Carlson, R.E. and Foley, T.A., 1991. Radial basis interpolation methods on track data. Lawrence Livermore National Laboratory, UCRL-JC-1074238. Choi, G., Kim, G. and Ahn, S., 1993. On the applicable ranges of kinematic and diffusion models in open channels. Proceedings of the A.S.C.E. Hydraulic Engineering 1993 Conference, 803-808, San Francisco, CA, 25-30 July, 1993. Chow, V.T., Maidment, D.R. and Mays, L.W., 1998. Applied hydrology. McGraw-Hill, Inc., New York, NY. Chow, V.T., 1959. Open channel hydraulics. McGraw-Hill, Inc., Boston, Mass. 162 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Christison, K.J., Steffler, P.M. and Katopodis, C , 1999. Two-dimensional hydrodynamic modeling of a braided river - a case study. Proceedings of the 3rd International Conference on Ecohydraulics, Salt Lake City, Utah, July 1999. Chu, Y., Metcalf, M., Smith, J.B., Puckette, T. and Nersesian, G.K., 1993. Acoustic Doppler current profiler for inlet current measurements. Proceedings of the A.S.C.E. Hydraulic Engineering 1993 Conference, 1671-1676, San Francisco, CA, 25-30 July, 1993. Church, M., Rempel, L. and Rice, Stephen, 2000. Morphological and habitat classification of the lower Fraser River gravel-bed reach. Prepared for Fisheries and Oceans Canada and the Fraser Basin Council, Department of Geography, University of British Columbia, Vancouver, BC, Church, M., Hassan, M.A. and Wolcott, W.F., 1998. Stabilizing self-organized structures in gravel bed stream channels. Water Resources Research, 34: 3169-3179. Church, M. and McLean, D.G., 1994. Sedimentation in lower Fraser River, British Columbia: implications for management, in The Variability of Large Alluvial Rivers: 221-241, Schumm, S.A. and Winkley, B.R. (eds.), A . S . C . E . , Monograph. Clifford, N.J., Robert, A. and Richards, K.S., 1992. Estimation of flow resistance in gravel-bedded rivers: a physical explanation of the multiplier of roughness length, Earth Surface Processes and Landforms, 17: 111-126. DHI Water and Environment, 2001. Mike 11 - river modeling software. On-line. Internet. 17 May 2001, Available: http://www.dhisoftware.com/mike11/index.htm. Davis, J.M., 1976. The finite element method: an alternative subdomain method for modelling unsteady flow in coastal waters and lakes. Proceedings of the International Symposium on Unsteady Flows in Open Channels, University of Newcastle-Upon-Tyne, England, 12-15 April, 1976. Demuren, A.O., 1993. A numerical model for flow in meandering channels with natural bed topography. Water Resources Research, 29(4): 1269-1277. Engrob, H.G. and von Lany, P.H., 1994. An application of 2-D mathematical modelling on the Brahmaputra River. Proceedings of the 2 n d International Conference on River Flood Hydraulics, York, England, 22-25 March, 1994. 163 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Environment Canada, 2000. Selected water management mathematical models. On-line. Internet. 16 March 2000, Available: http://www.ec.gc.ca/water/index.htm. Fassnacht, S.R., 1997. A multi-channel suspended sediment transport model for the Mackenzie Delta, Northwest Territories. Journal Of Hydrology, 197: 128-145. Fischer, H.B., List, E.J., Koh, R.C.Y., Imberger J . and Brooks, N.H., 1979. Mixing in inland and coastal waters. Academic Press, New York. Fraser Basin Council, 2001. Draft interim Fraser River management plan: Hope to Mission. Vancouver, BC. French, R.H., 1985. Open-channel hydraulics. McGraw-Hill, Inc., New York, NY. Ghanem, A., Steffler, P.M., Hicks, F.E., and Katopodis, C, 1995a. Two-dimensional finite element modeling of flow in aquatic habitats. Water Resources Engineering Report 95-S1, Dept. Of Civil Engineering, University of Alberta, Edmonton, Alta., 189 pp. Ghanem, A., Steffler, P., Hicks, F., and Katopodis, C, 1995b. Dry area treatment for two-dimensional finite element shallow flow modeling. Proceedings of the 12th Canadian Hydrotechnical Conference, Ottawa, ON, 1995. Ghanem, A., Steffler, P., Hicks, F. and Katopodis, C , 1996. Two-dimensional hydraulic simulation of physical habitat conditions in flowing streams. Regulated Rivers: Research and Management, 12: 185-200. Golden Software, Inc., 1999. Surfer 7.0 user's manual. Golden, Colorado. Graf, W.H., 1998. Fluvial hydraulics: flow.and transport processes in channels of simple geometry. John Wiley & Sons, Chichester, England. Henderson, F.M., 1963. Flood waves in prismatic channels. Journal of the Hydraulics Division, A.S.C.E., 89(HY4): 39-67. Henderson, F.M., 1966. Open channel flow. Prentice-Hall, Inc., Upper Saddle River, NJ. Hey, R.D., 1988. Bar form resistance in gravel-bed rivers. Journal of Hydraulic Engineering, 114(12): 1498-1508. Hicks, F.E. and Steffler, P.M., 1992. Characteristic dissipative Galerkin scheme for open-channel flow. Journal of Hydraulic Engineering, 118(2): 337-352. 164 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach HR Wallingford, 1998. Telecmac - a modelling system for hydrodynamics, sediment transport and water quality in the natural environment. On-Line. Internet. 17 March 2000. Available: http://www.hrwallingford.co.uk/software/telemac. Isaaks, E.H. and Srivastava, R.M., 1989. An introduction to applied geostatistics. Oxford University Press, New York. Journel, A .G . , 1989. Fundamentals of geostatistics in five lessons. American Geophyiscal Union, Washington D.C. Lacey, R.W.J., 2001. Application of a 2-dimensional hydrodynamic model for the assessment and design of instream channel restoration works. M.A.Sc. Thesis, Department of Civil Engineering, University of British Columbia. Vancouver, BC. Lane, S.N., 1998. Hydraulic modeling in hydrology and geomorphology: a review of high resolution approaches. Hydrological Processes, 12: 1131-1150. Lane S.N. and Richards, K.S., 1998. High resolution, two-dimensional spatial modelling of flow processes in a multi-thread channel. Hydrological Processes, 12: 1279-1298. McLean, D.G., Church, M. and Tassone, B., 1999. Sediment transport along lower Fraser River: 1. Measurements and hydraulic computations. Water Resources Research, 35(8): 2533-2548. McLean, D.G. and Church, M., 1999. Sediment transport along lower Fraser River: 2. Estimates based on the long-term gravel budget. Water Resources Research, 35(8): 2549-2559. Maghsoudi, N. and Simons, D.B., 1993. Important sources of errors in computational hydraulics. Proceedings of the A.S.C.E. Hydraulic Engineering 1993 Conference, San Francisco, CA, 25-30 July, 1993. Meselhe E.A., Sotiropolous, F. and Holly Jr., F.M., 1994. Numerical simulation of one dimensional transcritical flow. Proceedings of the A.S.C.E Hydraulic Engineering 1994 Conference, 512-516, Buffalo, NY, 1-5 August, 1994. Millar, R.G., 1999. Grain and form resistance in gravel-bed rivers. Journal of Hydraulic Research, 37(3): 303-312. Millar, R.G. and Barua, D.K., 1999. Hydrodynamic model selection study for the lower Fraser River between Laidlaw and Mission. Submitted to the District Of Chilliwack, Chilliwack, BC. 165 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Morse, B., Townsend, R.D., and Sydor, M., 1991. Mathematical modelling of riverbed dynamics - a Canadian case study. Canadian Journal Of Civil Engineering, 18: 772-780. Nikuradse, J . , 1933. Strbmungsgesetze in rauhen Rohren. Forschung auf dem Gebeite des Ingenieur-wesens, Ausgabe B Band 4. (English Translation: Laws of flow in rough pipes. Tech. Memo. 1292, Natl. Adv. Comm. for Aeron., Washinton D .C , 1950). Northwest Hydraulics Consultants Ltd., 1999. Fraser River - bank erosion downstream from Agassiz bridge. Prepared for Cheam Indian Band. North Vancouver, BC. Northwest Hydraulics Consultants Ltd., 2001. Appendix C Fraser River future bank erosion at Island #22. Prepared for City of Chilliwack Engineering Department, Chilliwack, BC. Oleson, K.W., 1992. 2-D mathematical modelling of morphological processes in the Brahmaputra River. Proc. Nile 2000, Cairo, 3-5 February, 1992. Pretegaard, K.L., 1983. Bar resistance in gravel bed streams at bankfull discharge. Water Resources Research, 19(2): 472-476. RD Instruments, 1996. Acoustic Doppler current profiler principles of operation - a practical primer, by Gordon R. Lee. Rennie, C D . , Millar, R.G. and Church, M.A., 2001. Measurement of bedload velocity using an acoustic Doppler current profiler. Submitted to A.S.C.E. Journal of Hydraulic Engineering, January 2001. Rodi, W., 1984. Turbulence models and their application in hydraulics - a state of the art review. 2 n d Ed., IAHR, Delft. Scott, R.D., Williamson, D.C,and Lord, W.S.A., 1999. Hydrodynamic modelling of the Fraser River. 1999 Canadian Coastal Conference, Royal Roads University, Victoria, BC, 19-22 May, 1999. Shrestha, P.L., DeVries, J .J . and Krone, R.B., 1993. Hydrodynamic modeling for channel barrier design. Proceedings of the A.S.C.E. Hydraulic Engineering 1993 Conference, 1079-1084, San Francisco, CA, 25-30 July, 1993. Steffler, P.M. and Jin, Y., 1993. Depth averaged and moment equations for moderately shallow free surface flow. Journal of Hydraulic Research, 31(1): 5-17. 166 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Steffler, P., 2000. River2D: Introduction to depth averaged modeling and user's manual. University Of Alberta, Edmonton, Alta. Strickler, A., 1923. Beitrazo zur Frage der Gerschwindigheits formel und der Rauhigkeirszahlen fur Strome Kanale und Geschlossene Leitungen. Mitteilungen des Eidgenossischer Amtes fur Wasserwirtschaft, Bern, Switzerland. Thompson, D.B., 1993. Numerical methods 101 - convergence of numerical methods. Proceedings of the A.S.C.E. Hydraulic Engineering 1993. Conference, San Francisco, CA, 25-30 July, 1993. UMA Engineering Ltd., 2000. Fraser River gravel reach hydraulic modeling study. Submitted to the City of Chilliwack, Chilliwack, BC. UMA Engineering Ltd with Ward & Associates Ltd. and Quick, M.C., 2000a. Draft report: Fraser and Harrison Rivers hydrologic and hydraulic investigations. Submitted to the City of Chilliwack, Chilliwack, BC. UMA Engineering Ltd., 2001. Assessment of Fraser River erosion at Island 22. Submitted to The City of Chilliwack Engineering Department, Chilliwack, BC. US Army Corps of Engineers, 1991. HEC-2 water surface profiles user's manual. The Hydrologic Engineering Centre, Davis, Ca. US Army Corps of Engineers, 1997. H E C - R A S hydraulic reference manual. The Hydrologic Engineering Centre, Davis, Ca. Viessman Jr., W. and Lewis, G.L., 1996. Introduction to hydrology. 4 t h Ed., HarperCollins College Publishers, New York, NY. Vreugdenhil, C.B., 1989. Computational hydraulics: an introduction. Springer-Verlag, Germany. Waddle, T., Steffler, P., Ghanem, A., Katopodis, C. and Locke, A., 2000. Comparison of one and two-dimensional open channel flow models for a small habitat stream. Rivers, 7(3): 205-220. Walter, R.A. and Casulli V., 1998. A robust, finite element model for hydrostatic surface water flows. Comm. Num. Meth. Engrg., 14: 931-940. Water Management Consultants, 2001. Harrison Bar pilot channel preliminary investigation. Submitted to the City of Chilliwack, Chilliwack, BC. 167 Application of a 2D Hydrodynamic Model to the Fraser River Gravel Reach Water Modelling Section, Environment Canada, 1988. One-dimensional hydrodynamic model: computer manual. Water Planning and Management Branch, Environment Canada, Ottawa, ON. Weiyan, T., 1992. Shallow water hydrodynamics. Elsevier Oceanography Series, 55, Elsevier, Amsterdam. Williamson, J . , 1951. The laws of flow in rough pipes. La Houille Blanche, 6(5): p.738, September-October 1951. Wood, W.L., 1993. Introduction to numerical methods for water resources. Oxford University Press Inc., New York. 168 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0063721/manifest

Comment

Related Items