MODELLING GEOMORPHOLOGY IN LANDSCAPE EVOLUTION by YVONNE MARTIN B A . (Hons.), The University of Western Ontario, 1989 M.Sc., The University of British Columbia, 1991 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Department of Geography) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April 1998 © Yvonne Martin, 1998 In presenting degree freely at this the thesis in partial fulfilment University of British Columbia, available for reference copying of department publication this or thesis by for his of this thesis and or her (jjCQCjr^t The University of British Vancouver, Canada Date DE-6 (2/88) I further Columbia requirements that agree may be It is representatives. for financial the I agree scholarly purposes permission. Department of study. of gain shall not that the an advanced Library shall make it permission for extensive granted by the understood be for allowed that head of my copying or without my written 11 ABSTRACT Many landscape evolution models have considered the interaction of exogenic and endogenic processes. However, geomorphological processes have not been successfully incorporated in landscape evolution models. The thesis begins with a critical analysis of methodologies for the study of large-scale geomorphological processes. A framework based on a generalization of the relevant processes is recommended. Hillslope and channel submodels, which are based on typical processes operating in coastal regions of British Columbia, are introduced. The following hillslope processes are considered: (i) slow, quasi-continuous mass movements; (ii) fast, episodic mass movements; and (iii) weathering. The transport relation for fast, episodic mass movements was found to be nonlinear. Fluvial transport in both low and high-gradient channels and debris flow transport are considered in the channel submodel. A bed load transport equation, which is a revised version of the Bagnold stream power formula, is derived. Suspended load is calculated using a suspended load/contributing area correlation. Connections between hillslope and channel processes are considered to ensure adequate representation in the model. The hillslope and channel submodels are explored in one-dimensional and surface model runs for small drainage basins in the Queen Charlotte Islands, British Columbia. Tests of the fluvial submodel demonstrate the robustness of the bed load equation used in this study. A conceptualization of the landscape into unstable and stable regimes is introduced. Results of surface model runs emphasize the key role of low-order channels in transferring sediment from hillslopes to main channels. The exercise of constructing and running the model highlighted major gaps in our present understanding of geomorphological process operation and sediment routing. Suggestions for future research are extensive and are outlined in the concluding chapter of the thesis. iii T A B L E OF CONTENTS ABSTRACT ii LIST OF TABLES ix LIST OF FIGURES xi ACKNOWLEDGEMENTS xvii 1. INTRODUCTION 1 2. NUMERICAL MODELLING OF LANDSCAPE EVOLUTION: RECONCILING PROCESS AND HISTORY IN GEOMORPHOLOGY 5 2.1 Introduction 5 2.2 Themes in geomorphological research 6 2.3 Reconciling historical and process-based studies 8 2.3.1 History and immanence 8 2.3.2 Scale and process in geomorphology 9 2.4 Scale 11 2.4.1 Spatial scale 11 2.4.2 Temporal scale 15 2.5 Process specification 17 2.6 Drainage initiation versus landscape evolution 19 2.7 Initial surface for subsequent landscape evolution 20 2.8 The role of modelling in landscape evolution studies 21 2.9 History of approaches to landscape modelling 24 2.9.1 Introduction 24 2.9.2 Classic models of landscape evolution 25 2.9.3 Mechanistic models of landscape evolution 27 IV 2.9.4 Generalized physics models of landscape evolution 29 2.9.5 Discussion 33 3. HILLSLOPE SUBMODEL 35 3.1 Introduction 35 3.2 Hillslope processes in past numerical landscape evolution models 36 3.2.1 Models with individual process specification 36 3.2.2 Generalized hillslope transport 37 3.2.3 Discussion of past models 40 3.3 Diffusion in the hillslope submodel 43 3.3.1 Introduction 43 3.3.2 Diffusivity for slow quasi-continuous mass movements 44 3.3.2.1 Examination of past creep data 44 3.3.2.2 Diffusivity value for creep 46 3.3.2.3 Nonlinear relation? 46 3.3.3 Diffusivity for rapid episodic mass movements: a case study of the Queen Charlotte Islands 49 3.3.3.1 Introduction 49 3.3.3.2 The approach 51 3.3.3.3 Diffusivity for the linear relation 52 3.3.3.4 Nonlinear transport relation 57 3.3.4 Comparison of diffusivities with other studies 61 3.3.5 Summary 66 3.4 Weathering processes 68 4. CHANNEL SUBMODEL 4.1 Fluvial transport relations in landscape evolution models 71 71 4.1.1 Introduction 71 4.1.2 Drainage initiation 72 4.1.3 Linear channel elements in a spatial model 73 4.1.4 Stream power transport relations 75 4.1.5 Net change in storage 77 4.1.6 Summary critique 78 4.2 Channel network submodel 79 4.3 Bed load transport 82 4.3.1 Bagnold-type relation 82 4.3.2 Variables required in the Bagnold relation 84 4.3.2.1 Discharge 86 4.3.2.2 River gradient 88 4.3.2.3 Channel pattern 88 4.3.2.4 Hydraulic geometry relation (width/depth) 91 4.3.2.5 Grain size 96 4.4 Bed load transport in steep fluvial reaches 98 4.5 Suspended load 103 4.6 Bedrock erosion 104 4.7 Debris flow activity 106 4.8 Sediment budget framework 111 5. VALLEY FLAT 5.1 Introduction 5.2 Causes of long-term aggradation and degradation 5.3 Valley filling 5.4 Valley incision through fill 5.5 Incorporating valley processes in the present model 5.5.1 "Valley rules" for the one-dimensional model 5.5.2 "Valley rules" for the surface model 5.5.3 Hillslope/fluvial coupling 5.6 Bedrock erosion 6. 1-DIMENSIONAL MODEL RUNS 6.1 Model outline 6.2 Study drainage basin 6.3 Profile development of main channel hillslope 6.3.1 Linear diffusion 6.3.1.1 Basic model runs 6.3.1.2 Change in space step 6.3.2 Nonlinear transport 6.3.2.1 Basic model runs 6.3.2.2 Change in space step 6.4 Valley evolution by processes in the main channel 6.4.1 Basic model runs 6.4.2 Change in space step 6.5 Profile evolution of hillslopes feeding into gully Vll 6.6 Weathering 155 6.6.1 Bedrock lowering rate with no transport 156 6.6.2 Profile evolution 156 6.6.3 Thickness of sediment layer 160 6.6.4 Gully recharge rates for model runs with weathering 167 7. TESTING THE CHANNEL SUBMODEL 7.1 Vedder River 168 168 7.1.1 Channel variables 169 7.1.2 Bed load transport 173 7.1.3 A further test 176 7.1.4 Suspended load 176 7.2 Fraser River 177 7.2.1 Channel variables 179 7.2.2 Bed load transport 179 7.2.3 A further test 181 7.2.4 Suspended load 182 7.3 Mosquito Creek Tributary Basin 182 7.3.1 Channel variables 184 7.3.2 Bed load transport 187 7.3.3 Suspended load 189 7.3.4 Volume changes 189 7.3.5 Debris flow activity 191 Vlll 8. SURFACE MODEL RUNS 194 8.1 Introduction 194 8.2 Study drainage basin 194 8.3 Supply-limited removal along the channel network 196 8.4 Hillslope processes (no weathering) 198 8.5 Hillslope processes (weathering) 202 8.6 Channel processes in the surface model 202 8.7 Conclusions 205 9. CONCLUSIONS 208 9.1 Introduction 208 9.2 Model construction 209 9.3 Model runs 210 9.4 "Unknowns" in the modelling of geomorphology 211 REFERENCES 224 APPENDIX I: REGRESSION ANALYSIS FOR NONLINEAR RELATIONS 239 APPENDIX II: ANALYSIS OF THE BAGNOLD BED LOAD FORMULA 240 II-A Relation between excess stream power and adjusted sediment transport 240 II-B Bagnold's depth and grain size adjustments 243 II-C Re-examination of depth and grain size adjustments 245 II-D Various stream power correlations 247 APPENDIX III: LIST OF SYMBOLS 252 ix LIST OF TABLES Number Title Page Number 2.1 Length and time scales for study units 13 3.1 Creep transport rates 45 3.2 Number of landslides and area for gradient classes, Queen Charlotte Islands, British Columbia 50 3.3 Diffusivities for linear diffusion 58 3.4 Transport/gradient categories and dominant rock types 62 3.5 (a) Comparison with diffusion coefficients derived in scarp studies 62 (b) Comparison with diffusion coefficients implemented in landscape evolution studies 62 4.1 Regression analysis for hydraulic geometry relations 94 4.2 Mean values of coefficients and exponents for sand and gravel rivers 94 4.3 Debris flow activity in unlogged basins in the Queen Charlotte Islands, British Columbia 108 6.1 Input parameters for initial model runs 133 6.2 Gully recharge rates (Oden, 1994) 154 6.3 Gully recharge rates in present study (m /yr) 154 7.1 (a) Volumetric transport rates for the Vedder River (m /yr) 174 Volume changes for the Vedder River (m /yr) 174 Discharge and gradient values and Fraser River cross-sections 178 (b) 7.2 (a) (b) 3 3 Daily mass and volumetric transport rates at Fraser River cross-sections using calculated variables 178 X (c) Daily mass and volumetric transport rates at Fraser River cross-sections using observed variables 178 7.3 Variables in the Mosquito Creek Tributary model runs 183 7.4 Debris flows in the Mosquito Creek Tributary model run 193 8.1 (a) Elevations (hillslope transport only) 199 (b) Changes in elevation (hillslope transport only) 199 (c) Gradients (hillslope transport only) 199 Elevations (hillslope transport and weathering) 203 (b) Changes in elevation (hillslope transport and weathering) 203 (c) Gradients (hillslope transport and weathering) 203 II-1 Functional analysis of it,* vs. excess stream power 242 II-2 Sediment transport vs. depth 246 II-3 Sediment transport vs. grain size 248 II-4 Functional analyses between excess stream power and bed load 8.2 (a) transport adjusted for various depth and grain size factors 250 xi LIST OF FIGURES Number Title Page Number 3.1 (a) Relation between gradient and creep transport rate (Finlayson, 1981) 47 Relation between gradient and creep transport rate (Owens, 1969) 47 3.2 Frequency distribution of volumetric creep rates 48 3.3 Relation between gradient and creep transport rate (VanAsch et al., (b) 1989) 3.4 48 Relation between volume and travel distance for landslides in the Queen Charlotte Islands, British Columbia 53 3.5 (a) Landslide transport rate vs. gradient for Group A drainage basins 54 (b) Landslide transport rate vs. gradient for Group B drainage basins 55 3.6 Frequency distribution of landslide transport rates at unit gradient 59 3.7 Sample plot of nonlinear relation used in analysis 59 3.8 (a) Relation between gradient and transport rate for Group A drainage basins (b) 63 Relation between gradient and transport rate for Group B drainage basins 63 3.9 (a) Transport relation for Group A drainage basins 64 (b) Transport relation for Group B drainage basins 64 3.10 Weathering rate vs. depth 69 4.1 The Bagnold relation 85 4.2 (a) Channel pattern criterion for sand rivers 90 Channel pattern criterion for gravel rivers 90 (b) Xll 4.3 (a) (b) Relation between discharge and width 93 Relation between discharge and depth... 93 4.4 Observed vs. calculated grain size 99 4.5 Wedge gradient vs. channel gradient 4.6 Initiation angle for debris flows in unlogged drainage basins, Queen ; Charlotte Islands, British Columbia 4.7 102 109 Depth of scour in the transport zone for debris flows in unlogged drainage basins, Queen Charlotte Islands, British Columbia 109 4.8 Sediment budget for a fluvial reach 112 5.1 Incision rates below dams 121 5.2 Valley filling in the channel submodel 122 5.3 Fluvial incision in the channel submodel 122 6.1 Flowchart of 1 -dimensional version of model 6.2 Mosquito Creek Tributary drainage basin 131 6.3 (a) Change in profile (diffusivity = 0.1 m /yr) 134 129 (b) Change in profile (diffusivity = 1 m /yr) (c) Change in profile (diffusivity = 10 m /yr) 135 (d) Change in profile (diffusivity =100m/yr) 135 6.4 (a) Comparison of profiles for dx = 50 m and dx = 100 m for initial profile (b) 134 138 Comparison of profiles for dx = 50 m and dx = 100 m at 100 000 years (diffusivity = 0.1 m /yr) 138 6.5 (a) Change in profile (Group A transport relation) 140 (b) Change in profile (Group B transport relation) 140 2 6.6 Comparison of linear versus nonlinear transport relations at 100 000 years 6.7 (a) 142 Comparison of profiles for dx = 50 m and dx = 100 m at 100 000 years (Group A transport relation) (b) 145 Comparison of profiles for dx = 50 m and dx = 100 m at 100 000 years (Group B transport relation) 145 6.8 (a) Fluvial aggradation in the main channel valley 147 (b) Fluvial degradation in the main channel valley 147 (c) Alternating episodes of aggradation and degradation in main channel valley 6.9 (a) 148 Comparison of profiles for dx = 50 m and dx = 100 m at 100 000 years (fluvial aggradation) (b) 150 Comparison of profiles for dx = 50 m and dx = 100 m at 100 000 years (fluvial degradation) 6.10 (a) 150 Change in gully cross-section profile (diffusivity = 0.1 m /yr) 152 (b) Change in gully cross-section profile (diffusivity =1 m /yr) 152 (c) Change in gully cross-section profile (Group A transport relation) 153 (d) Change in gully cross-section profile (Group B transport relation) 153 6.11 (a) 2 Change in gully cross-section profile with weathering (diffusivity = 0.1 .m /yr) 2 (b) 158 Comparison of gully cross-section profiles with and without weathering (diffusivity = 0.1 m /yr) 2 158 6.12 (a) Change in gully cross-section profile with weathering (diffusivity = 1 m /yr) 159 2 (b) Comparison of gully cross-section profiles with and without weathering (diffusivity = 1 m /yr) 159 2 6.13 (a) Change in gully cross-section profile with weathering (Group A transport relation) (b) Comparison of gully cross-section profiles with and without weathering (Group A transport relation) 6.14 (a) 2 2 2 163 Thickness of sediment layer for gully cross-section profile (Group A transport relation) (e) 163 Thickness of sediment layer for gully cross-section profile (diffusivity = 1 m /yr) (d) 163 Thickness of sediment layer for gully cross-section profile (diffusivity = 0.1 m /yr) (c) 162 Thickness of sediment layer for gully cross-section profile (diffusivity = 0.0002 m /yr) (b) 162 Comparison of gully cross-section profiles with and without weathering (Group B transport relation) 6.15 (a) 161 Change in gully cross-section profile with weathering (Group B transport relation) (b) 161 164 Thickness of sediment layer for gully cross-section profile (Group B transport relation)) 164 XV 7.1 (a) Calculated vs. observed channel widths for the Vedder River 171 (b) Calculated vs. observed channel depths for the Vedder River 171 (c) Calculated vs. observed channel grain sizes for the Vedder River 172 (d) Calculated vs. observed channel elevations for the Vedder River 172 7.2 (a) Calculated vs. observed channel widths for the Fraser River 180 (b) Calculated vs. observed channel depths for the Fraser River 180 (c) Calculated vs. observed channel grain sizes for the Fraser River 180 Discharge along main channel of Mosquito Creek Tributary 185 (b) Width along main channel of Mosquito Creek Tributary 185 (c) Depth along main channel of Mosquito Creek Tributary 185 (d) Gradient along main channel of Mosquito Creek Tributary 186 (e) Grain size along main channel of Mosquito Creek Tributary 186 7.4 (a) Unit mass transport rate for bed load along main channel of 7.3 (a) Mosquito Creek Tributary (b) Volumetric transport rate for bed load along main channel of Mosquito Creek Tributary (c) 188 188 Volumetric transport rate for suspended load along main channel of Mosquito Creek Tributary 188 7.5 Volume changes along main channel of Mosquito Creek Tributary 190 8.1 Landrick Creek drainage basin 195 8.2 Initial surface for Landrick Creek drainage basin 197 8.3 Changes in elevation after 100 000 years (hillslope transport only) 200 xvi II-1 Relation between stream power and ib* for Bagnold's bed load transport formula 242 II-2 Relation between depth and ib at constant excess stream power 246 II-3 Relation between grain size and ib at constant excess stream power 248 II-4 Relation between excess stream power and bed load transport 249 II-5 Relation between stream power and ib* for Bagnold-type formula derived in the present study 250 XVll ACKNOWLEDGEMENTS I want to express my most sincere appreciation to my supervisor, Michael Church, for the wise and wonderful guidance he has provided over the years. He has an uncanny ability to find the unique and interesting issues underlying whatever problem one may bring to him. This quality makes him an inspiring supervisor and researcher. He has left an indelible impression on how I view geomorphology, science and academia. I cannot express my gratitude enough. Thank you. I thank Olav Slaymaker for his solid support and enthusiasm for this research. It came at some much needed times and I will always be grateful. A big thank you is extended to Garry Clarke who showed undying patience as we attempted to transform me into a numerical modeller. I would also like to express thanks to Brian Klinkenberg for his friendly assistance as I ventured into the computer lab. Ken Rood provided helpful assistance and materials for the documentation of landsliding in the Queen Charlotte Islands. This work was supported by the Natural Sciences and Engineering Research Council of Canada. I would also like to thank the people who provided friendship and support over the years. A big thanks to Selina Tribe for being an inspirational friend in so many ways. Thank you also to Martin Evans, who day after day listened to my endless tales of woe with this project. Katrina Richards was a true source of support in the final few months when we were both trying desperately to finish our theses. Thank you to Matthias Roth, Darren Ham, Wendy Hales and Magdalena Rucker for the many conversations we have had over the years. A big "maholo" to the gang at dance class whom provided much needed ventures out into the "real world" away from the computers and books. I also want to thank Troy for always believing that this research really mattered and for being my friend for all of these years. Your help during the defense preparation will be remembered for a long time. And last, but certainly not least, a special thank you to my family for always being so wonderfully unique and fantastic - you are an important source of happiness and strength in my life. 1 CHAPTER 1: INTRODUCTION "The ball is now in the geomorphologists' court." -Anderson and Humphrey (1989, p. 350) Anderson and Humphrey made this remark in considering the requirements of modellers studying the combined interactions of tectonics and geomorphology at large scales. The adoption of a numerical modelling approach for the study of combined exogenic and endogenic processes represented a relatively recent development at the time. Researchers were turning to the geomorphological literature to extract information about the operation of surface transport processes at large scales in order to model the response of the lithosphere to erosion and sedimentation. Anderson and Humphrey found that the existing geomorphological research did not meet the requirements necessary for successful implementation in landscape evolution models. Much progress has been made in the modelling of landscape evolution. Geomorphological processes have been incorporated into several landscape models. However, there still remain many unanswered questions about the physical representation of basic transport processes. A methodological framework for large-scale process studies in geomorphology, which is defined clearly, must be adopted by the geomorphological community in order for progress to occur. This would contribute to an increased understanding of transport operation and interactions at large scales, and to an improvement in the representation of geomorphological processes in regional numerical models of landscape evolution. By addressing these issues, this thesis attempts to take on the challenge presented to geomorphologists by Anderson and Humphrey. The study of landscape evolution represents a formidable task for several reasons. The present-day morphology of the earth is a legacy of past geomorphological and tectonic processes, 2 and their interactions. The configuration of any landscape at a particular point in time affects its subsequent development. Because of the historical nature of landscape evolution, past landscapes are not directly observable. Therefore, it is not possible to ascertain past configurations of the landscape, which influence the sequence of following events. Moreover, the operation of processes and their interactions are very complex. W.M. Davis' "cycle of erosion" epitomizes early approaches to landscape evolution (Davis, 1899). He considered developmental sequences of landscapes, which are governed by both geomorphological and tectonic processes, in a conceptual manner. However, as the discipline of geomorphology evolved throughout the twentieth century, geomorphologists began to focus increasingly on the mechanics of process operation at smaller scales. Numerical representations of processes, which are most often based on a mechanically-oriented approach, are much easier to consider at smaller scales. There is greater control over variables, and transport data, which are necessary for the calibration of equations, are easier to procure. In recent years there has been a renewed interest in the study of landscape evolution. Most of the models adopt an approach in which tectonic and geomorphological processes are represented numerically. Unfortunately, the geomorphological rules adopted in these models must rest on a weak foundation. Geomorphologists have not explored adequately the operation of transport processes at large scales. An attempt is made to repair this shortcoming in this thesis. Although models must ultimately strive to integrate geomorphology and tectonics, this thesis focusses on geomorphological processes. For this reason, the operation of geomorphological processes is studied independently of tectonic processes. The landscape configurations resulting from model runs, therefore, cannot truly reflect the diversity of processes that are involved in landscape formation. When the primary focus of a model is the 3 tectonic component, formulations of geomorphological principles can be more generalized versions o f those incorporated into the present model. Landscape changes resulting from both hillslope and channel processes are considered in this thesis. This combination is of fundamental significance in driving landscape change nearly everywhere on the terrestrial surface of the planet. Landscape changes resulting from glacial processes are not considered, despite their obvious importance i n many world regions. Erosional and depositional processes occurring during large-scale glaciations are complex and require much further study before they can be incorporated effectively into models. The overriding objective of this thesis is to construct a model o f landscape evolution in which the resolved geomorphology is constrained by field observations. Small drainage basins in coastal British Columbia represent the "prototype" landscape used i n model development and in subsequent model runs. The suite o f processes considered in the model typifies this region. Calibrations of transport equations are based, when possible, on data collected in coastal regions of British Columbia. Progress i n numerical landscape modelling has been very rapid. For this reason, it prudent at this time to "step back" and explore its foundations. A major objective o f this thesis, which is addressed i n Chapter 2, is to define a framework for the quantitative study of large-scale geomorphological processes within the context o f landscape evolution. In order to achieve this goal, methodological issues surrounding the implementation o f geomorphological processes in landscape evolution models are explored thoroughly. In addition, the strengths arid weaknesses of past approaches are assessed to guide the development of the present model. A second objective o f this thesis is to define a set of equations describing the operation of hillslope and channel processes, and their interactions at the valley bottom, over large spatial and temporal scales. Such equations must reflect adequately the large-scale operation of processes, while remaining computationally feasible. Unlike most previous research, the equations used in 4 the present study are calibrated using field data. Hillslope, channel and valley processes are discussed in Chapters 3, 4 and 5. A third objective of this thesis is to incorporate the hillslope and channel equations into an initial model constrained by observed conditions on real landscapes. The results of model runs, which are based on small drainage basins in the Queen Charlotte Islands, British Columbia are presented in Chapters 6 through 8. A final objective of this thesis is to define clearly further research requirements for largescale geomorphological research. The very act of defining process equations and implementing them in the numerical model has served to highlight many aspects of geomorphological research which require further investigation. These suggestions are outlined in the concluding chapter of the thesis. 5 CHAPTER 2: NUMERICAL MODELLING OF LANDSCAPE EVOLUTION: RECONCILING PROCESS AND HISTORY IN GEOMORPHOLOGY 2.1 INTRODUCTION Methodologies adopted to study landscape evolution have evolved throughout this century in response to dominant themes in geomorphological research and reasoning. In the early 20 century, geomorphological studies became entrenched in geography departments and th as a result earlier concern to consider the effects of processes acting in the earth's interior on geomorphology decreased (Merritts and Ellis, 1994). In recent years, there has been an increasing geophysical interest in landscape modelling and a rebirth of interest in the role of tectonics in landscape evolution (Rind, 1992). This development can be associated with the success of the plate tectonics paradigm. Much recent research was inspired by influential papers by Adams (1980), Molnar and England (1990) and England and Molnar (1990), all of whom investigated the connections between uplift and erosion (Merritts and Ellis, 1994). From a geomorphological perspective, there remain many unanswered questions regarding the nature and rates of geomorphological processes at large scales. Many geomorphological relations should be considered at best tentative because they have been subjected to no rigorous evaluation. For example, a linear relation between gradient and soil creep has been posited and incorporated into several landscape evolution models in the form of a diffusion equation (e.g., Anderson, 1994). This relation has not been demonstrated at the landscape scale in the geomorphological literature - indeed the scanty available evidence appears to contradict it (Kirkby, 1967; Martin and Church, 1997). This type of problem arises frequently in geomorphology because of the difficulties associated with design of detailed observations in 6 the environmental sciences. There still remains significant uncertainty regarding the nature and rates of geomorphological processes at even small and medium scales. Before significant progress can be made in improving geomorphological components of landscape evolution models, it is critical to examine carefully the methodology surrounding the representation of landscape evolution and to ascertain reasonable approaches for its study. An attempt is made herein to outline the major methodological issues involved in the development of a numerical landscape evolution model. A central theme that is developed in this chapter is the necessity of unambiguous scale definition and the appropriate specification of geomorphological processes for the chosen scale. 2.2 T H E M E S IN G E O M O R P H O L O G I C A L R E S E A R C H Geomorphological phenomena have been studied throughout the period during which natural physical mechanisms have been invoked to explain the development of the earth and its associated features. W.M. Davis' (1899) "cycle of erosion" set the stage for the first dominating research theme in the emerging discipline of geomorphology. Starting at the end of the 19 th century and up until about 1950, geomorphological research focussed primarily on large-scale physiographic and historical studies. Many of the researchers from this era whose names are most recognizable in the present day, such as W.M. Davis and W. Penck, were conceptual modellers of landscape evolution. Theories of landscape formation were placed in a Darwinian framework and "evolutionary" sequences of landscape formation were put forward. A notable exception was G.K. Gilbert, whose research has been particularly influential for modern process studies. His work introduced a research paradigm based on Newtonian mechanics. The work of the conceptual modellers provided the template for geomorphological research during this era, which sought to explain the history of landforms (both general and 7 specific) and to elucidate relations between topography and the controlling variables, such as climate and geology, over large temporal and spatial scales. As such, the perspective of geomorphological research was predominantly historical in context. Historical studies focus on the formation of particular landscapes. On the other hand, the conceptual models were attempts to formulate somewhat more general sequences of landscape evolution. However, despite their seeming generality, they were designed to explain patterns of development in particular regions. They did not incorporate the flexibility required to render them generalized landscape evolution models. The models have been subject to frequent criticism for their inability to reflect adequately the nature of landscape evolution in different regions. Herein lies a major dilemma facing large-scale studies. Research conducted at large scales has usually examined in detail the formation of particular landforms, which are not obviously generalizable. Yet a major objective in science is the determination of "generalized" laws. A widely acknowledged shift in the nature of geomorphological investigations began ca. 1950. Modern process studies began to flourish, perhaps partly as a backlash against the "armwaving" geomorphological inquiry that existed up until this time. In addition, this quantitative movement may have occurred in part as a response to the positivist movement in science, which was flourishing at this time. The positivists demanded that science not rest itself on the shaky foundation of theories but on observable entities and "real" correlations. In an attempt to increase the rigor of the discipline and make it more "scientific", there was an increase in process quantification. Even those researchers still conducting research at the large scales during this latter era tended to reject the historical and evolutionary mode of thinking. For example, Hack (1960) proposed that landscape morphology is determined by a modern functional dependence between topography and controlling variables; his is essentially an ahistorical view of landscapes. This focus on process studies was coincident with a decrease in the dominant scales of 8 geomorphological research, which more readily led to quantification, and a new wealth of modern-day process knowledge was obtained. Schumm and Lichty (1965) attempted to reconcile the historical and process-oriented views of geomorphology by formulating the distinctions between studies at different scales. This thesis provides a framework which recognizes both the uniqueness and linkages among studies at all scales. There has been a renewed interest in large-scale geomorphological studies in recent years. Papers by M. Summerfield have documented this trend (e.g., Summerfield and Thomas, 1987; Thomas and Summerfield, 1987). Within geomorphological circles, there has been concern regarding the role of large-scale studies in the context of existing geomorphological research. For example, the particular issue of linking modern processes studies with large-scale landscape evolution was the theme of a recent issue of Earth Surface Processes and Landforms (March, 1997). However, the combined outcome of the papers suggests that there is, as of yet, no generally acknowledged consensus by geomorphologists regarding the resolution of this issue. Indeed, an answer to the issue of scale linkage in geomorphology has been developing, but not within traditional geomorphological circles where the initial discussions occurred. The numerical modelling of landscape evolution, which has been undertaken by researchers in earth science and geology departments (primarily in the U.S.A., but also in Canada and France), may represent the vital key for the resolution of the dilemmas encountered in the study of large-scale geomorphological phenomena. 2.3 RECONCILING HISTORICAL AND PROCESS-BASED STUDIES 2.3.1 HISTORY AND IMMANENCE The historical component of geomorphology does matter at large scales. Evidence of the importance of past events can be detected in many landforms and processes in the present day 9 (e.g., valley fill deposits that represent paraglacial sedimentation; Church and Slaymaker, 1989). It is proposed herein that numerical landscape evolution modelling can be viewed as an attempt to reconcile historical and process studies by "explaining" historical geomorphological landscapes in terms of physically-based processes. At the crux of large-scale modelling are the problems of how to: (i) structure the "laws" governing large-scale processes and (ii) connect the results obtained from the "laws" with specific (historical) landscapes. The historical and process studies, which have dominated geomorphological research respectively, coincide closely with two categories of phenomena discussed by Simpson (1963). He contrasted the character of geological phenomena with those in other physical sciences by assessing the relative roles of: (i) configuration and (ii) immanent processes. The former refers to the changing nature of realised phenomena over long time scales, which prevents definition of basic laws for historical processes. The latter refers to unchanging physical principles as described by the "laws" of nature, such as Newtonian mechanics. An examination of the fundamental nature of generic versus historical phenomena highlights the need to adopt appropriate methodologies at different scales. The determination of immanent laws entails general key relations that exist, all else being equal. Factors other than those under direct consideration must be controlled adequately in order to establish laws. Studies at smaller scales allow greater experimental precision, which increases the possibility for adequate control. Alternative approaches must be found for the study of large-scale phenomena. 2.3.2 SCALE AND PROCESS IN GEOMORPHOLOGY Strict experimental control is rarely obtained in the environmental sciences (Church, 1984). Even at smaller scales of geomorphological inquiry, the variability in controlling variables found in nature often confounds data and relations remain difficult to establish. 10. Despite the limitations for experimental control in the environmental sciences, adequate control of variables to the degree that is possible has not been achieved in geomorphological research. Many studies that claim to be experiments are no more than field measurement programs (Church, 1984). While such studies can provide important insights into the behaviour of a particular phenomenon, they do not lend themselves to the rigorous development of governing process equations. In order to achieve this result, a greater knowledge is required of the relations between controlling variables and the dependent variable. The issue then becomes how to approach the study of geomorphological processes at large scales for which experimental manipulation is effectively impossible. As scale increases, the complexity in the variation of controlling variables increases and only summary measures (e.g., sediment yield) are possible. The results of smaller scale studies cannot be transferred directly to the large scale as the details found in smaller scale process equations are not resolvable at the landscape scale. The importance of configuration (or contingency) increases in the study of landscape evolution. Other things are never equal, which makes the determination of immanent laws problematic. Significant changes in variables such as climate, lithology and tectonics occur which cannot be ignored. As a result there has been minimal study of the actual equations governing transport processes at extended scales. Instead, large-scale studies have focussed on the estimation of erosion and deposition rates based on analyses of landforms and sedimentary deposits. This approach is faced with the difficulties presented by the very nature of transport processes and the accompanying erosion, which are such that the original surfaces are eliminated and only the deposits remain for observation. A critical step in understanding the nature of long-term erosion and deposition is to link deposition sites to particular erosion or source sites. 11 Parameterizations of process at large scales must remain resolvable and yet maintain the essential nature of landscape evolution. In order to meet this objective, it is suggested that a greater emphasis should be placed on defining: (i) the structure of such equations and (ii) the methodology necessary to elucidate the nature and rates of processes at large scales. Exploration of methods for the estimation of transport rates over extended scales could lead to an improved knowledge of large-scale transport behaviour. In the past several years, exciting advances in cosmogenic isotope dating have occurred that increase the potential to estimate long-term process rates (Bierman, 1994). Further exploration of these methods may prove to be a key step in establishing rates of process operation in landscape evolution, which is critical for the calibration of process equations at large scales. 2.4 S C A L E Scales are a set of natural measures that are intrinsic to the system. Careful consideration of scale is necessary as the chosen scale guides the appropriate specification of processes. At smaller scales, greater detail can be resolved in process equations whereas larger scales require an approach in which "generalized" equations are defined. Careful consideration has to be given to scale in order to achieve appropriate representation of processes. 2.4.1 S P A T I A L S C A L E When studying landscape evolution it seems intuitively clear that we are referring to something greater in size than small-scale features such as river bars or local topographic irregularities on hillslopes. Summerfield (1991, p.13) states: "...while it is appropriate to consider a small section of a stream channel in terms of static, steady-state, and perhaps even dynamic equilibrium, it is inappropriate to discuss long-term landscape evolution and the 12 attainment of decay equilibrium in terms of such a restricted spatial scale." In this section, a rigorous approach is adopted in order to define appropriate scales of study for landscape evolution. If scales are based on natural measures in a system, then some defining criterion for a natural measure must be selected. One basis for the delineation of a natural measure is to locate natural breaks in some relevant attribute of the system. In landscape evolution studies, the primary attribute of concern is the morphological change of major landforms which define the landscape. Therefore, natural morphological breaks in the landscape can be used to guide the selection of appropriate scales for its study. Three alternative choices for the scale of landscape evolution studies are explored in this section (table 2.1): (i) tectonic units; (ii) drainage basins; and (iii) hillslopes. Tectonic Units The largest landscape unit for which there may be greater morphological variability between units than within a unit is the "tectonic unit". In many textbooks, the largest scale of landscape categorization is based on tectonically-driven topographic patterns (e.g., Strahler and Strahler, 1988; Summerfield, 1991). Although there are a number of different types of tectonic units, such as orogens located at subduction zones, continental/continental collision zones and rift zones, landscapes that fall under one such category display a reasonable degree of similarity to one another in surface form. Factors such as geology, climate, and length of time since an episode of tectonism was initiated determine the individual characteristics of a particular landscape falling under one category. However, there are gross similarities in both processes and resulting morphology, which allow landscapes to be grouped into categories 13 Table 2.1 Length and time scales for study units. Study Unit Ranges of diameter for study unit (km) Lower: 10 1 Upper: 10 3 Tectonic Unit * Ranges of virtual velocity (km/yr) Time scale (yr) IO'-IO 7 10 -10 9 10°-10 6 10 -10 9 10' - 10° 6 3 Lower: 10° Drainage Basin 10" -10° 6 Upper: 10 3 3 Lower: 10" 10" - 10 1 Hillslope 1 5 10" -10° 6 Upper: 10 1 10 - 10 1 7 * Virtual velocities cover a broad range due to the difficulties associated with the estimation of this parameter. Therefore, the time scales show a great range. In order to ensure that significant landscape changes are observed, the upper ranges of time scales may be the most appropriate for landscape evolution studies. 14 based on tectonism. The tectonic unit represents the fundamental unit for the definition of the balance between uplift and downwearing that exists in evolving landscapes. Several recent numerical models use the tectonic unit as the scale of study. For example, Anderson (1994) examined the evolution of the Santa Cruz Mountains, California. The evolution of rift zones was explored in the models of Tucker and Slingerland (1994) and Kooi and Beaumont (1994). Drainage Basins There is a basic similarity in morphology among drainage basins, making them an appropriate study unit for landscape evolution. Drainage basins are spatial units containing areal and linear pathways for sediment movement within structural/tectonic units. The particular characteristics of a drainage basin vary depending on the order of the basin and the geological and climatic characteristics. However, all drainage basins have the similarity that gradients are arranged so that as water and sediment are routed through the system, all paths are focussed on the drainage outlet. The drainage basin has long been recognized as an appropriate study unit for hydrological research and analyses of sediment delivery rates (e.g., Chorley, 1969). The drainage basin has been a study unit for research which examines the development of the fluvial system over large scales (e.g., Davis, 1899; Schumm, 1977). It has also been the focus of studies in which drainage initiation is investigated (e.g., Dunne, 1980). Hillslopes The smallest scale at which landscape evolution reasonably can be studied is the hillslope scale. All hillslopes have a general morphological similarity by their very definition. They are bounded by a slope base and a crest at the top of the hillslope. Hillslopes are planar or quasiplanar features which when assembled together make up the drainage basin surface. Sediment is 15 transported from the upper portions of the hillslope and is generally deposited along the slope base. Hillslopes have often been the study unit when considering topographic profile , development (Penck, 1953 (translation); Culling, 1960; Kirkby, 1978). Geomorphological studies have made frequent use of the drainage basin as a study unit. This is justifiable if the primary focus of a model is the erosional/sedimentary development of landscapes under the influence of exogenic forces. However, if the tectonic component of a model is relatively more important, then the tectonic unit may be a more appropriate study unit. Most terrestrial geomorphological processes occur within drainage basin boundaries (with the exception of perhaps glacial and aeolian processes), but endogenic processes know no such boundaries. 2.4.2 TEMPORAL SCALE The selection of spatial and temporal scales for a study cannot be made in isolation. As spatial scale increases, the detection of what is considered to be a "resolvable" amount of change requires significantly longer time periods of observation. An approach is herein explored which provides a method for evaluating complementary temporal scales for a particular spatial scale. Time and length are connected to one another via the measure of velocity. Temporal and spatial scales can also be connected to one another in a similar maimer whereby: . . . length scale virtual velocity = time scale (2.1) Virtual velocity denotes the amount of time it takes sediment, the transport of which induces morphological changes in the landscape, to move a certain distance (in this case the length scale). However, as sediment moves through the system, it is not in motion all of the time. In fact, the particle finds itself at rest during most of the journey, undergoing only intermittent periods of 16 actual transport. Virtual velocity refers to the apparent (or average) rate of movement through the system, including the time spent in storage. Equation (2.1) can be rewritten explicitly for time as: , length scale time scale = virtual velocity (2.2) The time scale defined by this equation can be thought of as representing a "cycle" time of sediment through the system. The cycle time is defined here as the characteristic (average) time that it takes sediment to move through the system. This provides some benchmark for defining time scales for a study. The observation of significant morphological changes in the landscape (note that what is considered "significant" depends on the spatial scale of a study) may require the passing of several cycle times of sediment. In particular, when tectonics are incorporated in a model, appropriate time scales may exceed the time scale suggested by equation 2.2. Ranges of time scales that are appropriate for the study of geomorphological processes at each of the spatially defined study units of landscape evolution introduced in the preceding section are presented in table 2.1. This approach for the definition of temporal scales requires a knowledge of virtual velocity. This, in turn, requires an understanding of sediment storage times in the hillslope and fluvial systems. Unfortunately, sediment storage has been a neglected topic in modern process studies. Perhaps this is because the actual movement of sediment is more obvious and "interesting" to study than a state of non-motion However, given the extreme rapidity of many transport processes when they actually occur (such as landsliding, debris flows and fluvial transport during significant flood events), it is the time spent in storage, in between transport events, that ultimately determines how quickly material is evacuated from the system. A major component of historical studies in geomorphology has been the attempt to date sedimentary deposits. Such information can be used to deduce long-term storage times of sediment. 17 However, because of the division that has existed between historically-based and modern process studies, a reconciliation of the relative roles of movement and storage has not been forthcoming. Storage times must be evaluated for sediment residing on hillslopes, the valley flat and in the active channel. Shreve (1979, p. 168) recognized that: "...the characteristic time scales of slope processes typically are orders of magnitude longer than those of channel processes." Creep processes are very slow and landsliding is episodic in space and time. Therefore, sediment may reside on hillslopes for long periods before it is transported to the valley flat. If the sediment then enters the active channel (and remains a part of the active sediment load), it moves quickly through the fluvial system. However, if sediment entering the valley flat does not enter the active fluvial system, it may enter long-term storage. Sediment that enters the active system, but is then deposited during an aggradation phase, may also go into long-term storage. Significant lateral or vertical erosion by the river into the valley fill is required to entrain such material into the active channel system. The issue of sediment storage has been receiving increasing attention in recent years. Sediment budget studies, which consider changes in storage, have provided a framework for the study of sediment storage (e.g., Dietrich and Dunne, 1978). Further research in this direction, with a particular focus on the evaluation of long-term residence times of stored material, is required in order to improve understanding of sediment routing and its associated time scales. 2.5 PROCESS SPECIFICATION Appropriate process specification requires that the level of detail incorporated into process equations is suitable for the particular scale of a study. As the study scale increases, it is not realistic to expect to resolve the same degree of morphological detail as at smaller scales. At the hillslope scale, relatively small changes in morphology are resolvable over time and space. 18 The frequency of sampling or grid points for the evaluation of elevation changes in numerical models can be relatively dense and relatively small changes in elevation may indeed be resolvable. Therefore, the parameterization of processes may include some mechanistic parameters, such as shear stress or stream power. Climate, geology and vegetation may be treated as independent parameters which are constant over time. Some spatial differentiation may be possible for these parameters because of the relatively large sampling frequency that is possible. As the scale of study is increased to that of drainage basins and tectonic units, local irregularities in morphology are no longer important. A high degree of detail in process equations is not suitable to this scale of study as detail focusses attention on small-scale features and events that are relevant only in some integral, or summary, fashion. Furthermore, any attempt to include such details will not be successful as it is impossible to achieve such detailed knowledge of controlling variables over large time and space scales because of difficulties associated with error propagation and the sparse sampling associated with large-scale models. At large scales, the major climate, vegetation and geology changes should be considered, although the large space and time steps in such models necessitate the inclusion of only relatively significant changes. Local variability of key parameters should be averaged out using appropriate techniques. Schumm and Lichty (1965) considered climate to be an independent variable at cyclic time scales. However, more recent research has suggested that there may be very complex feedbacks and interactions between landscape evolution and climate change (Molnar and England, 1990). Furthermore, variables which are considered to be independent at the hillslope scale may now have to be considered dependent variables. For example, at the hillslope scale vegetation may be an independent variable, which may be appropriate for a 19 process study bounded within the slope. As scales increase and the possibility of long-term climate change is introduced, vegetation varies in response to the climatic forcing. The overall objective of process studies at large scales is to ensure that parameterizations of processes replicate the essential character of landscapes as they evolve, without focussing on unresolvable detail. This objective is critical to this thesis and underpins the model that is presented herein. The concept has also been developed in recent numerical models of landscape evolution (e.g., Anderson, 1994; Kooi and Beaumont, 1994; Tucker and Slingerland, 1994). We have to move beyond the restrictions imposed by our own limited range of observation and perceptions of space and time, and consider carefully which phenomena are significant to the formation of landscapes over thousands or millions of years and how to parameterize them appropriately. 2.6 DRAINAGE INITIATION VERSUS LANDSCAPE EVOLUTION The creation of a landscape evolution model requires an initial surface upon which the chosen processes operate. This initial surface must be defined clearly as different specifications of the initial surface lead to studies which address different major questions. Two categories of numerical models of landscape evolution are explored in this section: (i) drainage initiation and (ii) subsequent landscape evolution. The central issue of concern to drainage initiation studies is the determination of where and when channels will begin to emerge (i.e., Ahnert, 1976; Willgoose et al., 1991a,b). The standard approach is to devise a set of quasi-mechanistic relations that define critical conditions required for incision by flow. These equations act on some initial landscape that exists before drainage is initiated. Typical initial landscapes for such studies are planar surfaces and fractal surfaces (Willgoose et al., 1991a,b; Kirkby, 1986). 20 Drainage initiation is supposed to begin with rill development and therefore requires knowledge of the relevant processes at reasonably small scales. In the early stages of rill evolution, a considerably detailed knowledge of hydrological variables and substrate properties must be included in the equations. However, as rills develop and the channel system enlarges, the inclusion of detail which was required to calculate rill incision becomes impractical. The latter stage of landscape evolution is herein referred to as "subsequent landscape evolution". Subsequent landscape evolution refers to changes in the landscape morphology that occur after the fluvial system has developed into an integrated drainage network which is at least approximately stable. Instead of focussing on details associated with rill incision, which necessitates smaller scales of study, attention must now be focussed on the issue of how to parameterize equations that are resolvable at extended scales. It is a weakness of some landscape evolution studies that both of these stages of landscape evolution are incorporated into one model (Ahnert, 1976; Willgoose et al., 1991; Rinaldo et al., 1995). The contrasting scales of drainage initiation mechanisms and subsequent landscape evolution suggest that such studies should be approached as two separate phenomena. Attempting to incorporate these issues into one model will lead to parameterizations of processes which may not be appropriate at either the early or later stages of landscape evolution. 2.7 INITIAL SURFACE FOR SUBSEQUENT LANDSCAPE EVOLUTION What is a suitable "beginning" point for the study of subsequent landscape evolution? Summerfield (1991) states that, except in a few rare situations, drainage development does not occur on "pure" surfaces. Pre-existing drainage and topography almost always exist for any point in time at which we choose to enter the system for study. Summerfield (1991, p. 412) summarizes this point eloquently: ".. .drainage systems have a heritage rather than an origin." 21 For the purposes of the current research it is not necessary to define a surface on which to initiate drainage. Since there is no one correct beginning point of subsequent landscape evolution then the modeller must make a decision as to what initial configuration lends itself best to the study of the particular aspects of landscape evolution that are being studied. The decision should be justified in relation to the aims of the research. For landscape evolution studies it may be preferable to set the initial surface so that it approximates a real landscape as opposed to some artificial surface. A real landscape does not necessarily have to be defined in the strict sense of the term. An alternative to using a real landscape (data for which can be obtained from a DEM) as the initial surface is to create a landscape that displays essential characteristics of real landscapes, as perhaps summarized in measures such as hypsometric integrals, fractal dimensions and drainage density. 2.8 T H E R O L E OF MODELLING IN LANDSCAPE EVOLUTION STUDIES Process can be studied in a quantitative manner at any scale if the process specification is adapted to suit the particular study. Past approaches to geomorphology have not generally been concerned with parameterizing processes at the large scales dictated by studies of landscape evolution. Although a few pioneering attempts were made to establish process relations at large scales over the past several decades (e.g., Culling, 1960; Scheidegger, 1970; Hirano, 1975), such work has not been in the mainstream of geomorphological research. The calculation of process equations at the large scales of landscape evolution studies requires a computing capability that was not possible until relatively recently. The numerical modelling approach has shown itself to be a flexible and useful framework around which to focus the study of landscape evolution. historical studies. This approach promises to reconcile process and 22 How can numerical models be used to study landscape evolution? The objective of a modelling exercise should not be to attempt to replicate exactly the details of development of a particular landscape. Any attempt to do so is bound to meet with failure as the boundary conditions and variability of controlling variables cannot practically be reconstructed. In order to define the role of numerical modelling in landscape evolution studies, the issues underlying model verification, validation, calibration and confirmation as discussed by Oreskes et al. (1994) are explored: 1. Model verification is the process of determining the "truth" of the model and its reliability. Earth science models are generally considered to be "open" due to an incomplete knowledge of input parameters. The models are effectively underspecified. This occurs because of: (i) spatial averaging of input parameters found in models; (ii) nonadditive properties of input parameters; and (iii) inferences and assumptions underlying model construction. This "incompleteness of information" means that model verification is not possible in open systems. 2. Model validation (a term that is often confused with verification) requires that a model contains no known or detectable flaws and is internally consistent. Earth science models should not be considered as definitive representations of physical reality due to the complexity of the systems being studied. 3. Model calibration involves the manipulation of parameters to improve the degree of correspondence between the simulated and observed results. However, obtaining consistent results does not imply a model's representation of reality is "correct". Oreskes et al. (1994) refer to the calibration procedure as "forced empirical adequacy". Moreover, just because a model has been forced to fit a certain set of data, it may not perform adequately for data collected at another time or location. 23 4. Model confirmation can be established by examining the ability of a model to match prediction with observation. However, the model performance can never be definitively confirmed because of our inability to know that our model will always predict the correct results - future testing may prove the model to be inadequate. Furthermore, more than one model may predict the correct results. Greene (1997) warns of the dangers of accepting constraints in model construction as legitimate constraints on theory. If models, by their very nature, are going to be imperfect representations of reality what, then, can numerical modelling contribute to the study of landscape evolution? Even if modelling exercises cannot be used to absolutely confirm or validate the ideas contained within, then they can nevertheless be useful tools for the investigation of scientific phenomena. Landscape evolution models can be used to support or more thoroughly explore ideas and hypotheses that have been partly established in other ways (Oreskes et al., 1994). Transport relations, which may have been shown to provide reasonable results, either in the field or experimentally in the laboratory, can be explored more fully in the model. Various controlling variables can be held constant, while others are allowed to vary and the implications of such manipulations over long time scales can be provisionally assessed. Perhaps the most important role of models is as a tool for the exploration of various "what-if' questions (Oreskes et al., 1994). Sensitivity analyses can be performed by changing the nature or intensity of various processes and observing the effects on the morphological evolution of landscapes. In addition, numerical landscape models can be used to explore theoretical ideas and conceptual models about which there is much conjecture, but little quantitative research. For example, Kooi and Beaumont (1996) recently explored the ideas of Davis and other classic modellers using their numerical model of landscape evolution. 24 Numerical modelling permits the exploration of joint or sequential action by several processes (geomorphological and/or tectonic) in a distributed field. This exercise usually confounds our analytical abilities and often our intuition. Hence, a modelling exercise provides an approach for assessing the plausibility of ideas about historical landscapes. This is not to say that landscape evolution studies should be restricted to numerical modelling exercises. The ideas contained within such models should be based on ideas, theories and equations that are derived from real-world observations. In order for significant progress to be made in landscape evolution studies both numerical modelling and field approaches should be integrated and be used to provide new insights for one another. 2.9 HISTORY OF APPROACHES TO LANDSCAPE MODELLING 2.9.1 INTRODUCTION There are four main factors that influence the nature of landscape evolution models created at a certain time period: (i) dominant themes of research in the discipline; (ii) available knowledge (i.e., existing research); (iii) available technology; and (iv) background and beliefs of the researcher. Three main groups of landscape evolution studies are now considered: (i) classic models; (ii) mechanistic models; and (iii) models incorporating generalized physics . Until the mid-20 century geomorphological research was largely conducted at larger th scales and a focus of much research was the development of conceptual landscape evolution models. W.M. Davis, W. Penck and L.C. King were, in effect, landscape modellers and therefore a review of landscape evolution modelling returns us back in time about a century. It is interesting to note that although there has been a significant amount of research conducted at the large scale, very few actual models of landscape evolution have been proposed. This may be because most research of large-scale earth history has been undertaken within a 25 purely historical paradigm. This tendency may be further exacerbated by the paucity of direct observations, thereby often necessitating space/time substitutions, and by the complexities encountered in large-scale geomorphological studies. The classic conceptual models represent the initial attempts by geomorphologists to explain how landscapes evolve. Throughout the latter half of the twentieth century geomorphologists made use of the new quantitative techniques and technology available to them and as a result several quasi-mechanistic and generalized physics models were created. Numerical modelling is the timely attempt to reconcile classic landscape evolution and historical studies with modern process studies. 2.9.2 CLASSIC MODELS OF LANDSCAPE EVOLUTION The classic landscape evolution models were created during the period 1890 to 1950 (Davis, 1899; Penck, 1953 (translation); King, 1962; Budel, 1982 (translation)). During this time most research focussed on historical and regional studies, making this a fruitful period for the creation of large-scale landscape evolution models. The large scale of study allowed for the incorporation of internal earth processes (e.g., uplift) into most of the models. Due to the fact that the entire modern approach to earth sciences was in its early stages, minimal exogenic process data and research existed, although such studies increased in number after the turn of the century. Endogenic knowledge consisted of such hypotheses as isostatic principles and the geosynclinal theory of sedimentary basin and mountain evolution. Technology posed a significant constraint. Equipment and techniques available at this time were considerably less advanced than those available in more recent years. The early geomorphologists were all acute "observers" of landscapes. In effect, they were studying the "observable" end products of the long-term processes involved in landscape 26 evolution. These researchers turned to space-time substitution in an attempt to make conjectures about the unobserved stages of landscape development. The most well-known of all landscape evolution models is that of W.M. Davis (1899). He proposed that, after a period of brief and episodic uplift, landscapes underwent downwearing and passed through a series of predictable stages referred to as the cycle of erosion. In later writings, Davis discussed various factors, such as climate change and renewed uplift, which might complicate this simple model. Penck (1953, translation) rejected the notion of disparate uplift and erosional events and instead focussed on their continuous interaction. Penck advocated the concept of slope replacement whereby the steep part of a slope retreats rapidly and leaves behind a lower angle debris pile at its base. King (1962), like Davis, believed that uplift is episodic. He proposed that slopes undergo parallel retreat and leave behind concave pediments at their base which eventually coalesce to form a pediplain. Budel (1982, translation) introduced the concept of etchplanation, which described the relation between the weathering mantle and removal of material. He asserted that tectonic and climatic stability result in the equality of weathering and denudation rates and the stability of weathering mantle depth. This stability can be offset by changes in controlling variables which strip the weathering mantle, thereby creating a situation of disequilibrium. An important observation to make about these classic models is that they are all qualitatively conceptual in nature and that, with the exception of Budel, both exogenic and endogenic processes are considered. Perhaps the most frequent criticism of these models is that exogenic processes are treated in a superficial and non-quantitative manner. This criticism is made in light of our increased understanding of process in recent years. These models also show a progressive tendency over time to focus on the hillslope system at the expense of the fluvial system. This perhaps reflects the increasing realization that hillslopes are the fundamental area! 27 unit of landscapes. In addition, interest in hillslope evolution may be further stimulated by the conundrum that although hillslope operation rates are slow, these processes seemingly lead to major morphological changes over the long term. The resolution of this enigma poses an interesting research challenge. 2.9.3 MECHANISTIC MODELS OF LANDSCAPE EVOLUTION The landscape evolution models of Ahnert (1976 ) and Kirkby (1976) are based on quasimechanistic process equations. The dominant themes of research during the period when these models were created and revised had changed enormously from those encountered by the classic landscape modellers: (i) the dominant scale of study decreased to small/medium scales; (ii) process quantification was a critical component of much geomorphological research; (iii) the degree of specialization of researchers increased; and (iv) endogenic processes came to be largely ignored by geomorphologists (this development is related to (i) and (iii)). These changes coincided with the advent of the so-called "quantitative revolution" in geomorphology. The exogenic process information available to landscape modellers during this era increased at a rapid rate whilst plate tectonics theory was changing the nature of endogenic studies. Unfortunately, at the time that this new wealth of endogenic process information became available, geomorphologists' interests focussed on smaller scale studies in which endogenic processes were largely ignored. Ahnert (1976, p.31) described the basic outline of his model: "...the model starts with identification of the initial surface, followed by bedrock weathering (i.e., waste production) base level lowering, if any, denudation processes, and finally recomputation of the resulting surface and its various parameters (slope etc.). This resulting surface then becomes the "initial surface" for the next passage through the sequence of processes." The processes included in this model 28 are weathering of bedrock and the movement of sediment by rainsplash, viscous flow, plastic flow, and wash. A computer program was created for the model runs, increasing the efficiency of process calculations. Uplift is not directly modelled; it is indirectly incorporated into the model through its effect on base level lowering. The model can be run in profile or surface form. The calibration of process equations was not discussed, suggesting that results of the basic model cannot be compared directly to real data. The model of Kirkby (1971) involves the calculation of the sediment transport rate for a variety of processes (slope wash, rainsplash etc.). These results are inserted into the continuity equation. Soil profile evolution and rapid mass wasting events are explored in his later modelling efforts. Endogenic processes are not included in this model. Kirkby kept the mathematics in the model tractable by focussing on the evolution of landscape profiles. Once again, the results have arbitrary units. The models of Ahnert and Kirkby are: (i) quantitative in nature; (ii) focus on exogenic processes; and (iii) have inadequate scale definition. These models introduced a quantitative sophistication that was not possible during the earlier period in which the classic models were created. Despite the obvious attraction of their quantitative rigor, it should not be overlooked that both of the models demonstrate a lack of adequate scale definition. The detail of the process equations and the overall neglect of endogenic processes show that these models are not suitable for very long time scales. Despite the detailed process specification found in his model, Ahnert applied the model to a variety of scales, ranging from individual point locations through to entire mountain ranges. It is not reasonable to apply this one representation of physics to such a large range of scales, as the parameterization must vary with scale. However, Kirkby (1971, p. 15-16) did demonstrate an awareness of the scale issue: "In considering the evolution of slope profiles through time, we are necessarily considering a system in 'cyclic time'...our specification of 29 variation in process must be in terms appropriate to cyclic time; that is to say, not in terms of hydraulic variables but in terms of relief variables only, even if the process is hydraulic in nature." Nonetheless, his process equations show an amount of detail which is not suitable for the cyclic time scale as defined by Schumm and Lichty (1965). 2.9.4 GENERALIZED PHYSICS MODELS OF LANDSCAPE EVOLUTION Generalized physics models have been created over a period starting about 1960 and lasting up until the present day. These models represent an attempt to overcome the constraints of mechanistic modelling. The earlier modellers of the genre (Culling, 1960; Hirano, 1975; Smith and Bremerton, 1972) were subject to many of the same intellectual and technological constraints as the mechanistic modellers. As a result they share some similarities, such as their focus on exogenic processes. Around the mid-1980's there was a renewed interest in large-scale studies accompanied by an interest in re-examining how processes are dealt with at these scales (Anderson and Humphrey, 1989). In addition, there has been an increasing recognition of the importance of endogenic and exogenic process interaction (e.g., Molnar and England, 1990). The increasing geophysical interest in landscape evolution may be an outgrowth of the maturing of plate tectonic theory. A notable characteristic of the most recent modellers in this category is that many of them are working in earth science departments in North America and, hence, this may restrict interaction between these modellers and European geomorphologists, who often are working in geography departments. A recurring theme in the description of the models that follows is the frequent use of diffusion to simulate hillslope processes. The diffusion equation is derived from a transport/gradient relation, which is inserted into a basic mass continuity equation. This equation constitutes a generalized physics approach in the sense that processes, which may actually be 30 quite variable in rate or occurrence, are assumed to act "continuously" over appropriately long time intervals. In addition, it is assumed that the details of the actual processes, such as exact initiation and deposition sites and the resulting morphometry of the depositional landforms produced by such events, "average" out in a manner that can be represented by one simple relation at large scales of space and time. Linear diffusion has been used in many landscape evolution models. In this form, a linear relation between sediment transport rate and gradient is assumed. In the less commonly used nonlinear equation, which is somewhat more complicated in structure, the rate of activity significantly increases at steeper gradients. Culling (1960, 1963, 1965) and Hirano (1975) recognized the potential of using a generalized transport relation, based on a linear diffusion model, to simulate transport processes at large scales. Culling (1960), working at the very beginning of the era of scientific computing, had to solve the equations analytically, making them cumbersome to manipulate and limiting the complexity that could be incorporated into his model. The model of Smith and Bretherton (1972) is based on a modified version of the diffusion equation. The particular strength of this model is the recognition that drainage initiation and subsequent landscape evolution represent two distinct stages in landscape evolution. Flemings and Jordan (1989) made use of linear diffusion in their model of foreland basin development. In this model, the primary concern is to establish realistic rates of erosion over long time scales in order to simulate the development of basin geometry and stratigraphy. Crustal shortening and lithospheric adjustment are also incorporated into the model. Anderson and Humphrey (1989), in response to the work of Flemings and Jordan (1989) and several other similar models, made a plea for geomorphologists to consider methods for assessing sediment delivery rates, and hence transport processes, at large scales. Anderson and Humphrey (1989) assessed the roles of linear diffusion and weathering in landscape evolution in their model runs. 31 During this same year, Koons (1989) introduced a surface model that also used the linear diffusion to simulate hillslope processes. Significant to his research are the lateral and vertical variations in the diffusion coefficient that are introduced to reflect changes in process rates over space and time that occur in response to changes in controlling variables. Stream elevations are calculated using an exponential function that is fitted to the hillslope morphology. Uplift can vary over both space and time in the model. The model of Willgoose et al. (1991a,b) includes hillslope, fluvial and endogenic processes. In particular, the model was used to simulate the initiation and development of drainage networks over time (see section 2.6). Hillslope and fluvial processes were simulated using linear diffusion and a stream power relation respectively. According to Chase (1992), models based on diffusion processes alone are insufficient to explain landscape evolution. He asserted that diffusion smoothes landscapes, whereas roughening of the land surface also occurs. Chase believed that, although individual processes are nonlinear and complicated, their synoptic effect can be reduced to simple, nearly linear laws. The following exogenic processes are included in the model: (i) diffusive smoothing weathering, slope wash, soil creep etc.; (ii) fluvial erosion - suspended and/or bed load; and (iii) deposition. Tectonic processes are simulated by the vertical uplift of grid cells. Several landscape evolution models were published in a special issue of Journal of Geophysical Research (1994) dedicated to the topic of tectonics and topography. These models, as well as several other publications, represent the state-of-the-art for landscape evolution modelling. The model of Howard (1994) includes both hillslope and fluvial processes. Slow mass movement and failures are modelled using linear diffusion and a threshold-driven equation 32 respectively. Both non-alluvial and alluvial river processes are incorporated into the model using complex variants of basic stream power relations. Rigon et al. (1994) studied the self-organized nature of channel networks by examining processes leading to drainage initiation. Linear diffusion is used to simulate hillslope processes while fluvial erosion occurs when a threshold shear stress is exceeded. The model of Anderson (1994) includes both tectonic and surface transport processes and is used to explore the evolution of the Santa Cruz Mountains, California. Slow mass movement is simulated using linear diffusion and landsliding is simulated using a modified equation that includes a critical angle for failure. The fluvial model consists of bedrock incision, which is modelled using a stream power function. Sensitivity analyses were performed to explore various relations and model performance. Rosenbloom and Anderson (1994) presented a similar version of this model in which soil creep is the dominant mass movement process. The numerical landscape evolution model of Kooi and Beaumont (1994) was used to explore the evolution of escarpments on rifted margins. Hillslope processes (including soil creep, rainsplash, earth flows, slides and rockfall) are simulated using linear diffusion. Stream power relations are used to simulate fluvial transport and redeposition along the river. Flexural isostasy is also included in the model. In a subsequent paper (Kooi and Beaumont, 1996), the model was used to explore some of the classical conceptual models of landscape evolution (Davis, 1899; King; 1962; Penck, 1953 (translation)). Tucker and Slingerland (1994) used linear diffusion to model slow mass movement and a threshold-driven algorithm to simulate fast mass movement. A stream power relation was used to simulate fluvial transport. Interactions between flexural isostatic uplift and geomorphological processes in escarpment retreat were examined. 33 Most recently, Avouac and Burov (1996) created a landscape evolution model to examine the role of erosion in driving intracontinental mountain growth. They hypothesized that removal of material from areas of high elevation to forelands opposes the spreading of the crustal root that would otherwise lead to the eventual collapse of the mountain range. Of particular interest is the fact that hillslope erosion is simulated using both linear and nonlinear diffusion. 2.9.5 DISCUSSION The approach to process specification followed in the most recent suite of landscape evolution models demonstrates an appreciation of the need to generalize process equations at large scales. In particular, many of these models use some variant of: (i) linear diffusion to simulate slow and/or fast mass movement; (ii) a threshold-based equation to simulate fast mass movement; and (iii) a stream power relation to simulate fluvial processes. The widespread use of these generalized equations is an important step towards establishing a sound methodological framework for landscape evolution studies. However, there are two concerns regarding this most recent set of numerical models: (i) geomorphological relations used in the models are often not evaluated rigorously in comparison with field data and (ii) the critical roles of sediment storage (particularly in valleys) and river incision through this fill are not recognized. However, it may be appropriate, when considering tectonic histories over very large time scales, to ignore valley storage. The diffusion analogy, which was often used to simulate various hillslope processes in past models, was not compared or tested rigorously against field data in order to evaluate its performance. Without this step in the analysis, the suitability of the diffusion concept to model hillslope processes remains questionable. Furthermore, several potential variants of diffusion have been proposed and research is needed to see which is most appropriate. Likewise, stream 34 power relations used in the various models have not been tested against data to evaluate whether reasonable results are produced. The final form of the geomorphological equations must be shown to behave in a manner which is reasonably consistent with the available field evidence. The importance of the connections between hillslope and fluvial processes is not recognized in existing landscape evolution models. The storage of sediment in valleys and the incision of channels through these sediments and the underlying bedrock are all key elements in determining the morphology of landscapes at time scales over which geomorphological processes are important (about < 10 years). The features produced by these processes represent the most 6 interesting elements of landscapes of geomorphological significance at intermediate to large time scales. These topics have been neglected in the geomorphological literature and in landscape evolution modelling and must be addressed in order to simulate realistically the routing of sediment through drainage basins. 35 CHAPTER 3: HILLSLOPE SUBMODEL 3.1 INTRODUCTION The modelling of hillslope processes at landscape scales requires resolution over relatively large temporal and spatial units. When model time scales are order of magnitude 10 4 to 10 years, the corresponding model time steps usually range from decades to thousands of 6 years. The modelling of large regions from 10 to 10 km requires model grid cell dimensions 1 3 2 ranging from several tens of m through to several km . Process parameterization at large scales 2 2 must lead to a computationally feasible model while still maintaining critical characteristics of evolving landscapes. The diffusion analogy was introduced into geomorphological reasoning for this purpose (see Nash, 1980ab for early references). Linear diffusion has been used subsequently to model the development of scarps (Nash, 1980a,b; Colman and Watson, 1984; Hanks et al., 1984). In recent years, diffusion has been employed in large-scale landscape evolution models to simulate slope evolution over long periods. The linear diffusion equation is derived from a statement of sediment continuity: 8h dt dx (3.1) dy and a sediment transport relation: -k\ q =-k' v dh (3.2a) dx dh (3.2b) dy wherein h is height, t is time, x andy are spatial dimensions, q and q are the x and>> components x y of the volumetric transport rate ( L ^ T " ) and k is & diffusion coefficient (L T ). 1 equations are combined to form the diffusion equation: 2 _1 These 36 dh dh 2 dh 2 (3.3) ~dt The equation assumes transport-limited removal of material from the slope (which is a standard assumption in many long-term models of landscape evolution). The diffusion analogy is attractive for modelling slope development in landscape evolution models because it eliminates mechanistic details that are resolvable only at smaller scales. It is flexible because the diffusion coefficient can be modified to reflect changes in space and/or time of the controlling variables. However, it remains to be determined if critical features of landscapes are simulated when using this equation. 3.2 HILLSLOPE PROCESSES IN PAST NUMERICAL LANDSCAPE EVOLUTION MODELS 3.2.1 MODELS WITH INDIVIDUAL PROCESS SPECIFICATION The early numerical landscape models of Ahnert (1976, 1977, 1987a, 1987b, 1988) and Kirkby (1971, 1976) specify individual equations for different hillslope processes. The model of Ahnert simulates morphological changes resulting from weathering, rainsplash, overland flow, plastic flow, viscous flow and debris slides. Transport rate is dependent on gradient in all of the equations. However, the relation for each process is defined by a unique assemblage of coefficients, exponents and other governing parameters. The empirical relations, which are based on theoretical-mechanistic principles, are considerably more complicated in structure than the diffusion equation. The model of Kirkby (1971, 1976) examines the effects of various transport processes on hillslope profile evolution. Hillslope processes included in earlier versions of his model are soil 37 creep, rainsplash and soil wash. The general form of the transport relation for each of these processes is: CocA S" m (3.4) where C is transporting capacity of the process, A is contributing area (surrogate for water discharge), S is slope, and m and n are exponents that depend on the process being modelled. For creep and rainsplash the value of m is 0 whereas it is greater than 1 for soil wash. Transport rates are inserted into a continuity equation to obtain changes in elevation along the profile. Other factors considered by Kirkby in various papers over the years include the interactions between the soil profile and the slope profile (Kirkby, 1977) and the evolution of hillslopes by mass movements (Kirkby, 1987). Considerable detail is introduced in both of these later additions, thus limiting their applicability for modelling long-term hillslope evolution. In the foregoing studies, a dependence on slope is incorporated into the hillslope transport relations. The equations are then evaluated within a continuity framework. 3.2.2 GENERALIZED HILLSLOPE TRANSPORT A diffusion analogy for modelling hillslope processes was used by Culling (1960) and Hirano (1975). Diffusion modelling was expanded in some recent models to include explicitly episodic processes. Diffusion was now used to model the combined effects of several hillslope processes including slow, quasi-continuous processes (e.g., creep) and fast, episodic processes (e.g., landsliding, rock slides) (Koons, 1989; Willgoose et al., 1991a,b; Chase, 1992; Kooi and Beaumont, 1994, 1996; Rigon et al., 1994). Rigon et al. (1994) discussed "hillslope processes" in a general manner, and did not explicitly consider the individual processes that are simulated using the diffusion analogy in their model. However, it seems reasonable to suppose that diffusion is representing both slow and rapid mass movements. 38 Kooi and Beaumont (1994, 1996) state explicitly that although in the short term these processes behave differently because of variations in threshold angles and transport intensity, in the long term these processes all transport material relatively short distances at a rate that is dependent on local slope. Assumptions in their model are that (i) changes in porosity and (ii) the removal of material by solution are both negligible. A notable feature of the models of Koons (1989) and Kooi and Beaumont (1994) is that the value of the diffusion coefficient can be changed to reflect variations in climate over space and time. Kooi and Beaumont (1994) allowed the strength of diffusion to be modified to simulate transport occurring on either loose sediment or bedrock. Chase included slope wash transport among the processes modelled by diffusion, whereas Willgoose introduced a slope dependency into the basic equation for slope wash. Anderson and Humphrey (1989) and Tucker and Slingerland (1994) studied the combined effects of weathering and sediment transport on hillslope profile evolution. Weathering rates are calculated using a relation whereby the weathering rate declines with increasing soil thickness. Diffusion is used to simulate slow, quasi-continuous hillslope transport processes in both of these models. Sediment flux cannot exceed sediment supply. In the model of Tucker and Slingerland (1994), rapid mass movements (including shallow debris landsliding and rock mass failure) occur when slopes are steeper than some specified threshold angle which is defined uniquely for each of sediment and bedrock. The material moves in a downslope direction until a stable gradient is established. oversteepened slopes. This procedure is repeated during a time step until there are no Tucker and Slingerland stated explicitly that their weathering model involves an assumption that weathering is isovolumetric, such that any change in density occurring when rock is converted to sediment is compensated for by the removal of mass in solution. 39 Rosenbloom and Anderson (1994) used a diffusion equation to model slow mass movements operating on marine terraces near Santa Cruz, California. Field observations indicate that the burrowing activity of ground squirrels and moles is the primary transport process operating on these grassy slopes in the present day. Slow and fast mass movements are treated as two additive terms in the model of Anderson (1994). Anderson (1994) expanded the model of Rosenbloom and Anderson (1994) to include the effects of slow, quasi-continuous processes (e.g., creep), in addition to more discrete processes such as tree throw and rodent burrowing. Moreover, Anderson (1994) also incorporated a term to account for the effects of landsliding on hillslope morphology: (3.5) exp where q mass is mass sediment discharge per unit width of slope, S is local slope, S is critical c slope, S* is a slope scale which determines the amount of difference between the actual slope and critical slope which is required for a significant increase in failure, kd is the normal diffusivity (for slower processes) and ki is the maximum effective diffusivity for landslides. s Howard (1994) also included a hillslope transport equation that consists of two additive terms: q K G(S) + K S f [\-K (abs(S)")] (3.6) x where q is the rate of movement of hillslope material, S is the local slope gradient, G(S) is an increasing function of slope gradient and s is the unit vector in the direction of S. The constants K , Kf, K and the exponent a are spatially and temporally invariant. Creep and rainsplash s x processes, which are represented in the first term, are modelled using linear diffusion. Mass 40 movement rates resulting from larger failures, represented in the second term, increase without limit as a threshold gradient is approached. In the numerical model of Avouac and Burov (1996), which is used to examine the effects of geomorphological processes on intracontinental growth, both linear and nonlinear diffusion are used to model hillslope processes. However, they pointed out in the study that their form of the nonlinear equation does not conserve mass and to do so would require an additional term. This shortcoming presents a significant problem as the equations used in a landscape model should conserve mass. 3.2.3 DISCUSSION OF PAST MODELS If such a mechanistic approach is chosen for a smaller hillslope study, consideration must be given to the relative efficacy of the various processes. Some processes operate at significantly faster rates than others, thereby eclipsing the effects of slower processes. The slower processes can be eliminated from the model in this event. Moreover, some processes may not even operate in some regions. A thorough examination of the relative importance of these processes in different landscapes is an essential step in the application of this model. While mechanistic approaches (e.g., Kirkby, 1971; Ahnert, 1977) may be reasonable for studies at the scale of relatively small hillslopes, this approach cannot be applied practically to large-scale studies. As scale increases the net effect of these various processes may be combined into a more generalized gradient-driven equation. Further examination is required to determine at what scales generalized transport models may replace more mechanistic process models. The sequence of studies outlined in the previous section demonstrates a progression in broadening the range of processes modelled using diffusion and recognizing that adjustments are needed when applying this equation to fast mass movement processes. 41 i The diffusion concept is used to simulate slow mass movements exclusively in some cases (e.g., Anderson, 1994; Howard, 1994; Tucker and Slingerland, 1994), while in other models this equation is used to simulate both slow and rapid mass movements (e.g., Koons, 1989; Kooi and Beaumont, 1994). The summary effects of various processes over long time periods may be such that they can be treated uniformly. All of the processes are assumed to be dependent on gradient. However, the use of a generalized equation requires study over minimum time scales for which both the quasi-continuous processes (e.g., creep) and episodic processes (e.g., landsliding) can be considered continuous. It may take only decades to homogenize the effects of slow processes, while time scales of 10 years may be required to amalgamate the 3 effects of fast, episodic processes. Furthermore, if fast, episodic processes can be considered as continuous for the purposes of representation in models, then diffusion can be used to simulate the combined effects of both fast and slow mass movements. This representation of fast, episodic processes may be appropriate at large scales for which detailed individual landslide erosional and depositional features are unresolvable. Further examination of field data is required in order to determine which transport processes can be modelled effectively using a diffusion analogy. Some modified version of the diffusion equation (e.g., one which incorporates a threshold angle for failure, or a nonlinear form of the relation) may be more appropriate for the simulation of fast, episodic mass movements (e.g., Anderson, 1994; Howard, 1994; Martin and Church, 1997). Some researchers have expressed concern that diffusive processes alone would smooth landscapes (e.g., Chase, 1989). At the grid scales of most numerical landscape evolution models, local irregularities cannot be resolved. It is true that diffusion alone acts to smooth out whatever larger-scale irregularities exist in the initial topography. However, if non-diffusive processes such as the channel system are also modelled then roughening effects in the landscape can be introduced. In addition, the 42 inclusion of infrequent, but nevertheless very important, processes such as deep-seated landsliding may continue to re-introduce large-scale topographic irregularities as diffusion acts to smooth the landscape. Several of the models include weathering in their assemblage of processes (e.g., Ahnert, 1976; Anderson and Humphrey, 1989; Tucker and Slingerland, 1994). In such models, the thickness of the sediment layer is tabulated in order to determine the availability of sediment for transport. The thickness of the sediment layer influences overall hillslope gradients as bedrock slopes remain stable at steeper angles than sediment-covered slopes. In view of the rich range of possibilities described here, and the specific spatial and temporal constraints that accompany many of them, it is surprising that there has been very little testing and calibration of transport relations with field data. Most of the modellers have not undertaken independent studies to obtain their diffusivity values. It appears that the choices are either based on past studies (i.e., the scarp studies of Nash, 1980a,b; Colman and Watson, 1984; Hanks et al., 1984) or are obtained by running the model and seeing which values produce "realistic results". In the former case, a problem arises as the scarp studies are limited to profile sections and hence are not necessarily representative of the larger areas simulated in landscape evolution models. Moreover, the diffusivity values for the scarp studies may be unique to the particular region in which they were derived. In the latter case, there is no real knowledge regarding long-term rates of hillslope processes. Such studies must be undertaken if we are to obtain an understanding of what constitutes realistic model results (Oreskes et al., 1984). 43 3.3 DIFFUSION IN T H E HILLSLOPE SUBMODEL 3.3.1 INTRODUCTION Equation 3.3 is rewritten with two diffusion terms in order to incorporate both slow and rapid mass wasting processes (Martin and Church, 1997): dh dh 7+ T dx by 2 dt 1 y i 2 2 (3.7) 2 wherein a is the diffusivity for slow, quasi-continuous mass movements (e.g., creep processes) and p is the diffusivity for rapid, episodic mass movements (e.g., shallow landslides). This approach is analogous with the treatment of molecular and eddy viscosity in fluid mechanics. The incorporation of p entails the assumption that, over the long-term, rapid mass movement affects the landscape everywhere that is above some threshold gradient. This approach is appropriate for landscape regions in which creep and/or landsliding constitute the principal processes operating on hillslopes. Depending on the characteristics of a particular region either slow or fast mass movements may dominate. The present model focusses on landscape evolution in coastal drainage basins in British Columbia, which is a mountainous, humid region. Other processes, such as earth flows and slope wash, which may be important in specific environments (e.g., badlands), are assumed to be important only at a local scale in the present study. Although it is recognized that such processes may require consideration when modelling landscape evolution in certain environments, they are not included in this version of the present model. However, such processes could readily be incorporated into the model at a later date. Debris flow processes are often considered as hillslope processes and represent an important mechanism delivering hillslope material to the fluvial system. However, they are considered to be channel processes in this study and, as such, are considered in the following chapter. 44 The basic fixed coefficient form of the diffusion equation assumes a linear dependence of transport on slope tangent. These questions are considered by writing the diffusion equation in a form that take advantage of available transport data: dh dt z= d dh m — dx dx — d +— dh m — (3.8) The diffusion coefficient is now represented as a standard mass transfer rate, m (ML" T ) at unit gradient divided by sediment bulk density (pt). This reveals that the transport/gradient relation can be assessed using either volumetric or mass transport data. 3.3.2 DIFFUSIVITY FOR SLOW QUASI-CONTINUOUS MASS MOVEMENTS 3.3.2.1 EXAMINATION OF PAST CREEP DATA In order to assess the appropriateness of a diffusion analogy, and to calibrate such an equation, if it indeed appears to be reasonable, requires long-term transport estimates. But longterm estimates of total sediment transport on hillslopes are not generally available and what measurements there are, have not generally been subject to adequate control. Nonetheless, the calibration of the diffusion equation for slow and fast processes is based on what transport data are available. To appraise diffusion coefficients for slow, quasi-continuous mass movements, the volumetric creep data presented in Young (1974), Saunders and Young (1983) and several more recent studies were examined (table 3.1) (Martin and Church, 1997). These two compilations, from which most data are extracted for the present study, were not originally assessed for reliability of measurements. Therefore, their adoption into the present study involves only a generalized analysis of the data. Some results (Owens, 1969; Finlayson, 1981) exhibit a weak, and statistically non- 45 Table 3.1 Creep transport rates. Study* Average creep rate •a Leopold and Emmett (1972) Chandler and Pook (1971) Young (1960,1963) Young (1978) Anderson (1977) Finlayson(1981) Carson and Kirkby (1972) Lewis (1975) Williams (1973) Day (1977) Kirkby (1964,1967) Slaymaker(1972) Williams (1973) Owens (1969) Williams (1973) Everett (1963) Dedkov and Duglav (1967) Lewis (1974, 1976) EylesandHo(1970) Barr and Swanson (1970) McKeanetal. (1993) Location (cm /cm/yr) 0.2 Washington D . C , U.S.A. 0.3 0.5 0.6 0.8 1.2 1.3 1.5 1.9 2.0 2.1 2.7 3.2 3.3 5.9 6.0 7.1 Central England Northern England Derbyshire, U.K. Weardale, U.K. Mendips, U.K. Maryland, U.S.A. Puerto Rico N.S.W., Australia Wales Scotland Wales N.S.W., Australia New Zealand Northern Territory, Australia Ohio, U.S.A. Tatar, U.S.S.R. 8.0 12.4 15.0 67 Puerto Rico Malaya Southern Alaska California, U.S.A. 46 significant, relation between gradient and creep transport (figure 3.1). Others (Williams, 1973) exhibit no correlation at all. Most of the studies are short-term and it is possible that a relation is detectable only at longer time scales when short-term variability of controlling factors, such as soil moisture and soil cohesion, may average out. It may also be the case that such properties vary sufficiently over both space and time to disguise the effects of gradient at all scales. 3.3.2.2 DIFFUSIVITY V A L U E FOR CREEP The diffusion coefficient is numerically equivalent to the volumetric transport rate at unit gradient (45°). All of the observations were made on considerably lower gradients but, since functional dependence is observed to be weak and in the absence of more definitive evidence, the median value from several studies of 0.0002 m yf (figure 3.2) is assumed to represent the 2 l correct order of magnitude value for the creep diffusion coefficient (Martin and Church, 1997). The data show an extremal distribution with lower values of creep transport rates occurring with a much greater frequency than high values. In a stochastic approach, the variation of diffusivity would be represented by an extremally distributed variate. There remains a need for additional study in order to evaluate the factors responsible for variation in creep rates. 3.3.2.3 NONLINEAR RELATION? Results of a field study by Schumm (1964) and a controlled experiment by Van Asch et al. (1989) suggest that there may be a nonlinear relation between creep rate and gradient (figure 3.3). This finding is supported by Andrews and Bucknam (1987), who proposed a nonlinear relation between transport and gradient on scarps. Confirmation of nonlinearity would be an important result because it would necessitate the formulation of a more suitable creep transport equation. 47 Figure 3.1 Relation between gradient (as slope tangent) and creep transport rate. a) Data are from Finlayson (1981). The negative values represent upslope transport. b) Data are from Owens (1969). a) 15 ~ 10 E E 5 CU t o o 4 CL W 4 Q_ CD o -10 -15 0.1 0.2 0.3 0.4 0.5 0.6 Gradient b) E o m E o, d) 2 -d o CL in c CO L_ -*—> CL CD CD o 0.1 0.2 0.3 0.4 Gradient 0.5 0.6 0.7 0.8 48 Figure 3.2 Frequency distribution of volumetric creep rates. Note outlier to the right of scale break. This creep estimate was obtained by McKean et al. (1993) using an isotope mass balance model whereas other studies are direct measurements. CNi r- n f m <D i i CM CO i i W N CO O ) O ^ CD N- CO CD O ( N CO T— CM o o coo) o"o co oo T-T-T-T"1--<-CM-<fr o i- m i i CD N i I i I I 7 CO Volumetric creep rate n^fincDNcoroo (cmVcm yr) Figure 3.3 Relation between gradient and creep transport rate. Gradient is represented as slope tangent. Experimental data are from VanAsch et al. (1989). Data sets represent different slope conditions simulated in the laboratory. 140 ^ trees, groundwater depth of 0 m 120 f • bare, groundwater depth of 0 m A trees, groundwater depth of 1 m 100 O bare, groundwater depth of 1 m 1 • E B 80 tr o Q_ w 60 40 20 n 0.2 0.25 0.3 0.35 6 Q Q 0.4 Gradient 0.45 0.5 0.55 0.6 49 3.3.3 DIFFUSIVITY FOR RAPID EPISODIC MASS MOVEMENTS: A CASE STUDY OF T H E QUEEN C H A R L O T T E ISLANDS 3.3.3.1 INTRODUCTION In order to approach rapid, episodic sediment transfer at the landscape scale data sets are required that incorporate many events. This approach is illustrated by using a data set from the Queen Charlotte Islands (Rood, 1984; 1990), which are located about 160 km off the coast of British Columbia. Landsliding rates are high due to the high annual rainfall (ranging from about 1500 to 5000 mm yr" ), rapidly weathered volcanic and sedimentary rocks, and glacially 1 oversteepened slopes. The landslide inventory was completed by identifying landslides on aerial photographs. The data set is believed to cover approximately a 40 year period as landslides older than this were not clearly visible on the photographs. In relation to landscape evolution, this period is very short. Nevertheless, the data constitute a comprehensive, extended record of landsliding activity. Forested portions of 23 watersheds (clearcut logging occurred in many basins) were analysed in order to obtain natural landsliding rates. The mean area of forested terrain in each basin is 12 km with a basin average over this area of about 30 landslides (table 3.2). An underlying assumption when applying the diffusion analogy to simulate fast, episodic mass movements is that such events can be considered to be continuous over the long time scales of a landscape evolution study. Rood (1984) estimates the total number of hectares per km of 2 forested steepland (> 20°) affected by landsliding per year and obtains a value of 0.03 ha/km /yr. 2 At this rate, about 3 000 years are needed for a failure to occur everywhere in the forested, steepland portion of study basins in the Queen Charlotte Islands. Therefore, optimum time scales for the simulation of landsliding using diffusion are greater than order of magnitude 10 years. 50 8 CU CO ccj CD CU fS S s CU U on T3 0) 2 00 2 vo 58 o o a CD CD — CD a of CD CO CO C I- cu c Cd CD o a CD CCj S-l bO ,o T3 a. o OD Ml -a <D 1.2 « HW CD 0) o CO 25 CD CO T3 •s «4H »H CD I I £* cu CD « H CD a CD cn .22 cs -a <D cd O CN 9 J3 o cd CN CN UH LH CO 0J) etT CD O cd g CCDD IQ =3 > CD cd -a cd is g 51 3.3.3.2 THE APPROACH The areas of each watershed that fall under particular slope gradient classes are calculated in order to evaluate the relation between gradient and transport rate. The following calculation is performed for each gradient class for each basin (Martin and Church, 1997): (3.9) 1=1 Cl, wherein qj is volumetric transport rate (m m per 40 yr), v is landslide volume, / is transport 3 _1 distance (typically the full length of the slope), a is area, n is the number of landslides and j is gradient class. This result is divided by the 40 years that the landslide inventory covers. Gradient distributions for forested regions in each basin are obtained from TRIM maps (scale 1:20 000) of the Queen Charlotte Islands. An overlay composed of 5mm circles (equal to 100 m on the map) was placed over the map sheet for each basin and numbers of contour lines in each cell were counted. The proportion of the forested region in each basin falling under a particular gradient class was calculated. These data were converted into areal values by multiplying each proportional value by the total forested area in each basin (refer back to table 3.2) Rood (1984, 1990) compiled data for debris slides and debris flows. Debris flow activity is included in the channel submodel in this landscape evolution model and is therefore eliminated from the hillslope analysis. The debris slide data include shallow landsliding events that are initiated on either open slopes or in gullies. Rood (1984) calculates landslide volumes by the summation of contributing volumes from the initiation and transport zones. The landslides are placed into the appropriate gradient class based on the slope angle in the initiation zone. Transport distance is required in order to obtain dynamic transport rates. Rood (1984) provided transport distances for approximately half of the landslides. In order to 52 determine transport distances for the missing values, correlations between transport distance and other relevant variables that were consistently measured by Rood (1984, 1990) were examined. A nonlinear relation between landslide volume and transport distance provided the strongest relation (figure 3.4). After a logarithmic transformation, the best-fit relation has an R value of 2 0.50 and a standard error of estimate of 0.33. The relation between transport distance and landslide volume with the correction factor for transformation bias (Miller, 1984; Appendix I) included is: / = 0.40v 077 (3.10) where / is transport distance and v is landslide volume. Finally, sediment transport rates for each gradient class for each of the 23 drainage basins are calculated according to equation 3.9 (p. 51). The plots of transport rate versus gradient for each drainage basins are shown in figure 3.5 (see section 3.3.3.4 for discussion of the two groups identified in the data). A threshold of about 30° is found for landslide activity. At higher gradients two interpretations are possible. Some basins appear to exhibit a distinctly nonlinear relation between transport and slope angle, whilst others show a step increase to a finite transport for slopes in the range 30°-45°. On slope angles greater than about 40° results are variable; this is not surprising since some cohesion mechanisms must come into play. A step-change is consistent with the belief that slopes in different gradient classes are dominated by distinct processes. Theory suggests that mass transport rates may be related to the sine of the slope angle rather than the tangent. However, this adjustment does not remove the nonlinearities from the data. 3.3.3.3 DIFFUSIVITY FOR T H E LINEAR RELATION Several landscape evolution models include both creep and landsliding among the 53 Figure 3.4 Relation between volume and travel distance for landslides in the Queen Charlotte Islands, British Columbia. Distance was measured to the nearest 5 meters. Data are from Rood (pers. comm). 54 Figure 3.5 (a) Landslide transport rate vs. gradient for Group A drainage basins. Gradient is represented as slope tangent. BURNABY ISLAND 0.8 0.6 0.4 0.2 • 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 • 0.8 1 1.2 1.4 1.6 1.8 2 MATHESON HEAD GOVERNMENT 0.8 0.8 0.6 0.4 0.4 0.2 o I • 0 • — • — • 0.2 0.4 0.6 t t 0.8 1 1.2 1.4 1.6 1.8 • 2 • 0.2 0.4 0.6 0.8 1 1.2 MOSQUITO CREEK TRIBUTARY E 1.4 1.6 1.8 2 POWRIVCO 0.8 0.8 0.6 0.4 0.4 0.2 0.2 o Q. W C CO 0 0.2 • 0.4 •-- 0.6 0.8 -•—•1 1.2 1.4 1.6 1.8 • 2 • 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 SECURITY SACHS 0.8 0.8 0.6 0.6 0.4 0.4 0.2 • 0 4 - * - 4 0 • 0.2 0.4 • 1 1 1 1 1 1 0.6 0.8 1 1.2 1.4 1.6 1 — 1.8 2 0.2 0.4 0.6 0.8 1 1.2 TALUNKWAN ISLAND 1.4 1.6 1.8 2 TARUNDL 0.8 0.8 0.6 0.6 0.4 0.2 0.2 ol • 0 •—•—•—•—• 0.2 0.4 0.6 0.8 1 oJ-*-# 1.2 Gradient 1.4 1.6 1.8 2 0 0.2 • • • • 0.4 0.6 0.8 1 1.2 Gradient 1.4 1.6 1.8 2 55 Figure 3.5 (b) Landslide transport rate vs. gradient for Group B drainage basins (Bonanza - Marshall Head). Gradient is represented as slope tangent. BONANZA 0.8 0.6 0.6 0.4 0.4 0.2 0.2 o J - 0 * 0.2 - * — 0.4 • — • 0.6 — * 0.8 - 1 1.2 1.4 1.6 1.8 GREGORY 0.8 2 0 I • • 0 0.2 • 0.4 • 0.6 1 1 —I 1 1 0.8 1 1.2 1.4 1.6 1.8 2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 1 — 1 HANGOVER 0.8 0.6 0.4 E 0.2 • 0 CD 2 •c o C L CO £= CD • 0.2 • 0.4 0.6 0.8 1 1.2 1.4 1.6 1 1.8 2 0 0.2 1 JASON 0.8 LANDRICK 0.8 0.6 • • 0.6 0.4 0.4 0.2 0 0.2 0.4 • t 0.6 0.8 H 1 1.2 1.4 1 H 1.6 1.8 2 1 - 0 •—•—• 0.4 0.6 0.2 + 0.8 1 1 1 1 1- 1.2 1.4 1.6 1.8 2 1 MACMILLAN MARSHALL HEAD 0.8 • 0.8 0.6 0.6 0.4 0.4 0.2 • • 0.2 0.4 0.6 0.8 1 1.2 Gradient 1.4 1.6 1.8 2 0 I • • 0 0.2 • • • 1 1— 0.4 0.6 0.8 1 1.2 Gradient 1.4 1.6 1.8 2 56 Figure 3.5 (b) Group B drainage basins cont'd (Mountain - Windy Bay). Gradient is represented as slope tangent. 1 1 MOUNTAIN E 0 RILEY O.S 0.8 0.6 + 0.6 0.4 0.4 0.2 0.2 0 • • 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 o 4 - 4 - » — • — • — 0 0.2 0.4 0.6 0.8 2 1 1.2 1.4 1.6 1.8 2 -— t • 1 WINDY 0.8 BAY 0.6 0.4 0.2 o 4 - * * — • — • — # — ? — i 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 • SOUTH BAY 2.8 | 0.4 0.6 DUMP 2.6 2.4 2.2 2 1.8 1.6 j I 1.4 1.2 1 0.8 0.6 0.4 0.2 0 • •—•—• 0.2 0.4 0.6 1 0.8 1 1 1 1.2 Gradient 1 1.4 1 h~ 1.6 1.8 2 0.8 1 . 1.2 1.4 1 1— 1.6 1.8 2 57 hillslope processes simulated using diffusion (Koons, 1989; Kooi and Beaumont, 1994). A threshold gradient is not incorporated in the diffusion equation in these models. The Queen Charlotte Islands' data suggest that threshold gradients may be critical for the appropriate parameterization of landsliding. However, the analysis of the transport/gradient relation is initially undertaken without consideration for threshold gradients in order to make comparisons between results found in the present study and diffusivity values used in other landscape evolution models. The following approach is adopted in order to obtain a representative linear diffusion coefficient. Best-fit linear relations were determined in order to obtain the transport rate at unit gradient (45°) for each basin (table 3.3). When the results are tabulated, an extremal distribution is once again evident, as was the case for the compilation of creep transport rates (figure 3.6). Therefore, the median value of 0.1 m /yr is adopted as the best representation of the 2 value for linear diffusivity. 3.3.3.4 NONLINEAR TRANSPORT RELATION A nonlinear relation, strictly observed, requires a modification of the governing equation. Nonlinearity may appear in the transport/gradient relation for several reasons (Martin and Church, 1997). One is that transport distance itself depends on slope angle. Kirkby (1992) has elaborated such a model from his general formulation of hillslope evolution. However, a related possibility is that starting angle and slope length (which enter the transport estimates) are themselves confounded in the data. In fact, they are independent of each other. It finally must be acknowledged that the gradient/transport relations are based on slope angle in the detachment zone. But slopes in the Queen Charlotte Islands are typically relatively short (about several hundred meters to one kilometer) and nearly rectilinear, so this factor appears not to be a source of serious bias. It remains to explore more complex models for these data. Table 3.3 Diffusivities for linear diffusion. Drainage Basin Diffusivity (m /yr) 0.072 0.27 0.14 0.029 0.41 0.18 0.28 0.24 0.090 0.25 0.22 0.037 0.0062 0.16 0.022 0.11 0.14 0.044 1.9 0.0041 0.0043 0.13 0.28 2 Armentieres Bonanza Burnaby Island Government Gregory Hangover Inskip Jason Landrick Macmillan Marshall Head Matheson Head Mosquito Mountain Powrivco Riley Sachs Security South Bay Dump Talunkwan Island Tarundl Two Torrent Windy Bay Figure 3 . 6 Frequency distribution of landslide transport rates at unit gradient. The outli< at the right of the graph is the datum for South Bay Dump Creek. in o ro £ CD CO O CD E —s^r— LO q in co in CD T—' in" cn in GO Landslide transport rate at unit gradient (m /m yr) 3 Figure 3 . 7 Sample plot of nonlinear relation used in analysis. Gradient 60 Nonlinear best-fit models were initially tested for the transport/gradient relations for each watershed. However, these models did not capture the key characteristics of the relations. For this reason, it was decided instead to adopt a suitable nonlinear relation that embodies the critical patterns observed in the data. A hyperbolic tangent relation was found to replicate reasonably the transport/gradient relations for the Queen Charlotte landslide data. The present model incorporates the notion of a critical angle for landsliding. This relation shows zero or low transport rates at gentle gradients and then rises sharply (figure 3.7), thereafter stabilizing at a higher value. It was decided to include this latter feature of the curve as in some cases the transport rates keep rising (e.g., Inskip Creek; see figure 3.5) and in other cases there is a stepchange to a lower value (e.g., Hangover Creek; see figure 3.5). The option of maintaining a constant value at high slopes seems to express some sort of average for the types of behaviour that might occur. In reality, landsliding does not generally occur on slopes greater than about 5060° (Rood, 1995, pers. comm.), where weathering rates cannot keep up with transport rates. In this case, the hillslopes are depleted of sediment cover and bedrock slopes appear. An option to reduce landsliding transport rates to zero at slopes greater than 60° (approximately tangent slope 1.7) can also be included. The hyperbolic tangent relation used in this analysis is of the form: q = K + [0.5 * (K, - K ) * (l + tanh(value))] 0 0 (3.11) where: value = abs V<3xy ^x scale and q is transport rate, tanh is the hyperbolic tangent, Ko and K; are diffusivity parameters (m /yr), h o represents the location on the x-axis of the mid-point of the rise for the tanh function 2 x 61 (m/m) and h x ie sca is the width parameter for the rise (which determines the steepness of the rise) (m/m). Examination of the data reveals that there are two major patterns of transport/gradient relations among the 23 drainage basins (table 3.4; figure 3.8). Each group shows similar threshold gradients and transport rates. The drainage basins in Group A exhibit characteristically low sediment transport rates and threshold gradients ranging from tangent slope 0.6 to 1.0. Sediment transport rates remain relatively low at unit gradient and peak at tangent slopes of about 1.2 and 1.4. In other cases, transport rates show a gradual increase in transport rates over tangent slopes of 1.0 and 1.2. The threshold transport tangent slope for Group B watersheds is about 0.8. Transport rates increase abruptly at this threshold. The median transport values for each gradient class are determined for the two groups. These values are used to define hyperbolic tangent relations based on equation 3.11 (figure 3.9). Geological formations found in the drainage basins were examined in order to determine if geological differences might be responsible for the two distinctive patterns of transport/gradient relations (refer back to table 3.4). The less active drainage basins of Group A show a larger proportion of basins composed of hard volcanic and granites. Group B, on the other hand, has a larger proportion of basins composed of soft volcanic and sedimentary rocks. Therefore, there appears to be a reasonable basis for the suggestion that geology is a major determinant of the landsliding transport rates in this study. 3.3.4 COMPARISON OF DIFFUSIVITIES WITH OTHER STUDIES Previous attempts to estimate linear sediment diffusivities have been based on scarp erosion studies (table 3.5a). These are local studies in which a diffusion coefficient was fitted to slope profile development. The primary process in operation on some slopes is slow mass 62 Table 3.4 Transport/gradient categories and dominant rock types. Group A Armentieres (Hard Volcanics) Burnaby Island (Granites) Government (Hard Volcanics) Matheson Head (Hard Volcanics) Mosquito Creek (Hard Volcanics) Tributary Powrivco (Soft Volcanics) Sachs (Clastic Sedimentary) Security (Hard Volcanics) Talunkwan Island (Clastic Sedimentary) Tarundl (Clastic Sedimentary) Bonanza Gregory Hangover Inskip Jason Landrick MacMillan Marshall Head Mountain Riley South Bay Dump Two Torrent Windy Bay Group B (Soft Volcanics) (Soft Volcanics) (Soft Volcanics) (Hard Volcanics) (Hard Volcanics) (Soft Volcanics) (Clastic Sedimentary) (Hard Volcanics) (Granites) (Soft Volcanics) (Clastic Sedimentary) (Clastic Sedimentary) (Soft Volcanics) Table 3.5 (a) Comparison with diffusion coefficients derived from scarp studies. Study Present study Diffusion Coefficient (m /yr) 2 — i_j Landslides: 2x10" Creep: 2x 10" 1.2 x 10" 4.4 x 10" 9x 10" 1.1 x 10" 1.1 x 10" 1.6 x 10" 4 Nash (1980a) Nash (1980b) Colman and Watson (1984) Hanks etal. (1984) 2 4 4 3 2 2 (b) Comparison with diffusion coefficients implemented in landscape evolution models. Study Present Study Diffusion Coefficient (m /yr) 2 i Landslides: 2x10" Creep: 2x 10" 0.001, 10 1 x 10 to5x 10 1.5 x 10" to 1.5 x 10 Order of 10" 4 Anderson and Humphrey (1989) Flemings and Jordan (1989) Koons(1989) Anderson (1994), Rosenbloom and Anderson (1994) Howard (1994) Kooi and Beaumont (1994,1996) 1 2 3 1 1 2 0.004-7.0 2x 10" to 1 x 10 2 2 63 Figure 3.8 (a) Relation between gradient and transport rate for Group A drainage basins. Gradient is represented as slope tangent. Note scale change on y-axis in part (b). 0.16i 0.14~ >. E 0.120.10+ = o 0.08 + CO £ 0.06 • Mosquito Creek Tributary • Government A Matheson Head s Burnaby Island • Tarundl • Security Talunkwan Island • Sachs A Powrivco 0 0.04 + 0.02 + A a 0 —• 0.2 H 0.4 H 0.6 -*1 B— 0.8 —i 1.2 1 1 1.4 1.6 1— 1.8 1.8 2 Gradient (b) Relation between gradient and transport rate for Group B basins. Gradient is represented as slope tangent. 1.0" E 0 •4—« 2 t o CL CO c CD 0.8 1 Gradient 1.8 Figure 3 . 9 a) Transport relation for Group A drainage basins, b) Transport relation for Group B drainage basins. Gradient is represented as slope tangent. a) 0.6 Variable Group A 0.5 K 0 K, 0.027 h 0.7 * *x sca/e 0.1 0 J= 0.4 M 0.3 .2 0.2 0.1 + —I— 0.5 1 1.5 Gradient 0.6 Variable Group B 0.5 K 0 K, 0.43 0 0.4 + • f t 1 1 1.06 ~ ^ 0.3 x 1 1 0.19 sca/e 0.2 / 0.1 -+- — / / / • # / / / / / / / / / / / / /• • • " I — » " 1 0.5 1 1.5 Gradient 65 movement such as creep and/or transport due to basal undercutting (Colman and Watson, 1984; Nash, 1980ab). Hanks et al. (1984) stated that diffusion: "...accounts, correctly or incorrectly, for all the physical processes that contribute to erosional mass transport on the slopes of interest in this study." The diffusivity for creep derived in this study is at the lower end of the range of scarp values and the landslide diffusivity derived here is at least an order of magnitude greater than these values. From the perspective of landscape modelling, data derived from an ensemble of many possibly representative hillslopes, as in the present study, are apt to be more informative than results from individual scarp slopes. However, the scarp slopes are relatively steep. This may possibly be indicative of diffusivities for a combination of processes operating near unit gradient, which is what the diffusion coefficient actually represents. The linear diffusivity estimates (non-threshold) for the present study are compared with diffusivities adopted in landscape development models (table 3.5b). The diffusivities employed by Anderson (1994), Rosenbloom and Anderson (1994) and Howard (1994) represent values for slow mass movements but are greater than the creep diffusivity found in the present study. Anderson (1994) and Howard (1994) included an additional term to simulate mass movement activity. The diffusivities implemented by Koons (1989) and Kooi and Beaumont (1994, 1996) were assumed to represent all slope processes including fast mass movements. The values at the higher end of the ranges are used when modelling humid regions. Kooi and Beaumont (1994) increased diffusivities when modelling landscapes covered in loose sediment and not bedrock. The values used in these models are higher than the landslide diffusivity estimated in the present study (and higher even than an upper extreme). The Queen Charlotte Islands are known to have high landsliding rates; it is unlikely that diffusivities for landsliding would exceed those found in 66 the present study by the three orders of magnitude suggested by the upper estimate of Kooi and Beaumont (1994). The lower diffusivity value of Anderson and Humphrey (1989) was used for the modelling of creep processes on scarps and moraines. The higher value of diffusivity used by Anderson and Humphrey (1989) and the range of values used by Flemings and J o r d a n (1989) exceed the diffusivities found in the present study. In the former case, the diffusion coefficient is derived from an estimate of debris flow delivery rates to fans. The diffusivity values of Flemings and J o r d a n (1989) are based on estimates of mean regional gradients in mountain belts and basin fill rates; both hillslope and fluvial processes may be reflected in these high values. 3.3.5 SUMMARY On balance, relations between transport rate and gradient appear to be nonlinear. If this is the case, then the hillslope process equations used in landscape evolution models must be selected to incorporate these nonlinearities. This has not, heretofore, been the rule. However, significant variations in the transport/gradient relations among the Queen Charlotte Islands' study basins were found. It is suggested that geology is a significant factor influencing landsliding activity. The linear diffusivity for landsliding found in this study is about 3 orders of magnitude greater than the creep diffusivity. This result suggests that in susceptible terrain only the latter need be considered in long-term landscape evolution models (Martin and Church, 1997). However, if a threshold gradient landslide model is used then creep is the dominant mass wasting process operating on lower gradient hillslopes. The low value of creep diffusivity suggests that the morphological evolution of hillslopes by creep occurs at very slow rates. In such cases, channel processes may have a particularly critical role in landscape formation. In steep basins 67 showing rectilinear slopes such as those on the Queen Charlotte Islands, significant net hillslope degradation by landsliding is likely to occur only on the upper slopes, whereas significant deposition of hillslope material may be focussed at the slope base. When gradients, geology and vegetation are approximately constant along the slope length, erosion and deposition may be approximately balanced in the mid-sections of hillslopes. To confirm current appearances, field data must cover longer time periods and must be controlled to identify sources of variability. The creep data are confounded by inadequately controlled measurement programs. The time period of such measurements, often only several years, is usually too short for the purposes of long-term modelling. Landsliding activity, on the other hand, is visible on aerial photographs, making it possible to obtain medium-term inventories of landsliding rates. Landslide inventories collected for the purpose of estimating transport rates must include key variables such as initiation angle and transport distance, in addition to landslide volumes. There are severe limitations to the time scales of transport rates that are available by direct measurement. Techniques based on absolute dating methods should be subject to much further investigation. The study of McKean et al. (1993) uses cosmogenic isotopes for the estimation of creep rates. The values found in this study are significantly higher than values obtained by traditional field methods. Further studies are required in order to assess more thoroughly the ability of isotope dating methods to estimate creep. The ability of such techniques to estimate transport rates of more rapid, episodic transport processes should also be investigated. Although laboratory experiments do not provide transport estimates that are representative of those found in nature, further experimentation (particularly relevant for creep) might help us to elucidate the nature of the transport/gradient relation. Recently, Densmore et al. (1997) used a simple physical model (consisting of beans!) to simulate patterns of bedrock 68 landsliding across hillslopes. Such experimentation can provide critical insights into landscape evolution that can be incorporated into existing models. 3.4 W E A T H E R I N G P R O C E S S E S Transport-limited and weathering-limited behaviours are widely acknowledged concepts in the geomorphological literature. In order to determine conditions necessary for the occurrence of one of these particular types of transport behaviour, the relative efficacy of erosion/deposition rates and soil production rates must be known. Despite the relevance of soil production to geomorphological processes, there have been very few field studies which document its rates. It is generally assumed that soil production rates due to bedrock weathering reach a maximum at either the surface or some intermediate depth of soil thickness (Gilbert, 1877; Dietrich et al., 1995). Production rates thereafter decline (an exponential decrease is usually assumed) until such a depth below which the value is zero. In the several more recent studies which consider this phenomenon, an exponentially declining model has been utilized (Anderson and Humphrey, 1989; Tucker and Slingerland, 1994; Dietrich et al, 1995; Heimsath et al., 1997) (figure 3.10): = 1 x 10" exp [-10 (h - h )] Anderson and Humphrey (1989) (3.12a) = 5 x 10" exp [-10 (h - h )] Tucker and Slingerland (1994) (3.12b) = 1.9 x 10" exp [- 5 (h - h )] Dietrich et al. (1995) (3.12c) 4 dt dh, dt dh, dt dh b dt s b 5 s b 4 s b = 7.7 x 10" exp [- 2.3 (h - h )] Heimsath et al. (1991) 5 s b (3.12d) where h is surface height, hb is the height of the soil/bedrock interface, and the change in height s is given in m/yr. Anderson and Humphrey (1989) did not discuss the source of the constants used in their relations. Tucker and Slingerland (1994) followed the approach taken by Anderson 69 Figure 3.10 Weathering rate vs. depth 0.0002 Soil depth (m) 70 and Humphrey (1989) when constructing their weathering equation. A value of 10, which was used by Anderson and Humphrey (1989), was inserted into their equation for the constant in front of the term for sediment layer depth. These two relations plot below the relation of Dietrich et al. (1995), and show only negligible soil production rates when soil thicknesses exceed about 0.5 m. Dietrich et al. (1995) derived their relation from soil production estimates in the field at soil thicknesses of 1.5 m and 0.30 m. The production rate is 0 m/yr in the former case. In the latter case, the value of 4.2x10" m/yr is based on an assumption that soil thickness at the chosen site 5 remained constant over the study period of sediment accumulation in the associated valley. Accumulation rates are assumed to also represent erosion rates, and in turn soil production rates. Heimsath et al. (1997) used measurements of in situ produced Be and A l concentrations in 10 26 bedrock under soils of different thicknesses. The study was conducted in Tennessee Valley, Marin County, California. The average rainfall in the area is 760 mm and the vegetation is grassland and scrub. Results of the study suggest that a model in which soil production rates decline with increasing depth provides a good fit to the data. The basic exponentially declining soil production/depth relation of Heimsath et al. (1997), which is based on cosmogenic data, is adopted for use in the present landscape evolution model. The study is field-based and makes use of a long-term dating technique. Elevations of the surface and bedrock/soil boundary are determined at each time step. When the soil layer is depleted, diffusivity falls to 0 m /yr. Therefore, an additional dependency enters the transport equation whereby diffusivity is a step-function of soil layer depth. Soil production rate is dependent on soil layer depth, which is simply the difference between the elevations of the surface and soil/bedrock interface. The equations should be solved simultaneously in order to obtain the best results when incorporating weathering in the model runs. 71 C H A P T E R 4: C H A N N E L S U B M O D E L 4.1 F L U V I A L T R A N S P O R T R E L A T I O N S IN L A N D S C A P E M O D E L S 4.1.1 I N T R O D U C T I O N The movement of sediment through the channel system is an integral component of landscape evolution. The channel system determines the local lowest elevations in the landscape, which provide the base level and the potential energy gradient, as defined by local relief, for hillslopes. The approach used to describe channel processes varies considerably among landscape evolution models. Several of the early modellers (e.g., Penck, King, Budel) emphasized hillslope processes and excluded detailed consideration of channel processes. However, unlike the other classic modellers, Davis (1899) considered carefully both hillslope and channel processes in his cycle of erosion. Several aspects of fluvial processes, which are important to landscape evolution, were considered by Davis: (i) sediment delivery from hillslopes to channels; (ii) quantity and coarseness of load; (iii) fluvial aggradation; (iv) the concept of grade; and (v) river meandering processes. Numerical landscape evolution models which incorporate fluvial activity can be grouped into three major categories based on overall model structure and fluvial process specification: (i) Certain general models which can be adjusted to simulate either hillslope or stream profile development by varying the values of some key coefficients and exponents (Culling, 1960; Kirkby, 1971). Such an approach does not allow for the integration of both hillslope and fluvial processes into one comprehensive model. (ii) Some later models explore the initiation of channel networks, with a particular emphasis on processes leading to channelization (e.g., Ahnert 1976; Willgoose et al., 1991a,b; Rinaldo et al., 72 1994). (iii) Many of the more recent models include fluvial processes as a significant component of an integrated model of landscape evolution (e.g., Kooi and Beaumont, 1994; Tucker and Slingerland, 1994). Transport in smaller debris flow channels is not considered in most cases. 4.1.2 DRAINAGE INITIATION Several existing landscape evolution models focus on the initiation of channel networks and not the later stages of landscape evolution. In these models, transporting processes driven by non-channelized water flow, such as slope wash, operate over the entire model domain. Ahnert (1976) assessed drainage initiation by resolving the net effects of non-channelized transport processes over the landscape. Willgoose et al. (1991a,b) included a channel initiation function, which is nonlinearly dependent on slope and discharge, in order to determine if channelization occurs. The model of Chase (1992) invokes rules for the transport of sediment by water flow which allow the model to eventually form its own channels. The model of Rigon et al. (1994) involves an exploration of the self-organized nature of channel networks by considering the processes leading to channel initiation. Fluvial erosion occurs when shear stress exceeds some threshold value. Smith and Bremerton (1972) recognized the difference between drainage initiation and subsequent landscape evolution and created a distinct model for each phenomenon. Several landscape models also do not explicitly focus on the issue of drainage initiation (Tucker and Slingerland, 1994; Kooi and Beaumont, 1994; Koons 1989; Rosenbloom and Anderson, 1994; Anderson, 1994). The location of the drainage network is defined at the beginning of a model run. The different objectives of these two classes of models and the implications this has for 73 model creation have not otherwise been carefully considered in the literature. Because of the different aims and scales of these types of studies, process parameterization in each case should distinctly reflect the requirements of the particular approach. The primary objective for drainage initiation studies is to determine the factors influencing the creation of drainage networks and to evaluate the resultant patterns. Because of the importance of processes such as rill formation at the beginning phases of channel initiation, model scales have to be reasonably small in the early stages of network development. Once the drainage network is established, the scale of study should become significantly larger as sediment transfers over significant distances are considered. 4.1.3 LINEAR CHANNEL ELEMENTS IN A SPATIAL MODEL In many models, hillslopes and streams are not treated as unique physical entities, despite the specification of unique transport equations for each of these classes of process (Chase, 1992; Howard, 1994; Kooi and Beaumont, 1994; Tucker and Slingerland, 1994). Hillslope and channel equations are applied across the entire model domain. The rationale behind this approach has been based on the contrasting spatial scales of hillslopes and fluvial processes. Typically grid elements in a model are significantly larger than lateral channel dimensions (channel width). Therefore, when grid cells are large (in many cases >1 km) an assumption of models has been that each cell contains channel elements in addition to the hillslope elements. The morphological changes occurring in the channel are effectively "smeared" across the entire grid cell. Chase (1992, p. 41) described his rationale for following this approach: "...in a model intended to work at grid resolutions up to km, there is no point in attempting to distinguish between hillslope and channel elements of the topography. Each grid cell represents topography containing both. Therefore, no a priori distinctions are made between slope and channel, and the 74 rules are uniform for all grid cells." Kooi and Beaumont (1994, 1996) also justified their decision to use this approach on the grounds that the large spatial scale and low grid resolution of the study do not allow distinctions between hillslope and channels to be made. The adoption of expedient approaches, such as the smearing of morphological changes and the operation of all processes across every cell, are necessary because of the difficulties encountered when incorporating the linear channel system into a numerical model structure based on grid cells. These assumptions allow for simplicity in the model structure. The question remains as to whether such approaches reasonably simulate the significant features of evolving landscapes. The ability of such approaches to provide reasonable results depends on the intended resolution of a study. However, because of the complexity of many landscape models, computational expediency may often require a certain degree of generalization. Other models distinctly define the location of channel cells in the model surface. Smith and Bretherton (1972), in the second model presented in their paper, applied one set of equations to hillslopes and treat the channel as a line of zero thickness that exists at the intersection of the slopes. Koons (1989) defined the position of the main channel although he did not explicitly calculate fluvial transport rates in his model. Instead, stream elevations are obtained by fitting a logarithmic function between sea level and the dynamic elevation at the main divide. Willgoose et al. (1991a,b) applied overland flow and hillslope processes across the entire grid. However, fully-developed channelled flow occurs only in those grid cells defined as channel grid elements according to a channel initiation function. Rosenbloom and Anderson (1994) and Anderson (1994) set the initial location of the stream network, which consists of only lightly etched channels, at the start of the model run. The channel network is assumed to remain in the same planform position throughout the course of a model run. 75 4.1.4 S T R E A M P O W E R T R A N S P O R T R E L A T I O N S When simulating fluvial transport in landscape evolution models, the chosen relation should include variables that are resolved in the model. From this viewpoint, an attractive formulation is a stream power relation as the principal variates of width, depth and slope are resolvable at landscape scales. Derivations of this basic relation have been used in many landscape evolution models. Stream power per unit length of channel is defined (Bagnold, 1977, 1980, 1986) as: n = p„QS (4.1) where Q is stream power per unit length, pw is density of water, Q is discharge and S is slope. The right-hand side of the equation should also include the acceleration of gravity, but Bagnold eliminated gravity from the definition of stream power in his analysis. Stream power per unit width, co, is obtained by dividing n by channel width. Bagnold's stream power definitions specify stream power per unit length (and per unit width if co is used). Because Bagnold's research represents perhaps the most influential statement concerning stream power/sediment transport relations his definitions of these terms are followed. The most basic form of the sediment transport/stream power specification is related to equation 4.1 such that (see equation 3.4, p. 37): q oc Q S" m (4.2) where q is sediment transport rate, and m and n are exponents that vary depending on the particular formulation of the equation. Linear and nonlinear forms of the transport/stream power relation have been implemented in landscape evolution models by manipulating the exponents m and n. Area is often used as a surrogate for discharge because of the well-demonstrated relation between contributing area and discharge: Q = aA b (4.3) 76 where A is contributing area and the exponent b is typically less than 1 due to storage of water throughout the basin. Kirkby (1971, 1976) adjusted equation 4.2 to simulate sediment transport rates for slope wash and fluvial activity. The values of the exponents m and n vary depending on the particular process being modelled. For sediment transport in fully developed channels m is assigned a value between 2 and 3 and n is assigned a value of 3 (these values were derived from Leopold and Maddock, 1953). The variable discharge is replaced with area draining to a point in order to define the models in terms of morphological variables only. Smith and Bretherton (1972), Kooi and Beaumont (1994) and Tucker and Slingerland (1994) used the basic stream power relation of equation 4.2 (i.e., m and n are both set equal to 1). Howard (1994) introduced a more complex variation of the basic stream power equation. Willgoose et al. (1991) inserted an additional term into the right-hand side of the basic equation that adjusts the transport rate to reflect movement by either overland or channelized flow. Sediment transport equations for channelled flow usually incorporate some threshold condition before transport occurs (Gomez and Church, 1989). The incorporation of such a threshold has not been generally adopted in most landscape models. However, in the drainage initiation study of Rigon et al. (1994) the exceedance of a critical shear stress serves this purpose. The channel initiation function in the model of Willgoose (1991) also serves this function. The use of a threshold for transport is well established in gravels but has not been as closely adhered to for sand rivers (e.g., Einstein, 1950). Practically, a threshold condition is one way to "adjust" a model to account for the effects on sediment transport rates of surface properties such as grain size and armouring. The origin of the final form of the transport equation used in most landscape evolution models, including the exponents and coefficients, has not been adequately discussed and no 77 mention is made of a calibration procedure using field data. Furthermore, most of the models do not explicitly address the distinct nature of channel processes in low-order mountain streams and debris flow channels, which differ significantly from processes occurring in low-gradient streams. These channels are often important mechanisms for the delivery of hillslope material to higher-order channels. Stream power relations have been shown in the literature to provide reasonable estimates of transport rates for streams of relatively gentle gradient. However, the ability of such equations to describe processes occurring in low-order channels has not been evaluated. 4.1.5 NET CHANGE IN STORAGE Once sediment transport rates are calculated, they are often inserted into continuity equations within the model (Willgoose et al., 1991a,b; Kooi and Beaumont, 1994, Tucker and Slingerland, 1994). The most basic form of a continuity equation defined for a channel reach is: A Storage = q -q in out (4.4) However, in more recent models, the channel exists within the hillslope grid cells. The translation of calculated storage change into elevation changes is an important step in the numerical coding that is considered explicitly only by Tucker and Slingerland (1994). Channel elevations cannot be calculated independently of hillslope elevations, as the channel provides the base level for hillslope processes. Independent assignment of channel and hillslope elevation could lead to the absurd situation of a river being perched above its surrounding valley flat and/or hillslopes if channel aggradation exceeds the rate of aggradation of the cell in which it resides. The approach that has been adopted by Tucker and Slingerland (1994) is to obtain the change in elevation by dividing the change in storage by the grid cell dimensions such that: 78 d h dt = ^ Z ^ L ( 4 > 5 ) Ax Ay Changes i n channel storage are effectively spread out across the entire cell. When grid cells are significantly greater i n size than river and/or valley flat dimensions, this approach may significantly reduce the apparent elevation changes. A s such, it compromises an explicit representation o f channel morphology. The implications of this approach on landscape evolution should be assessed. 4.1.6 SUMMARY CRITIQUE Through a careful review o f modelling strategies several key issues emerge regarding the simulation o f fluvial processes in landscape models: (i) It is imperative that both hillslope and fluvial processes be integrated into comprehensive models o f landscape evolution which intend to include meaningful geomorphology. Interaction between these classes o f processes affects their respective rates of operation. (ii) A clear distinction must be made between models investigating drainage initiation and those emphasizing later stages of landscape evolution. In the former case, a criterion for channel initiation must be incorporated into the model, whereas the latter requires careful consideration of the long-term effects o f prolonged aggradation or degradation episodes on landscape evolution. (iii) Careful consideration should be given to the selection of resolved variates i n the sediment' transport equations. These variables drive the operation of processes and determine transport rates. (iv) Most o f the equations used in previous models do not recognize the use o f a threshold stream power for transport, although this is the standard for gravel transport equations. 79 (v) In most cases, there has been no solid justification given for the particular choice of the final form of the fluvial transport equation used in the model. There is minimal discussion about the criteria used in the selection of the coefficients and exponents in the equations. More generally, the equations used to model fluvial transport must be subjected to much further testing in order to assess their effectiveness in predicting transport rates. Critical evaluations of performance in comparison to real data are necessary for the evaluation of the proposed relations. (vi) Most of the models do not include sediment transport processes occurring in high-gradient streams or debris flow channels. The high activity rates occurring in debris flow channels provide an important mechanism for the delivery of sediment from hillslopes to the main stream channel in many steepland environments. (vii) Further research is required in order to determine the most effective methods for incorporating both fluvial and hillslope processes into the grid frameworks used in standard numerical models. The primary difficulty arises because of the fact that channel widths are considerably smaller than the usual sizes of grid cells which are implemented in numerical models. 4.2 T H E CHANNEL NETWORK SUBMODEL The actual channel network occupies only a small portion of basin area in comparison to hillslopes. Nonetheless, valleys are the location which experience the greatest elevation change during the early course of landscape evolution. This is because of high activity rates in channels. Furthermore, the delivery of sediment to the valley flat represents a critical stage in the transport of material through the basin. From this point, sediment has the potential to be transported very efficiently to the basin outlet. An objective of the present study is to develop a channel submodel that accounts for the 80 significance of channel processes in the geomorphological development of landscapes. This channel submodel will be constructed so that typical patterns of sediment movement through the channel system based on current configuration can be studied independently of evolving hillslopes (although both hillslope and channel processes can be integrated in the complete model). The channel submodel will incorporate: (i) fluvial transport occurring in streams of relatively low gradient; (ii) fluvial transport occurring in high-gradient mountain streams (transport behaviour differs from streams in the first category due to the particular characteristics of these channels); and (iii) debris flow transport. Each of these transporting regimes occurs in small, steep drainage basins in coastal British Columbia, the "prototype" region for submodel development and subsequent model runs. In addition to simulating patterns of sediment redistribution along the channel network, an objective of the present study is to simulate the aggradation and degradation in valleys which can, over the long-term, lead to the formation of significant valley fill deposits, canyons and terraces. Such features are significant as they affect hillslope base levels and regional gradients. In addition, these features document recent sedimentary history. However, prior models have not explicitly considered the aggradation and degradation of valleys. Both bed load and suspended load are considered in the channel submodel. A generalized bed load transport formula is adopted which requires the incorporation of several key parameters such as channel dimensions and gradient. A generalized approach based on sediment yield patterns is taken to estimate suspended load evacuation rates. Before further discussing the equations used in the channel submodel, a distinction must be made between sediment quantities that are measured in the field which are used in the development of a transport equation for the model, and those that are important in terms of 81 morphological change. Bed load and suspended load can be measured in the field. The erosion and/or deposition of sediments which comprise the bed and lower banks of a river are responsible for significant channel changes. Such material, herein referred to as bed material, can be moved as bed load and/or suspended load. Transport rates obtained by direct measurement at the bed involve only the fraction of bed material moving close to the channel floor. Bed load transport formulae are calibrated using such measurements and therefore exclude consideration of the finer portion of bed material. The bed load data sets compiled by Gomez and Church (1988), which are used to develop the transport equation in this study, have D o 5 values ranging from 0.9 mm to 32 mm. Therefore, the equation can be considered applicable for both gravel-bed (> 2 mm) and coarse-sand systems, although it does not cover the full range of bed material grain sizes. The erosion and deposition of material moved as bed load contributes most significantly to morphological changes in gravel-bed rivers. The finer portion of bed material load consists of material moved as bed load (coarse sands) and suspended load. Such material travels greater distances than coarse bed load and is often deposited as interstitial fill in gravel-bed rivers, and therefore does not significantly affect channel morphology in these circumstances. Although data often indicate that suspended load accounts for a considerably greater proportion of the total load, it is not necessarily always the most important from a channelforming perspective. Therefore, transport of both bed load and suspended load must be considered when modelling sediment transport in gravel-bed rivers. In sand-bed rivers, bed material can consist of the finer portion of bed load (coarse sands) and suspended load. Both of these categories of material are also important in defining channel morphology in this case. Finer material which is not re-deposited in the channel, and therefore moves through the 82 system very quickly, is referred to as wash material. Wash material can be important from a channel-forming perspective for both sand and gravel rivers. This material usually comprises a significant portion of material involved in bank erosion. In addition, it can be deposited on the floodplain. 4.3 BED LOAD TRANSPORT 4.3.1 BAGNOLD-TYPE RELATION In order to simulate realistically the development of the fluvial system over large temporal scales, it is essential to adopt relations that provide reasonable estimates of transport rates, yet remain manageable from a computational perspective. There are a number of physically-based equations available for the estimation of bed load transport rates. Most of these equations incorporate a level of detail that is not resolvable at larger scales. Moreover, a thorough investigation of the ability of such equations to simulate sediment transport rates over significantly longer time scales than individual flood events has not been undertaken. Gomez and Church (1989) evaluated a number of bed load transport relations over individual flood events and found that the Bagnold formula performed well in comparison to other formulae. The Bagnold stream power correlation is an appropriate formulation for landscape evolution modelling as it incorporates parameters that are resolvable at the landscape scale. It was decided to pursue the Bagnold form of the stream power relation, as opposed to the more general form of equation 4.2, as the former has been much more widely tested than the latter. The Bagnold bed load formula was developed in a series of papers (Bagnold, 1977, 1980, 1986), undergoing significant modifications along the way. equation is (Bagnold, 1980; Church and Gomez, 1989): The final form of the Bagnold 83 3/2 Ys CO-Cd (d/dJ "(D/D y" 0 2 b ref l Ys ~Y ref (4.6) ref wherein ib is specific bed load transport rate (dry weight), y is the specific gravity of fluid, y is s the specific gravity of sediment, co is specific stream power, a>o is critical specific stream power, d is depth, D is characteristic particle size and the subscript ref refers to some reference value obtained from a reliable data set (Bagnold used data from Williams, 1970). Ys I (.Ys ~ Y) w a The term introduced by Gomez and Church (1989) for the conversion of immersed bed s load weight to dry weight. The latter measure is the standard for most transport formulae. The following variables also require definition: co = p dSu w = p QS/w w (pw is density of fluid, S is slope, u is velocity, Q is discharge, w is river width) co = 5.75{0.04(^ -r)p f (g/pJ D /2 0 /2 logfyld/D) 3/2 w h f = 0.1; (co-oh\ef = 0.5; d f = 0.1; D re re ref (g is the acceleration of gravity) = 0.0011 For the purposes of the present study, the objective is to derive a simplified form of the Bagnold relation that eliminates many of the details found in the constants of the original equation, instead collapsing them into a scaling coefficient. Such a relation would, in essence, represent a scale relation between bed load and discharge which can be applied to rivers of various magnitudes. Width, depth and slope are all correlates of discharge and, therefore, form the basis to expect such a scale relation. The equation also includes grain size, which is a key factor in providing resistance to transport by the material. The general form of the relation is: ib =f (discharge, slope, width, depth, grain size) (4.7) The notion that a Bagnold-type relation may be a scale relation of the system, that does not necessarily reflect the fundamental physics of entrainment and transport, does not present a concern for its use in landscape evolution models. The time scales of such models are 84 sufficiently large that the details of small-scale physics are not resolvable. Therefore, it is sufficient that the effect of the phenomenon be appropriately preserved. A seemingly remarkable graph, which suggests the strength of the Bagnold bed load correlation, is presented in Bagnold (1986) (modified here as figure 4.1). In this graph, excess stream power (CO-COD) is plotted against which is the bed load transport rate adjusted to a common flow depth and grain size according to equation 4.6: f i \ h* = h \ ref d J 2/3/ xl/2 ' \ D ^ D ref (4.8) J The data collapse onto a straight line, which is the expected result if relations amongst these variables hold. Bagnold (1986) reassures us that : "In compiling figure 1 there has been no discarding of awkward evidence, and all the data for it were published before the idea of a conversion factor appeared in Bagnold (1980)". Based on the apparent strength of the Bagnold bed load correlation, an investigation of its development was undertaken. This investigation and the subsequent re-analysis of the formula are given in Appendix II. The re-analysis of the formula is based on an extensive data compilation consisting of transport data for flumes and rivers covering a range of discharges between orders of magnitude 10" m /s and 10 m /s. The final form of the relation chosen for 3 3 3 3 use in the present model is: i„ = 0.005{co-coJ d- Dl2 x U2 (4.9) The dimensions of this equation are not balanced, making this a scale correlation, rather than a physically meaningful relation. 4.3.2 VARIABLES REQUIRED IN T H E BAGNOLD RELATION The Bagnold-type relation chosen for the present channel submodel is such that 85 10 6 10" 41 T -icr 4 1 1 10' 2 1 1 10° 1 h 10 2 Excess stream power (kg/m s) Figure 4.1 The Bagnold relation. Modified from Bagnold (1986). 86 discharge, slope, width, depth and grain size must be defined for the calculation of sediment transport. The determination of these variables in the model is discussed in the following sections. Their definition in a numerical model requires either: (i) appropriate equations for their calculation, for which all terms in the equation can be resolved or (ii) some justification for the assumption that a particular variable can be treated as an independent variable. Discharge and slope calculations are based on a combination of topographic configuration and suitable equations in the channel submodel. Grain size is determined on the basis of a simple relation. These three variables together determine the channel pattern type for a particular set of conditions (i.e., single-thread vs. braided). Channel pattern type is required for the selection of appropriate hydraulic geometry equations for the calculation of width and depth. Once all of these variables are defined, the bed load transport rate can be calculated. 4.3.2.1 DISCHARGE Discharge values are required for the calculation of stream power, channel pattern, depth and width along the stream channel, making it a critical variable in the channel submodel. Some hydrologically based reference discharge is needed for calculations of discharge values along regional channel networks. It is expected that the mean annual flood should transport significant volumes of sediment. As channels decrease in size the most significant flows become rarer, and comparably more powerful. At the limit, it is about the order-of-magnitude 100 year flood in debris flow channels. For simplicity the 2-year flood is chosen as the reference discharge for the present study. Further elaborations of this basic approach should include the option of varying the return flows used in the model over time and space. In the absence of a model to simulate the hydrology of drainage basins, a discharge/contributing area relation may provide an adequate estimate of discharge for landscape 87 evolution models. Contributing area can be resolved in the model based on topographic configuration and an assumption that drainage occurs in the direction of steepest descent. The standard form of the scale relation between discharge and contributing drainage area is: Q =aA (4.10) b ref where Q / is some reference discharge and A is contributing area. The exponent b has a value re less than 1 in most cases. This suggests that during transmission of the reference flow to downstream portions of the basin, some amount goes into storage along the route or runoff production is not uniformly intense. The coefficient, a, reflects the hydroclimatological and geological conditions found in the particular study region. An analysis of regional variation in the discharge/area relation across British Columbia was undertaken by Church (1997). Discharge data were analyzed for a period of approximately constant climatic conditions from 1965 to 1984. The relation between mean annual flood and drainage area for British Columbia was found to be: Q =kA°- (4.11) 6S maf where Q f is the mean annual flood and k, which is the unit area runoff, varies with climate and ma geology across the province. Other studies have found that the value of the exponent is somewhat higher, with values ranging from 0.75 to 0.80 (Dooge, 1986). The value may be somewhat lower for British Columbia because of the steep landscapes for which the greatest runoff generation occurs in mountain headwaters. An exponent of 0.70 is chosen for the present model. Church (1997) found that coastal regions of British Columbia, which experience higher rainfall rates, show considerably higher values of k than the interior regions of the province. The analysis shows that the coefficient k varies by 2 orders of magnitude (lC'-lO ) with a median 1 value of about 1. The chosen k value for humid climates should be at the upper end of this range. 88 For particularly arid regions the designated k values should be lower than the values found in Church's analysis. 4.3.2.2 RIVER GRADIENT Gradient is calculated on the basis of grid elevations in the model. The locations of channel points do not necessarily coincide with grid points defined by the basic grid network. Therefore, each channel point is assigned the elevation of the nearest model grid point. A best-fit exponential curve is fit to assigned channel elevations. This approach ensures a smooth transition from lower to higher gradients. Unrealistically abrupt changes in gradient can adversely affect bed load transport rates calculated using the Bagnold equation. 4.3.2.3 CHANNEL PATTERN For a given discharge, a braided river has a wider and shallower channel than a singlethread channel. Therefore, the relations between discharge and width/depth vary according to channel pattern. In order to select appropriate hydraulic geometry equations for the model, it is necessary to first define the channel pattern. The variation is expressed in distinct coefficients in the empirical equations of hydraulic geometry (Leopold and Maddock, 1953) (see section 4.3.2.4 for further details). There is a continuum of channel patterns found in nature, of which multiplechannel and single-thread are just two possible end-members. However, differences among hydraulic geometry relations for different channel types has not been studied systematically. Therefore, the two general categories, single-thread and multiple-channel, will be retained in the present study. Rivers of both patterns can be composed of either sand or gravel. It is expected that the criterion will be different for sand and gravel rivers. A braided channel (multiple thread) is herein defined (Church, 1996, pers. comm.) as a "channel which, at low and moderate flow, 89 divides about channel bars which become drowned at high flow. The bars are unvegetated and more or less transient, so that alignments of individual channels shift frequently and rapidly." It has been found in previous studies that a discharge/slope criterion can be used to distinguish between channel pattern types. Leopold and Wolman (1957) found a transition from meandering to braided channels based on the equation: S = 0.0125Q /' (4.12) 044 b where Qbf is bankfull discharge. Discharge and slope, when combined, are two of the important terms in the calculation of stream power. Therefore, this equation suggests that at some threshold stream power, a river undergoes a transition from a single-thread to a braided pattern. In this analysis, data compiled by Church (1996, pers. comm.) and the meandering river data of Neill (1973) are analyzed in order to define a channel pattern criterion. The reference discharges accepted in the study of Church (1996, pers. comm.) were (i) mean annual flood; (ii) bankfull discharge; and (iii) the highest observed flow in relatively short records (used in only a few cases). Grain size data are based on D50 or D90 sizes of surface samples. The criteria for both sand and gravel rivers show a fairly sharp transition between regimes. A common slope of -0.5 is assigned for both types of river. A slope of -0.5 is not optimum for either data set. However, considering the restricted data base, it is sufficiently close in both cases and is therefore chosen as the common exponent. The criteria are defined by the equations (figure 4.2): Sand rivers: S = 0.01£T = O . K T 05 Gravel rivers: S = 0.07 Q~° = 4.9Q 5 (4.13a) 10 -10 (4.13b) In this case, the relation between slope and stream power is derived from the relation: S ocQ' 2 1 (4.13c) Figure 4.2 (a) Channel pattern criterion for sand rivers. Gradient is represented as slope tangent. 0.01 • Sand braided o Sand meandering m 0.001 S=.01Q-° • 5 c co CD CD 0.0001 0.00001 10 100 1 000 10 000 100 000 Discharge (m /s) 3 (b) Channel pattern criterion for gravel rivers. Gradient is represented as slope tangent. 0.01 0.001 0.0001 10 100 Discharge (m7s) 1 000 10 000 100 000 91 4.3.2.4 HYDRAULIC G E O M E T R Y RELATIONS (WIDTH/DEPTH) Meandering/Wandering Rivers Hydraulic geometry relations are an empirical set of equations that link changes in the river width, depth and velocity with changes in some reference discharge (Leopold and Maddock, 1953). Equations can be established either for changes that occur at-a-station or downstream. It is the latter case which is of interest in the present study. The relations have the form: w = aQ (4.14a) d = cQ f (4.14b) v = kQ m (4.14c) h Continuity necessitates that the product of the coefficients must equal 1, and the sum of the exponents must also equal 1. These equations imply that the river can adjust and scale itself to the discharge regime. Therefore, these equations ideally should be applied only to alluvial rivers, and not bedrock rivers. A distinct set of relations should be expected only in situations where hydrology and alluvial sediment characteristics are the major variables controlling the hydraulic geometry and other factors such as geology and physiography are kept approximately constant (Church, 1980). Within a setting, only discharge and grain size change. Church (1980) compiled discharge, width and depth data for rivers, canals and controlled experiments for single-thread channels. A re-analysis of this data is undertaken for the present study. Church (1980) selected data points only when flows were sufficiently regular or flow patterns were stable. This requirement is necessary for a chosen reference discharge to provide an adequate flow index. For river data, reported reference discharges are either the 2-year flood or bankfull flood. Church (1980) recognized that width and depth can be measured with greater precision than discharge. This violates the critical assumption in regression analysis whereby all error is attributed to the dependent variables. Therefore, the approach adopted by Church (1980) 92 is followed in which the functional analysis is performed. The error is instead allocated to the variable discharge in this case. For data sets in which the Froude number or other potentially significant controls on channel development changed, regression analyses were performed for several data subsets. This step ensures the elimination of anomalous data from the analysis. In addition, a correction factor is calculated to account for the bias introduced in the double log transformations. Plots of width and depth versus discharge for all data points are presented in figure 4.3. The results of the regression analysis for each data set show fairly strong relations with R values ranging from 0.67 2 to 0.995 (table 4.1). The mean values of the width and depth coefficients and exponents for both the sand and gravel rivers were calculated (table 4.2). An independent t-test was performed to see if there are significant differences between the width and depth exponents and coefficients for sand and gravel channels. There is no significant difference at the 0.05 significance level between the width and depth exponents for the two types of channels. While there is no significant difference between the depth coefficients, there is a significant difference found for the width coefficients. For the final set of relations it was decided to adopt the mean values of width/depth coefficients and exponents (refer back to table 4.2): Sand rivers: width = 7 . 4 0 ° (4.15a) 52 depth = 0.340° 42 Gravel rivers: width = 3.1(2° (4.15b) (4.15c) 53 depth = 0.22(2° 40 (4.15d) The data of Bray (1972, 1979) are based on Alberta rivers, mostly on the east slope of the Rockies, which tend to show higher than average widths for a given discharge in comparison to other rivers. When applying the channel submodel to British Columbia basins, it is most appropriate to use the hydraulic geometry relations of Bray. 93 o Stebbings (sand)* • Ackers (sand) Figure 4.3 (a) Relation between discharge and width. A Wolman and Brush (sand) o Simons and Bender (sand) • Lane and Carlson (gravel) Sand channels are open symbols and gravel channels are filled symbols. a Bray (gravel) A Charlton (gravel) a Emmett (gravel) * flume studies are in italics 1 000 0.01 0.00001 0.0001 0.001 0.01 0.1 1 10 100 1 000 10 000 Discharge (m /s) 3 (b) Relation between discharge and depth. 10- 0.01 0.00001 0.0001 0.001 0.01 0.1 1 Discharge (mVs) 10 100 1 000 10 000 94 Table 4.1 Regression analysis for hydraulic geometry relations. Study Stebbings (1963) Ackers (1964) Wolman and Brush (1961) Simons/Bender data in Simons and Albertson (1963) Lane and Carlson (1953) Bray (1972;1979) Charlton et al. (1978) Emmett (1975) Relation w= 11.7 Q d = 0.124 Q w = 6.08 Q d = 0.435 Q ° W-8.10Q d = 0.353 Q 0 3 5 0 5 0 4 4 0 5 1 0 4 7 w = 3.89 Q d = 0.448 Q ° 0 5 4 w = 3.22 Q d = 0.283 Q° w = 4.29 Q d = 0.191 Q w = 2.50 Q ° d = 0.197 Q w = 2.60 Q ° d = 0.201 Q R 0.995 0.98 0.98 0.94 0.83 0.67 Comments • Flume/sand 2 0 5 2 4 2 0 5 1 3 4 0 5 4 0 3 8 52 0 4 7 5 6 0 4 0 • flume/sand • flume/sand 0.95 0.96 • canal/sand 0.96 0.93 0.97 .0.92 0.87 0.81 0.89 0.81 • canal/gravel • river/gravel • river/gravel • river/gravel Table 4.2 Mean values of coefficients and exponents for sand and gravel rivers. Sand Gravel p-value Width Exponent 0.52 0.53 0.31 n = 4 for each category Depth Exponent Width Coefficient 0.42 7.4 0.40 3.1 0.55 0.046 Depth Coefficient .34 .22 0.22 95 Braided Rivers Drage and Carlson (1977) assessed hydraulic geometry relations for major sub-branches of braided streams in the Brooks Range in Alaska and the Yukon Territory: w = 4.66C/ (4.16a) 47 (4.16b) d = 0A3Q 0M The different coefficients for braided rivers than for gravel rivers emphasizes that the subchannels of these rivers are characteristically wider and shallower (Drage and Carlson, 1977). However, the width and depth can be calculated only when the discharge of an individual anabranch is known. It is first necessary to appropriately divide the known discharges among the appropriate number of anabranches for an individual river. Howard et al. (1970) studied a range of rivers in the United States and reported a relation between the number of sub-channels, discharge and gradient of the form: E = 0.675°Q ° M (4.17) 29 maf where E is the average number of sub-channels and Q j is the mean annual flood in m /s (it is 3 ma assumed in the present model that the mean annual flood is approximately equivalent to the 2year flood used in the calculation of other variables). It is also assumed in the present model that a fractional value, for example 2.4, implies that the number of sub-channels varies between 2 and 3. The discharge is then divided by the lower number of sub-channels (in this case 2). This discharge value is inserted into equations 4.16ab to obtain values of width and depth. The value for width is multiplied by the number of sub-channels to obtain the total active width. Following these calculations, the discharge is then divided by the greater number of sub-channels (in our example, the value is 3) and a similar procedure is followed to obtain width and depth. In order to calculate the proportion of reach length represented by the lower number of 96 sub-channels, the following equation is used: Proportion = 1 - (number lowermlu - integer(number subchannels )) lowersubchamels (4.18) The percentage of reach length represented by the higher number of sub-channels is determined by: Proportion =(l-Proportion uppervalue lower ) value (4.19) These weights are then used to calculated the weighted averages of width and depth for the lower and upper number of sub-channels, which represent the final values used in the model. 4.3.2.5 GRAIN SIZE Grain size can affect sediment transport due to its role in determining flow properties near the bed and by providing grain inertia. It has generally been supposed that there is a decrease in grain size along rivers in the downstream direction. Fining occurs in response to selective transport and/or abrasion of material. The change in grain size along gravel-bed rivers has often been found to fit an exponential relation (Church and Kellerhals, 1978): D = D exp(-a x) 0 D (4.20) where D is grain size, Do is a characteristic particle size at a downstream distance x=0 and a is D the diminution coefficient. In smaller channels which are strongly coupled with hillslopes, colluvial inputs and the existence of obstructions such as large organic debris can interrupt the pattern of systematic downstream fining of channel-bed gravels (Rice and Church, 1995). Rice and Church (1995) found that, due to random occurrence of log jams, grain sizes varied significantly and unpredictably over relatively short distances. Rice (1996) looked at overall fining over distances of about 100 km. He found negligible fining over a distance of about 100 km due to tributary resetting of grain size at this scale of 97 study. Rice (1996) found that the complexity of grain size structure could be explained by: (i) tributary inputs (ii) non-alluvial sediment inputs and (iii) the legacy of Pleistocene glaciation. However, within links delimited by significant lateral sediment sources, a downstream fining structure was detected. The importance of tributary inputs in resetting the grain size distribution along rivers is also recognized by Church and Kellerhals (1978). Research at larger scales suggests that geology may significantly affect the breakdown of particles in the river. For example, Shaw and Kellerhals (1982) studied reaches of order 10 km 3 in Alberta. The rivers studied by Shaw and Kellerhals (1982) display the following overall characteristics. Grain size changes in mountain reaches are highly variable with some suggestion of downstream coarsening. The central reaches show exponentially decreasing grain sizes and, finally, lower reaches are predominantly sand bedded. In addition, they found that quartzites break down into larger particles than limestones and, hence, larger material was found in river reaches with a high quartzite content. The study of downstream changes in particle size in natural rivers, which are subject to sediment inputs from various sources, is in its relatively early stages. Furthermore, such research has been conducted in order to determine spatial patterns at a particular point in time. At present there is no basis on which to model numerically the changes in sediment size over time along the river system. Therefore, one option for defining grain sizes in the model is to assume that grain size is an independent variable in the model that remains approximately constant over time at a location. When there are no real grain size values available for guidance in the setting of model grain size values, some mild downstream fining following a negative exponential model is an option. Grain size can be increased at locations where a tributary enters the channel. Consideration should also be given to rock type within the basin when defining grain sizes. Alternatively, sediment transport theory implies that there is some competence limit, 98 reflected in the largest grain size able to be moved, in a channel which is dependent on stream power. Therefore, grain size may be computable from relations amongst grain size, discharge, gradients and channel width suggested by several researchers (Griffiths, 1981; Parker, 1979; Henderson, 1963; Kellerhals, 1967). The relation, in each case, has the form: w where c is a coefficient. The ability of this equation to predict channel gradient is tested using the data set of Gomez and Church (1988). A coefficient of 1.26 was obtained by optimizing the value of c such that the mean of the calculated/observed values is equal to unity (figure 4.4). The relation between calculated and observed grain sizes shows that although the bias is small, the precision is very poor. In this case, the low precision might reasonably reflect stochastic effects which influence grain size. 4.4 BED LOAD TRANSPORT IN STEEP FLUVIAL REACHES The data used in the derivation of the bed load equation in this study have maximum slopes of about 0.01 (approximately 0.6°). Therefore, the bed load formula can only be applied confidently to channels with relatively low gradients. As gradients in natural stream channels increase above about 1° or 2°, channel characteristics begin to change. High-gradient streams are made up of step-pool sequences, which are alternating reaches of relatively steep and gentle gradients (Wohl et al., 1997). Steps, which represent sudden drops in elevation, are composed of clasts which form imbricate clusters or line, or they may form as a result of the accumulation of woody debris in the channel. Steps have a high resistance to transport due to these structural constraints and are generally immobile. Pool reaches in these channels have much lower gradients and no structural constraints and, therefore, are more apt to be locations of sediment Figure 4.4 Observed versus calculated grain size. Data shown are individual points for data sets in Gomez and Church (1988). Data appear grouped as grain size for each data set is constant. 0 10 20 30 40 50 Observed grain size (mm) 60 70 100 transport. Sediment transport processes in high-gradient streams are not understood. To date, there has been only minimal research investigating these phenomena. Wohl et al. (1997) explored sediment transport in step-pool reaches using a revised version of the Meyer-Peter and Muller equation for bed load revised for high-gradient streams (Smart, 1984). However, this equation was not designed for the complicated hydraulics of step-pool channels and hence results should be considered exploratory in nature. Transport rates in step reaches averaged about 55% of rates calculated for pools. High-gradient streams in coastal British Columbia drainage basins, the location of study basins modelled in this thesis, are dominated by log jams. These log jams create a series of sediment wedges behind them which, like the pool locations in step-pool streams, have lower than average channel gradients. In this initial model, it is assumed that fluvial transport in steep reaches is dominated by sediment wedges. Therefore, the effective transport gradients are reduced to values below average reach gradients. Given the exploratory nature of this model component and the lack of previous research, this simple assumption seems reasonable as a first approximation. Appropriate adjustments must be made in the model to account for this transport behaviour. In the present study, it is assumed that bed load transport operates normally when gradients are below 1°; that is, the effective transport slope is equal to the average slope. Above this gradient, high-gradient channel characteristics are assumed to emerge. This criterion is supported by a study of bar development in gravel rivers by Church and Jones (1982), in which they estimated the range of gradients above which bars do not develop under ordinary circumstances. The value at the lower end of this range was about 1°. Wedge features in coastal British Columbia develop in streams with average gradients as low as 1°. Above this criterion it 101 is assumed that sediment wedge gradients represent the hydraulically important channel gradient for sediment transport. Hogan (1997, pers. comm.) measured average channel gradients and associated gradients of sediment wedges upstream of log jams in several basins in the Queen Charlotte Islands, British Columbia (figure 4.5). A nonlinear functional model was fit to the data and the resulting equation used in the model is: Gradient = IXAAverage gradient ^ 1 wedge (4.22) with an R value of 0.86 and the standard error of the estimate is about 0.16 log units. The latter 2 results indicates that in unlogged units results can be expected to be between 0.7 and 1.45 times the estimate of wedge gradient. The wedge is steepening faster than the average gradient. However, the wedge gradient equals the average gradient when the latter reaches a value of about 11°. This value is greater than the cutoff for fluvial transport, hence this does not present a difficulty in the channel submodel. This equation is used to calculate the effective transport gradient in high-gradient reaches. In addition, the hydraulic geometry of high-gradient streams is expected to change from that of lower-gradient rivers. Day (1969) studied the hydraulic geometry of steep channels in the Coast Mountains and found a width relation of the form: w = 0.95 Q° 76 (4.23) The equation given by Day to calculate cross-sectional area is re-arranged in order to obtain an equation for depth: d = mQ" Day then provides relations for the calculation of m and n: (4.24) 102 Figure 4 . 5 Wedge vs. channel gradient. Gradient is represented as slope tangent. 0.01 0.02 0.03 0.04 Channel gradient 0.05 0.06 0.07 0.08 0.09 0.1 103 ,0.47 m = 0.9W (4.25a) n =0.3w',0.17 (4.25b) Depth, in this case, is a conceptual depth as real depth will vary dramatically from pool to step. 4.5 S U S P E N D E D L O A D Suspended load often comprises between 90% and 97% of total annual load (Schumm, 1977). There is little basis on which to derive a physically-based equation for suspended load transport as suspended sediment measurements include wash sediment, which is supply determined. Hence, a more empirical approach based on regional sediment yield data is adopted (Slaymaker, 1987; Church et al., 1989; Church and Slaymaker, 1991). The sediment yield study of Church et al. (1989) for drainage basins in British Columbia was examined. They assessed changes in suspended load with increasing contributing drainage area. Suspended load data were compiled for 63 stations within the period 1966-1985. It was found that there is an increase in specific sediment yield until a drainage area of 30 000 km . Thereafter, specific sediment yield decreases. This finding contrasts with the conventional model in which sediment yield decreases downstream due to deposition. Church et al. (1989) suggest that the downstream increase in sediment yield is due to the erosion of Quaternary sediments along stream banks and valley sides. The main trend line (until 30 000 km ) has the form: 0.6 Specific Sediment Yield = aA (4.26) wherein specific sediment yield has the units Mg/km /day. The value of the coefficient is 0.003 when only undisturbed basins are considered. Once drainage area exceeds 30 000 km the main trend for undisturbed basins can be approximated by the equation: -0.6 Specific Sediment Yield = aA (4.27) 104 where a is equal to about 700. The sediment yield pattern is for naturally disturbed basins in British Columbia and provides a way to consider regional disturbance in the channel submodel. Specific sediment yield is calculated at points along the channel length. The specific sediment yield is then multiplied by the contributing drainage area to get the total sediment yield transported past a certain river point. Data are then converted into annual volumetric transport rates. A bulk •a density of 1800 kg/m is assumed. The difference in suspended sediment transport rates between successive points defines the net change in suspended load for a reach. Based on this equation, there is expected to be net erosion of material moved as suspended load along the river up to a contributing drainage area of 30 000 km . 2 4.6 B E D R O C K E R O S I O N The importance of bedrock erosion lies in its ability to carve the bedrock valley thereby altering the underlying configuration of the landscape upon which weathering and redistribution of sediments occur. Tucker and Slingerland (1994), in their landscape evolution model, calculated the rate of lowering by bedrock erosion for a grid cell using a stream power/erosion relation: EBR = (k /Ax) QS (4.28) b where kb is the proportionality constant for bedrock channel erosion (L ) (kb values range from -1 10" to 10" m" in their initial model runs) and A x is the space step. Therefore, the lowering of a 5 4 1 cell is proportional to stream power scaled by the grid cell size. Rosenbloom and Anderson (1994) calculate bedrock incision using the equation: 105 (4.29) where x is upstream distance, y/ is a constant describing incision efficiency (U'T" ) and K is a 1 c diffusion coefficient. In this case, bedrock erosion is assumed to be proportional to stream power; in addition, a diffusion component is incorporated into the equation. Area is held constant in the application of this equation to the model; the effects of changing area contributing to discharge are assumed to be minimal. The values of a and K are determined from modelling c of the stream profiles in the marine terraced landscape in Santa Cruz, California and are assumed to be spatially and temporally constant. The model was found to be insensitive to channel diffusivities in the order 10-500 m kyr" . The incision efficiency constant was found to be about 5to7x 10- m-'kyr" . 7 1 Anderson (1994) calculated bedrock incision using a basic stream power/erosion relation: (4.30) where c is an empirically determined constant that reflects the proportion of stream power expended in channel incision and R is the average precipitation rate in the vicinity. The equation can be adjusted to account for debris flow activity in the upper channel network by making the change in elevation dependent on S , 2 Seidl et al. (1994) proposed that bedrock erosion is proportional to stream power. Assuming that discharge of peak runoff events scales directly with contributing drainage area then the relation, which is related to the basic stream power/transport given in equation 4.2, is: n (4.31) Seidl et al. (1994) assessed bedrock erosion rates for Hawaiian channels. They assumed that the exponents m and n were both equal to 1. Linear regression analysis was performed to 106 determine an equation for bedrock erosion: dh - — = 6.1xl0- +6.3xl0-"AS dt (4.32) 6 where the erosion rate is given in m/yr and drainage area in m . 2 Seidl et al. (1994) re-analyzed the bedrock erosion data of Hack (1965) collected in Michigan and found the relation: dh = 9.8xl0- +3.5xl0" AS 4 9 (4.33) dt It was concluded that the lower coefficients of the Hawaiian channels indicate that the rivers are downcutting into material of greater resistance (basalt) than the Michigan rivers (glacial till and bedrock; sandstone and shale). Therefore, in the present landscape evolution model the appropriate equation is selected on the basis of the overall resistance of bedrock in the chosen study area (equation 4.32 for Group A drainage basins and equation 4.33 for Group B drainage basins). 4.7 D E B R I S F L O W A C T I V I T Y Debris flow activity is important in the coastal regions of British Columbia as it provides an efficient mechanism whereby material is transferred from hillslopes to the river system. Material that collects in gullies is periodically flushed into channel zones of lower gradient. Moreover, gullies are an important roughness element in such landscapes. They provide strong downslope, topographically linear perturbations in the landscape. Rood (1984, 1990) compiled debris flow events for his study basins in the Queen Charlotte Islands, British Columbia. The present study is concerned with rates of natural landscape evolution and, therefore, debris flow rates were assessed for the 8 basins in which no logging occurred. The number of debris flows in each basin was recorded and divided by the 107 basin area in order to obtain the number of flows per km over the 40 year period that the 2 inventory covers (table 4.3). These values are then divided by this 40 year interval in order to determine the annual rate of debris flow occurrence. The average of these values, 0.01 flows/km yr, is assumed to be representative of natural debris.flow rates in the Queen Charlotte Islands. The number of debris flows in a particular drainage basin is obtained by multiplying this value by both basin area and the number of years in a model time step. These values are probably in the upper region of values expected in coastal regions of British Columbia because of the very high rainfall rates. A range of values from 0.001 flows/km yr to 0.01 flows/km yr is 2 2 probably appropriate for coastal British Columbia Based on initiation angles reported in this study, a value of 25° is selected as the threshold channel gradient for debris flow initiation (figure 4.6). The initiation of debris flows in the model is treated as a random occurrence, as the factors that make one site more prone to failure than another site are not resolvable at the scales of this study. Moreover, over the significant time periods considered in the present model (i.e., 10 -10 years), most gullies above 3 5 this threshold will experience debris flow activity, thereby making the actual ordering of trigger initiation sites insignificant. Debris flows are assumed to evacuate material across their width and to a depth of about 0.1-1.0 m (typical of debris flows in the Queen Charlotte Islands). Rood assigned a scour depth of 0.5 m in the initiation zone for most of the debris flows in his analysis. Scour depths for debris flows are slightly greater in the transport zone than in the initiation zone (figure 4.7). In order to simulate the variability in scour depths, a probability function is used to select a scour depth for each debris flow event. approximated from the function: The probability of a particular scour depth can be Table 4.3 Debris flow activity in unlogged basins in the Queen Charlotte Islands, B.C. Basin Area (km ) #flows / 40yr Hangover Government Jason Inskip Windy Bay Marshall Head Matheson Head Burnaby Island 19.4 15.2 14.4 13.4 18.3 8.6 10.0 11.0 2 #flows/km yr 3 6 10 3 1 6 0 #flows/km 40y r 0.15 0.39 0.69 0.22 0.055 0.70 0 3 0.27 0.0068 2 2 0.0038 0.0098 0.017 0.0055 0.0014 0.018 0 average: 0.0078 109 Figure 4.6 Initiation angle for debris flows in unlogged drainage basins, Queen Charlotte Islands, British Columbia. Data are from Rood (pers. comm.). 11 10 CO 9 c g 8 CD 7 '•4—* cu CO o 6 5 cu 4 n E 3 3 2 1 0 15-20 20-25 25-30 30-35 35-40 40-45 45-50 50-55 Initiation angle (degrees) Figure 4.7 Depth of scour in the transport zone for debris flows in unlogged drainage basins, Queen Charlotte Islands, British Columbia. Data are from Rood (pers. comm.). 10 9 8 c g 7 erv CO 6 ro CO _o o o «*— 5 4 l_ CD JZl E 3 z 2 .2-.3 .3-.4 .4-.5 .5-6 .6-7 .7-.8 .8-.9 Scour depth (m) .9-1.0 1.0-1.1 1.1-1.2 1.2-1.3 110 sum of depths where P is the probability of occurrence for the i scour category. There are assumed to be 8 th t scour categories in this case. The scour categories are numbered from 3 to 10 in order to represent scour depths ranging from 0.25 to 0.95 in the numerator of equation 4.34. The sum of the scour depths is obtained from the equation: 10 sum of depths = ]T(0.h'-0.05) (4.35) i=3 It is assumed that gully width increases as contributing drainage area increases, although no specific study of this relation has been undertaken for debris flow channels. Therefore, the width relation defined by Day (1969) for high-gradient channels is applied to debris flow channels in the present model. Once a debris flow is initiated material continues to be entrained until it reaches a zone where the slope falls below 8° (typical of coastal regions of British Columbia). The total volume of the debris flow event is calculated by multiplying scour depth by gully width along the length of the entrainment zone (defined as the distance from the initiation point to the start of the deposition zone): volume = j scour depth x width ) t i (4.36) 1=1 where j is the total number of reaches which are subject to debris flow erosion. The material is assumed to deposit over a distance that is related to the magnitude of the event. Typical lengths of deposition range from about 0.1 to 1 km (Church, 1997, pers. comm.). Unfortunately, Rood did not analyze deposition lengths of debris flows in his study. Based on the magnitudes of the debris flows in Rood's study, a relation is suggested which may provide reasonable estimates of Ill deposition length: deposition length = debrisflowvolume/10 (4.37) The material begins to deposit downstream of the point at which the gradient is lower than the debris flow stopping angle. The volume of deposition is calculated for each successive reach in a downstream direction from the stopping point. The deposition volumes are calculated according to the following equation until such a point that the cumulative distance between the stopping point and the lowermost reach boundary exceeds the deposition length: reach length volume deposited in reach = * debris flow volume deposition length (4.38) At this point the remaining volume is deposited in the lowermost reach such that the total deposition length equals the value calculated in equation 4.37. 4.8 S E D I M E N T B U D G E T F R A M E W O R K The channel system is divided according to a gradient criterion based on the predominant channel process in operation. Channel points having slopes greater than 8° are defined as debris flow channels and points below this criterion are defined as stream channels. When a reach is defined by two stream points in the model, it is assumed to be acted on by fluvial processes only, except in the event that a debris flow is deposited in this reach. Likewise, when a reach is defined by two channel points with gradients above the threshold, this reach is subject to debris flow activity. An appropriate framework for the channel system submodel is a sediment budget approach. The sediment budget for the case of fluvial transport processes is defined by the equation (figure 4.8): bed loadi -bed load +susp load -susp load +tributary t+debris flow = A storage n out in out inpu in (4.39) 112 Figure 4.8 Sediment budget for a fluvial channel reach. Any debris flow deposition is added to the change in storage. V (if tributary) in Bed load and suspended load calculations Bed load and suspended load calculations V (if tributary) in 113 which is evaluated for each reach along the length of the river. The Bagnold equation provides the fluvial transport rates for incorporation into this equation. Any sediment transport that enters the river channel from its tributary must be treated as an input to the reach at that location. The mass transport rate per unit width is obtained from the bed load equation. Therefore, in order to obtain the total mass transported at a location this value must be multiplied by channel width. In addition, this rate is calculated per second and therefore must be integrated over the total length of time the river spends in flood (in our case the mean annual flood was the chosen reference discharge) per model time step. The stream channel length between each pair of river points represents a channel reach. The calculated transport rates at either end of a reach represent the input and output transport values for incorporation into the sediment budget. The net change is the mass change per reach. In order to obtain the net volume change per reach, the mass value is divided by sediment bulk density. The sediment bulk density of 1800 kg/m used in this study 3 assumes a porosity of about 33% (the average of typical porosity values ranging from 25% to 40%). The incorporation of debris flow activity into the model framework differs somewhat from the approach for fluvial transport, as the actual transport rates are not calculated. Instead, the actual volumes of material eroded and deposited in a reach are known. Therefore, the changes in volume for a reach are immediately known; the intermediate step of obtaining this value by the difference in transport into and out of the reach is not necessary. 114 CHAPTER 5: THE V A L L E Y FLAT 5.1 INTRODUCTION The valley flat represents the interface between the hillslope and channel submodels. The fluvial system, which is located in this zone, is one of the most active locations in the landscape. The valley flat undergoes significant change during the course of landscape evolution. Moreover, valley fill and subsequent incision through this fill provide sedimentary records and evidence of erosional activity which contain critical information for the interpretation of landscape history. The valley flat is an important location in terms of sediment routing for several reasons: (i) The characteristics of the valley flat (e.g., valley width, relative position of the channel in the valley flat) determine the degree of interaction between hillslopes and channels. When sediment delivered from hillslopes to the valley flat enters the active channel relatively quickly, there is a . high degree of coupling between these systems. A lower degree of coupling results in longerterm storage of sediment before it is incorporated into the channel system. In the latter case, a wider valley flat will result in deposition of material along the foot of the slope. This material goes into storage, perhaps for long time periods. Stored material is eventually entrained by the river as the channel migrates across the valley flat and/or incises through the valley fill. (ii) Long-term episodes of fluvial aggradation and degradation involve significant changes in the elevation of the valley flat. As the elevation of the valley flat changes, so too does the base level for hillslope processes. In addition, the local relative relief and hence local gradients, are affected by these elevation changes. These factors are important in determining rates of hillslope erosion. 115 (iii) The fluvial system, which is located within the valley flat, contains short-term (active channel zone), long-term (floodplain) and very long-term (alluvial fans, terraces) sediment stores. The patterns of storage and re-entrainment of material in these two locations determines the transport rates of sediment through the fluvial system. The purpose of this chapter is to define "typical" patterns for the long-term erosion and deposition of material in the valley. These aggradational and degradational behaviours must then be translated into appropriately generalized rules for implementation in the landscape model. Appropriate simulation of these processes over long time scales allows relative increases or decreases in the width of the valley flat to be estimated. In addition, reasonable approximations of the elevation along the valley floor are important as they determine hillslope gradients and effectively increase or decrease the length of the hillslope from which material may be eroded. Few field studies have examined the actual processes and patterns of change that occur in the valley flat. This represents a significant shortcoming in geomorphological research as a greater understanding of this critical interface is required in order to improve our knowledge of long-term movements of sediment through the drainage system. Landscape evolution models have not explicitly invoked rules for the simulation of long-term aggradation or degradation in the valley. If tectonic effects on landscapes are the focus of a study then it may be appropriate to ignore the valley flat. More specifically if the time scale of interest is much greater than the basin diameter divided than virtual velocity, then details such as valley flat storage may no longer be important. Anderson (1994) considers only bedrock fluvial processes in his model, thereby not allowing for long-term aggradation or incision through these deposits, which may be appropriate given the importance of tectonics in his study. Other models (e.g., Tucker and Slingerland, 1994) allow for re-deposition of sediment by the river, but given the large size of grid cells the details of deposition are not resolvable. 116 5.2 CAUSES OF L O N G - T E R M AGGRADATION AND DEGRADATION Material which is deposited in the valley flat is derived from upstream channel reaches, tributaries, and laterally from the hillslopes. The erosion of material in channels along the valley bottom is initiated by either fluvial or debris flow processes. Net storage changes along the valley are determined by assessing the differences between erosion and deposition rates which occur in response to changes in transport capacity along the channel. Rivers may experience alternating episodes of aggradation and degradation over short and medium time scales. However, such episodes are "superimposed" on underlying longer-term trends (Schumm and Lichty, 1965). From the perspective of landscape evolution modelling, it is these latter, longerterm trends that are of particular interest. Topographic configuration, geology and climate affect long-term patterns of erosion and deposition in the valley flat. Given the relation that has been shown to exist between gradient and transport rate for some processes, the topographic configuration represents an important determinant of relative process rates occurring throughout the basin. Rock type is expected to affect the efficacy of process operation through: (i) its influence on weathering rates and (ii) its role in establishing sediment characteristics, which in turn affect transport rates (e.g., cohesive properties and grain size). Changes in climate over space and time modify the hydrological characteristics and/or vegetation types found within a basin. Fluvial transport rates are enhanced when a change in climate increases flood activity. Hillslope processes may respond to a wetter climate in one of two ways: (i) greater transport rates if vegetation type remains unchanged (greater groundwater pore pressures) or (ii) reduced transport rates if the amount of vegetation increases significantly (stabilization of the surface). Past phases of increased erosion may influence the routing of sediment through river basins for many years. For example, sediment yield patterns suggest that material eroded during 117 the last glaciation in British Columbia has resulted in large "pulses" of sediment which are still working their way through drainage systems (Church and Slaymaker, 1989). Increases or decreases in transport rates within a basin cannot result in phases of either basin-wide aggradation or degradation, except in two specific circumstances. If the basin is very small, such as headward tributaries, then there can be basin-wide degradation. The basin diameter is small in this case, making the time for evacuation of material from the basin relatively short (see equation 2.2, p. 16). Deposition in one location necessitates corresponding erosion at some other location. The inclusion of valley filling and channel incision into landscape models may contribute to an increased understanding of the relative importance of factors responsible for sustained periods of aggradation and degradation in valleys by providing an appropriate framework within which to formulate suitable research questions and perform sensitivity analyses. 5.3 V A L L E Y FILLING Material which comprises valley fill was originally deposited in one of several locations (Schumm, 1977): (i) the valley margin (hillslope deposits), (ii) the channel bottom (material at rest in the bed) (iii) the channel margin (point bar deposits) and (iv) the floodplain (overbank fine-grained deposits). In the present model, which represents a first modelling effort, the various fluvial processes which deposit material are not differentiated. As confidence is gained in the model's performance, details regarding fluvial deposition may be incorporated into the model. In this model, any material which is not defined as "bed load" (in the present study the bed load equation was calibrated on data > 1 mm) is assumed to be transported directly through the system once it is entrained. Therefore, it cannot be involved in the valley filling process. 118 However, in reality some proportion of this material is either incorporated into the valley fill as overbank floodplain deposits or as interstitial fill. For simplicity, this material is ignored in the channel submodel. The width of the valley flat increases as a channel aggrades. However, when material is deposited by rivers having a width significantly smaller than the valley flat, the sediment cannot simply pile up in the channel. If this approach is followed in a model then the eventual result would be the creation of mounds of sediment many meters high in the middle of the valley flat. The significant issue then becomes the consideration of how material is spread across the valley flat so that deposits are approximately level. The lateral migration of the river as it swings across the width of the floodplain may contribute to the even distribution of river deposits across the valley flat. In addition to distributing material across the valley flat, this process also mobilizes previously inactive sediment. An assumption is made in this model that any material which is deposited by channel processes is evenly distributed across the valley flat. It is essential to investigate recorded rates of movement in order to understand better the lateral migration of rivers. Factors which may affect bank erosion rates are the texture of the valley fill through which the river is migrating, discharge and incoming sediment load. Lateral migration rates provide information about the time scales over which the deposition of sediment occurs as the river moves across the valley flat. Hooke (1980) compiled bank erosion data for meandering rivers and documented lateral migration rates ranging from several cm through to several hundreds of meters per year. Most of the bank erosion rates in the data are restricted to rivers in temperate or continental regions in the Northern Hemisphere. Hooke (1980) also found that there is a nonlinear relation between bank erosion rate and drainage area. His data show that for a given drainage area, there is generally a 2 order-of-magnitude variation in the bank erosion rate. 119 The length of time required for a river to traverse its valley flat is dependent on both the rate of lateral movement and the width of the valley flat. The data of Hooke suggest that it may take anywhere from decades to thousands of years for rivers to complete a traverse of the valley flat. In any event, these time scales are shorter than the usual time scales associated with landscape evolution models. This suggests that the approach whereby sediment is "spread" across the width of the valley flat may be reasonable at the large scales of landscape models. Bridge and Leeder (1979) provided a useful compilation of floodplain accretion rates over time scales ranging from individual flood events through several thousands of years. They found rates of floodplain accretion ranging from 0.1 mm to 10 mm/year. 5.4 V A L L E Y INCISION THROUGH FILL Long-term fluvial incision through valley fill leaves its imprint on landscapes in the form of terraces. Incision may be episodic in nature, with depositional phases interspersed within longer periods of downcutting. Vertical incision may be accompanied by a lateral component of erosion. However, the relative roles of vertical and lateral erosion have not been studied in a rigorous manner and therefore an assumption is made in this model that the river maintains the same planform position as it erodes vertically downwards. When a river incises through fill (as opposed to bedrock) the adjacent slopes cannot maintain an angle greater than some maximum angle of stability. Nearby slopes erode backwards in order to maintain a gradient below this value and, as a consequence, the eroded material makes its way to the valley flat. Long-term degradation rates can be estimated from dated terrace deposits in conjunction with heights of the terraces above the present river level. The terrace data of Drozdowski and Bergland (1976) and Personius et al. (1993) indicate incision rates of about 0.5 mm/yr. 120 Degradation rates which are measured below dams over medium time scales may provide an extreme upper limit of incision rates (Williams and Wolman, 1984) (figure 5.1). The median value of the median degradation rates calculated for the study rivers is about 25 mm/yr. Natural rates are expected to be considerably lower than this value. 5.5 INCORPORATING V A L L E Y PROCESSES IN T H E PRESENT M O D E L 5.5.1 " V A L L E Y RULES" FOR T H E ONE-DIMENSIONAL M O D E L In the present study, an attempt is made to generalize appropriately the processes leading to valley filling. The width of the valley flat must first be determined as the deposit is "spread" across this width in the model. It is necessary to define a criterion which represents the maximum angle that can exist between two points for them to be considered a part of the valley flat. A gradient of 1° is chosen to represent this value in the model. The gradients between the river and its adjacent points are evaluated (figure 5.2). If the gradient between the river grid cell and its neighbouring grid point is below the criterion angle, then the gradient is considered negligible and the neighbour point is assigned to the valley flat. The procedure is repeated for any newly-assigned valley flat grid cells and their neighbours. The procedure continues until the maximum lateral extent of the valley flat is determined. The net cross-sectional area of the material to be deposited is then evenly distributed across the width of the valley flat. The calculation of sediment transfers along the channel requires the implementation of the channel submodel within a surface model run. Therefore, in one-dimensional model runs a representative value of channel changes must be inputted directly during model runs. An implicit assumption in this model is that compaction of the underlying sediment does not occur. It seems reasonable to speculate that any changes due to compaction that would occur are smaller than the error bounds of this study, thereby justifying their omission in the model. However, it may be 121 Figure 5.1 Incision rates below dams. Data are from Williams and Wolman (1984). ro cu CO X) o cu XI E 0-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80 80-90 90-100 Erosion rate (mm/yr) 122 Figure 5.2 Valley filling in the channel submodel. Gradient < 1 0 Gradient < 1° Lateral extent of valley flat Figure 5.3 Fluvial incision in the channel submodel. Time step: a is less than 35° so river incises at a width of 1 grid cell Time step: a is greater than 35° so slopes are unstable and when next incision occurs the incision width increases. Time step: Channel incision width increases and slopes stabilize. 123 beneficial to consider compaction in further studies which focus on the filling of valleys with sediment and, perhaps, in later versions of the model. The rules for valley incision require the definition of a gradient which represents the maximum angle of stability for slopes adjacent to the incising river. A criterion of 35° is suitable for a river which is incising through valley fill (see section 5.6 for the case of bedrock channels). Gradients between the river point and its adjacent points are evaluated to determine if this angle is exceeded. If the calculated gradient is below 35°, then river incision is initiated (figure 5.3). In the event that the criterion is exceeded, the angle between the river and the next successive point is evaluated. This procedure is repeated until the incision width is great enough for stability to be achieved. The cross-sectional distribution of net erosion is then determined according to the newly-defined incision width. 5.5.2 " V A L L E Y RULES" FOR T H E SURFACE MODEL Net channel changes are calculated in the channel submodel for surface (2-dimensional) model runs. The general principles for the "smearing" of fluvial deposits are similar to the onedimensional case, except that the algorithm to define the portion of the valley flat across which deposition occurs is more complex. In the algorithm, each river grid point is associated with its nearest neighbour on the grid network; this point effectively becomes the new river grid point. The gradients between each of the new river grid points and its eight neighbours are evaluated. If the absolute value of the gradient is less than the defined threshold for valley aggradation (in this model it is defined as 1°), then that point is considered to be a part of the valley flat across which the river will eventually migrate and deposit sediment. Gradients between the eight neighbours of each new point and the original river point are evaluated in the next step. If gradients are lower than the criterion then they, too, are assigned as 124 part of the valley flat. If a particular point has already been defined as part of the valley flat for this particular river point, then it is skipped. The procedure continues until all gradients exceed the criterion. This condition signals that the valley flat associated with a particular river point has been defined. The procedure is repeated for the next river grid point. A particular grid point may be assigned to several river points as long as the gradient criterion is met in each case. This situation implies that, given enough time, the channel length associated with each of these several river points migrates towards a particular location in the valley flat and deposits material. A similar algorithm is applied for the case of fluvial degradation in the surface model as for the case of aggradation. However, in this case a particular area contributes to river erosion when the threshold gradient defined for slope stability (35° in our case) is exceeded. 5.5.3 HILLSLOPE/FLUVIAL COUPLING Depending on the degree of coupling between hillslopes and channels, colluvial material that reaches the slope base may: (i) be deposited in the active channel and transported downstream relatively quickly or (ii) go into long-term storage along the valley sides. The length of time before this material is entrained by the channel system depends on the particular characteristics of the valley flat, such as valley width, river location in the valley flat and the rate of lateral migration by the river. Despite the importance of these processes in the routing of sediment, there is a paucity of research that explicitly addresses these interactions. In general, the width of the valley flat is expected to increase as channel order increases. The width of the valley flat is an important determinant of whether or not hillslope inputs directly enter the channel, or go into storage. 125 Whiting and Bradley (1993) made the following observations about the connectivity of hillslope and channel systems. In low-order channels, colluvial inputs clog the valley flat, which is usually narrow at these scales. Water flow is generally incapable of moving material in loworder channels due to the large grain sizes and relatively low discharge rates found in these channels. Debris flows are the primary transport mechanism in low-order mountainous channels. As order increases, the transporting capacity of the flow increases and flood benches may form. Landslide and debris flow deposits entering the valley flat move variable distances across its width, in some cases reaching the active channel zone. Increasing discharges are more capable of eroding slope base and slope base deposits. However, valley widths continue.to increase as channel order increases. In consequence, the potential for material derived from colluvial inputs to go into longer-term storage also increases. Because the rules in the present model allow the width of the valley flat to vary, the strength of hillslope and fluvial coupling can also vary. The gradient between the lower valley cell and the grid cell at the edge of the valley flat induces erosion of sediment according to the hillslope rules introduced in chapter 3. When the valley flat is narrow, this "valley edge cell" is equivalent to the river cell and, hence, material is transferred directly into the river. This immediate transfer of material into the river cell represents a situation of strong coupling. When the valley flat is wide, perhaps several model grid cells in width, the materialfromthe hillslopes is deposited into the cell at either side of the valley flat, and not the river cell. In this case, there is minimal deposition of hillslope material in the mid-valley zone, which indicates a low degree of coupling. 126 5.6 B E D R O C K E R O S I O N Bedrock erosion is an important process in many regions over the long time scales of landscape evolution. Bedrock erosion is an important process in badlands (Howard and Kerby, 1984) and in mountainous regions (Seidl et al., 1994). Furthermore, many reaches that are generally considered to be dominated by alluvial processes often show indications of significant bedrock control. A commonly held assumption regarding bedrock river incision is that erosion occurs primarily in a downward direction which results in very steep adjacent canyon walls. However, several studies have suggested that bedrock valley floors underlying valley fills in higher-order streams are often relatively flat or somewhat irregular in shape and not "V-shaped" (Burrin and Jones, 1991). Low order streams are more likely to display the characteristic "V-shaped" bedrock cross-section. The occurrence of bedrock strath terraces also indicates the potentially significant lateral component to bedrock erosion (Burbank et al., 1996). Crickmay (1974) draws attention to the potential significance of lateral, vertical and oblique corrasion by rivers through underlying bedrock, although his ideas have not been subject to rigorous evaluation. There is an unfortunate lack of knowledge regarding the relative roles of vertical and horizontal bedrock erosion by rivers. Therefore, it was decided to incorporate the simple rule into the model that bedrock erosion occurs only in a downward direction (see equations 4.32 and 4.33 for bedrock erosion on p. 106). As knowledge of bedrock erosion processes increases, consideration should be given to both the horizontal and vertical components of bedrock erosion. 127 CHAPTER 6: 1-DIMENSIONAL MODEL RUNS 6.1 M O D E L OUTLINE The hillslope and valley submodels presented in this thesis are investigated initially using a 1-dimensional (profile) version of the program. Diffusive evolution of the landscape is examined using the linear diffusivity found in the present study (0.1 m /yr) and a range of 2 diffusivities implemented in previous landscape evolution models (1-100 m /yr). The nonlinear 2 transport functions introduced in Chapter 3 are also examined. An optional sink/source term is incorporated into the model to simulate either the deposition or erosion of material by the river and associated valley filling or incision. Typical heights of aggradation or degradation occurring in channels are read into the model in this suite of runs. In later surface model runs, channel transport rates are determined within the model itself. The channel is assumed to be located at the lowermost point in the valley. Typical river widths expected for the profiles used in this study are smaller than the space step (the grid cell dimension). The morphological effects of erosion and deposition cannot be resolved for lengths smaller than the space step. The procedures for simulating aggradation and degradation in the model were outlined in detail in chapter 5. When there is net channel aggradation over a time step the material is effectively "smeared" across the width of the valley flat. When there is net channel degradation, the incision width increases until such a point that the gradient of the valley sidewalls falls below the threshold gradient. Partial differential equations are solved in the model using the "method of lines". An implicit procedure is used to solve the first term in the right-hand side of the linear diffusion equation: 128 dh_ &h = k^i + <f> k dt ~ dx 2 (6.1) wherein 0 is the river sink/source term. The term <f> is solved explicitly. In this method the second-order derivative is replaced by a series of algebraic equations and the time derivative is initially preserved in its original form. This results in a system of differential equations in time (t) and since there is now only one independent variable, the equation is an ordinary differential equation (ODE). An implicit Runge-Kutta scheme, in this case a subroutine called ddriv2.f created by Kahaner et al. (1989), is used to perform the integration and the new value of h is calculated. A flowchart of the model is given in figure 6.1. After the input data are read (including initial elevations and river input values), the input <p value for river activity is used in the aggradation/degradation subroutines and changes in elevation are calculated. These heights are then incorporated into the differential equations for the relevant grid points. The subroutine DERIVS is used to solve the right-hand side of the differential equation. At this point the major do-loop for the time steps is entered. Within each time step, a call is made to the integrator ddriv2.f, which in turn continues to call the subroutine DERIVS until the solution converges and a value for h is obtained. A new phi-value is entered at each time iteration and is incorporated into the appropriate differential equations. In all cases, model results were checked to confirm that continuity was preserved. 6.2 STUDY DRAINAGE BASIN It is preferable for model runs to be of limited spatial extent when initially evaluating model performance. In addition, an ideal study basin exhibits a simple topographic pattern, which enhances the interpretability of results. The model can be applied to more complex situations as R e a d in model variables, input elevation data and input river data t =0 Call "report" Writes initial data to output file Initial call to "derivs" "derivs" calculates the right- hand side of the differential equation for use in "ddriv2" | Initial call to "phi_calc" "phi_calc" distributes> input values of fluvial aggradation and degrad ation across the valley flat and calcula e values of 0 tjo = t t hi = t + dt Call "ddriv2" "ddriv2" makes internal calls to "derivs" as necessary and calculates model elevations t = t_hi Call "report" Call "phi_calc" calculates value of 0 for next iteration Figure 6.1 Flowchart of 1-dimensional version of model. 130 confidence is increased in its ability to reasonably simulate landscape evolution. The model was run for a drainage basin in the Queen Charlotte Islands, British Columbia for several reasons: (i) the hillslope analysis in this study is based on Queen Charlotte Islands' data; (ii) the channel submodel was created with coastal British Columbia drainage basins in mind, which are typically steep and humid as are basins found in the Queen Charlotte Islands; and (iii) there is a good regional knowledge of geomorphological processes operating in the Queen Charlotte Islands, which is supported by a substantial literature. Mosquito Creek Tributary was chosen as the study basin as it is a small watershed of area 5.5 km and exhibits a 2 reasonably simple topographic structure (figure 6.2). Slopes are steep and approximately rectilinear in this basin. A major channel runs along the length of the basin with a series of simple gully channels feeding into this main channel. Moreover, Mosquito Creek Tributary has been a research site in several other studies (Roberts, 1984; Rood, 1984, 1990) and, therefore, many aspects of this basin have been previously analyzed. The initial hillslope profiles (time = 0 years) were created using 1:20 000 TRIM (Terrain Resource Information Management) map data for the Queen Charlotte Islands, British Columbia. The two profiles investigated in this analysis are (see figure 6.2 for locations): (i) a major hillslope feeding into the main channel and (ii) a hillslope feeding into a typical gully. In both cases, the profile actually represents the elevations along the path of greatest descent which converge to a common point in the channel (not necessarily a cross-sectional profile perpendicular directly across the channel). In each case the profile is extended beyond the divide so that model behaviour at the boundaries does not adversely affect the results. Figure 6.2 Mosquito Creek Tributary drainage basin. 132 6.3 PROFILE DEVELOPMENT OF MAIN CHANNEL HILLSLOPES 6.3.1 LINEAR DIFFUSION 6.3.1.1 BASIC M O D E L RUNS The development of the hillslope profile feeding into the main channel is first assessed for the case of basic linear diffusion. In these model runs, fluvial processes are not considered. Initial values of model variables and diffusivities used in this analysis are given in table 6.1. The model is run using the creep diffusivity. Although a linear relation may not be strictly appropriate for creep processes, there is no systematic study of creep available on which to base an alternative relation. The hillslope analysis in chapter 4 indicated that linear diffusion is not an appropriate model for rapid landsliding in the Queen Charlotte Islands. Nevertheless, the profile model is run using the linear diffusivity found in this study and values of several other studies in order to allow comparison with the preferred nonlinear transport equation. The resulting profiles for model runs for each diffusivity value are plotted after (figure 6.3): (i) 0 years (initial profile) (i) 3 000 years (ii) 10 000 years (Holocene time-scale) (iii) 30 000 years and (iv) 100 000 years. The results immediately suggest that creep is a relatively ineffective process, even at long time scales. The morphological changes are not even visible at the resolution of the graphs and, therefore, are not illustrated. After 100 000 years the landscape shows only very minor changes in elevation with maximum erosion rates at the divides of about 0.2-0.3 m and a deposition height in the river channel cell of 0.12 m. These values represent average rates of change of order 1 mm/1000 yr and suggest that in landscapes for which creep is the dominant process, rates of change are very slow. The relative roles of lateral and vertical river activity in landscape evolution are enhanced in such regions. However, when the linear diffusivity value for landsliding obtained in the present study (0.1 m /yr) is implemented, there is much greater activity (figure 6.3a). Significant erosion on 2 Table 6.1 Input parameters for initial model runs. Variable Number grid points (nx) Space step (dx) Time step (dt) Final time Diffusivities (k) Value 37 100 m 100 yr 100 0 0 0 yr 0.0002 m /yr: Present study (creep) •y 0.1 m /yr: Present study (landsliding) 1 m /yr, 10 m /yr, 100 m /yr: Values cover range of other landscape evolution studies. 2 Elevation change due to channel processes (c/5) 0 2 2 (ill) U0!1BA9|3 (ill) U0I}BA3|3 (Ui) U0I}BA9|3 (Ol) U0!1BA8|3 136 the upper slopes and deposition on the lower slopes are observed. After 10 000 years the divides show decreases in elevation ranging from 8-13 m (800-1300 mm/lOOOyr). After 100 000 years, the two divides have eroded about 65 m and 81 m (650-810 mm/1000 yr). In these initial model runs the sediment supply is assumed to be unlimited (the weathering of bedrock is not considered), which results in relatively high rates of erosion at the divide. The lower slopes show aggradation rates slightly lower than the central valley points, while upper slopes show some minor erosion. The middle zones of the hillslopes represent, in effect, transport slopes. Material is expected to travel through these areas relatively efficiently, until it reaches the lower valley slopes where it may enter longer-term storage. This is the expected behaviour for essentially rectilinear slopes, such as those found in the Mosquito Creek Tributary profiles used in this study. There is net aggradation of about 6 m (600 mm/1000 yr) in the valley after 10 000 years and 50 m in 100 000 years (500 mm/1000 yr). In this model, the hillslope material deposited in the valley flat becomes incorporated into the valley fill. The valley filling rate associated with this deposition rate is 0.6 mm/yr. This value is at the lower end of floodplain aggradation rates given in Bridge and Leeder (1979). The relatively high aggradation rate suggests that significant amounts of hillslope material are transferred to the valley flat. In the real world, this material may go into long-term storage or be evacuated relatively quickly, depending on the transporting capability of the river and the characteristics of the valley flat (e.g., width). The profile is still recognizable, even after 100 000 years. The only major topographic differences are the changes in elevation for the upper and lower slopes, and the associated decrease in gradients along the slope length. However, as diffusivities increase to values ranging from 1-100 m /yr (these values cover the range of diffusivities used in other landscape evolution models; figure 6.3b-d), erosion and 137 deposition rates are considerably greater than changes expected to occur in reality. For a diffusivity of 1 m /yr , the elevations of the crests have decreased by about 65 and 81 m and 2 valleys have increased by about 50 m in 10 000 years (order 10 mm/lOOOyr). For diffusivities of 3 100 m /yr the entire landscape has been reduced to a relatively flat surface in only 10 000 years. The diffusivity value adopted in the present study appears to be more viable than values adopted in some other studies. 6.3.1.2 C H A N G E IN S P A C E S T E P The effect on model results of a reduction in space step from 100 m to 50 m was studied using the landsliding diffusivity of 0.1 m /yr obtained in the present study. The differences 2 between the initial profile for each of the space steps are shown in figure 6.4a. Some minor topographic irregularities are captured for the 50 m space step profile at time 0 (initial profile) which are not observed for the greater space step (figure 6.4a). In particular, there is a high point on the right-hand drainage divide which is not visible for the 100 m spacing. However, after the model has run for 100 000 years, these initial differences disappear and only negligible differences exist between the profiles for most points (figure 6.4b). The difference between the elevations for two model runs at corresponding points along the slopes are of order 10° m after 100 000 years. The difference in valley aggradation after 100 000 years is only 1.5 m. The differences between elevations for each space step are considerable for one of the divides (nearly 30 m), while the difference is 1.5 m for the other divide. These results suggest that diffusive hillslope processes are relatively insensitive when space steps are changed by a factor of 2. Any large-scale irregularities in the landscape which are captured at the 50 m space step and not at the 100 m space step are efficiently smoothed out in less than 10 000 years when a diffusivity of 0.1 m /yr is used. Thereafter, the landscape evolves in a similar manner. (oi) uoueAaig (uu) uo!ieA9|3 139 6.3.2 NONLINEAR TRANSPORT 6.3.2.1 BASIC M O D E L RUNS The nonlinear transport functions defined in Chapter 3 are implemented in the next set of model runs. Morphological changes for the Group A basins (resistant geology) show that landscape changes are significantly reduced from those for the linear diffusivity value of 0.1 m /yr (see figure 6.5a). The transport relation for Group A basins shows an upper limit transport rate of about 0.04 m Ira. yr, which restricts sediment transport rates to low values across all gradients. Elevation changes at the main divides are about 2.5 m after 10 000 years (about 0.25 mm/yr). Changes at the divides range from 20 to 24 m after 100 000 years, about one-third of the values obtained when the linear diffusivity of 0.1 m /yr is implemented. Due to the relatively 2 gentle gradients on the lower slopes adjacent to the grid cell containing the river (for which transport rates are about Om/m yr), sediment is deposited on these slopes and is not transported to the river grid cell. The central valley aggrades only 3 mm and 200 mm in 10 000 years and 100 000 years respectively. The lower slopes adjacent to the river grid cell aggraded about 2 m and 16 m after 10 000 years and 100 000 years. Nonlinear diffusivity for the Group B basins (non-resistant geology) shows somewhat greater activity (figure 6.5b) than for Group A basins, but changes in elevation still remain considerably lower than for linear diffusivity in the present study. Erosion values of about 8 m and 37 m occurred at the divides after 10 000 years and 100 000 years respectively. The central valley shows aggradation values of only about 4 cm and 2 m for these same two times, while aggradation values at points adjacent to the river were about 2.5 m and 25 m. The very high transport rates that the function defines at high gradients for Group B basins might lead to the expectation that elevation changes for this function would exceed those for linear diffusivity. This is not the case. The reason for this lies in the range of gradients found (w) uorjeA3|3 (w) U O H B A 3 | 3 141 in the Mosquito Creek Tributary hillslope profile and the associated diffusivities for both the linear and nonlinear case over this particular range of gradients (figure 6.6). If there are steep gradients in a basin, they are worked upon very quickly when the nonlinear transport model is implemented. Once these exceptionally steep slopes are eradicated, the gradients are no longer in the range of steep gradients for which high transport rates are expected. Throughout most of the 100 000 years of landscape evolution, the gradients lie in the range where the diffusivities for the linear equation defined in this model exceed the values for the Group B nonlinear function. The nonlinear function provided a considerably stronger fit, as shown in Chapter 3, for the Queen Charlotte Islands' landsliding data set. The general form of the nonlinear function provides the basis for a reasonable scenario of changes in transport activity in landscapes as they evolve. It is only after landscapes are subject to significant events which lead to increased gradients and an "unstable" landscape that hillslope erosion will proceed at a rapid rate. The operation of hillslope processes removes the large-scale regional instabilities. Anything that can increase gradients can lead to the "re-activation" of significant transport activity. Some possible causes of steepening are: (i) uplift events which lead to increased incision by the river and associated increases in hillslope gradients; (ii) deleveling uplift events which cause localized steepening of certain slopes by tilting; (iii) glaciation; and (iv) lateral river activity. In addition, changes in climate and geology may also lead to increased landscape activity by changing the nature of the transport relation and the threshold gradient of stability. This may involve a change in : (i) the resistance of landscapes (e.g., change in vegetation due to climate change; change in rock resistance due to exposure of a new surface) or (ii) the hydrological regime, which is a major factor determining transport activity. It can be argued that some locally steep areas will remain in a landscape, even after significant time periods have lapsed. In addition, steep slopes may be recreated by gullying or by 143 deep-seated landslides, neither of which is included in the present model. Because the Queen Charlotte Islands have been subject to recent glaciation and tectonic activity, and are subject to an exceptionally moist climatic regime, they are notably unstable and thus transport rates are very high. Because the transport functions in this study are based on a specific study region, they are not directly transferable to other regions. However, it seems likely that nonlinear relations may be representative of transport functions in other regions where there are steep, rough upper slopes and gentler, smoother lower slopes. Many regions of the world have been subject to recent glaciation and uplift events, so upper slopes may be operating in the "active" transport zones defined by nonlinear relations. Vertical and lateral erosion by rivers may destabilize lower slopes and contribute to further instability and transport. The interplay between hillslope and channel processes is particularly interesting in active landscapes as the slopes undergo periods of relative stability and instability. Instabilities in the landscape may create episodes of high transport rates. After enough time has elapsed for unstable slopes to be eradicated, the landscape relaxes into a state of "quasistability" and geomorphological activity proceeds at a reduced rate. Relatively ineffective creep processes become the dominant geomorphological transport mechanism on hillslopes once gradients are below the threshold for landsliding (or other rapid transport processes of significance in a particular region). But in comparison to the overall mass of the mountain range, transfers by creep processes appear not to result in obvious configurational changes according to transport rates estimated in this model. Furthermore, the persistence of Holocene-spanning colluvial footslopes and alluvial fans suggests that processes operating at low gradients are not effective transporting agents. As another example consider the gentle landscapes found in southern England or prairie regions of North America. Over most of their area, these landscapes appear to be relatively inactive in the present day. This situation can be expected to continue for 144 some time into the future. Therefore, lateral fluvial activity appears to be the dominant process in low-gradient landscapes. Further research is required to assess whether other processes, such as creep and surface wash (depending on the region being studied) can, in particular circumstances, be significant at lower gradients. 6.3.2.2 CHANGE IN SPACE STEP The effects of a change in space step to 50 m on the nonlinear transport model runs are now examined (figure 6.7). The differences in elevations at the divides for the Group A nonlinear transport model are about 5 m in both cases after 100 000 years while the elevation difference is only about 2 m in the valley bottom. The differences at the divides for the Group B nonlinear transport model are also about 5 m and the difference at the river grid cell is about 1 m. Such discrepancies are likely within the error range of this study. The nonlinear transport model runs are relatively insensitive to a change in the space step of 2 times. 6.4 V A L L E Y EVOLUTION BY PROCESSES IN T H E MAIN CHANNEL 6.4.1 BASIC M O D E L RUNS In the next set of model runs, hillslope transport was turned off in order to study the "valley rules" in isolation. The model was run initially with a space step of 100 m and a time step of 100 years, over a time period of 100 000 years. The input files for the phi values covered situations of: (i) aggradation over 100 000 years (ii) degradation over 100 000 years and (iii) alternating periods of aggradation and degradation (60 000 years aggradation; 20 000 years degradation; 20 000 years aggradation). The assumed magnitude of annual elevation change for all cases (aggradation and degradation) was 5 mm and the channel width was defined as 25 m (the latter value is based on channel widths for Mosquito Creek Tributary given in Roberts, 145 o o CD o o CN CO lath o CD ran tr o Q. CO o o 00 CN < CL ye; rs(( 3 O i_ rn o o CN CD o o o o o O O O -*—' CO CN CD O £Z E o o — 2 II CO X O b TJ O TJ CD C CO T— E o dx o o LO II CN co Q) «= o o 00 O CL o c o CO 'L. CO CL o o E o O 3s to CD i_ 3 (Ul) U0I}BA9|3 (Ol) U0!iBA3|3 146 1984). This defines a cross-sectional area of sediment accumulation of 0.125 m . In order to 2 conserve mass this value is effectively "spread" over the 100 m width of the cell containing the river (the space step in the model). Direct comparison with degradation rates would be misleading in this case, as the actual area of sediment accumulation is input into the model. The purpose of these model runs is to evaluate how the aggradation and degradation algorithms, which distribute actual channel changes across the grid according to particular "rules", affect the final erosion and deposition rates calculated in the model. In order to calculate river changes in the model which can be compared to real data, the 2-dimensional version of the program must be run. Results for continuous aggradation in the valley are plotted after 3 000, 10 000, 30 000 and 100 000 years of evolution (figure 6.8a). Although it is not likely that continuous aggradation (uninterrupted by erosion) would occur over such long periods, the performance of the valley rules can be evaluated by testing the model for this simple case. After 10 000 years the valley flat still consists only of the grid cell containing the river. The aggradation at this point is 12.5 m, which translates into an annual deposition rate across this cell of 1.25 mm/yr. After 100 000 years, the valley flat is defined by 3 grid cells and the annual rate of aggradation over the entire 100 000 year period is 0.6 mm/yr. The annual aggradation rates reported here are at the lower end of floodplain aggradation rates found by Bridge and Leeder (1979). After fluvial erosion has continued for 10 000 years, the river incision width is one grid cell and the incision depth across the cell is about 12.5 m, which represents an annual incision rate of about 1.25 mm/yr (figure 6.8b). This value is within the range of incision values (0.5-25 mm/yr) presented in section 5.4. By 100 000 years, the incision width has increased in order to maintain side slopes of less than 35° (maximum angle of stability). However, in order to achieve this situation, the river now cuts into the base of the main hillslope and sets off progressive (w) uoueA3|3 (LU) U0!JBA9|3 (W) U 0 U B A 3 B 149 upslope erosion. After 100 000 years the entire length of the right-hand main slope has been affected by "basal undercutting" due to channel incision. The incision rate over the 100 000 year time period is about 0.5 mm/yr. This value is lower than the rate over the first 10 000 years because, as time progresses, the incision is distributed across a greater area. In the third model run aggradation continues for 60 000 years, followed by 20 000 years of degradation, and a final 20 000 years of aggradation. In order to most effectively demonstrate patterns of infilling and incision, the results are now plotted at 0 years, 60 000 years, 80 000 years and 90 000 years (figure 6.8c). For the first 60 000 years the river grid cell aggrades about 46 m and the valley width extends across 3 grid cells. Thereafter, the river erodes at a width of one grid cell. The minimum width of incision possible in the model is the grid cell width. This, in many cases, will result in incision widths which exceed values expected in the field. The remaining outside two points of the old floodplain constitute, in effect a terrace. In all diagrams the lines between calculated grid points are plotted directly from point to point, thereby precluding direct observation of the terrace in figure 6.8c. Finally, after 80 000 years, the river starts to aggrade once again. The aggradation is spread across the width of one grid cell, and river sediments fill in the earlier incision. 6.4.2 CHANGE IN SPACE STEP The model runs in section 6.4.1 are repeated using a space step of 50 m. The results are compared to results for the earlier model runs, which used a space step of 100 m. At the final time step of 100 000 years the decrease in space step significantly increases the height of aggradation in the valley, increasing the elevation of the grid cell containing the river by 39 m (figure 6.9a). The same amount of sediment is effectively spread out over a smaller width (50 m instead of 100 m) and, hence, there is a greater increase in elevation. In the model runs for 150 151 continuous degradation over 100 000 years, the incision depth for the 50 m space step is slightly greater than for the model run with a 100 m space step (figure 6.9b). The depth of incision is greater by about 15 m in the former case. These results show that the valley submodel is more sensitive to changes in space steps than the hillslope submodel. There is a considerable computational cost in reducing the space step which favours the selection of a larger space step. Moreover, the increased area over which river processes operate when space steps exceed the river width can be thought of as adding in a factor to account for the lateral activity of rivers that is involved when both aggradation and degradation occur. At worst, it introduces a consistent negative bias into changes in elevation. 6.5 PROFILE EVOLUTION OF HILLSLOPES FEEDING INTO G U L L Y Model runs were conducted using the gully hillslope profile as the initial landscape in order to evaluate the deposition of material in the gully axis. The model was run for the case of linear diffusion (0.1 and 1 m /yr ) and the 2 nonlinear functions (Groups A and B) (figure 6.10). The rate of deposition in the gully channel represents what is usually referred to as the "recharge rate" of debris flow channels. The valley at the bottom of this hillslope contains a gully channel, and hence accumulation rates in the model can be compared with debris flow recharge rates for gullies in the Queen Charlotte Islands. Oden (1994) estimated recharge rate for 13 debris flow gullies in several basins in the Queen Charlotte Islands (table 6.2). She estimated recharge rates by calculating the volume of material stored in the gully and dividing this value by the age of the deposit (based on an estimate of the date of occurrence of the last debris flow). The data values of Oden (1994) were converted into cross-section recharge rates by dividing the volume of material stored in the gully by the length of the gully. These values were then adjusted to obtain recharge rates over a 100 year period. 153 154 Table 6.2 Gully recharge rates (Oden, 1994). Gully Identifier Gl G2 G3 G4 G5 G6 G9 G10 Gil G12 R4 R5 Slog Volume (m ) 388 576 339 385 432 744 231 309 468 51 228 121 50 3 Length (m) 218 341 232 235 247 512 410 180 280 40 352 296 118 Age (years) 14 14 14 14 75 57 14 14 14 14 6 10 8 Recharge rate(m /100yr) 12.7 12.1 10.4 11.7 2.33 2.54 4 12.3 11.9 7.29 10.8 4.09 5.3 2 Table 6.3 Gully recharge rates in present study. Transport rule Recharge rate: No weathering (m /100yr) dx = 100 m 15 Recharge rate No weathering (m /100yr) dx = 50 m 11 Recharge rate: Weathering (m /100yr) dx =100 m 6.5 144 107 7.2 3.3 0.35 2.9 2.5 0.45 2.6 2 Linear diffusion (0.1 m /yr) Linear diffusion (1 m /yr) Nonlinear function (Group A) Nonlinear function (Group B) 2 2 2 2 155 Recharge rates are determined in the present study by running the model for a 100 year time period. The accumulation depths obtained for each transport rule are multiplied by the grid cell width to obtain areal increases in sediment storage in the gully axis (table 6.3). The cross-sectional recharge rates based on the data of Oden range from 2.33-12.7 m /yr 2 with a mean value of 8.27 m /yr and a median value of 10.4 m /yr. The values found in the model results compare favourably with these values in all cases, except for the diffusivity value of 1 m /yr, which produces a recharge rate about an order of magnitude greater than those found 2 by Oden. When a larger space step is used the actual sediment accumulation rates are more consistent with Oden's results. However, it should be kept in mind that, although the net areal change is consistent with Oden's results, the heights of sediment accumulation differ in each case. In real gullies, the sediment accumulation is spread over a much smaller area with a greater depth. This is in contrast to the wide distribution areas which result from the model space step, and the accompanying decreases in elevation change. However, it is encouraging that values of the change in cross-sectional area are preserved, even if the actual distribution in the landscape is not. 6.6 W E A T H E R I N G The gully profile for Mosquito Creek Tributary is used in this suite of model runs. The space step is 100 m in all cases. The thickness of the sediment layer at the start of the model runs is set to 0 m. The weathering equation of Heimsath et al. (1997) (equation 3.12d, p. 68) is selected to simulate the weathering of bedrock as it is derived from field evidence. Bedrock conversion into sediment and transport of material interact to create a loose sediment layer. The heights of the surface and bedrock are evaluated at each time step. 156 The bedrock weathering equation is incorporated into the "DERIVS" subroutine. During each call to "DERIVS" the depth of the sediment profile, which is equal to the difference between the surface and bedrock heights, is calculated. Appropriate restrictions are placed on transport rates such that, when the depth of the sediment profile approaches 0, the transport rate is reduced to 0 m /m yr. In this manner, the rate of bedrock weathering limits the total amount of 3 hillslope transport. 6.6.1 BEDROCK LOWERING R A T E WITH NO TRANSPORT The version of the 1-dimensional model which includes weathering was initially run for a zero transport rate in order to evaluate changes in sediment thickness which occur when there is no hillslope transport. The sediment thickness is set to 0 m at the start of the model run. Weathering initially proceeds rapidly, achieving a depth of about 0.19 m in the first 3 000 years. Thereafter, sediment thickness increases at a decreasing rate, achieving values of about 0.44 m after 10 000 years and 1.3 m after 100 000 years. The increase is logarithmic because of the negative exponential expression. 6.6.2 PROFILE EVOLUTION Both transport and weathering are incorporated into the following model runs. An interesting set of relations is introduced into the model when weathering, in addition to transport, is considered. Bedrock lowering is expected to proceed at a very slow rate in areas subject to significant deposition. The increase in the thickness of the sediment layer which occurs is due primarily to the deposition of material, and not bedrock weathering, in this case. 157 Erosion acts to reduce the thickness of the sediment layer. However, as the thickness is reduced, the bedrock weathering rate increases. In these situations, the complex interplay between these processes determines the thickness of the sediment layer at a point in time. Linear Transportfor Creep Results of the model run for creep transport with weathering indicate that creep transport occurs at a slower rate than bedrock weathering. During the entire 100 000 years of the model run, the height of the surface remains well above that of the bedrock surface. Under the rules in the model, the evolution of the land surface is not materially affected by weathering in regions dominated by creep processes. Linear Transportfor Landslides The introduction of weathering into the model significantly reduces transport activity for the linear diffusivity value of 0.1 m /yr. (figure 6.11). Erosion rates at the hillslope crests are about one-half of their original values (non-weathering limited) after 10 000 and 100 000 years. Deposition on the lower slopes is reduced, although accumulation rates in the gully axis remain approximately the same. The overall decreases in delivery are a result of the reduction in delivery rates. Hillslope transport rates are considerably lower for the case with weathering than for the original results when the linear diffusivity is increased to 1 m /yr (figure 6.12). The rate of bedrock conversion into sediment provides a constraint which considerably lowers transport activity. 158 o o o o o o o o o o o o o o o o o o o c n o o i ^ - C D L n - ^ c o c M T (Ui) U0JJEA9I3 o o o o o o o o o o o o o o o o o o o D c o N t o m f r o N T - C (Ol) U0HBA9|3 159 o o o o o o o o o o O O O O O O O (Ul) U0|}BA3|3 O O o o o o o o o o o o O O O O O O O (LU) U0|iBA9|3 O O 160 Nonlinear Transportfor Landslides The Group A and B nonlinear transport equations are applied to model runs which incorporate bedrock weathering (figures 6.13 and 6.14). The morphological differences between this set of results and the original model runs are minimal, particularly for the Group A results for which the differences are not readily observable on the plots. Transport rates are operating relatively slowly and are not significantly constrained by the weathering rate. The results for non-linear transport for landsliding do not show a significant change when weathering is introduced. This fact suggests that the gradients across most of the landscape are below the critical threshold for landsliding. Discussion The results of this analysis suggest that for the high linear diffusivity of 1 m /yr, the 2 introduction of weathering into the model produces results which appear to be more realistic. The weathering sufficiently constrains transport so that unusually high transport rates no longer occur. This line of reasoning seems to suggest that the higher diffusivities implemented by other researchers may be acceptable as long as weathering is accounted for in the model. However, at this point it becomes critical to examine the accompanying patterns of sediment layer thicknesses. These data provide further information on which to evaluate the performance of the different transport functions. 6.6.3 THICKNESS OF SEDIMENT L A Y E R Linear Transportfor Creep The resulting thicknesses of the sediment layer across the gully profile for the creep diffusivity of 0.0002 m /yr are presented in figure 6.15a. Sediment depth is equal to about 0.45 2 163 900-] 800~ 700• § 600+ .1 500- 1.6- Gully cross-section profile "Note: The legend shown in (b) applies to all graphs. The x-axis is the same for all graphs. 1.41.21.0Q. CD Q 0.80.60.40.20- Figure 6.15 (a) Thickness of sediment layer for gully cross-section profile (diffusivity = 0.0002 rrrVyr). 70 Initial profile 60 3 000 years 50 Q. CD Q 10 000 years 40 30 000 years 30 100 000 years 20 i I I •. i/ \ ! 10 / i 0 (b) Thickness of sediment layer for gully cross-section profile (diffusivity = 0.1 m7yr) 70 60 50 40 CL CD Q 30 20 10 0 \ V . X \ :.' ^ \ \\ I i -f* • 1500 1800 2100 2400 2700 Distance (m) (c) Thickness of sediment layer for gully cross-section profile (diffusivity = 1 mVyr). 300 600 900 1200 3000 164 900 800 700 600 500+ 400 300200100 0 Gully cross-section profile *Note: The legend shown at bottom of page applies to both graphs. The x-axis is the same for all graphs. 18 16 14 12 10 8 6 4 2 0 Initial profile 3 000 years 10 000 years 30 000 years II 1 \ 100 000 years . i • ' s' . rUfrf- - -- X— liE (d) Thickness of sediment layer for gully cross-section profile (Group A transport relation). 18 16 j 14 12 10 8 6 4 « A ". / \ i i \ i \ \ 2 i--=-^'----.V-. 0 ••^'•-'i 300 600 900 .xavty.^. 1200 1500 1800 2100 2400 2700 3000 Distance (m) (e) Thickness of sediment layer for gully cross-section profile (Group B transport relation). 165 m after 10 000 years and ranges from 1.2 to 1.4 m after 100 000 years. These values are very similar to sediment thicknesses developed in the earlier model run in which there was no transport activity. The soil depths computed by the model are of the same order of magnitude as those on the Queen Charlotte Islands, which are typically about 0.5 to 1.0 m on steep slopes and somewhat deeper in hollows (Church, 1997, pers. comm.). Linear Transportfor Landslides The results of model runs using a linear diffusivity of 0.1 m /yr for landsliding are given 2 in figure 6.15b. In this case, there is no development of a sediment layer at the hillslope crest and only minimal development along the slopes. The negligible thicknesses of the sediment layer on the slopes suggests that: (i) material is eroded as quickly as it is weathered and (ii) there is minimal deposition of material along the slopes. Trees can be found throughout the entire area of most basins, even on the upper slopes, in the Queen Charlotte Islands. However, soil development is often restricted to a very shallow organic layer over rock with trees often anchored in joints. Sometimes, there may be a till soil, which is not produced by weathering. Therefore, this result is not necessarily unreasonable given observations on the Queen Charlotte Islands. The depth of the sediment layer is about 8 m in the gully axis after 10 000 years and increases to about 60 m after 100 000 years. The bedrock surface is expected to be relatively close to the surface along the gully axis. In order for this rate of accumulation to be considered realistic, high debris flow rates are necessary to counteract these rates of accumulation in the gully bottom. Assuming that about 0.5 m of erosion is associated with a debris flow event, 16 debris flows must occur every 10 000 years (equal on average to about 1 every 600 years) to counteract this accumulation rate. This seems to be a reasonable expectation. In this way the 166 characteristic thin veneer of sediment in these locations could be explained while still accepting this diffusivity value. When the model is run using a higher landsliding diffusivity of 1 m /yr there is no 2 sediment layer development on the lower, middle and upper slopes after 100 000 years (6.17c). This result is not supported by observations of typical hillslopes in the Queen Charlotte Islands, which suggests that this diffusivity value is too high. Sediment accumulation in the central axis is about 8 m after 10 000 years and 70 m after 100 000 years. These accumulation rates are of the same order of magnitude as for the case of the lower diffusivity. Therefore, it is possible for these rates to be offset by debris flow activity. Nonlinear Transportfor Landslides Sediment thicknesses on the upper slopes are about 0.4 m after 10 000 years for Group A and Group B transport relations (figure 6.15e). After 100 000 years the values are the same order of magnitude, although the variability is somewhat increased. There is somewhat less development on the middle slopes than the upper slopes after 10 000 years, which indicates that this area is a transport zone with minimal deposition. These values seem reasonable in light of the relatively shallow soil depths found on middle and upper slopes in the Queen Charlotte Islands. Sediment thicknesses on the lower slopes range from negligible values to about 1.6 m after 10 000 years and 3-10 m after 100 000 years. In the gully axis, sediment thicknesses for both groups of basins range from several meters to about 14 m after 10 000 and 100 000 years respectively. These low sediment accumulation rates do not require exceptionally high debris flow activity in order to evacuate this material. The low sediment thicknesses in the valley flat reflect the overall reduced activity rates for the nonlinear transport functions. The values of 167 sediment thickness produced by the model for Groups A and B basins seem reasonable in comparison to observations in the Queen Charlotte Islands. 6.6.4 G U L L Y R E C H A R G E RATES FOR MODEL RUNS WITH WEATHERING Model results for two linear diffusivities and two nonlinear diffusivities are examined to assess the impact of weathering on gully recharge rates (table 6.3). Recharge rates for the linear diffusivity of 0.1 m /yr and the two nonlinear models were lowered only slightly by this 2 modification. Therefore, the new results do not affect the earlier interpretation of recharge rates for these functions in comparison to Oden's data. In the case of the linear diffusivity of 1 m /yr, 2 the gully recharge rate was reduced from an unrealistically high value of 144 m /100yr to a more 2 reasonable value of about 7 m /100yr. 2 168 CHAPTER 7: TESTING T H E CHANNEL SUBMODEL The channel submodel is tested in isolation from the hillslope submodel in order to evaluate its performance without the added effects of hillslope coupling. The channel submodel is tested for channel systems covering a broad range of scales in order to examine the suitability of its application to rivers of varying magnitudes. Fluvial transport rates are estimated for the Vedder and Fraser Rivers in British Columbia using the fluvial submodel. The latter river is significantly larger than any river used in the calibration of the bed load transport equation. Successful application at this scale would suggest the robustness of the revised bed load relation. Values of relevant channel variables and bed load transport are compared to the results of prior studies of these rivers (McLean, 1990; Martin, 1991; Martin and Church, 1995). In the final section of this chapter, channel changes resulting from fluvial and debris flow activity are calculated for Mosquito Creek Tributary basin. In addition, estimates of annual suspended load are contrasted with estimates of bed load for each study basin. 7.1 VEDDER RIVER The Vedder River is located 160 km east of Vancouver and drains an area of 1230 km . 2 The "Vedder River", the name given to the most distal reach of the Chilliwack River, flows across an alluvial fan. Significant aggradation occurs along the Vedder River due to reduced gradients on the fan. The Vedder and Chilliwack Rivers are particularly active in the present day, possibly because of increased sediment loads resulting from the last Neoglaciation in this area (1500 to 1900 A.D.) (Church, 1997, pers. comm.). The 2-year rainfall generated flood for the Vedder River has a discharge of about 450 m /s (Church, 1997, pers. comm.). 3 This value is input directly into the model for later calculations. Elevations along the Vedder River are obtained from prior studies and are input 169 directly into the program (Martin, 1991; Church, 1997, pers. comm.). The 11 cross-sections used in this analysis are the bounding cross-sections of the 10 condensed study reaches defined in Martin (1991). The average length of the study reaches is about 800 m and together they cover a total distance of just over 8 km. Model results are compared to transport rates and volume changes determined from repeated cross-section surveys (Martin, 1991; Martin and Church, 1995). Channel changes are inserted into a sediment budget framework in order to obtain transport rates. 7.1.1 C H A N N E L V A R I A B L E S Channel pattern type must be defined in order to select the appropriate equations for the calculation of width and depth. The algorithm for channel pattern requires a knowledge of whether sediment (the D o in this case) at the channel points consists of sand or gravel. Grain 5 size measurements in Martin (1991) show that at all locations along the Vedder River the D o is 5 > 2 mm. This information is fed directly into the program. The channel is classified as wandering along the entire river, except for the uppermost cross-section which is defined as braided. These results are reasonable in comparison with observed channel patterns along the Vedder River, which show increased braiding in the upstream reaches. However, the uppermost several reaches are braided while the model predicts braiding only in the uppermost reach. The hydraulic geometry equations of Bray (1972, 1979), which form the most appropriate set for British Columbia rivers, are used in this model run. Model results are compared to width and depth estimates made by Church (1997, pers. comm.) on the basis of the cross-section survey data. While in reality, width and depth vary along the length of the Vedder River, the fixed discharge used in this model run yields only one value of each parameter for each type of pattern. Moreover, the width and depth estimates in the model do not account for changes in controlling 170 variables along the channel or other apparently stochastic perturbations. The effects of "smoothed" widths and depths eliminate locally-determined erosional and depositional sites, thereby providing only a broad, generalized pattern of channel change. From the perspective of landscape evolution, such generalizations may be appropriate. The calculated data show a relatively insignificant change in width and a somewhat more significant change in depth at the location where the transition from a wandering to braided pattern occurs (figure 7.1a,b). Calculated widths provide reasonable approximations to the observed widths in the wandering reaches. The model predicts widths of about 116 m for the wandering sections of the river while observed widths at channel cross-sections range from about 60 m to 180 m. Braided reaches are generally expected to be wider and shallower than wandering reaches. The model estimates approximately the same width for the braided reach as for wandering reaches. Width estimates in these braided reaches are lower than the observed values of Church (1997, pers. comm.). The bank-to-bank width of a braided river is considerably greater than the channel width, the latter being important from a transporting perspective. Church (1997, pers. comm.) attempted to account for this factor by eliminating the higher bar surfaces from his estimates of width. The approach adopted in the model to calculate channel width for braided reaches (section 4.3.2.4) also attempts to include only the active channel. However, a discrepancy still remains between model results and observed values of width. It is possible that the truly active portion of a braided channel may be smaller than the values obtained by eliminating the high bar tops. In addition, the approach adopted to estimate the active width in the model is based on two individual studies. There remains a need for additional study as the active widths of braided reaches have not been well documented. Depths of 1.9 m and 1.2 m are calculated for wandering and braided reaches respectively. Observed depths range from less than 1 m to greater than 3 m. Calculated depths deviate 171 300 - • Calculated width a Observed width 250 200 •g u 150 o 100 • Dp H • • • e • • • i • 50 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 Distance (m) Figure 7.1 (a) Calculated vs. observed channel widths for the Vedder River. 3.5 U II • Calculated depth a Observed depth 3 2.5 1 Q. CD Q 2 1.5 • H • 7000 8000 B 1 0.5 0 0 1000 2000 3000 4000 5000 6000 Distance (m) (b) Calculated vs. observed channel depths for the Vedder River. 9000 172 • 0.04 0.035 • Calculated grain size • Observed grain size U • 0.03 f 0025 m I (3 0 0 2 H 0.015 m 0.01 o • • » • • • • N m 0.005 0 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 Distance (m) (c) Calculated vs. observed grain sizes for the Vedder River. 35 •Calculated elevation m Observed elevation 30 25 c o 20 a> LU 15 10 5i 0 -I 0 • 1 1 1 1 1 1 1 1 1 1000 2000 3000 4000 5000 6000 7000 8000 9000 Distance (m) (d) Calculated vs. Observed channel elevations for the Vedder River. 173 significantly from individual observed values (refer to figure 7.1b). Nonetheless, they represent a reasonable average of observed depths along the Vedder River. Although the pattern of measured grain sizes shows more overall variability than model results, grain sizes calculated in the model show a pattern of increasing grain size with distance upstream (figure 7.1c). The values are somewhat lower than observed values, although they provide a reasonable lower-bound estimate. Lower-bound estimates of grain size imply a reduced threshold for transport and higher estimates of transport rates. An exponential model is fit to the channel elevation data. It provides a good approximation of elevations, and hence gradients, along the river (figure 7.Id). 7.1.2 B E D L O A D T R A N S P O R T Unit mass transport rates calculated in the model range from 0.05 kg/m s to 0.62 kg/m s. The results fall within typical ranges of transport rates for rivers in flood (Gomez and Church, 1988). Transport decreases in a downstream direction, which indicates an aggradational trend in the data. Unit mass transport rates are converted into annual volumetric transport rates (using a sediment bulk density of 1 800 kg/m ) (table 7.1a). Discharge data were inspected to determine "typical" flood events on the Vedder River. Floods typically last for several days, with a 4-day duration being an "average" value. This value is used in model runs. The model estimates that about 11 400 m /yr are transported into the uppermost reach of the Vedder River. Transport decreases to a value of 1 000 m /yr in the lowermost reach. The 2-year discharge value used in the present calculations (450 m /s) falls within the range of floods occurring over the period from 1982-1987 in the study of Martin and Church (1997), which allows for comparison of results. Transport rates reported in Martin and Church (1995) are several times higher than transport rates calculated by the model (table 7.1a). The 174 Table 7.1 (a) Volumetric transport rates for the Vedder River (m /yr). Data are from Martin and Church (1995). 3 Year: 1981-82 (142 m /s) 3 Reach 1 2 3 4 5 6 7 8 9 10 11 1982-83 1983-84 (313 m /s) (524 m /s) 3 3 1985-87 (344m /s; 349m /s) 1 000 0 1 000 4 000 8 000 10 000 15 000 14 000 20 000 28 000 29 000 3 3 8 000 8 000 8 000 1 000 5 000 4 000 5 000 15 000 25 000 20 000 20 000 0 12 000 18 000 40 000 40 000 48 000 42 000 72 000 73 000 62 000 67 000 0 3 000 5 000 12 000 11 000 12 000 20 000 23 000 25 000 34 000 33 000 Model Model Results Results (450 m /s) (450 m /s) 1 000 654 1 220 388 1 420 910 1 720 0 1 990 2 320 2 350 977 2 870 0 3 530 0 4 180 0 5 140 2 630 11 400 9 920 8 3 b 3 (b) Volume changes for the Vedder River (m /yr). Data are from Martin and Church (1995). 3 Year: 1981-82 (142 m /s) 3 Reach 1 2 3 4 5 6 7 8 9 10 1982-83 1983-84 (313 m /s) (524 m /s) 3 3 1985-87 (344m /s; 349m /s) -2 000 +1 000 +7 000 +8 000 +4 000 +9 000 -2 000 +5 000 +12 000 +3 000 3 3 0 -2 000 -10 000 +3 000 -2 000 +3 000 +21 000 +11 000 -6 000 0 +21 000 +8 000 +39 000 +2 000 +12 000 -8 000 +45 000 +3 000 -20 000 +6 000 +4 000 +7 000 +9 000 -3 000 +2 000 +13 000 +7 000 +3 000 +14 000 -2 000 "Results calculated using calculated values for input variables. Results calculated using observed values for input variables. b Model Model Results Results (450 m /s) (450 m /s) +214 -267 +205 . +522 +298 -910 +270 +2 320 +356 -1 340 +518 -977 +665 0 +652 0 +959 +2 630 +6 300 +7 290 8 3 b 3 175 model predicts aggradation along the length of the river, as more material is transported into each reach than is evacuated. The results of Martin (1991) also show an aggradational trend along the river. Volumetric changes are calculated by incorporating transport rates into a sediment budget framework. Model results estimate volume changes ranging from 200 to 6 300 m /yr (table 3 7.1b). These values are lower than those reported by Martin and Church (1995) (table 7.1b) but they are well within the correct order of magnitude. There are several possible explanations for the somewhat low transport rates calculated by the model: (i) Errors associated with the input variables of width, depth and grain size may affect transport calculations in the model. Local variability in these quantities is not captured in the model calculations and this may affect transport patterns. At a constant discharge, increases in width are expected to be associated with decreases in depth. When depth is decreased the transport rate is expected to increase (according to the threshold stream power relation and the final transport relation; see section 4.3.1); an increase in depth would decrease transport rates. Likewise, a change in the grain size results in the opposite change in the transport relation. (ii) The entire calcuation is keyed to a particular return flow. Lesser competent flows may affect the results as they may return with a greater frequency. In addition, exceptionally large flows may lead to particularly significant channel changes, as was shown to be the case in the study of the Vedder River by Martin and Church (1997). (iii) There is expected to be some error associated with the use of the Bagnold-type equation derived in this study. Despite the high value of R found in the calibration of this equation, there was noticeable scatter about the best-fit line, which was documented in Appendix II. There is a large, nonsystematic variation in the observations. 176 7.1.3 A F U R T H E R T E S T The concern that inexact estimates of channel parameters may contribute to the low values of transport calculated by the model was addressed by re-calculating transport rates, this time using observed values of channel variables given in the summary data of Church (1997, pers. comm.). The unit mass transport rate for the most upstream cross-section increases to a value of 0.99 m /yr from its previous value of 0.62 m /yr . Values in the lower reaches are significantly lower than calculated in the initial model run. Furthermore, several cross-sections have a 0 transport rate. Volumetric transport rates range from 0 to 9 900 m /yr, while volume changes range 3 from -1 340 m /yr through to a maximum of 7 290 m /yr (table 7.1a,b). The overall transport 3 3 pattern along the river is more variable than for the initial model run. The introduction of observed values for channel variables does not result in increased transport rates. In addition, it is problematic that transport rates of 0 were calculated at several cross-sections. The somewhat low transport estimates do fall within the 3-times standard error for the Bagnold-type equation used in the model (see Appendix II). Furthermore, the threshold stream power calculation of Bagnold overestimates the threshold for transport as indicated by the negative values of excess stream power found in the analysis in Appendix II. If threshold stream power were reduced in the model, then there would be an accompanying increase in transport estimates. 7.1.4 S U S P E N D E D L O A D The suspended load transport rate is calculated using equation 4.26 (p. 103). A result of 0.21 Mg/km day is obtained. Church et al. (1989) found that the daily suspended load (based on 2 a 10 year record) measured for the Chilliwack basin at Vedder Crossing was 0.29 Mg/km day. 177 These suspended load transport rates translate into annual mass evacuation rates of 9.4 x 10 7 Q kg/yr and 1.3 x 10 kg/yr respectively. Martin (1991) found that typical annual transport rates of bed material along the Vedder River were about 5.4 x 10 kg/yr. Based on these data, bed load 7 comprises about one-third of the total load in this basin. This value is significantly higher than estimates of about 3-10% usually given for bed load rivers (Schumm, 1977). The suspended load rates calculated above are compared to the bed load transport rate calculated at Vedder Crossing in the model (section 7.1.2). Bed load makes up about 14 to 18% of the total load in these cases. This value remains somewhat high, but it is very close to the observed record. 7.2 FRASER RIVER The Fraser River is about 1 360 km in length and drains an area of about 232 000 km . 2 McLean (1990) studied sediment transport along a reach of the lower Fraser River extending from about 85 km below the Fraser Canyon (near Hope) to Mission. The 5 cross-sections used in this study correspond to the boundaries of 4 major study reaches in McLean (1990); Sumas, Chilliwack, Rosedale and Cheam reaches. The total length of these reaches is about 60 km. Values for the 2-year flood and gradients used in the channel submodel are based on data from McLean (1990) (table 7.2a). In addition, width, depth and grain size values are extracted from McLean (1990) for comparison with model results. Transport measurements reported in McLean (1990) are derived from several sources. The values given by McLean which are used in this analysis are: (i) bed load measurements made at the Agassiz gauging station and (ii) sediment transport rates estimated using a sediment budget approach in combination with measurements of morphological change along the river. In the latter case, measurements are based on bathymetric survey data. Transport estimates are for grain sizes > 2 mm whereas the bed load equation used in the model is calibrated for bed load > 178 Table 7.2 (a) Discharge and gradient values at Fraser River cross-sections (data extracted from McLean, 1990). Cross-section 1 2 3 4 5 Discharge (2-year flood) (m /s) 10 500 10 300 9 350 8 550 8 530 3 Gradient 0.000085 0.000132 0.000325 0.000495 0.000535 (b) Daily mass and volumetric transport rates at Fraser River cross-sections using calculated variables. Crosssection Unit mass transport rate (kg/m s) 1 2 3 4 5 0.0092 0.0195 0.0465 0.0691 0.0745 Mass transport rate (tonnes/day) 440 1 100 2 400 3 400 3 700 Volumetric transport rate (m^/day) 240 590 1 300 1 900 2 000 Volumetric transport rate (m /year) 18 300 44 400 100 000 142 000 153 000 (c) Daily mass and volumetric transport rates at Fraser River cross-sections using observed variables. Crosssection 1 2 3 4 5 Unit mass transport rate (kg/m s) 0.028 0.0022 0.038 0.6004 0.0875 Mass transport rate (tonnes/day) 2 100 210 2 300 13 000 3 800 Volumetric transport rate (m7day) 1 180 116 1 280 7 200 2 100 Volumetric transport rate (m /yr) 88 600 8 670 95 600 540 000 157 000 3 179 1 mm. The bias is expected to be minimal since a small fraction of bed material falls in 1-2 mm category. 7.2.1 CHANNEL VARIABLES Initial grain size estimates used in the determination of channel pattern are based on data from McLean (1990). The model predicts that the lowermost cross-section is a sand meandering channel and upstream cross-sections are gravel wandering. These data correspond to the classification given by McLean (1990). Once again, the hydraulic geometry relations of Bray (1972, 1979) are used in this test of the river submodel. Calculated widths are significantly lower than observed values at several cross-sections and at other locations the model overestimates width (figure 7.2a). The model provides "smoothed" values of width which do not account for the natural variability found along real channels. However, the results do provide a reasonable average of widths along the study reach. Estimates of depth are consistently too high (figure 7.2b). Calculated grain sizes compare well to measured values, except at the upstream limit of the study reach (figure 7.2c). 7.2.2 BED LOAD TRANSPORT Mass transport rates per unit width calculated in the model range from 0.0092 kg/m s to 0.075 kg/m s (table 7.2b). These values fall within typical ranges of unit transport rates compiled by Gomez and Church (1988). Total daily mass transport rates for the 2-year flood range from 440 tonnes/day through to about 3 700 tonnes/day along the study reaches (table 7.2b). These values compare favourably to bed load mass transport rates measured at the Agassiz gauging station of several thousand tonnes/day for the 2-year discharge (McLean, 1990). The gauging station is located close to cross-section 4 as defined in the present study. McLean (1990) states that discharges of about 5 000 m /s are required before significant 3 180 1200 u • Calculated width H Observed width 1000 11 ~ I 5 800 600 o • • w • > • 400 11 200 0 10 20 30 40 50 60 70 Distance (km) 7.2 (a) Calculated vs. observed channel widths for the Fraser River. 12 1 8 1 • Calculated depth • Observed depth 10 o rip * n 6 4 • 1 iiiji I1 2 0 10 20 30 40 50 60 70 Distance (km) (b) Calculated vs. observed channel depths for the Fraser River. 0.03 • Calculated grain size i i Observed grain size 0.025 • r Lr 0.02 1 (IH— 0.015 • 0.01 0.005 o 0 • il 10 20 30 40 Distance (km) (c) Calculated vs. observed grain sizes for the Fraser River. 50 60 70 181 gravel transport occurs along this study reach of the Fraser River. Discharge records for relevant stations were assessed for years experiencing a maximum annual flood about equivalent to the 2year flood. Discharges exceed the threshold for gravel transport (5 000 m /s) over extended time 3 periods during the spring and summer snowmelt season, with an average flood period of about 75 days. Therefore, the daily transport rate is multiplied by this value in order to obtain annual transport rates (table 7.2b). Values range from about 150 000 rnVyr in the upper study reaches to 18 000 m /yr in the lower reaches. These values compare favourably with the gravel transport rates estimated by McLean (1990). 7.2.3 A FURTHER TEST Transport rates were re-calculated using observed values of channel variables. Depths and active widths for the 2-year flood and grain sizes along the study reach were extracted from various tables and diagrams presented in McLean (1990). Mass transport rates per unit width calculated by the model range from 0.0022 kg/m s to 0.60 kg/m s (table 7.2c). These rates exceed the results of the initial model run for some crosssections, while at other locations the re-calculated values are lower. Daily mass transport rates and annual volumetric transport rates are significantly greater than the results of the initial model run in several cases (table 7.2c). The incorporation of observed values of channel variables in the model run reduces the compatibility of model results with the transport estimates of McLean (1990). It is possible that difficulties associated with the estimation of threshold stream power for transport may be at least partially responsible for the poor transport estimates. In natural rivers, there is significant natural variability in resistance to transport due to properties of bed material and its structure. Such variability is not captured in the equation used to calculate threshold stream power. 182 7.2.4 SUSPENDED LOAD A suspended load transport rate of 0.44 Mg/km day is calculated using equation 4.27 (p. 2 103). Church et al. (1989) report suspended transport rates for the Fraser River near Agassiz of 0.2 Mg/km day. The associated annual mass transport rates are 3.5 x 10 kg/year and 1.59 x 2 10 kg/year respectively. 10 10 McLean (1990) estimated a gravel transport rate of about 125 000 m /year near Agassiz, which converts to an annual mass transport rate of 2.25 x 10 kg/year. The 3 8 gravel load comprises about 0.64% to 1.4% of the total load in these cases. Bed load finer than 2 mm is not considered here, which means that this value is a lower-bound estimate. A volumetric transport rate near Agassiz of about 2.75 x 10 kg/yr was calculated by the 8 model. Bed load accounts for about 0.8% to 1.7% of the total load in this case. 7.3 MOSQUITO C R E E K TRIBUTARY BASIN The complete channel submodel, which includes the possibility for transport by debris flow mechanisms, is tested for the Mosquito Creek Tributary basin. This basin, which was used in the hillslope profile model runs, has a drainage area of 5.5 km . The channel network in this 2 basin consists of both fluvial and debris flow channels. The hillslope transport function is set to 0 in this model run in order to observe changes in the channel system in isolation, without the effects of hillslope coupling. A data file of grid cell elevations for the study basin was created based on 1:20 000 TRIM map data. The outline of the channel network was overlaid on the elevation grid and channel points were defined where the channel intersects the grid network. The drainage network consists of a main channel running along the length of the basin, into which 23 major gullies feed. A further 6 gullies feed into several of these major gullies. Channel points are assigned the elevation of the nearest grid point. Model variables are defined in table 7.3. In this model run, variables are calculated Table 7.3 Variables in the Mosquito Creek Tributary model run.* Variable nx (# grid points in x-direction) ny (# grid points in y-direction) dx, dy (space step) Q=aA (discharge/area relation) Suspended load b Value 40 53 100 m a=6.0 (from Church, 1997) 0.01 Mg/km day (from Church et al., 1989) 2 * Coefficients and/or exponents in equations which are common to all model runs are not given here. Refer back to Chapter 4 for details on equations. 184 according to the complete set of rules given in chapter 4. 7.3.1 CHANNEL VARIABLES The model defines about 4 km of the main channel length as a fluvial channel. The upper 730 m of the main channel are greater than 8° and thus defined as a debris flow channel. Gradients in most gullies exceed 8° along their entire length. However, gradients for the lowermost channel points of 6 gullies are lower than 8°. The 2-year flood discharge is less than 1 m /s in upper gullies and reaches a maximum of 3 20 m /s at the basin outlet. Discharges along the main channel increase as contributing area to the channel increases (figure 7.3a). Discharge increases in a step-like manner at channel points where significant gullies join the main channel. Gradients range from a low of 0.01 at the basin outlet to a maximum of 0.56 at the extreme upper region of debris flow channels. Gradients along the main channel show a steadily increasing trend and reach a maximum of about 0.22. Calculated widths are greatest in the downstream portions of the main channel and decrease through high-gradient fluvial reaches and debris flow reaches (figure 7.3b). In an initial model run, which is not reported in this thesis, there was an abrupt change in widths at the boundary between braided reaches and high-gradient (>1°) reaches immediately upstream. Abrupt changes in variables are found to adversely affect the sediment budget by causing unrealistically high erosion or deposition rates. Therefore, it was decided to smooth the transition between the braided and high-gradient channel zone by defining a zone of intermediate width. The width is defined as the average of the widths calculated for a braided reach and a high-gradient reach. This intermediate zone was found to most effectively eliminate abrupt changes in transport by defining it as the region with gradients between 1° and 2°. Roberts (1984) reported widths in the lower braided reaches of the main channel of around 20-25 m, 185 20 18 16 < 14 £ CD E> co 5 Q 12 10 8 6 4 2 0 • 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Distance (m) Figure 7.3 (a) Discharge along main channel of Mosquito Creek Tributary. 25 20 IT 15 * 'WW 9 i • < TJ § 10 • 5 ••• • • • • * 0 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 4500 5000 Distance (m) (b) Width along main channel of Mosquito Creek Tributary. 1.2 • • 1 0.8 •4-* Q_ CD '••>•> • < • •••• •••••>• * • « • • < • • - 0.6 Q 0.4 0.2 0 0 500 1000 1500 2000 2500 3000 Distance (m) (c) Depth along main channel of Mosquito Creek Tributary. 3500 4000 186 • Gradient • Adjusted Gradient » •• i >•••• 3 minn ( 0 500 1000 1500 2000 • B 2500 m• a H 3000 3500 4000 4500 5000 4500 5000 Distance (m) (d) Gradient along main channel of Mosquito Creek Tributary. 4 — - • A ^ 0 A • • 500 1000 1500 2000 2500 3000 Distance (m) (e) Grain size along main channel of Mosquito Creek Tributary. 3500 4000 187 which suggests that model approximations are reasonable. Calculated depths of flow in the main channel are of order 10° m (figure 7.3c). No field data are available for comparison with these values. An exponential line was fitted to the elevations assigned to the river points. The resulting gradients which were calculated along the river are plotted in figure 7.3d. There is a steady increase in slope in an upstream direction. The hydraulically effective slopes, used for the calculations of bed load transport in the high-gradient, fluvial reaches (1° < gradient < 8°), are also shown on this graph. Grain sizes calculated by the model range from 45 mm at the basin outlet through to 400 mm in the upper fluvial reaches in the main channel (figure 7.3e). These values are reasonable when compared to the field data reported in the study of Roberts (1984). Grain sizes of 17-32 mm were found in the lower 2 km of the main channel. 7.3.2 BED LOAD TRANSPORT Unit mass transport rates are 0.7 kg/m s at the basin outlet and reach a maximum of about 3.4 kg/m s in the upper fluvial reaches of the main channel (figure 7.4a). Values are in the upper limit of rates reported in Gomez and Church (1988). The values translate into an annual volumetric transport rate of 500 m /yr at the basin outlet (invoking a 2-year flood of duration 24 hours) (figure 7.4b). Transport increases in an upstream direction starting at the basin outlet and moving upstream to the start of the transition zone (gradient>l°) at about 800 m. Thereafter, transport decreases to zero transport for several km. At about 2 700 m transport rates increase up until the end of the fluvial reach at about 4 100 m. The transport rates in the fluvial reach are somewhat higher than would normally be expected for a basin of this size. In fact, the Mosquito Creek Tributary basin is highly disturbed 188 500 2500 3000 3500 4000 4500 5000 Distance (m) Figure 7.4 (a) Unit mass transport rate for bed load along main channel of Mosquito Creek Tributary. 1600 CO 1500 2000 rl 1400 O 1000 1200 1000 ** o 800 600 -* 400 200 0 2500 3000 3500 4000 4500 Distance (m) (b) Volumetric transport rate for bed load along main channel of Mosquito Creek Tributary. 10 9 8 7 t l 6 5 SS -o 4 O TJ 3 2 1 a 0 0 500 1000 1500 500 1000 1500 2000 5000 II > 2000 2500 3000 3500 4000 4500 5000 Distance (m) (c) Volumetric transport rate for suspended load along main channel of Mosquito Creek Tributary. 189 due to extensive logging in the basin. However, in this model run the effects of increased erosion on hillslopes are not considered and therefore do not account for the high fluvial transport rates. Volumetric transport rates of 15 m /km yr were estimated for Carnation Creek, 3 2 British Columbia, a coastal basin with an drainage area of about 10 km (Tassone, 1987). When 2 the bed load transport rate at the basin outlet in the present study is translated into a comparable 3 2 form, a value of 90 m /km yr is obtained. This value exceeds the Carnation Creek datum by a significant amount, suggesting that the model is overestimating transport rates. 7.3.3 SUSPENDED LOAD The volumetric transport rate for suspended load at the basin outlet of Mosquito Creek -a Tributary is about 9.3 m /yr. In contrast, the volumetric transport rate for bed load at the same location is about 500 m /yr (figure 7.4c). Bed load comprises about 98% of the total load. A comparison of values along the fluvial portion of the main channel shows that in most cases bed load is about 98-99%) of the total load. However, the low values of bed load as percentage of total load given in Schumm (1977) may not be representative of small, mountainous drainage basins, for which bed load transport may be relatively more important. However, these results lend further support to the suggestion that bed load transport rates for the Mosquito Creek Tributary are overestimated by the model. 7.3.4 V O L U M E CHANGES Examination of the sediment budget constructed along the main channel shows that there is an aggradational trend in the lower reaches of the main channel (see figure 7.5). In the lower portion of the high-gradient channel zone, transport decreases to zero as a result of the increasing importance of steps in reducing effective transport gradients. However, in the upward portions of the high-gradient channel zone, the configuration of channel parameters leads to increasing Figure 7.5 Volume changes along main channel of Mosquito Creek Tributary. 191 transport. An anomalously high value of volume change occurs in the reach that defines the boundary between wandering and braided rivers. The relatively significant differences in the width and depth parameters result in large differences between inputs and outputs from this reach and, hence, high volume changes. The problem of unusually high volume changes is representative of the types of problems that arose frequently in model runs (including the final results illustrated in this chapter). Whenever channel variables change significantly, such as at channel pattern boundaries, accompanying significant changes in the transport results occurred. Therefore, a significant research problem is to further address changes in channel pattern. Of particular significance is an examination of characteristics of these changes (for example, the abruptness of changes over space and time). Along the zone of zero transport that is defined in the lower portion of the high-gradient fluvial reach, aggradation occurs in locations where significant amounts of sediment are transported from tributaries to the main channel and are deposited. The approach used in the later combined hillslope/channel surface model (see chapter 8), whereby the volume change is effectively "smeared" across a model grid cell, will reduce anomalously high volume changes, thereby decreasing their impact on model results. 7.3.5 DEBRIS FLOW ACTIVITY The debris flow activity rate of 0.01 events/km yr for the Queen Charlotte Islands 2 translates into an estimate of 6 debris flow events every 100 years in Mosquito Creek Tributary -y basin (drainage area is 5.5 km ). Initiation sites are determined using a random number generator. The program is set up to choose randomly six locations as initiation sites from among channel locations with gradients of at least 25°. Total volumes of debris flow erosion for each event range from 206 m through to 1 658 192 m (table 7.4). The values are on the low side in comparison to the average value of about 5 000 m based on the data of Rood (1984, 1990) for forested basins. The lengths of the erosional zones are of order 10 m. The debris flow data for forested basins of Rood (1984, 1990) have an 2 average erosion length of about 400 m, which is consistent with the present study. The lengths of the deposits calculated in the model run are generally of order lOMo 2 m. The somewhat low deposition lengths occur because the erosion volumes are somewhat on the low side. The calculated widths for debris flow channels, which are estimated from hydraulic geometry equations for steep mountain streams, are significantly lower than the values found in Rood's data. The debris flow data of Rood indicate widths in the initiation zone of order 10 m. 1 Past studies have shown that channel changes resulting from debris flows, when converted to "annual" rates of change, are similar in magnitude to fluvial changes (Church, 1997, pers. comm). However, because of the apparently high fluvial transport rates calculated in the model, annual changes due to fluvial processes significantly exceed those due to debris flow processes. 193 Table 7.4 Debris flows in the Mosquito Creek Tributary model run. Debris flow event 1 2 3 4 5 6 Channel 10 22 19 13 23 27 River Point 5 7 6 8 15 4 Total erosion (m /100yr) 262 416 780 262 1 658 206 3 194 CHAPTER 8: SURFACE MODEL RUNS 8.1 INTRODUCTION A two-dimensional (surface) version of the model is now investigated. The nonlinear transport model for landsliding is used in these final model runs as it provided a considerably better fit to the landsliding data for the Queen Charlotte Islands than the linear model. Cell dimensions in surface model runs are 100 m x 100 m (the same length as for the one-dimensional model runs). The model is run for 100 000 years, which is about ten times greater than the current interglacial. The numerical code is based on the translation of the partial differential equation used in the profile model runs (see equation 6.1, p. 128) into a form for two spatial dimensions. Differential equations are once again solved using the implicit Runge-Kutta subroutine "ddriv2.f' created by Kahaner et al. (1989). The number of equations to be solved in the 2-dimensional model increases exponentially as surface area increases. Hence, the running time of the numerical programming code also increases significantly as area increases. When channel processes are incorporated into the surface model the value of <p for each time step is based on volume changes calculated in the channel submodel, which incorporates both changes in transport along the channel and lateral inputs (in the profile runs the value was input directly into the model). The volume changes calculated along the channel are assigned to the nearest grid cell and distributed across the area of that cell. 8.2 T H E STUDY DRAINAGE BASIN The surface model is run for a small drainage basin in the Queen Charlotte Islands, Landrick Creek (figure 8.1). The input elevation data and river locations were obtained from 1:20 000 TRIM data for the Queen Charlotte Islands. The basin has a drainage area of about 2.1 km and the channel network is simple, consisting of a main channel with several tributary Figure 8.1 Landrick Creek drainage basin. 196 branches. The transport rule for Group B basins in implemented in these model runs as the landsliding data for Landrick Creek falls under the Group B transport rule category. The initial landscape used in the surface model run for Landrick Creek is shown in figure 8.2. The model surface extends beyond the drainage divide so that the boundary conditions specified in the model do not adversely affect results within the actual drainage basin. Although the model grid is only 24 x 24 grid points, the data are converted for plotting to 48 x 48 points using the Kriging geostatistical gridding method in order to improve the readability of the surface diagrams. 8.3 S U P P L Y - L I M I T E D R E M O V A L A L O N G T H E C H A N N E L N E T W O R K The incorporation of weathering into the calculations for hillslope transport is similar to that discussed for the case of one-dimensional model runs (see section 6.6). Supply-limited transport along the channel network, which has not heretofore been discussed, is dealt with in the following manner. The sediment budget is initially calculated assuming unlimited sediment supply. Next, sediment availability is calculated for each channel reach by multiplying the channel area (average width multiplied by reach length) by the depth of the sediment layer (heightsurface - heightbedrock) calculated for the grid cell associated with the particular channel reach. When net erosion is calculated for a reach, the length of time for available sediment in the channel reach to be depleted is calculated. If this time is less than the model time step, then the sediment supply is exhausted during this time step. Clastic sediment transport then ceases and bedrock erosion is initiated. The appropriate bedrock erosion equation (see equation 4.32 and 4.33, p. 106) is applied to the channel reach for the remaining length of time in the model time step. Figure 8.2 Initial surface for Landrick Creek Drainage Basin 198 8.4 H I L L S L O P E P R O C E S S E S (NO W E A T H E R I N G ) T h e surface m o d e l is i n i t i a l l y r u n f o r the case o f h i l l s l o p e transport w i t h u n l i m i t e d sediment supply. T h e m e a n , m a x i m u m and m i n i m u m changes i n e l e v a t i o n after 3 000, 10 000, 30 000 a n d 100 0 0 0 years are s h o w n i n table 8.1a ( a l l e l e v a t i o n and gradient statistics reported i n m o d e l results e x c l u d e p o r t i o n s o f the m o d e l surface w h i c h are not a part o f the L a n d r i c k C r e e k drainage basin). T h e m a x i m u m e l e v a t i o n w i t h i n L a n d r i c k C r e e k drainage b a s i n w a s i n i t i a l l y 500 m. This v a l u e w a s r e d u c e d b y o n l y about 3 m after 100 0 0 0 years. T h e m i n i m u m e l e v a t i o n increased b y about 3 c m . T h e m e a n change i n height w a s about - 3 c m after 10 0 0 0 years a n d about - 1 5 c m after 100 0 0 0 years (table 8.1b). T h e m a x i m u m and m i n i m u m changes i n height reached magnitudes o f about 15 m and - 1 7 m r e s p e c t i v e l y after 100 0 0 0 years. W i t h o u t the r e m o v a l o f mass f r o m the b a s i n b y f l u v i a l processes, there is not a s i g n i f i c a n t a m o u n t o f m a t e r i a l r e m o v e d f r o m the b a s i n . T h e decrease i n e l e v a t i o n due to e r o s i o n is g e n e r a l l y c o m p e n s a t e d for b y aggradation at other l o c a t i o n s i n the basin. H o w e v e r , some o f the m a t e r i a l e r o d e d at drainage d i v i d e s w a s transferred to the adjacent b a s i n , r e s u l t i n g i n a net loss o f m a t e r i a l (and, hence, net l o w e r i n g o f m e a n elevation) f r o m the L a n d r i c k C r e e k drainage basin. In order to assess changes i n gradient distributions across the landscape, gradients are c a l c u l a t e d u s i n g an a l g o r i t h m b y R i t t e r (1987), w h i c h is based o n the f o u r nearest neighbour elevations a n d is r e c o m m e n d e d b y H o d g s o n (1995). T h e average gradient o f the i n i t i a l landscape is about 2 2 . 6 ° , w i t h a m a x i m u m v a l u e o f 39.6° a n d a m i n i m u m v a l u e o f 1.4° (table 8.1c). A f t e r 100 0 0 0 years the average gradient is r e d u c e d o n l y m a r g i n a l l y to a v a l u e o f 2 2 . 2 ° , w i t h a m a x i m u m v a l u e o f 36.4° a n d a m i n i m u m v a l u e o f 1.5°. T h e changes i n e l e v a t i o n across the b a s i n after 100 0 0 0 years are s h o w n i n figure 8.3. clear pattern o f s e d i m e n t transfer w i t h i n the b a s i n emerges f r o m this d i a g r a m . A T h e greatest 199 Table 8.1 (a) Elevation (hillslope transport). Mean elevation (m) Standard deviation (m) Maximum elevation (m) Minimum elevation (m) Initial 244.199 3 000 years 244.189 10 000 years 244.171 30 000 years 244.136 100 000 years 244.051 120.740 120.690 120.588 120.357 119.775 500 499.958 499.847 499.439 497.228 15.0 15.001 15.003 15.010 15.033 (b) Changes in elevation (hillslope transport). After: Mean change in elevation (m) Maximum change in elevation (m) Minimum change in elevation (m) 3 000 years -0.010 10 000 years -0.028 30 000 years -0.063 100 000 years -0.148 1.092 3.189 7.379 15.365 -1.771 -4.694 -9.459 -16.692 Aggradation is represented by positive values and degradation is represented by negative values (c) Gradients (hillslope transport). Initial Tangent: 0.43 Mean gradient 22.6° Standard Tangent: 0.18 deviation 8.6° Maximum Tangent: 0.83 gradient 39.6° Minimum Tangent: 0.025 gradient 1.4° 3 000 years Tangent: 0.43 22.6° Tangent: 0.18 8.6° Tangent: 0.82 39.4° Tangent: 0.025 1.4° 10 000 years Tangent: 0.43 22.5° Tangent: 0.18 8.5° Tangent: 0.81 39.1° Tangent: 0.025 1.4° . • 30 000 years 100 000 years Tangent: 0.42 Tangent: 0.42 22.5° 22.2° Tangent: 0.17 Tangent: 0.16 8.4° 8.1° Tangent: 0.79 Tangent: 0.74 36.4° 38.3° Tangent: 0.025 Tangent: 0.026 1.5° 1.5° Figure 8.3 Changes in elevation after 100 000 years (hillslope transport only) 201 erosion occurred at the divides in the upper parts of the basin. The most significant aggradation occurred in the valleys of low-order tributary channels, and not the main channel valley. The steep and most active hillslopes deposit their eroded material into these steep channel reaches. Rather than moving directly down slopes and into the main channel, most material eroded on basin slopes goes into temporary storage in these low-order valleys before incorporation into the main channel, primarily by episodic debris flows. This result demonstrates clearly the importance of incorporating these steep channels into models of drainage basin evolution. The linear perturbations in the landscape represented by these channels affect the local rates and directions of sediment movement. Their unique transport and storage characteristics also alter the timing and rates of sediment transfers. High-gradient, low-order channels are the primary deposition zone for material eroded on hillslopes, and they represent the key link between hillslope processes and fluvial transport in the main valley. The critical role of these channels in sediment routing suggests the importance of advancing our knowledge of high-gradient channel transport processes. The gradients along these channel reaches are sufficiently low such that hillslope processes alone will be ineffective in moving this material. Further movement of material deposited in these channels within the model requires the successful representation of high-gradient channel processes. Adequate representation of fluvial transport and debris flow activity in these steep reaches is critical, yet these phenomena remain among the least well understood and least recognized aspects of sediment movement through the channel network. Most landscape evolution models have sufficiently coarse grids that these low-order channels are not resolved. This conceptualization of landscape configuration does not lend itself to adequate realization of the suite of processes that actually transfer sediment from hillslopes to the main channel. Low-order channels were resolved in the present model as the drainage basin was sufficiently small to allow for relatively small grid cell dimensions. 202 8.5 H I L L S L O P E P R O C E S S E S ( W E A T H E R I N G ) The adoption of weathering into model runs considerably slowed the running time of the model as there was a constant interplay between the transport of material and the rate of conversion of bedrock into loose sediments. The depth of the soil layer was initially set to 0.5 m across the model surface. The introduction of weathering into the hillslope models is expected to lead to decreased sediment transfers due to restrictions on sediment supply. The elevations and changes in elevation calculated in this model run are given in table 8.2. Although there was an overall reduction in the amount of erosion occurring in the basin due to limited sediment supply, the patterns of change are the same as those found in the model run with no weathering (see figure 8.3). The mean change in elevation of the landscape was about -10 cm, in contrast to the value of -15 cm for the case without weathering. The maximum deposition rate after 100 000 years remained approximately the same, about 15 m, in the model runs with and without weathering. However, the maximum erosion rate increased from -17 m for the case without weathering to about -8 m when weathering was introduced. The change in slope distributions across the landscape did not change significantly when weathering was introduced into the model (table 8.2c) 8.6 C H A N N E L P R O C E S S E S IN T H E S U R F A C E M O D E L Hillslope, weathering and channel processes were combined in the next model run. However, after several thousand years of model run time several locations of very high sediment accumulations developed at the junction of tributaries and the main channel and in the lower reaches of some of tributaries. These accumulations were often many tens of meters high and exceeded a hundred meters in several exceptional cases. This peculiarity precluded successful 203 Table 8.2 (a) Elevations (hillslope transport and weathering). Mean Elevation (m) Standard deviation (m) Maximum elevation (m) Minimum elevation (m) Initial 244.199 3 000 years 244.194 10 000 years 244.187 30 000 years 244.168 100 000 years 244.097 120.740 120.697 120.622 120.428 119.866 500 499.958 499.854 499.537 498.122 15.0 15.001 15.003 15.010 15.033 (b) Changes in elevation (hillslope transport and weathering). After: Mean change in elevation (m) Maximum change in elevation (m) Minimum change in elevation (m) 3 000 years -0.005 10 000 years -0.012 30 000 years -0.031 100 000 years -0.102 1.090 3.183 7.341 15.281 -0.696 -1.225 -2.737 -8.024 Aggradation is represented by positive values and degradation is represented by negative values (c) Gradients (hillslope transport and weathering). Mean gradient Standard deviation Maximum gradient Minimum gradient Initial Tangent: 0.43 22.6° Tangent: 0.18 8.6° Tangent: 0.83 39.6° Tangent: 0.025 1.4° 3 000 years Tangent: 0.43 22.6° Tangent: 0.18 8.6° Tangent: 0.82 39.4° Tangent: 0.025 1.4° 10 000 years Tangent: 0.43 22.5° Tangent: 0.18 8.6° Tangent: 0.81 39.0° Tangent: 0.025 1.4° 30 000 years Tangent: 0.42 22.5° Tangent: 0.17 8.4° Tangent: 0.79 38.2° Tangent: 0.025 1.5° 100 000 years Tangent: 0.42 22.2° Tangent: 0.17 8.1° Tangent: 0.74 36.4° Tangent: 0.026 1.5° 204 completion of the model run. Mass transport rates per unit width along the fluvial portion of the main channel (< 8°) ranged from 0 to 1.2 kg/m s at the beginning of the model run. These are within the expected range of values as indicated by the data compilation of Gomez and Church (1988). Sediment transport occurred in the lower reaches of many of the tributaries and values once again were reasonable. Furthermore, sediment transport rates remained at acceptable magnitudes even as the large sediment accumulations continued to grow. Therefore, it was not high magnitudes of transport alone which were the source of difficulties in the model run. Sediment transport rates calculated in the channel submodel for the main channel and its tributaries suggest that the large anomaly noted above is a result of high aggradation rates which continue to persist for long time periods. Aggradation rates may be high at locations where an upstream reach boundary has a high sediment transport rate and the corresponding downstream reach boundary has a low transport rate. In particular, it was thought that this situation may occur where sudden changes in channel variables, which result from changes in channel pattern, occur. In fact, the anomalously high aggradation rates did not coincide with locations where channel pattern changed. High aggradation rates occurred in some tributary reaches as a result of contrasting transport rates at the lower and upper boundaries of a reach. The high-gradient transport equation requires further investigation to ensure reasonable transport patterns are calculated in the channel submodel. However, most of the very high aggradation rates occurred where tributaries enter the main channel. Sediment was transported out of the tributary channel and large amounts of sediment deposited in the main channel reach located at the junction. In both of the above cases, the coarse grid resolution of the model ultimately leads to the continuation of sediment deposition over long time periods. Once some sediment begins to accumulate at a particular location, then it continues to accumulate unabated, the reasons for which are now explored. Channel variables controlling sediment transport in the model remain 205 approximately constant over time and change only relatively slowly. Therefore, zones of high aggradation remain locations of high aggradation for long time periods. In the real world, once sediment begins to accumulate at a particular location in the channel then the local gradient increases. The increased gradient leads to higher sediment transport rates at this location and the sediment is transported away from this location. However, these steep local gradients are not resolved in the model because of the relatively coarse grid resolution. Channel points are assigned the elevation of the closest grid cell in the channel submodel (see section 4.3.2.2). During initial model testing, the results of which are not reported herein, these elevations were used directly in the calculation of channel gradient. However, the direct application of these elevations in the calculations resulted in rather abrupt changes in channel gradient. The abrupt changes in gradient, in turn, led to unrealistically large and abrupt changes in sediment transport and very high volume changes in channel reaches. Therefore, it was decided instead to adopt a best-fit exponential relation to calculate gradients along the stream and eliminate these abrupt changes in gradient from the model (see section 4.3.2.2). However, significant problems arise when using this approach to calculate gradient in the surface model over extended time periods. Local channel gradients, which are in large part responsible for the elimination of major irregularities in channel morphology, are not resolved. Gradients remain relatively unaffected by the rapidly developing deposit and sediment continues to accumulate. Calculations of gradient must be made tractably adjustable to situations of accumulation and degradation in order to avoid unrealistic aggradation in the model. 8.7 CONCLUSIONS The results of the model run which considered the movement of sediment by hillslope processes only revealed noticeable erosion and deposition across the landscape. In particular, aggradation occurred in low-order valleys while erosion occurred at the upper divides. The 206 overall character of the landscape did not, however, change appreciably in 100 000 years of model run time. It is suggested that rates of landscape change calculated in this model are probably low compared to sediment transfers that occur during glacial episodes and immediately thereafter when sediment supply is particularly high and slopes are oversteepened. It is imperative to consider the influence and timing of major landscape perturbations, such as glaciation, on sediment transfers when examining landscape evolution. The effects of glacial episodes are likely to be particularly exaggerated immediately following glaciation (Church and Ryder, 1972), although the movement of paraglacial sediments through drainage basins may continue to affect transport patterns for many thousands of years (Church and Slaymaker, 1991). The surface model runs highlighted several key factors that must be considered when modelling sediment transfers across the landscape. Low-order channels create linear perturbations in the landscape which are locations of significant sediment accumulation. Therefore, the removal of the material from these locations by transport processes occurring in high-gradient channels (debris flows/fluvial transport) is of key importance when simulating sediment routing through the landscape. The adoption of a coarse model grid eliminates the possibility to consider high-gradient channels in the landscape, which are an essential link in the transfer of sediment through the landscape. Most recent landscape evolution models consider only low-gradient fluvial transport processes and, moreover, grid resolution is often too coarse to resolve low-order channels. In addition, the coarse resolution used in this and many other landscape evolution models cannot adequately resolve certain channel characteristics, such as gradient, which are required for the successful incorporation of channel processes in an integrated model of landscape evolution. In order for computing to remain tractable for large drainage basins, the number of grid cells must be minimized. However, increased resolution at channel network locations remains desirable in order to resolve low-order channels. It is recommended that adaptive grid methods, 207 which would allow for higher grid resolution at channel locations and lower grid resolution at interfluve locations, be considered for use i n future landscape evolution models. Increased resolution at channel locations would allow for better estimates of channel gradients, which would act to counter the accumulation of unusually large sediment deposits. 208 CHAPTER 9: CONCLUSIONS 9.1 INTRODUCTION There has been no clearly identified unifying theme for geomorphology since the period starting about 1950, other than the very general notion that as geomorphologists we are all interested in the study of exogenic processes. However, geomorphological studies sometimes appear to be disparate pieces of research. There is often no attempt to place individual studies within the broader scopes and aims of the discipline. In addition, the interactions amongst the various traditional process categories (e.g. fluvial, hillslope, glacial, periglacial) have been largely ignored. This categorization of the discipline, although useful in many ways, can exacerbate the neglect of process interaction which lies at the very heart of the study of morphological changes in the landscape. The structuring of research within a general framework would serve to highlight to geomorphologists the basic commonality in their research. It is suggested that the relations between scale and process specification outlined in this thesis represent a possible means of unifying many heretofore apparently isolated studies. The approach adopted for the study of smaller-scale phenomena has been dominantly process-oriented in recent years whereas larger scale studies have traditionally been more sedimentological and qualitative in nature. Most of the quantitative elements of large-scale studies have focussed on correlations between variables or dating sequences of deposits. A conundrum, until recently, has been how to address large-scale phenomena in a rigorous, process-oriented manner. This thesis attempts to resolve the issues raised by this conundrum. The contributions of the thesis to the modelling of landscape evolution are discussed in the following sections. 209 9.2 MODEL CONSTRUCTION There is general agreement in the discipline that direct scaling of small-scale process studies to large scales is an inadequate approach for large-scale studies. The approach followed in this thesis requires the careful definition of process at appropriate scales. It is proposed that a key to solving the difficulties involved in process specification at large scales lies in the realization that it is acceptable to neglect details of process operation which become irrelevant at large scales. For example, this notion has long existed in fluid mechanics. "Average" flow properties are often considered in hydraulic studies, despite the known fluctuations in the real system. With this realization it can be seen that both descriptive studies and numerical process studies are possible at all scales. A significant distinction between small-scale and large-scale process studies is the historical nature of the latter. The historical element leads to difficulties in the calibration of transport equations, which requires the use of process data obtained from real landscapes. The historical component of large-scale phenomena diminishes the convenience with which such data can be obtained. However, this problem can begin to be addressed by the collection of large data bases (e.g. landslide inventories such as the data set of Rood, 1984; 1990) over sufficiently large spatial scales and the largest resolvable time scales. In addition, the development of absolute dating methods may provide a means to estimate transport rates over very long time periods. The hillslope analysis in this thesis is based on the concept that gradient is the principle controlling factor affecting transport rates. Hillslope processes were broken down into two main components: (i) slow, quasi-continuous mass movements and (ii) fast, episodic mass movements. Theoretical considerations suggest that a relation may exist between creep rate and gradient. However, such a relation was not identified on the basis of existing data. This may partly reflect the short-term nature of most studies and inadequate control of variables during measurement 210 programs. The analysis of fast, episodic mass movements was based on an extensive data set for landsliding in the Queen Charlotte Islands, British Columbia (Rood, 1984, 1990). Results of this study may be more representative than previous studies in which transport equations were not calibrated. A nonlinear relation was found to be more appropriate than linear diffusion for the modelling of fast, episodic mass movements. Hillslope characteristics do not change dramatically as model scales increase. Therefore, process description is not challenged as hillslope models are applied over a variety of scales. However, as scales of channel systems increase their characteristics change significantly. There must be some assurance that the chosen transport relation can be applied over a wide range of scales. A robust equation is required for the modelling of sediment transport. A revised form of the Bagnold equation was calibrated on the basis of field data for flumes and rivers at a variety of scales. A strong relation between excess stream power and transport rate was demonstrated. The large range of model scales over which the channel submodel was applied and found to perform reasonably demonstrate its robustness. The application of the revised Bagnold relation to steep channels required the adjustment of average channel gradients to account for reduced gradients found behind sediment wedges. The adjustment was made on the basis of field data from the Queen Charlotte Islands, British Columbia. Relations for suspended load and bedrock erosion are based on past field studies of these phenomena (Church et al., 1989; Seidl et al., 1994). Equations for debris flow transport were constructed empirically on the basis of data collected by Rood (1984, 1990) for the Queen Charlotte Islands, British Columbia. 9.3 M O D E L RUNS Results of 1-dimensional model runs demonstrated that diffusivities used by previous modellers overestimated hillslope transport rates. The landscapes became relatively flat in 211 unreasonably short time periods. The nonlinear transport rules which were defined in this thesis provided realistic rates of landscape change. The landscapes exhibited appropriate relief based on a qualitative interpretation of model results. Furthermore, recharge rates calculated in gullies compared well with field data. Nonlinear diffusion allows for rapid sediment transfers on steep slopes and lower transfer rates, as a result of creep processes, on gentler slopes. The distinction between these two transport regimes is preserved. Landscapes with steep gradients are relatively unstable and hillslope erosion is rapid. Unstable slopes are eventually eliminated and landscape change proceeds at a reduced rate. The incorporation of weathering into the model runs for the nonlinear transport relation produced minimal changes in overall rates of landscape change. The resulting depths of the sediment layer were reasonable in comparison to observations of soil depths made in the Queen Charlotte Islands, British Columbia. Results of the channel submodel showed that the model provided reasonable estimates of channel variables. The relation performed reasonably in model tests covering a wide range of scales. The surface model runs highlighted the key role of high-gradient channels in sediment transfers through drainage basins. Steep channels represent the primary link between hillslopes and the main channel. In addition, the difficulties encountered when channel processes were incorporated into the model runs highlighted the contrasting scales of resolution which are required when considering channel versus hillslope elements in a landscape evolution model. 9.4 "UNKNOWNS" IN T H E MODELLING OF GEOMORPHOLOGY The most intriguing finding of this thesis may lie in the realization of how much is not known about the operation of geomorphological processes and their interactions at large scales. This was noticed as result of the need for explicit information in order to develop reasonable 212 formulations for processes in the model. This "missing information" has proven to be a difficult, and often frustrating, aspect of this research project. Each unknown or unsubstantiated aspect of geomorphology identified in this thesis represents a major research project in itself. Moreover, this list of unknowns is long. In several instances, expedient approaches to certain dilemmas had to be adopted, despite the obvious shortcomings they represent. The identification of weaknesses in geomorphological research provides a basis for the definition of an agenda for the further study of large-scale geomorphological processes. What, then, are these shortcomings in our knowledge of geomorphological processes? The major concepts which require visitation (or re-visitation) by geomorphologists are outlined below (the topics are addressed in the order in which they first appeared in the thesis). Virtual Velocity (Chapter 2) The definition for time scales introduced in this thesis requires a knowledge of the "virtual velocity" of sediment movement through the system. The concept of virtual velocity requires consideration of both the transport and storage of sediment. Smaller-scale studies often focus on the actual movement of sediment, which occurs generally during relatively discrete events. However, over the longer time scales of landscape evolution, the relative importance of sediment storage increases and the discrete nature of transport events is highlighted. Hence, valley storage was highlighted in the model. The amount of time sediment is actually in transport (excluding slow mass movements such as creep and the movement of fine material by rivers) is negligible in comparison to the time that is spent at rest. A key step in estimating virtual velocity is an examination of residence times of stored sediment. Sedimentological studies often investigate long-term reservoirs of stored sediment (e.g. alluvial fans, fluvial deposits etc.). However, such features are seldom approached from the 213 perspective of associated residence times. Basic concepts of reservoir theory (Bolin and Rodhe, 1973) can be applied to sediment storage reservoirs, such as those listed above. In addition, a sediment budget framework, which has gained prominence in geomorphological research in recent years, focusses on changes in storage. Such studies provide an obvious entry point for studies considering residence times of stored sediment. Slow Mass Movement (Chapter 3) Results of modelling in this thesis showed that creep transport is insignificant relative to faster processes in high-gradient landscapes. However, creep and associated effects, such as animal processes, represent the primary hillslope process in operation on gentler slopes. Therefore, it is imperative to continue investigation of long-term creep rates and their controlling factors. Despite the significant number of creep studies that have been undertaken in the past several decades, there is a notable lack of knowledge regarding factors controlling creep transport rates. The time-scales of most studies are too short to understand adequately variations in process rates. Studies have generally not been adequately controlled and, hence, our knowledge of creep remains insufficient for modelling purposes. The research reads as a set of individual field observation studies, from which few generalizations can be made. Future studies of creep processes should be approached in a more rigorous manner. The laboratory study of VanAsch et al. (1989) is important for two reasons: (i) it demonstrates the potentially useful role of laboratory experimentation in geomorphology and, in particular, for creep studies (it is easier to control variables in a laboratory situation) and (ii) it suggests the possibility of a nonlinear creep/gradient relation (which contradicts the formulation of creep often used in landscape evolution models). 214 Rapid, Episodic Mass Movements (Chapter 4) Many rapid, episodic mass movement studies have been approached from the perspective of mechanics, with a particular emphasis on slope stability and hazard prediction. There has been minimal documentation of long-term landsliding rates in the literature. However, rapid mass movements are extremely important transport processes in steep terrain. Despite the importance of such events, few studies have documented quantitatively the long-term role of fast mass movements in landscape evolution. The landslide analysis in the present study represents only a preliminary effort. However, the calibration of transport equations undertaken for landsliding, which was based on a large data set, is exceptional in the landscape modelling literature. It represents a viable approach for improving large-scale models of rapid, episodic mass movements. In particular, the claim put forth in this thesis, that landslide occurrence is described as a nonlinear process, requires further examination. The nonlinearity implies that after events which lead to oversteepening in the landscape (e.g. glaciation, tectonic movements), landscape change will proceed at an accelerated rate until the very steep slopes are eliminated. In addition, large landslide data bases are required to further assess the role of geology and climate in determining rates of fast, episodic mass movements. Weathering (Chapter 3/Chapter 6) The one-dimensional model runs in Chapter 6 demonstrated the importance of limited sediment supply on erosion rates. This knowledge is critical in landscape evolution as weathering determines the rate at which sediment becomes available for transport. The importance of weathering was recognized by several other landscape modellers, several of whom included it in their models (e.g. Tucker and Slingerland, 1994; Anderson and Humphrey, 1989). Given the persistent usage of the terms "transport-limited" and "supply-limited", and the 215 potential impact of weathering on transport rates, the lack of knowledge which remains regarding long-term weathering rates (i.e. conversion rates of bedrock into sediment) represents a significant oversight on the part of the geomorphological community. The neglect may be due in part to the logistical difficulties involved in its study as the conversion to bedrock occurs at the base of the regolith column and is not readily observable. Nonetheless, carefully controlled field observations would contribute significantly to our understanding of this phenomenon. Bed Load Transport in High-gradient Streams (Chapter 4/Chapter 7) Most earlier landscape models made use of a stream power-type equation for the calculation of bed load transport. This thesis improves upon earlier modelling efforts by reexamining the Bagnold relation carefully and calibrating it using a large data set. The strong performance of the relation used in this thesis suggests that Bagnold-type relations should be subject to further investigation. This equation may be, in effect, a finely tuned scale correlation. However, the possibility that the equation may have some underlying mechanical basis has not been eliminated and should be explored. The Bagnold-type equation provided unrealistic transport estimates for steep channels when used in its original form. The poor performance of the relation when applied directly to high-gradient channels is not surprising, given that it was calibrated using the data of Gomez and Church (1988), for which gradients were always lower than 0.5°. The approach which was adopted in this study, whereby the "effective" transporting gradient became the gradient of the sediment wedge when average channel gradients exceeded 1°, represents an initial attempt to address this issue. The morphology of high-gradient streams has been the subject of recent research (e.g. Grant et al.,1990; Wohl et al., 1997). However, a detailed investigation of sediment transport in step-pool streams has not been undertaken in the geomorphological literature. These high-gradient fluvial reaches represent an important "connecting" zone for 216 sediment deposited by debris flow events in steep channels (which, in turn, arrived through mass movements) to eventually be incorporated into higher-order, low-gradient channels. For this reason, the study of transport through steep reaches is important in terms of continuity of transport through the channel system. Suspended Load (Chapter 4/Chapter 7) A distinction is made in this thesis between the transporting mechanisms which occur in the fluvial system, namely the suspended and bed load. Suspended load was calculated from a suspended load/contributing drainage area correlation (there is an assumed relation between contributing drainage area and discharge). However, the suspended load data of Church et al. (1989), which were used in this thesis, show considerable variation for drainage basins of approximately the same size. Further research is required to establish the factors which control suspended load transport rates. It was assumed that suspended sediment, once entrained, is evacuated directly through the system (i.e. it is considered to be a part of the wash load). This is a simplification of the true behaviour of suspended sediment, some of which may be re-deposited between high flows. Suspended load, in reality, can be a part of the wash or bed material load. This material may be evacuated directly through the system or go into storage. Transport patterns of suspended sediment along the fluvial system require further consideration. This is particularly important as in many cases suspended load comprises a large proportion of total sediment load. Variables Required for Fluvial Transport Calculations (Chapter 4/Chapter 7) The discharge/area relation adopted in this thesis was based on a study by Church (1997), in which regional variations in discharge/area relations across British Columbia were investigated. The results supported the exponent of 0.7 found in many other studies. The 217 coefficient was found to vary significantly across the province. discharge/area relations need to be undertaken in other locations. Regional studies of In particular, the relative importance of climate, geology and topography in determining the value of the coefficient should be investigated. The grain size relation used in this thesis is a reorganized threshold regime equation. This approach led to reasonable results in the model runs and provides a relatively straightforward method for evaluating grain size in future modelling efforts. However, the theoretical underpinnings of this relation require further research. In addition, as theory of grainsize changes along channel networks progresses, relations used in the model should become more physically based. In particular, abrupt grain size changes which may occur at confluences should eventually be incorporated into the channel submodel. Channel pattern is calculated in the model as this information is necessary for the selection of appropriate hydraulic geometry relations. However, further investigation of factors controlling channel pattern changes should be undertaken in order improve our physical understanding of this phenomenon and to account for the strength of the criterion found in this study. The use of hydraulic geometry equations to calculate channel width and depth needs further evaluation. The data used in such analyses should ideally range from small, high- gradient channels through to large, main river channels. Most studies have focussed on channels of intermediate scale. Of further significance is the fact that the model predicts only "gross" downstream changes in width and depth as a result of using these equations. The equations do not capture variations that typically exist within the "broad" downstream trends. It is possible that natural variability of width/depth may affect significantly the sediment transport regime. This notion should be investigated further to obtain a better understanding of the role of variability in channel morphology in determining transport regimes along the channel system. 218 High rates of sediment erosion and deposition were often calculated in the model at locations where channel variables changed abruptly. Abrupt changes in channel variables result from sudden switches in the defined channel type definition from "wandering" to "braided" to "debris flow". When channel pattern changes so does the channel morphology. This, in turn, causes relatively abrupt changes in transport rates which are dependent on channel morphology. Model results highlighted the importance of exploring the effects on results of continuity/discontinuity in the definition of model parameters. While there is a significant literature describing the hydraulic geometry of meandering and wandering channels, there has been little documentation of downstream changes in river variables for braided rivers. A simple set of rules was created for the calculation of width and depth for braided rivers, which provided reasonable results. However, an extended effort is required to more formally address the hydraulic geometry of braided rivers. Debris Flows (Chapter 4/Chapter 7) Many studies have focussed on the detailed mechanics of debris flow events or approached them from a hazard prediction perspective. In contrast, a model of long-term debris flow initiation, transport and deposition in drainage basins was created in this thesis. The model was created with reference to the large data base of Rood (1984), which provided an essential link to reality. Further data collection, of the type performed by Rood (1984, 1990) is necessary in order to define long-term transport rates and their controlling factors. The effects of climate and geology on debris flow occurrence can be examined by ensuring adequate control of variables. In addition, studies of downstream changes in gully characteristics, such as width, are required in order to estimate the amount of debris stored in gullies better. 219 The resolution of grid cells at gully locations should ideally be smaller than grid cells for hillslopes. This would allow for better resolution of gully recharge and scour rates, which are dependent on gully morphology and local depth to bedrock. This would lead to improved estimates of debris flow volumes. This issue is discussed further in this chapter in a section concerned with model resolution. Valley Flat (Chapter 5) The valley flat was examined in this thesis from the perspective that it represents the interface between hillslope and fluvial regimes. Hillslope/fluvial interactions have been largely neglected in geomorphological research. In this thesis, algorithms were created for filling and incision in the valley flat. Further field investigation is required in order to better understand the processes occurring during infilling and incising episodes in valleys and to obtain data for model calibration. Sediment derived from hillslopes can be evacuated from the drainage basin only after it is entrained by the river. It is important to assess the length of time that sediments are stored in the valley flat before they are entrained by the fluvial system. Time-scales of deposition and reentrainment of sediment after it has been incorporated into the fluvial system should also be investigated. These issues are relevant to the definition of virtual velocities that were discussed earlier in this chapter. Unstable vs. Stable Landscapes (Chapter 6) The concept of "unstable" and "stable" landscapes provides a basis for defining dominant styles of landscape evolution. The frequency of events which lead to "oversteepening" of slopes has direct implications on the rate of progress of geomorphological changes. In landscapes with "oversteepened" slopes, hillslope transport occurs at accelerated rates. Thereafter, when creep 220 processes dominate on hillslopes, fluvial transport becomes the primary transport mechanism in the landscape. In the latter situation, some sediment is released from long-term storage while other sediment goes into long-term storage as the river migrates laterally across the valley floor. The implications of these ideas on landscape evolution in different regions should be explored. Flow Regime for Sediment Transport (Chapter 7) The annual flow regime of a channel system is very important in determining annual transport rates. This was illustrated by the contrasting flow regimes of, for example, the Fraser River and the Vedder River in this study. In the former case, bed load transport occurs over several months each year whereas bed load transport in the Vedder River is restricted to several days a year. The duration of effective flows is expected to change downstream in a channel system. Significant events in steep debris flow channels are relatively rare, while sediment transport events in major rivers are common in comparison. The effective flows for sediment transport should be investigated more thoroughly in order to quantify the downstream changes in transporting events. This information could then be translated into patterns of effective discharge frequencies for modelling major channel systems. Model Resolution (Chapter 5/Chapter 8) The contrasting scales of grid cells required for convenient numerical modelling of hillslope and fluvial processes proved to be an obstacle during the conversion of transport rules into numerical modelling code. Model grid cells were restricted to a sufficiently small size in the present model so that the "smearing" of fluvial processes across the whole grid cell would not result in significant distortions in the patterns of morphological change. This approach is likely to meet with some limitations as drainage basin size increases. It would become necessary to 221 increase cell dimensions in order for the number of model grid cells to remain tractable for computing. Hence, the "smearing" approach would become problematic as grid cells would be significantly greater than channel dimensions. In addition, the adoption of a uniform grid cell size across the model domain led to particular difficulties in the modelling of fluvial processes. Fluvial transport is sensitive to local gradient. However, grid cells remain relatively large in the model and only regional river gradients are obtained from the elevation data in the present model. Therefore, when significant accumulations of sediment were deposited at locations where tributaries enter the river, local channel gradients, which would normally act to erode away some of the deposit, were not calculated in the model. Adaptive grid strategies, which allow for changes in cell dimensions (over space or time), should be explored. Cell density would decrease in hillslope regions and increase at channel locations. This would allow for greater resolution of channel variables which are important in determining sediment transport rates. Model Assessment (Whole thesis) A particular difficulty encountered in large-scale landscape evolution modelling is the issue of how to assess model performance. No acceptable approach has been put forward to resolve this issue and any such resolution is not expected to be a simple matter. The difficulties are not surprising given the historical nature of landscape evolution. To date the most common approach is to compare the resultant landscape with the region being modelled and to see if it appears to provide reasonable approximations to expected rates of landscape change. However, this approach is unsatisfactory as it relies on largely subjective judgements. The careful calibration of equations in the present model, which made use of real data, provides some basis for accepting the overall rates of landscape change as being reasonable. 222 One possible approach for model assessment is to compare model results with key landscape parameters of real landscapes (e.g. hypsometric integral, fractal dimension). However, such parameters may not be sensitive to the rates of landscape change found in the present study. Further exploration is required in order to identify landscape parameters which are sensitive to the overall rates of morphological change resulting from geomorphological processes. Process Generalization/Long-term Transport Rates (Whole thesis) A fundamental shortcoming in large-scale geomorphology, until recently, has been the lack of a framework for process specification at large scales. Process generalization in numerical modelling has shown that it has the potential to eliminate this void. However, this approach has not yet been adopted in geomorphological research as a whole. Its use remains confined primarily to modellers of landscape evolution, who often approach landscape modelling from a geodynamic/geophysical perspective. Assessment of current generalizations used in landscape models, most of which are based on a diffusion-type equation for hillslopes and a stream power equation for rivers, must continue. The relations should be compared to field evidence documenting transport rates. Absolute dating methods should be further examined as they may provide the key for the estimation of long-term transport rates. 9.5 FINAL THOUGHTS This thesis has attempted to illustrate possible methods of approach for the numerical representation of geomorphological processes at large scales. A framework for the study of large-scale process geomorphology, based on a generalized physics representation of processes, was outlined and assessed. Calibrations of transport equations were based on field data in order to ensure adequate representation of processes. The model runs provided new insights into the operation of transport processes and, most importantly, provided a stimulus for the framing of a 223 series of key issues that must be addressed in order to improve our understanding of large-scale geomorphological processes. Although the generalized physics approach has been adopted previously by modellers of landscape evolution, it has not yet had a large impact upon the geomorphological community as a whole. Many of the recent modellers have been concerned with both the endogenic and geomorphological components of landscape evolution, the latter which is the focus of this study. The incorporation of geomorphology into landscape evolution models in a meaningful way requires that shortcomings in our geomorphological knowledge be addressed. This study represents an initial attempt to address these issues in a critical manner. It is hoped that this thesis contributes to a reconsideration of how geomorphological processes operate at large scales. Maybe the time will soon arrive when Anderson and Humphrey can feel that geomorphologists are attempting to return the ball in a meaningful way. 224 REFERENCES Ackers, P. (1964) Experiments on small streams in alluvium, American Society of Civil Engineers, Proceedings, Journal of the Hydraulics Division, 90, 1-37 Adams, J. (1980) Contemporary uplift and erosion of the Southern Alps, New Zealand, Geological. Society of America Bulletin, 91, Part 1, 2-4 Ahnert, F. (1976) Brief description of a comprehensive three-dimensional process-response model of landform development, Zeitschrift fur Geomorphologie., 25, 29-49 Ahnert, F. (1977) Some comments on the quantitative formulation of geomorphological processes in a theoretical model, Earth Surface Processes, 2, 191-201 Ahnert, F. (1987a) Approaches to dynamic equilibrium in theoretical simulation of slope development, Earth Surface Processes and Landforms, 12, 3-15 Ahnert, F. (1987b) Process-response models of denudation at different spatial scales, In F. Ahnert, (ed.) Geomorphological models - theoretical and empirical aspects, Catena Supplement, 10, 31-50 Ahnert, F. (1988) Modelling landform change, In M.G. Anderson Geomorphological Systems, John Wiley and Sons, Chichester, 375-400 (ed.) Modelling Anderson, E.W. (1977) Soil Creep: An Assessment of Certain Controlling Factors with Special Reference to Upper Weardale, England, Ph.D. thesis, University of Durham (cited in Saunders and Young, 1983) Anderson, R.S. (1994) Evolution of the Santa Cruz Mountains, California, through tectonic growth and geomorphic decay, Journal of Geophysical Research, 99, 20161-20179 Anderson, R.S. and Humphrey, N.F. (1989) Interaction of weathering and transport processes in the evolution of arid landscapes, In T.A. Cross (ed.) Quantitative Dynamic Stratigraphy, Prentice Hall, Englewood Cliffs, N.J., 349-361 Andrews, D.J. and Bucknam, R.C. (1987) Fitting degradation of shoreline scarps by a nonlinear diffusion model, Journal of Geophysical Research, 92, 12857-12867 225 Avouac, J.P and Burov, E.B. (1996) Erosion as a driving mechanism of intracontinental mountain growth, Journal of Geophysical Research, 101, 17747-17769 Bagnold, R A . (1977) Bed load transport by natural rivers, Water Resources Research, 13, 303312 Bagnold, R.A. (1980) A n empirical correlation of bedload transport rates in flumes and natural rivers, Proceedings of the Royal Society of London, A , 372, 453-473 Bagnold, R A . (1986) Transport of solids by natural water flow: evidence for a worldwide correlation, Proceedings of the Royal Society of London, A, 405, 369-374 Barr, D.J. and Swanson, D.N. (1970) Measurement of creep in a shallow, slide-prone till soil, American Journal of Science, 269, 467-480 Bierman, P.R. (1994) Using in situ produced cosmogenic isotopes to estimate rates of landscape evolution: A review from the geomorphic perspective, Journal of Geophysical Research, 99, 13885-13896 Bolin, B. and Rodhe, H. (1973) A note on the concepts of age distribution and transit time in natural reservoirs, Tellus, X X V , 1, 58-62 Bray, D.I. (1972) Generalized regime-type analysis of Alberta rivers, Ph.D. thesis (unpublished), University of Alberta, Edmonton, 232 pp. Bray, D.I. (1979) Estimating average velocity in gravel-bed rivers, American Society of Civil Engineers, Journal of the Hydraulics Division, 105, 1103-1122 Bridge, J. and Leeder, M . (1979) A simulation model of alluvial stratigraphy, Sedimentology, 26, 617-644 Budel, J. (1982) Climatic geomorphology (translated by L. Fischer and D. Busche), Princeton University Press, Princeton, 443 pp. 226 Burbank, D., Leland, J., Fielding, E., Anderson, R., Brozovic, N., Reid, M. and Duncan, C. (1996) Bedrock incision, rock uplift and threshold hillslopes in the northwestern Himalayas, Nature, 379, 505-510 Burrin, P.J. and Jones, D.K.C (1991) Environmental processes and fluvial responses in a small temperate zone catchment: A case study of the Sussex Ouse Valley, Southeast England, In L.Starkel, K.J. Gregory and Thornes, J.B. (eds.) Temperate Palaeohydrology, John Wiley and Sons, Chichester, 217-251 Carson, M A . and Kirkby, M.J. (1972) Hillslope Form and Process, Cambridge University Press, Cambridge, 475 pp. Chandler, R.J. and Pook, M.J. (1971) Creep movement of low gradient clay slopes since the Late Glacial, Nature, 229, 399-400 Charlton, F.G., Brown, P.M. and Benson, R.W. (1978) The hydraulic geometry of some gravel rivers in Britain, Hydraulics Research Station, Wallingford, IT180, 48 pp. Plus tables and figures Chase, C.G. (1992) Fluvial landsculpting Geomorphology, 5, 39-57 and the fractal dimension of topography, Chorley, R.J. (1969) The drainage basin as the fundamental geomorphic unit, In R.J. Chorley (ed.) Water, Earth and Man, Methuen and Co. Ltd., London, 77-98 Church, M . (1980) On the equations of hydraulic geometry, Department of Geography, University of British Columbia unpublished report, 59 pp. Church, M. (1984) On experimental method in geomorphology, In T.P Burt and D.E. Walling (eds.) Catchment Experiments in Fluvial Geomorphology, Geo Books, Norwich, United Kingdom, 563-580 Church, M . (1996) Space, time and the mountain - how do we order what we see?, In B.L. Rhoads and C.E. Thorn (eds.) The Scientific Nature of Geomorphology, John Wiley, Chichester, 147-170 227 Church, M. (1997) Regionalised hydrological estimates for British Columbia: first approximation of scale effects, Report for Resources Inventory and Data Management Branch, British Columbia Ministry of Environment, Lands and Parks, Victoria, Contract DEF6156, 47 pp. Church, M. and Jones, D. (1982) Channel bars in gravel bed rivers, In R.D. Hey, J.C. Bathurst and CR. Thorne (eds.) Gravel-bed Rivers, John Wiley, Chischester, 291-338 Church, M. and Kellerhals, R. (1978) On the statistics of grain size variation along a gravel river, Canadian Journal of Earth Sciences, 15, 1151-1160 Church, M., Kellerhals, R. and Day, T.J. (1989) Regional clastic sediment yield in British Columbia, Canadian Journal of Earth Sciences, 26, 31-45 Church, M . and Ryder, J.M. (1972) Paraglacial sedimentation: a consideration of fluvial processes conditioned by glaciation, Geologica. Society ofAmerica Bulletin, 83, 3059-3071 Church, M. and Slaymaker, O. (1989) Disequilibrium of Holocene sediment yield in glaciated British Columbia, Nature, 337, 452-454 Colman, S.M. and Watson, K. (1984) Ages estimated from a diffusion equation model for scarp degradation, Science, 221, 263-265 Crickmay, C H . (1974) The Work of the River, Macmillan, London, 271 pp. Culling, W.E.H. (1960) Analytical theory of erosion, Journal of Geology, 68, 336-344 Culling, W.E.H. (1963) Soil creep and the development of hillside slopes, Journal of Geology, 71, 127-161 Culling, W.E.H. (1965) Theory of erosion on soil-covered slopes, Journal of Geology, 73, 230254 Davis, W.M. (1899) The geographical cycle, Geographical Journal, 14, 481-504 Day, T.J. (1969) The channel geometry of mountain streams, M A . thesis (unpublished), The University of British Columbia, 59 pp. 228 Day, M.J. (1977) Some Effects of Water on Slope Recession, Ph.D. thesis, University of East Anglia, Norwich (cited in Saunders and Young, 1983) Dedkov, A . V . and Duglav, V . A . (1972) Slow movements of soil masses on grassed slopes, National Lending Library Transl. Progr., R T S 7548 (Russian original Izv. Akad. Nauk SSSR ser. Geogr. (1967), 90-93) Densmore, A . L . , Anderson, R.S., McAdoo, B.G. and Ellis, M.A. (1997) Hillslope evolution by bedrock landslides, Science, 275, 369-372 Dietrich, W.E. and Dunne, T. (1978) Sediment budget for a small catchment in mountainous terrain, Zeitschrift fur Geomorphologie Supplementband, 29, 191-206 Dietrich, W.E., Reiss, T., Hsu, M.-L. and Montgomery, D.R. (1995) A process-based model for colluvial soil depth and shallow landsliding using digital elevation data, Hydrological Processes, 9, 383-400 Dooge, J.C.I. (1986) Looking for hydrologic laws, Water Resources Research, 22, 46S-58S Drage, B. and Carlson, R. (1977) Hydraulic geometry relationships for northern braided rivers, Proceedings, 3 National Hydrotechnical Conference, Canadian Society for Civil Engineering 1977, 235-251 r d Drozdowski, E. and Bergland B. (1976) Development and chronology of the lower Vistula River valley, North Poland, Boreas, 5, 95-107 Dunne, (1980) Formation and controls of channel networks, Progress in Physical Geography, 4, 211-239 Einstein, H.A. (1950) The bedload function for sediment transportation in open channel flows, U.S. Department ofAgriculture, Technical Bulletin, 1026, 70 pp. Emmett, W.W. (1975) The channels and waters of upper Salmon River area, Idaho, United States Geological Survey, Professional Paper, 870-A, 116 pp. 229 England, P.C. and Molnar, P. (1990) Surface uplift, uplift of rocks, and exhumation of rocks, Geology, 18, 1173-1177 Everett, K.R. (1963) Slope movement, Neotoma Valley, southern Ohio, Report for Institute of Polar Studies, Ohio State University, 6 (cited in Young, 1974) Eyles, R.J. and Ho, R. (1970) Soil creep on a humid tropical slope, Journal of Tropical Geography, 31, 40-42 Finlayson, B.L. (1981) Field measurements of soil creep, Earth Surface Processes and Landforms, 6, 35-48 Flemings, P.B. and Jordan, T.E. (1989) A synthetic stratigraphic model of foreland basin development, Journal of Geophysical Research, 94, 3851-3866 Gilbert, G.K. (1877) Report of the Geology of the Henry Mountains (Utah), US Geographical and Geological Survey of the Rocky Mountains Region, US Government Printing Office, Washington, D C , 169 pp. Gomez, B. and Church, M . (1988) A catalogue of equilibrium bedload transport data for coarse sand and gravel-bed channel, Department of Geography Report, University of British Columbia, 90 pp. Gomez, B. and Church, M . (1989) A n assessment of bed load sediment transport formulae for gravel bed rivers, Water Resources Research, 25, 1161-1186 Grant, G., Swanson, F. and Wolman, M.G. (1990) Pattern and origin of stepped-bed morphology in high-gradient streams, Western Cascades, Oregon, Geological Society of America Bulletin, 102, 340-352 Greene, M.T. (1997) What cannot be said in science, Nature, 388, 619-620 Griffiths, G.A. (1981) Stable-channel design in gravel-bed rivers, Journal of Hydrology, 52, 291305 Hack, J.T. (1960) Interpretation of erosional topography in humid temperate regions, American Journal of Science, 258A, 80-97 230 Hanks, T.C., Bucknam, R.C., Lajoie, K.R. and Wallace, R.E. (1984) Modification of wave-cut and faulting-controlled landforms, Journal of Geophysical Research, 89, 5771-5790 Heimsath, A . M . , Deitrich, W.E., Nishiizumi, K. and Finkel, R.C. (1997) The soil production function and landscape equilibrium, Nature, 388, 358-361 Henderson, F.M. (1963) Stability of alluvial channels, Transactions of the American Society of Civil Engineers, 128, 657-686 Hirano (1975) Simulation of developmental process of interfluvial slopes with reference to graded form, Journal of Geology, 83, 113-123 Hodgson, M.E. (1995) What cell size does the computed slope/aspect angle represent, Photogrammetric Engineering and Remote Sensing, 61, 513-517 Hooke, J.M. (1980) Magnitude and distribution of rates of river bank erosion, Earth Surface Processes, 5, 143-157 Howard, A . (1994) A detachment-limited model of drainage basin evolution, Water Resources Research, 30, 2261-2285 Howard, A.D., Keetch, M.E. and Vincent, C L . (1970) Topological and geometrical properties of braided rivers, Water Resources Research, 6, 1674-1688 Howard, A . D . and Kerby, G. (1984) Channel changes in badlands, Geological Society of America Bulletin., 94, 739-752 Kahaner, D., Moler, C. and Nash, S. (1989) Numerical Methods and Software, Prentice-Hall, 495 pp. Kellerhals, R. (1967) Stable channels with gravel paved beds, Journal of the Waterways, Harbours and Coastal Engineering Division, Proceeding of American Society of Civil Engineers, 93, 63-84 231 King, L. (1953) Canons of landscape evolution, Bulletin of the Geological Society of America, 64, 721-752 King, L. (1962) The Morphology of the Earth, Hafner, New York, 699 pp. Kirkby, M.J. (1964) in 'Slope profiles: a symposium', Geographic Journal, 130, 86 Kirkby, M.J. (1967) Measurement and theory of soil creep, Journal of Geology, 75, 359-378 Kirkby, M.J. (1971) Hillslope process-response models based on the continuity equation, In Brunsden, D. (ed.) Slope form and process, Institute of British Geographers Special Publication, 3, 15-30 Kirkby, M.J. (1976) Deterministic continuous slope models, Zeitschriftfur Geomorphologie, 25, 1-19 Kirkby, M.J. (1977) Soil development as a component of slope models, Earth Surface Processes, 2, 203-230 Kirkby, M.J. (1986) A two-dimensional simulation model for slope and stream evolution, In A.D. Abrahams (ed.) Hillslope Processes, 203-222 Kirkby, M.J. (1987) General models of long-term slope evolution through mass movement, In M.G. Anderson, and Burt, T.P. (eds.), Slope Stability (John Wiley, Chichester, 259-379 Kirkby, M.J. (1992) A n erosion-limited hillslope evolution model, Catena Supplement, 23, 157187 Kooi, H. and Beaumont, C. (1994) Escarpment evolution on high-elevation rifted margins: Insights derived from a surface processes model that combines diffusion, advection, and reaction, Journal of Geophysical Research, 99, B6, 12191-12209 Kooi, H. and Beaumont, C. (1996) Large-scale geomorphology: Classical concepts reconciled and integrated with contemporary ideas via a surface processes model, Journal of Geophysical Research, 101, 3361-3386 232 Koons, P.O. (1989) The topographic evolution of collisional mountain belts: a numerical look at the Southern Alps, New Zealand, American Journal of Science, 289, 1041-1069 Lane, E.W. and Carlson, E.J. (1953) Some factors affecting the stability of canals constructed in coarse granular materials, International Association for Hydraulic Research, Conference Proceedings, Minneapolis, Minnesota, September, 1953, 37-48 Leopold, L.B. and Emmett, W.W. (1972) Some rates of geomorphological processes, Geographia Polonica, 23, 27-35 Leopold, L.B. and Maddock, T. (1953) The hydraulic geometry of stream channels and some physiographic implications, United States Geological Survey Professional Paper, 252, 57 pp. Leopold, L.B. and Wolman, M.G. (1957) River channel patterns: braided, meandering and straight, United States Geological Survey Professional Paper, 282-B, 39-85 Lewis, L.A. (1974) Slow movement of earth under tropical rainforest conditions, Geology, 2, 910 Lewis, L.A. (1975) Slow slope movement in the dry tropics; La Paguera, Puerto Rico, Zeitschrift fur Geomorphologie, 19, 334-339 Lewis, L.A. (1976) Soil movement in the tropics Geomorphologie Supplementband, 25, 132-144 a general model, Zeitschrift fur Mark, D. and Church, M . (1977) On the misuse of regression in earth science, Mathematical Geology, 9, 63-75 Martin, Y . (1991) Sediment budget from morphology, M.Sc. thesis (unpublished), University of British Columbia, 184 pp. Martin, Y . and Church, M . (1995) Bed-material transport estimated from channel surveys: Vedder River, British Columbia, Earth Surface Processes and Landforms, 20, 347-363 Martin, Y . and Church, M . (1997) Diffusion in landscape development models: On the nature of basic transport relations, Earth Surface Processes and Landforms, 22, 273-279 233 McKean, J . A . et al. (1993) Quantification of soil production and downslope creep rates from cosmogemc Be accumulations on a hillslope profile, Geology, 21, 343-346 10 McLean, D. (1990) Channel instability of the lower Fraser River, Ph.D. thesis, University of British Columbia, Ph.D. thesis, 275 pp. Merritts, D. and Ellis, M . (1994) Introduction to special section of tectonics and topography, Journal of Geophysical Research, 99, 12135-12141 Miller, D.M. (1984) Reducing transformation bias in curve fitting, The American Statistician, 38, 124-126 Molnar, P. and England, P.C. (1990) Late Cenozoic uplift of mountain ranges and global climate change: chicken or egg? Nature, 346, 29-34 Nash, D.B. (1980a) Forms of bluffs degraded for different lengths of time in Emmet County, Michigan, U.S. A , Earth Surface Processes and Landforms, 5, 31-345 Nash, D.B. (1980b) Morphologic dating of degraded normal fault scarps, Journal of Geology, 88, 353-360 Neill, C R . (1973) Hydraulic geometry of sand rivers in Alberta, Proceedings of Hydrology Symposium, 9, Queen's Printer and Controller of Stationary, Ottawa, 453-461 Oden, M . (1994) Debris recharge rates in torrented gullies on the Queen Charlotte Islands M.Sc. thesis (unpublished), The University of British Columbia, 100 pp. Oreskes, N., Shrader-Frechette, K. and Belitz, K. (1994) Verification, validation, and confirmation of numerical models in the earth sciences, Science, 263, 641-646 Owens, LF. (1969) Causes and rates of soil creep in the Chilton Valley, Cass, New Zealand, Arctic and Alpine Research, 1, 213-220 Parker, G. (1979) Hydraulic geometry of active gravel rivers, Journal of the Hydraulics Division, Proceeding of the American Society of Civil Engineers, 105, 1185-1201 234 Penck, W. (1953) Morphological Analysis of Landforms (translated by H. Czech and K.C. Boswell), Macmillan, London, 429 pp. Personius, S., Kelsey, H. and Grabau, P. (1993) Evidence for regional stream aggradation in the Central Oregon Coast Range during the Pleistocene-Holocene Transition, Quaternary Research, 40, 297-308 Popper, K.R. (1959) The Logic of Scientific Discovery, BasicBooks, New York, 479 pp. Rice, S.P. and Church, M . (1995) Bed material texture in low order streams on the Queen Charlotte Islands, British Columbia, Earth Surface Processes and Landforms, 2 1 , 1-18 Rice, S.P. (1996) Bed material texture along gravel bed rivers with confluences, Ph.D. thesis (unpublished), The University of British Columbia, 238 pp. Rigon, R., Rinaldo, A . and Rodriguez-Iturbe, I. (1994) On landscape self-organization, Journal of Geophysical Research, 99, 11971-11993 Rinaldo, A . , Dietrich, W.E., Rigon, R., Vogel, G.K. and Rodriguez-Iturbe, I. Geomorphological signatures of varying climate, Nature, 374, 632-635 (1995) Rind, D. (1992) Ice ages: an uplifting experience, Nature, 360, 414-415 Ritter, P. (1987) A vector-based slope and aspect generation algorithm, Photogrammetric Engineering and Remote Sensing, 53, 1109-1 111 Roberts, R.G. (1984) Stream channel instability in the Queen Charlotte Islands: An examination of Schumm's fluvial model, M.Sc. thesis (unpublished), University of British Columbia, 198 pp. Rood, K. (1984) An aerial photography inventory of the frequency and yield of mass wasting on the Queen Charlotte Islands, British Columbia, B.C Ministry of Forests, Land Management Report, 34, 55 pp. Rood, K. (1990) Site characteristics and landsliding in forested and clearcut terrain, Queen Charlotte Islands, B.C., B.C Ministry of Forests, Land Management Report, 64, 46 pp. 235 Rosenbloom, N.A. and Anderson, R.S. (1994) Hillslope and channel evolution in a marine terraced landscape, Santa Cruz, California, Journal of Geophysical Research, 99, 1401314029 Saunders, I. And Young, A . (1983) Rates of surface processes on slopes, slope retreat and denudation, Earth Surface Proceses and Landforms, 8, 473-501 Scheidegger, A . D . (1970) Theoretical Geomorphology, Springer-Verlag, Berlin, 435 pp. Schumm, S.A. (1964) Seasonal variations of erosion rates and processes on hillslopes in western Colorado, Zeitschrift fur Geomorphologie Supplementband, 5, 215-238 Schumm, S.A. and Lichty, R.L. (1965) Time, space and causality in geomorphology, American Journal of Science, 263, 110-119 Schumm, S.A. (1977) The Fluvial System, John Wiley and Sons, New York, 338 pp. Seidl, M.A., Dietrich, W.E., Kirchner, J.W. (1994) Longitudinal profile development into bedrock: A n analysis of Hawaiian channels, Journal of Geology, 102, 457-474 Shaw, J., and Kellerhals, R. (1982) The composition of recent alluvial gravels in Alberta river beds, Alberta Research Council Bulletin, 41,151 pp. Shreve, R.L. (1979) Models for prediction in geomorphology, Journal of the International Association for Mathematical Geology, 11, 165-174 Simons, D.B. and Albertson, M.L. (1963) Uniform water conveyance channels in alluvial material, American Society of Civil Engineers, Transactions, 128, 65-167 Simpson, G.G. (1963) Historical Science, In C.C. Britton (ed.) The Fabric of Geology, Freeman, Cooper, Stanford, California, 24-48 Slaymaker, H.O. (1972) Patterns of present sub-aerial erosion and landforms in mid-Wales, Transactions of the Institute of British Geographers, 55, 47-68 236 Smart, G . M . (1984) Sediment transport formula for steep channels, Journal of Hydraulic Engineering, 110, 267-276 Smillie, G . M . and Koch, R.W. (1984) Bias in hydrologic prediction using log-log regression model, American Geophysical Union, Fall Meeting, San Francisco, December 3-7, Abstract in Eos, 65, 894 (No. H52B-04) Smith, T.R. and Bremerton, F.P. (1972) Stability and the conservation of mass in drainage basin evolution, Water Resources Research, 8, 1506-1529 Stebbings, J. (1963) The shapes of self-formed model alluvial channels, Institution of Civil Engineers Proceedings, 25, 485-510, Discussion, 26, 225-232 Strahler, A . H . and Strahler, A . N . (1992) Modern Physical Geography, John Wiley and Sons, New York, 638 pp. Summerfield, M.A. and Thomas, M.F. (1987) Long-term landform development: editorial introduction, In V . Gardiner, (ed.) International Geomorphology 1986 Part II, John Wiley and Sons, 927-933 Summerfield, M.A. (1991) Global Geomorphology, Longman Scientific and Technical, Essex, 537 pp. Tassone, B.L. (1987) Sediment loads from 1973 to 1984:08HB048 Carnation Creek at the mouth, British Columbia, In T.W. Chamberlin (ed.) Applying 15 years of Carnation Creek results, Proceedings of the workshop, January 13-15, 1987, Nanaimo, British Columbia, Carnation Creek Steering Committee, c/o Pacific Biological Station, Nanaimo, 46-58 Thomas, M.F and Summerfield, M.A. (1987) Long-term landform development: Key themes and research problems, In V . Gardiner (ed.) International Geomorphology 1986 Part II, John Wiley and Sons, 935-956 Tucker, G.E. and Slingerland, R.L. (1994) Erosional dynamics, flexural isostasy, and long-lived escarpments: A numerical modeling study, Journal of Geophysical Research, 99, 1222912243 237 Van Asch, Th.W.J., Demiel, M.S., Haak, W.J.C. and Simon, J. (1989) The viscous creep component in shallow clayey soil and the influence of tree load on creep rates, Earth Surface Processes and Landforms, 14, 557-564 Whiting, P. and Bradley, J. (1993) A process-based classification system for headwater streams, Earth Surface Processes and Landforms, 18, 603-612 Willgoose, G., Bras, R. and Rodriguez-Iturbe, I. (1991a) Results from a new model of river basin evolution, Earth Surface Processes and Landforms, 16, 237-254 Willgoose, G., Bras, R.L. and Rodriguez-Iturbe, I. (1991b) A coupled channel network growth and hillslope evolution model 1. Theory, Water Resources Research, 27, 1671-1684 Willgoose, G., Bras, R.L. and Rodriguez-Iturbe, I. (1991c) A coupled channel network growth and hillslope evolution model 2. Nondimensionalization and applications, Water Resources Research, 27, 1685-1696 Williams, G.P. (1970) Flume width and water depth effects in sediment transport experiments, United States Geological Survey Professional Paper, 562-H, 37 pp. Williams, M.A.J. (1973) The efficacy of creep and slopewash in tropical and temperate Australia, Australian Geographical Studies, 11, 62-78 Williams, G.P. and Wolman, M.G. (1984) Downstream effects of dams on alluvial rivers, United States Geological Survery Professional Paper, 1286, 83 pp. Wohl, E., Madsen, S. and MacDonald, L. (1997) Characteristics of log and clast bed-steps in step-pool streams of northwestern Montana, USA, Geomorphology, 20, 1-10 Wolman, M.G. and Brush, L.M. Jr. (1961) Factors controlling the size and shape of stream channels in coarse noncohesive sands, United States Geological Survey Professional Paper, 282-G, 183-210 Wolman, M.G. and Miller, J.P. (1960) Magnitude and frequency of forces in geomorphic processes, Journal of Geology, 68, 54-74 Young, A. (1960) Soil movement by denudational processes in slopes, Nature, 188, 120-122 238 Young, A. (1963) Soil movement on slopes, Nature, 200, 120-130 Young, A . (1974) The rate of slope retreat, Institute of British Geographers, Special Publication, 7, 65-78 Young, A . (1978) Slopes: 1970-75, In C. Embleton et al. (eds.) Geomorphology, Problems and Future Prospects, Oxford, Oxford University Press, 73-83 Present 239 APPENDIX I: REGRESSION ANALYSIS FOR NONLINEAR RELATIONS Miller (1984) documents the tendency to obtain a biased result when nonlinear relations are determined by transforming the data and then performing regression analysis. A reverse transformation is required to derive the best nonlinear equation for prediction of the dependent variable, in the sense that it gives the mean response. The regression model for data transformed using a base-10 logarithm is (Miller, 1984): logy - A + blogx + £ (1-1) where e is the prediction error of the model, A is the intercept and b is slope. The error term is assumed to be an independent, identically distributed random variate with zero mean (expected value=0) and variance &' . Transformation of the equation into its original form gives the e equation: y = a be x (1-2) z where z=2.303f for base-10 logarithms (the coefficient is eliminated for natural logarithms). However, the equation is unbiased if the expected value of e is equal to exactly 1. Smillie and z Koch (1984) found that for base-10 logarithms the expected value is: 5302a\ r E(e ) - exp z ^ (1-3) Because the above value is always positive, this means the predicted value of y is always biased negatively. The error variance required in this equation is evaluated using the following calculation: CT s = C7 \o y{\.-R \o% '\o%x) 2 1 2 i y (1-4) 240 APPENDIX II: ANALYSIS OF T H E BAGNOLD BED LOAD FORMULA II-A RELATION BETWEEN EXCESS STREAM POWER AND ADJUSTED SEDIMENT TRANSPORT An attempt was made to reconstruct the original Bagnold 1986 plot presented in figure 4.1 . In several cases the data points used by Bagnold were easily identified as they correspond to measurements presented in the data tables given in his earlier 1980 paper. In some cases wherein there were many measurements available for a particular data set, Bagnold included only one or several data points (e.g., Elbow River). Unfortunately, Bagnold does not discuss the basis for his selection of particular data sets and individual measurements used in his study. These points do not appear to be averages of the complete data sets and no rational reason for his particular choice of points could be deduced. A decision was made to reconstruct the Bagnold (1986) graph using a larger data set consisting of reliable bed load measurements. The bed load data compilation of Gomez and Church (1988) meets this requirement. Gomez and Church (1988) compiled 410 measurements made for various flumes and rivers, ranging in discharge from 10" to 10 m /s, which meet a 3 3 3 requirement of equilibrium transport conditions; i.e., flow and sediment conditions remain approximately constant throughout the measurement procedure. Under such conditions, sediment transport is supposed to be operating at capacity. Hence, there is expected to be a relation between flow parameters and transport rate. Such conditions are implicit in all computational attempts to estimate bed load transport (Gomez and Church, 1988). The data were obtained from measurements in gravel and coarse sand-bed channels (median grain size > 1 mm) in which no or minimal bedform development was observed. All river points were assumed to represent equilibrium transport conditions. In addition, Gomez and Church (1988) also included 240 points in their compilation that were obtained during the flume experiments 241 reported in their collection, but that did not meet the strict equilibrium requirements; these points are clearly flagged in their compilation. The reconstructed graph of ib* versus excess stream power is presented in figure II-1. Of the total (650) data points, 4 are eliminated as there was incomplete definition of key variables for the calculation of stream power and 169 are eliminated from the analysis as measured stream power is lower than calculated threshold stream power (although sediment transport actually occurred). This latter shortcoming applies primarily to the flume data, where 166 of a potential 516 points encountered this problem. The river data, in contrast, did not experience this difficulty to the same degree, with only 3 out of 130 data points being problematic. This shortcoming is particularly prevalent for the flume data, as opposed to the river data, because of the very low transport values measured in several of the studies (in particular those at the St. Anthony Falls Hydraulic Laboratory and Iowa Institute for Hydraulic Research). These data are eliminated from the analysis. Functional analysis is performed for three groups of data in order to evaluate the strength of transport/excess stream power relations and to determine if the slope of 1.5 found in Bagnold's bed load formula is preserved (see equation 4.6, p. 83) (table II-1): (i) river data (ii) flume data and (iii) all data points. The functional relation between the excess stream power and ib* is of greater concern when assessing the functional (and not predictive) relation that exists between these variables. The functional slope can be calculated from the equation (Mark and Church, 1977): - (b 2 b,=— f IR 2 + J(b 1R 2 2 - A) + AXb 2 2 — y (II-l) 2b r wherein b is the regression slope and X, which was calculated from the original data set, is the r ratio between the error variances for the x and y variables. The functional slope is calculated using the regression results between excess stream power and ib*. The slope for the river data is 242 Figure 11-1 Relation between stream power and i * for Bagnold's bed load transport formula. b Excess stream power (kg/m s) Table II-l Functional analysis of ib* vs. excess stream power. Data Set River Flume All data points # data points 127 350 477 R 0.94 0.86 0.90 J Slope 1.47 1.38 1.38 Confidence Limits of Slope 1.41-1.52 1.33-1.43 1.35-1.41 243 1.5 while the slope for the flume data is 1.4. A slope of 1.4 is obtained when river and flume data are analyzed together. The slope in each case is 1.5 (or very close to it), which indicates that the Bagnold correlation is functionally reasonable. The standard error of the estimate is about 0.48 log units, which indicates in unlogged units that model results can be expected to be between 0.3 and 3 times the estimate of ib* II-B BAGNOLD'S DEPTH AND GRAIN SIZE ADJUSTMENTS Attention is now focussed on the derivation of depth and grain size adjustments of the basic excess stream power/transport relation. To begin this examination, we must return to the form of the stream power relation presented in Bagnold (1977). This relation is similar to, but not identical with, the final form of Bagnold's equation (see equation 4.6, p. 83). Bagnold (1977) asserted that sediment transport is dependent on both excess stream power and relative roughness (d/D). In order to determine the nature of the relations, Bagnold assumed a best-fit line of 1.5 for the relation between threshold stream power and sediment transport for 3 river and 2 flume data sets. The relative roughness for each data set was calculated. From this information Bagnold extracted the relation: CO — COQ cc d_ (11-2) He then derived the relation between sediment transport and relative roughness by taking the inverse of this relation such that: h d_ 0 0 (II-3) In a subsequent 1980 paper, Bagnold explains that, in the above relation, the grain sizes (D50) of the various data sets were approximately equivalent and therefore grain size should not have been included in the resultant equation. He states that the correct result, instead, should be: 244 i oc d~ (II-4) in b Therefore, in this paper, Bagnold retains the d term but reassesses the relation between 2/3 sediment transport and grain size. Following a procedure similar to that followed to determine the depth relation, Bagnold (1980) adjusted several data sets, each having a different D 50 , to a standard depth. Excess stream power versus transport rate is plotted for each data set and it is concluded that: zr i oc b ,/2 (n-5) Before moving on it is necessary to first return to equations II-2 and II-3 in order to clear up an error in Bagnold's analysis. Close examination of Bagnold's 1977 paper shows that he made a simple mathematical error which, if corrected, changes his final result. The relation shown in equation II-2 still holds and the correct exponent is 2/3. The approach of taking the negative exponent (-2/3) (equation II-3) is correct only if the best-fit lines have a unit slope (tangent slope = 1). However, the relation between excess stream power and sediment transport in Bagnold's relation has a slope of 1.5. The correct mathematical manipulation is as follows. If the relation: i cc(a>-a> y n b o (11-6) is inserted into equation II-2 it transpires that: h d_ x (II-7) These results demonstrate that there is some confusion regarding the appropriate scaling for depth and grain size. It was decided to investigate this issue further. 245 II-C RE-EXAMINATION OF DEPTH AND GRAIN SIZE ADJUSTMENTS The relation between sediment transport and depth is first analyzed. An approach similar to that taken by Bagnold (1977; see previous section) was followed, once again using the large data set of Gomez and Church (1988). The data are separated into four groups of approximately equal grain size (1-2 mm, 2-5 mm, 5-8 mm and 22-32 mm). Sediment transport rates at a constant value of excess stream power for each data set in each grain size category are obtained by the following procedure. The average excess stream power and sediment transport rate are calculated for each data set. A best-fit line with a slope of 1.5 is assumed to hold in each case. The value of ib at constant excess stream power of 1.0 kg m'V (any other constant value would 1 suffice) for this relation is then calculated. The relation between this transport rate and depth is assessed by performing functional analyses (see figure II-2 and table II-2). In this case, it is reasonable to assign all of the error to the sediment transport term (depth and grain size can be measured with much greater accuracy), which is the dependent variable. This is the assumption underlying regression analysis, and hence a regression model can be fit to the data. The only individually significant relation is for the grain size category of 1 -2 mm with a slope of -1 ± 0.062. The outlier point in the graph represents the value for an equilibrium flume data set and as such it is preserved in the present analysis. When all of the data are analyzed together the slope of the line is -0.9 ± 0.2. On the basis of this analysis and the revised results of Bagnold, the best relation appears to be: i„ oc d~ (II-8) x Next the relation between sediment transport and grain size was analyzed. The data are separated into four groups of approximately equal depth (2-8 cm, 8-14 cm, 14-20 cm, 60-700 cm) (the large range of the last category is necessary in order to ensure there are enough data points for the regression analysis). Once again a best-fit line with a slope of 1.5 is fit to each data set. The value of ib at constant excess stream power of 1.0 kg m'V is calculated. 1 246 Figure 11-2 Relation between depth and i at constant excess stream power. The constant value of excess stream power used is 1.0 kg/m s. The outlier'( )' is the datum for the St. Anthony Falls Hydraulic Laboratory: 1 flume data. b 10 Grain s i z e O 1 - 2 mm A 2 - 5 mm E • 5 - 8 mm • 22 - 32 m m o Q. E CO CD GO GO 0.1 o X a; c ro -* —i GO o o ro 0.01- 0.001 0.01 0.1 10 1 Depth (m) T a b l e II-2 Sediment transport vs. depth. Grain Size Category 1-2 mm 2-5 mm 5-8 mm 22-32 mm All categories # data points 9 6 4 7 26 R 1 0.97 0.015 0.0036 0.019 0.47 Slope -0.97 -0.93 Standard Error of Slope 0.062 0.20 247 Regression analyses between this transport value and grain size are performed for each depth category (figure II-3 and table II-3). The relations show reasonable results with three of the depth categories having slopes of approximately -1.5. The depth category of 60-700 cm has a slope of about -0.5. When all of the data are analyzed together, a slope of -1.3 is obtained. On the basis of the calculated regression slopes and Bagnold's (1980) results, there appear to be two reasonable choices for the grain size/sediment transport relation: i„ oc Z T 1/2 (II-9a) i oc Z T 3/2 (II-9b) b II-D VARIOUS STREAM POWER CORRELATIONS In order to obtain the best stream power equation for use in the present channel submodel, relations between excess stream power and sediment transport (the latter is adjusted for various combinations of depth and grain size factors) were assessed by performing a functional analysis (table II-4) for a fixed slope of 1.5, which was shown to be reasonable in section II-A. Perhaps the most remarkable finding is that without any adjustment for depth and grain size the relation between excess stream power and sediment transport is already quite strong, having a best-fit functional relation with an R value of 0.74 (figure II-4). This immediately 2 suggests that we may be dealing with a simple scale relation in the system. The four combinations of depth and grain size adjustment factors are assessed. They all show similar R 2 values, ranging from 0.87 to 0.90. The relation for the adjustment factors of d and D ° has the 1 5 highest R value (figure II-5). Therefore, these are the adjustment factors chosen to include in 2 the bed load transport relation used in this study. According to the standard error of the estimate 248 Figure 11-3 Relation between grain size and ib at constant excess stream power. The constant value of excess stream power used is 1.0 kg/m s. 10 Depth Categories O 2-8cm • 8 - 1 4 cm A 1 4 - 2 0 cm • 6 0 - 7 0 0 cm ^4 GO E p 0 o Q. E co CD £ 0.1 CO CO CD o X CD -*-» c CO +-> CO c o o 0.01 co 0.001 0.0001 0.1 0.01 0.001 D (m) 50 Table II-3 Sediment transport vs. grain size. Grain Size Category 2-8 cm 8-14 cm 14-20 cm 60-700 cm All categories # data points 6 9 5 6 26 R 1 0.87 0.51 0.59 0.47 0.67 Slope -1.50 -1.29 -1.46 -0.52 -1.30 Standard Error of Slope 0.29 0.46 0.70 0.28 0.19 Figure 11-4 Relation between excess stream power and bed load transport. 250 Figure 11-5 Relation between stream power and i * for Bagnold-type formula derived in the present study. b Excess stream power (kg/m s) Table II-4 Functional analyses between excess stream power and bed load transport adjusted for various depth and grain size factors. A slope of 1.5 is assumed in all cases. Excess Stream power versus: i d V" 2y b i d v i dD2 b b _ Functional R 0.88 0.88 0.90 0.87 Confidence Limits of Slope 1.47-1.53 1.47-1.53 1.47-1.53 1.47-1.53 251 the c a l c u l a t e d transport rate c a n b e expected to v a r y ± 3 t i m e s the estimated v a l u e . When c o r r e c t i o n i s m a d e f o r the b i a s i n t r o d u c e d i n the d o u b l e - l o g transformations the f i n a l relation i s : i„ =0.005(co-co y d- D,2 o ] U2 (11-10) 252 APPENDIX III: LIST OF SYMBOLS List of symbols for equations used in the development of the present model and in the the final model (given in order of appearance in thesis). Symbols used in other models are given in text. h t X y k A S a P m Pb % j n V I Clj k 0 kj ho x hx scale h h ib h* n s b CO coo Pw Q d dref D D ref r Ys u w height time spatial dimension spatial dimension x component of transport rate y component of transport rate diffusion coefficient contributing drainage area slope diffusivity for slow, quasi-continuous mass movements diffusivity for rapid, episodic mass movements mass transfer rate sediment bulk density volumetric transport rate for j" gradient class gradient class number of landslides landslide volume transport distance for landslides basin area for a particular gradient class diffusivity parameter diffusivity parameter location on x-axis of midpoint of rise for nonlinear transport relation width parameter for rise in nonlinear transport relation height of the surface height of sediment/bedrock interface specific bed load transport rate bed load transport rate adjusted to a common flow depth and grain size stream power per unit length specific stream power critical specific stream power density of water discharge depth reference depth grain size reference grain size specific gravity of water specific gravity of sediment velocity width 1 253 g Q f Qbf E Do C7D Pi (j> s b b A, ma f r gravity mean annual discharge bankfull discharge average number of sub-channels in braided river characteristic particle size at downstream distance x=0 diminution coefficient Probability of a particular scour depth for debris flows river sink/source term prediction error functional slope regression slope ratio between error variances for x and y variables
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Modelling geomorphology in landscape evolution
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Modelling geomorphology in landscape evolution Martin, Yvonne. 1998
pdf
Page Metadata
Item Metadata
Title | Modelling geomorphology in landscape evolution |
Creator |
Martin, Yvonne. |
Date Issued | 1998 |
Description | Many landscape evolution models have considered the interaction of exogenic and endogenic processes. However, geomorphological processes have not been successfully incorporated in landscape evolution models. The thesis begins with a critical analysis of methodologies for the study of large-scale geomorphological processes. A framework based on a generalization of the relevant processes is recommended. Hillslope and channel submodels, which are based on typical processes operating in coastal regions of British Columbia, are introduced. The following hillslope processes are considered: (i) slow, quasi-continuous mass movements; (ii) fast, episodic mass movements; and (iii) weathering. The transport relation for fast, episodic mass movements was found to be nonlinear. Fluvial transport in both low and high-gradient channels and debris flow transport are considered in the channel submodel. A bed load transport equation, which is a revised version of the Bagnold stream power formula, is derived. Suspended load is calculated using a suspended load/contributing area correlation. Connections between hillslope and channel processes are considered to ensure adequate representation in the model. The hillslope and channel submodels are explored in one-dimensional and surface model runs for small drainage basins in the Queen Charlotte Islands, British Columbia. Tests of the fluvial submodel demonstrate the robustness of the bed load equation used in this study. A conceptualization of the landscape into unstable and stable regimes is introduced. Results of surface model runs emphasize the key role of low-order channels in transferring sediment from hillslopes to main channels. The exercise of constructing and running the model highlighted major gaps in our present understanding of geomorphological process operation and sediment routing. Suggestions for future research are extensive and are outlined in the concluding chapter of the thesis. |
Extent | 16809720 bytes |
Subject |
Geomorphology - Mathematical models |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-06-02 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0088766 |
URI | http://hdl.handle.net/2429/8596 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geography |
Affiliation |
Arts, Faculty of Geography, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1998-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1998-271986.pdf [ 16.03MB ]
- Metadata
- JSON: 831-1.0088766.json
- JSON-LD: 831-1.0088766-ld.json
- RDF/XML (Pretty): 831-1.0088766-rdf.xml
- RDF/JSON: 831-1.0088766-rdf.json
- Turtle: 831-1.0088766-turtle.txt
- N-Triples: 831-1.0088766-rdf-ntriples.txt
- Original Record: 831-1.0088766-source.json
- Full Text
- 831-1.0088766-fulltext.txt
- Citation
- 831-1.0088766.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0088766/manifest