Computer simulations of cobble structure on a gravel stream bed by Selina Tribe B.Sc, University of British Columbia, 1994 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E DEGREE OF Master of Science in T H E FACULTY OF GRADUATE STUDIES (Department of Geography) We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA September 1996 © Selina Tribe, 1996 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f t h e r e q u i r e m e n t s f o r a n a d v a n c e d d e g r e e at t h e U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e t h a t t h e L i b r a r y s h a l l m a k e it f r e e l y a v a i l a b l e f o r r e f e r e n c e a n d s t u d y . I f u r t h e r a g r e e t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s m a y b e g r a n t e d b y t h e h e a d o f m y d e p a r t m e n t o r b y h i s o r h e r r e p r e s e n t a t i v e s . It is u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l n o t b e a l l o w e d w i t h o u t m y w r i t t e n p e r m i s s i o n . D e p a r t m e n t T h e U n i v e r s i t y o f B r i t i s h C o l u m b i a V a n c o u v e r , C a n a d a DE-6 (2/88) Abstract In Harris Creek, a gravel bed river in southcentral British Columbia, the large stones protruding from the surface armour of the bed are arranged as curvilinear ridges oriented at varying angles to the flow, defining a net-like surface structure. A computer simulation, using various rules governing the entrainment and deposition of stones, was written to reproduce these structures. The working hypothesis for the simulation is that the patterns are the result of the mutual interaction of stones as they are mobilized into the flow and collide with downstream obstacles. An opposing hypothesis would be that the pattern being imposed on the stones by the flow. A statistical algorithm was developed to identify cobble ridges and their orientations and to allow a comparison between the simulation results and the real example. Simulations using different sized stones in a bed with 10% of its area covered by stones produce flow-parallel, linear clusters resembling cluster bedforms. Simulations using a more densely populated bed, up to 25% covered area, produce a variety of gravel bedforms including cluster bedforms, transverse ribs, some trending at oblique orientations, and arcuate ridges. These results indicate that the magnitude of the stone population influences the suite of bedforms that emerge. Cluster bedforms may be the result of obstacles to the downstream movement of stones in a sparsely populated bed. Furthermore, cluster bedforms may be nuclei from which oblique ridges develop out of excess sediment. Finally, simulations in which large stones are rarely entrained and travel shorter distances than small stones, produce structures very similar to those found in Harris Creek. ii Contents Abstract ii List of Tables v List of Figures vi Acknowledgements xi Chapter 1 Introduction 1 Chapter 2 The Cobble Patterns in Harris Creek 3 2.1 Description 3 2.2 Generating the map of Harris Creek 8 2.3 Comparison with Cluster Bedforms and Transverse Ribs 8 2.3.1 A Hydraulic Calculation to Test for Transverse Ribs 10 2.4 Digitizing Harris Creek 13 2.5 A Statistical Algorithm to Identify Ridge Structures 22 Chapter 3 Computer Simulations of Cobble Patterns 26 3.1 Facts about coarse sediment transport 27 3.2 Riverbed Computer Simulation Model 28 3.2.1 Initial and Boundary Conditions 30 3.2.2 Experimental Setup 33 3.2.3 Model Calibration and Testing 33 3.2.4 Programming Details 35 3.3 Simulations using equal sized stones 36 3.3.1 Discussion 37 3.4 Simulations using unequal sized stones 38 iii 3.4.1 Discussion 40 3.5 Simulations with reluctant large and mobile small stones 41 3.5.1 ' Discussion 42 Chapter 4 Conclusions 44 4.1 Future Work . . . 46 Appendix A Statistical Algorithm to Identify Ridge Structures 47 Appendix B Mathematica Script for Statistical Algorithm 53 Appendix C Computer Code for Riverbed 63 Bibliography 78 iv List of Tables 2.1 Statistics of the rose diagrams (figure 2.3) for the line segment approximation of the cobble patterns (figure 2.4) 6 2.2 Size distributions of identified stones in the digitized Harris Creek reach 20 2.3 Size distributions and area covered by stones in the digitized Harris Creek reach. Refer to figure 2.7 for location of boxes 21 v List of Figures 1.1 Vertical airphoto of the Harris Creek study reach showing cobbles protruding from the stream bed and defining a net-like structure. This photo shows the upstream third of section B as mapped in figure 2.1. The scale is visible as two perpendicular white 1.2 meter long sticks resting on some stones. A rope by which the suspended camera was towed along the reach is visible in the bottom right. . . 2 2.1 P lan view map of the cobble patterns in Harris Creek showing the outline of all protruding stones on the river bed. Section B continues approximately 5 m downstream of the end of section A . The three rectangular boxes outline representative examples of the cobbles patterns and are shown magnified in figure 2.2 4 2.2 Large scale view of protruding stones in three select areas of Harris Creek where the net-like pattern of cobble ridges is well displayed. The cobble patterns can be seen to incorporate linear, arcuate and ramifying ridges of stones as well as isolated stones, clusters and rib elements. A well-developed cluster bedform is evident in the center of bed at the downstream margin of Box 1 and several others are visible in the other boxes. A n example of a transverse rib can be found 3.5 m upstream of the cluster bedform in Box 1. Flow direction is to the right 5 2.3 Rose diagrams for each of the representative segments of Harris Creek's cobble patterns. The arcs represent the angular deviation in the population. Flow direction is toward 0° . The statistics for these roses are tabulated in Table 2.1 6 2.4 Line segment approximation of Harris Creek's cobble patterns drawn to the same scale as figure 2.1 and used to generate the data for the rose diagrams of figure 2.3. Each line represents the trend of a cobble ridge of 3 or more clasts in extent 7 2.5 Diagram of a single cluster bedform showing its three major components. The height of the obstacle clast is on order of 20 to 30 cm 9 vi 2.6 Sketch of transverse ribs. The side view includes a profile of an antidune wave illustrating that the increased water depth upstream of the antidune crests causes flow deceleration and the deposition of clasts which make up the ribs 10 2.7 Stones digitized as covering circles of five different radii. This figure is drawn to the same scale as figure 2.1. Careful comparison of the two figures will reveal a few places where the covering circle of a large stone completely overlaps a smaller adjacent stone. In these instances, computed areal coverage will be reduced slightly but stone population counts will be accurate since the smaller stone, though invisible, still remains in the digital file 14 2.8 Stones of radius less than 6.25 cm. These stones tend to cluster together into ridges that trace the cobble pattern 15 2.9 Stones of radius 6.25 cm to 12.5 cm. These stones tend to cluster together into ridges that in some places trace the cobble pattern 16 2.10 Stones of radius 12.5 cm to 18.75 cm. These stones are less clustered and more scattered throughout the reach than smaller stones 17 2.11 Stones of radius 18.75 to 25.0 cm. These stones appear randomly scattered throughout the reach and do not delineate the cobble patterns 18 2.12 Stones of radius 25.0 to 31.25 cm 19 2.13 Density plot (top) and PC2 graph (bottom) for Box 3 in Harris Creek along x. For the density plot, white pixels indicate regions with zero percentage of stones and black pixels indicate regions of high percentage of stones. The peaks in the PC2 graph correspond to angles at which the distribution of stones is least uniform. In the PC2 graph maxima appear at angles +15° and —7° and can be seen to be coincident with the angles at which there exist regions of alternating light and dark pixels in the density plot, that is, along vertical lines in the top graph. These angles correspond to ridges aligned at —75° and 83° in the Harris Creek map 23 2.14 Density plot (top) and PC2 graph (bottom) for Box 1 in Harris Creek along x. For the density plot, white pixels indicate regions with zero percentage of stones and black pixels indicate regions of high percentage of stones. The peaks in the PC2 graph at —13° and 19° correspond to ridge structures oriented at angles of 77° and —71° in the Harris Creek reach 24 2.15 Density plot (top) and PC2 graph (bottom) for Box 2 in Harris Creek along x. For the density plot, white pixels indicate regions with zero percentage of stones and black pixels indicate regions of high percentage of stones. The peaks in the PC2 graph and the angles at which there are regions of alternating light and dark pixels in the density plot occur at —15°, 10° and 16° corresponding to ridge structures oriented at angles of 75°, —80° and —74° in the Harris Creek reach 25 vii 3.1 Schematic of the rebound rule showing that an incoming stone whose rebound angle is determined to be greater than 90 degrees from the flow direction (left diagram) wi l l stop whereas a stone whose rebound angle is computed to be less than 90 degrees from the flow direction (right diagram) wi l l be deflected by the obstacle and continue on a new trajectory 29 3.2 Diagram showing how Riverbed calculates the impact position xnew, ynew and the rebound angle vtemp2 using knowledge of the stone's starting coordinates and trajectory, xtemp, ytemp and vtemp respectively, and the obstacles coordinates obx and oby. The equations used to calculate the impact position and rebound angle can be found in part 4 of the entrain() function in all versions of the Riverbed code. The grey filled stone is the obstacle, the black unfilled stone skims this obstacle and moves on as calculated by the black angles and vectors, and the grey unfilled stone stops behind the obstacle, as calculated with the grey angles and vectors. A l l angles and vectors have the same names as in the Riverbed code. Flow direction is along the positive x axis. The line segments labelled 10 are center-to-center lines between the obstacle and the impinging stone, both of radius 5.0. If different sized stones are used then this length must be adjusted accordingly 31 3.3 The init ial seeded beds for experiments with 10% (top), 20% (center), and 25% (bottom) coverage. The top bed has stones of equal radius and is the initial condition for experiments using Riverbeds 5.0 and 5.1. The center bed is the init ial condition for experiments using Riverbeds 6.1 and 6.2. The bottom bed is the init ial condition for an experiment using Riverbed 7.0 32 3.4 Diagram showing the river bed resulting from 20 randomly entrained stones and 20 newly injected stones in a short Riverbed test reach using the rebound rule on stones of equal size. Black circles represent the positions of stones at the initial conditions of the experiment. Entrained stones travel along vectors until they are deposited as a grey stone. Newly injected stones are transported along vectors originating on the upstream boundary (at x = 0) and are deposited in positions denoted by grey circles. Stones which travel out of the reach have vectors extending past the downstream boundary and are counted as throughput. Note that the upstream margin of the reach, at y = 350 -400, begins to clog up with stones. Even after only 40 entrainments the stones begin to cluster. . . 34 3.5 F ina l bed of Riverbed 5.0 using the rebound rule after 2000 random entrainments of 238 stones of radius 10 cm with F U L C R U M = 90 and S E E D F R A C T I O N = 0 . 1 . Flow is to the right 36 3.6 F ina l bed of Riverbed 5.1 using the rebound and neighbor rules after 2000 random entrainments of the 238 stones of radius 10 cm in the bed. F U L C R U M = 90 and S E E D F R A C T I O N = 0.1. Flow is to the right 36 vii i 3.7 Density plot (top) and PC2 graph (bottom) for the Riverbed 5.1 result along x. For the density plot, white pixels indicate regions with zero percentage of stones and black pixels indicate regions of high percentage of stones. The peaks in the PC2 graph correspond to angles at which the distribution of stones is least uniform. In the PC2 graph maxima appear at angles + 1 5 ° and —7° and can be seen to be coincident with the angles at which there exist regions of alternating light and dark pixels in the density plot, that is, along vertical lines in the top graph. These angles correspond to ridges aligned at —75° and 83° in the Harris Creek map 38 3.8 Final bed of Riverbed 6.1 using 477 unequal sized stones and the rebound rule after 3000 random entrainments. S E E D F R A C T I O N is 0.20 and F U L C R U M is set to 60. Flow is to the right 3 9 3.9 Final bed of Riverbed 6.2 using the obstacle and rebound rules on 477 unequal sized stones after 3000 random entrainments. S E E D F R A C T I O N is 0.20 and F U L C R U M is set to 60. Flow is to the right. Because of the primacy of the obstacle rule over the rebound rule, some overlap of stones occurs. . . 39 3 .10 Size distribution of stones in the initial bed (left), final bed (center) and of the stones that travel through the reach (right) for Figure 3.9. . 41 3.11 Final bed of Riverbed 7.0 using 596 unequal sized stones, the entrainment probability, travel distance and rebound rules after 3500 random entrainments. S E E D F R A C T I O N is 0.25 and F U L C R U M is set to 60. Flow is to the right 41 3.12 Density plot (top) and PC2 graph (bottom) along x for the Riverbed 7.0 simulation result in figure 3.11. The density plot exhibits alternate regions of light and dark pixels at an angle of —19° . The PC2 graph maxima appear at angles + 1 5 ° and —35° . These angles correspond to ridges aligned at 7 1 ° , 55° and - 7 5 ° in the simulated creek bed 42 3 .13 Final bed of Riverbed 7.0 using 238 unequal sized stones, the entrainment probability, travel distance and rebound rules, after 2500 random entrainments. S E E D F R A C T I O N is 0.10 and F U L C R U M is set to 60. Flow is to the right 43 A . l Schematic diagram illustrating how ridges oriented at different angles register in the x and y his-tograms. The transverse ridge on the right causes a maximum in the x histogram since, for this angle of 9, it is parallel to the bin strips along x. The longitudinal ridge on the left causes a peak in the y histogram. Flow direction is parallel to the box and in the positive x direction when the box is rotated to 0° 48 A . 2 A graph showing how the span of the x (top curve) and y (bottom curve) histograms change depending on the angle the bounding box is rotated to. The vertical axis is the length of the histogram in cm and the horizontal axis is 0 in degrees. The span of the x histogram is more constant than that of the y histogram 49 ix Given a clast that spans several bins in either the x or y dimension, it is necessary to determine which bin the leftmost (bottommost) point, the centre, and the rightmost (topmost) point are contained. The area of the stone in each bin is computed by the integrals Ix or Iy Acknowledgements Mike Church initially suggested this topic for study and the working hypothesis. He gave insightful comments, suggestions and encouragement throughout the course of this work. Birger Bergersen suggested the linked list approach and Rik Blok provided several tutorials on linked lists. Both Birger Bergersen and Brian Klinkenberg provided access to computational facilities. Jerry Maedel in the Forestry Department provided digitizing facilities and both he and Kurtis Halingten assisted in the digitizing. Kurtis Halingten assisted in developing the statistical algorithm, writing the Mathematica code, and solving numerous problems relating to conflicting file formats and non-compatible software. This work was done under the auspices of a University Graduate Fellowship. S E L I N A T R I B E The University of British Columbia September 1996 xi Chapter 1 Introduction When viewed from above, the cobbles protruding from the bed of Harris Creek are arranged into a net-like structure of curvilinear ridges (figure 1.1). How does this pattern form? The purpose of this study is to reproduce the cobble patterns in Harris Creek in a computer simulation. The design of the simulation is guided by the working hypothesis that the pattern is due to the mutual interaction of cobbles as they are entrained, transported and deposited, in contrast to the hypothesis that the pattern is imposed on the cobbles by the water flow. If the simulated structures resemble those in Harris Creek then inferences may be made about how the real structure is formed and possibly about coarse sediment transport in general. For those readers not versed in sediment transport, the following explanatory statements will help demystify the jargon used here. Although most of the terms cobble, gravel, clast and stone have strict size definitions in sedimentology, they are used interchangeably here. The size of a stone is denoted D50 if its diameter is that of the median size of the population. In a like manner, a stone denoted Dg5 has a diameter as large as the 95th percentile of stones in the population. The shape and orientation of irregular stones can be characterized as the 3 orthogonal axes, a, b and c of an approximating ellipsoid. The a axis defines the largest diameter of the stone and the b and c axes define the intermediate and smallest diameters respectively. The most hydrodynamically stable position for stones is to be oriented with the plane containing the ab axes dipping upstream because if they were positioned such that this plane was dipping downstream, the flow would catch the upturned end and flip the stone over. Stones in the stable orientation are called imbricated. The lee is the area immediately downstream of a stone while the stoss is the area immediately upstream of the stone. If a stone is positioned on the bed without any touching or near neighbours, it is in an open or isolated position. This is in contrast to stones touching or close to other stones which are said to be in closed or shielded positions. A stone is entrained at the precise moment it is mobilized into the flow and then it is transported along the river bed until finally it stops, hence becoming deposited. 1 Figure 1.1: Vertical airphoto of the Harris Creek study reach showing cobbles protruding from the stream bed and defining a net-like structure. This photo shows the upstream third of section B as mapped in figure 2.1. The scale is visible as two perpendicular white 1.2 meter long sticks resting on some stones. A rope by which the suspended camera was towed along the reach is visible in the bottom right. 2 Chapter 2 The Cobble Patterns in Harris Creek 2.1 Description The cobble patterns in Harris Creek comprise curvilinear ridges of stones protruding from the plane of the river bed defining a moderately well-developed, net-like structure (figure 2.1). The ridges are commonly 1 to 3 stones thick and 1 stone high and composed of the coarsest cobble and boulder fraction available in the creek. The ridges extend from 3 to 20 stone diameters and any one ridge can span approximately 25 to 75% of the uncovered width of the reach. The ridges are straight in some regions and in other regions are arcuate with the concave side facing either upstream or downstream. Despite the irregularity in their pattern, the ridges tend to form a net-like or rib-like structure whose segments are oriented oblique to the flow direction (figure 2.3) and spaced on order of 1 to 2 meters or 5 to 10 stone diameters apart. The cobble ridges are made up of sub-rounded to sub-angular cobbles and boulders of diverse litholo-gies with an average diameter of 20 centimeters (as measured from figure 2.1 and described in [9], riffle surface material) identifying them as the coarsest fraction of sediment available in the creek, comparable with the Dg5 of the bulk bed material. The cobble ridges have been observed to form during moderate floods capable of moving coarse stones and typical of snowmelt gravel bed streams. This arrangement of stones is also observed in Furry Creek, located in the southwestern Coast Mountains, British Columbia [10]. The study reach in which the cobble patterns are found is a 60 m long riffle portion of Harris Creek, a snowmelt flood regime creek draining a portion of the Okanagan highlands in southcentral British Columbia [9]. The reach is approximately 5 to 7 meters wide with an upstream boundary located approximately 6 m upstream of the first appearance of the structure and a downstream boundary located where the river narrows to half its average width resulting in a zone of convergent, super-critical flow. In the study reach, three domains of approximately 15 m long and 5 m wide exhibit representative examples of the cobble 3 Figure 2.1: Plan view map of the cobble patterns in Harris Creek showing the outline of all protruding stones on the river bed. Section B continues approximately 5 m downstream of the end of section A . The three rectangular boxes outline representative examples of the cobbles patterns and are shown magnified in figure 2.2. 4 Box 1 Box 2 Figure 2.2: Large scale view of protruding stones in three select areas of Harris Creek where the net-like pattern of cobble ridges is well displayed. The cobble patterns can be seen to incorporate linear, arcuate and ramifying ridges of stones as well as isolated stones, clusters and rib elements. A well-developed cluster bedform is evident in the center of bed at the downstream margin of Box 1 and several others are visible in the other boxes. An example of a transverse rib can be found 3.5 m upstream of the cluster bedform in Box 1. Flow direction is to the right. 5 flow direction Box 1 Box 2 Box 3 F i g u r e 2.3: Rose diagrams for each of the representative segments of Harris Creek's cobble patterns. The arcs represent the angular deviation in the population. Flow direction is toward 0 ° . The statistics for these roses are tabulated in Table 2.1. Tab le 2.1: Statistics of the rose diagrams (figure 2.3) for the line segment approximation of the cobble patterns (figure 2.4). statistic B o x 1 B o x 2 B o x 3 N a 37 39 24 class interval 5 ° 5 ° 5 ° m a x i m u m 6 5 5 3 m a x i m u m percentage"1 13.5 12.8 12.5 m o d e d 1 6 0 ° 1 6 5 ° 1 6 5 ° m e a n 6 3 4 2 . 7 ° 9 . 0 5 ° 1 5 . 2 5 ° angular deviation 3 7 . 2 5 ° 4 8 . 8 5 ° 4 5 . 8 ° "Total number of observations. ^Maximum number of observations in a single class. cMaximum percentage of total observations in a single class. dMode: class interval with maximum number of observations. eMean: azimuth of resultant vector. patterns and are outl ined on figures 2.1, 2.4 and 2.7. T h e verisimilitude of the computer simulations is judged by compar ing computed results with the areas within these three boxes. Harr i s Creek's cobble ridges are aligned at sub-transverse angles to the mean flow direct ion as shown by the bi-direct ional rose d iagrams 1 in figure 2.3. T h e associated statistics appear in Table 2.1. T o generate these diagrams, cobble ridges of three or more stones in extent were approximated by line segments. For each of the boxes outl ined in figures 2.1, 2.4 and 2.7, the angle that each line segment made wi th a line runn ing perpendicular to the box's long dimension, taken to be transverse to the mean flow direct ion, was tabulated. A r c u a t e ridges require several measurements. T h i s qualitative test for preferred orientation reveals that ridges tend to be aligned transverse rather than parallel to the flow. 1Rose diagrams made using Vector Rose 2.08, Orientation analysis and statistics for directional data, 1988-1992, software for the Macintosh written by P. Zippi, 1987. 6 Figure 2.4: Line segment approximation of Harris Creek's cobble patterns drawn to the same scale as figure 2.1 and used to generate the data for the rose diagrams of figure 2.3. Each line represents the trend of a cobble ridge of 3 or more clasts in extent. 7 2.2 Generating the map of Harris Creek In 1989, when the creek was at low flow, Harris Creek's cobble patterns were photographed with a gimbal-mounted camera suspended by helium-filled balloons. The camera and balloons were anchored and positioned by three people holding onto kevlar lines who slowly towed the apparatus above and along the reach. Every 5 to 10 meters the shutter was activated by remote control in order to obtain overlapping photographs for stereo viewing. The scale, in this case two 1.2 meter long sticks oriented perpendicular to each other, was laid on exposed rocks and repositioned every twenty meters or so. The resulting photographs provide an aerial, stereo view of the cobble structures on the river bed with a scale in almost every frame (figure 1.1). The plan-view map of the cobble structures was made by viewing each portion of the reach in stereo and outlining on a mylar overlay every clast observed to protrude from the river bed. Care was taken to avoid excess translation or rotation of each set of stereo photos used in tracing successive portions of the reach. The lateral boundaries of the study reach are the limit of overhanging trees for approximately one quarter of the length and the river bank for the remaining length. After the cobble outlines and reach boundaries were traced onto the mylar, the overlay was scanned and used as a digital template in Adobe Illustrator to produce figure 2.1. Cobbles protruding from the bed were identified by several criteria with varying degrees of ease and accuracy. Those that stuck out above the local water surface were the easiest to map, especially when leaves and other debris were trapped in their lee. Submerged cobbles often produced a disturbance of the water surface making them fairly easy to outline. Cobbles which were significantly submerged rarely disturbed the water surface; however, a careful observation of the colour, contrast and texture of the river bed through the water allowed clast outlines to be drawn. The submerged yet protruding stones were slightly more in focus with moderately more distinct outlines and a brighter colour than those stones locked into the basal armour of the creek bed. This outline can be thought of as representing roughly the plane defined by the oc or ab axes of the cobble. Some outlines may not be exact representations of the cobble sizes because of the variable orientation and uncertainties in the position and shape of the submerged clasts. In some case it is possible that submerged yet protruding clasts were not mapped because the water obscured them or that mapped clasts are not significantly proud of the bed and rightfully belong to the basal armour. 2.3 Comparison with Cluster Bedforms and Transverse Ribs Our understanding of the form and process of gravel bedforms lags behind our understanding of these matters in sandy sediment. This is because of the logistical problems inherent in dealing with gravel. For example, since coarse clasts range in diameter from centimeters to meters and require high shear stresses to entrain them, flumes capable of handling these conditions are expensive to make while field observations are difficult 8 to accomplish through flood-force, turbid water. Nonetheless, a suite of gravel bedforms has been established by various researchers over the years, usually by observations of outwash plains [23, 15] or by examining the bed configuration of a coarse sediment stream after a flood has passed [22, 14]. Harris Creek's cobble patterns resemble cluster bedforms and transverse ribs. The cellular aspect to the pattern is similar to the step-pool geometry of steep mountain streams [14], though at a smaller scale and in a stream with a lower gradient. ^ - a c c u m u l a l i o n - ^ o b s l a c l e ^ w a k e H plan view Figure 2.5: Diagram of a single cluster bedform showing its three major components. The height of the obstacle clast is on order of 20 to 30 cm. Cluster bedforms (figure 2.5), sometimes referred to as pebble clusters [12] and imbricate clusters [14], are perhaps the most common and best understood gravel bedform [12, 7, 27, 5, 13]. They consist of upstream-trending accumulations of coarse particles immediately upstream of a coarser particle. The cluster is typically several clast diameters long and one or two clast diameters wide. A single cluster bedform can be divided into three parts: the obstacle, the stoss accumulation and the wake accumulation [5]. The obstacle is usually the central and largest clast of the cluster behind which pile up clasts of comparable but usually smaller diameter. These piled-up clasts are often imbricated and are termed the stoss accumulation. Sand and gravel clasts of significantly smaller diameter than the obstacle are trapped in the lee of the obstacle and comprise the wake accumulation. Cluster bedforms are visible in the Harris Creek study reach. A good example can be found in the center of the bed at the downstream limit of Box 1 (figure 2.2) and less well developed examples elsewhere. A group of pebble clusters will define a pattern of points distributed in two dimensional space whereas Harris Creek's cobble patterns are more linear in two dimensional space. The linearity of the patterns and their preferred orientation sub-transverse to flow resembles transverse ribs. 9 flow direction side view Figure 2.6: Sketch of transverse ribs. The side view includes a profile of an antidune wave illustrating that the increased water depth upstream of the antidune crests causes flow deceleration and the deposition of clasts which make up the ribs. Transverse ribs (figure 2.6) are rows of gravel clasts typically 1 to 2 clasts thick, several clasts wide, oriented transverse to the mean flow direction and semi-regularly spaced along the stream at distances that vary proportionally with the diameter of the constituent clasts [15, 2]. Commonly ribs are roughly linear, though some are concave facing upstream [22], and the constituent clasts are imbricated. On gravel outwash plains, transverse ribs are observed to grade into a cellular arrangement of clasts usually 1 to 2 clasts thick with a diameter on order of 1 m called stone cells [15]. The stones constituting these cells are also imbricated. In Harris Creek, an example of a transverse rib can be seen 4.5 m upstream of Box l's downstream boundary (figure 2.2). The strong sub-transverse alignment of the cobble ridges resembles transverse ribs and the net-like, arcuate form traced by some the ridges resembles the description of stone cells [15]. Harris Creek's cobble patterns may be considered as the large scale structure defined by a continuum of gravel bedforms; an admixture of transverse ribs, stone cells, cluster bedforms and isolated stones. Moreover, cluster bedforms may be the articulation points of ridges. 2.3.1 A Hydraulic Calculation to Test for Transverse Ribs Harris Creek's cobble patterns more strongly resemble transverse ribs than cluster bedforms, prompting a closer investigation of their similarity. It is a widely accepted hypothesis that transverse ribs form in almost critical to supercritical flow regimes. That they form in supercritical flows is supported by the fact that ribs form preferentially on gravel bars or riffle sections in gravel bed streams, braided streams and outwash plains [23, 15] where a shallower water depth causes high Froude number flows. Also, flume experiments in gravel 10 show that clasts tend to accumulate slightly upstream of a hydraulic jump [24]. That ribs form in almost critical flows is supported by observations that floods rarely became supercritical in regions where transverse ribs were known to form. Koster [19] suggests that the increase in flow depth at an antidune crest causes a decrease in flow velocity and subsequent deposition of coarse particles underneath the crest (figure 2.6). Rib formation during almost critical flows is also supported by observations of rooster tail and checkerboard antidunes, whose pattern of crests closely resembles the linear pattern of ribs and the more circular pattern of stone cells respectively, in regions where ribs and cells are known to form [15]. Assuming ribs require flows with a Froude number approaching or exceeding unity, if calculations show that Harris Creek has near critical flows then its cobble patterns could be transverse ribs. Alternatively, if the calculated Froude numbers are well below 1.0 this can be interpreted to mean that either the cobble patterns are distinct from transverse ribs or that transverse ribs have a formative mechanism that does not require critical flows. The latter interpretation would suggest that the widely held hypothesis for the origin of transverse ribs is in error. The starting point for the calculation is the critical Shields' number, 6, which is related to the value of shear stress T felt by a particle of diameter D at the moment of entrainment, where ps = 2.7g/cm3 and p — 1.Op/cm 3 are the densities of sediment and water respectively, and g = 9.8 m/s2 is gravitational acceleration. This formula was developed for well-sorted natural gravels and for individual particles of varying sizes occupying open positions on a flume bed. Shields found that as the particle diameter increased, 6 attained a constant value of 0.06. Later investigators found that 6 can be as low as 0.03 for large, isolated stones and as high as 0.09 for shielded stones [10]. After obtaining a characteristic clast diameter for Harris Creek, Equation (3.1) can be used to determine the shear stress r of the flow required to entrain that clast. The value of shear stress resulting from assuming the Shields relation holds for a particular particle size can then be equated to the shear stress exerted by water of average depth d flowing uniformly down a bed of slope 5, r = pgdS (3.2) Given the slope S of the reach, equation (3.2) will give the depth d of the flow at the moment of entrainment. Using discharge-depth and discharge-velocity relations obtained from gauging station records near the study reach [30], the depth can be used to find the discharge of the flow which in turn is used to find the velocity v. The final step is to determine the Froude number, 11 where crit ical flows have Fr = 1.0. T h e above calculation was performed for a bed slope S of 0.013 and D50 and DS4 of bed material on the riffle surface of Harr is Creek, approximately 75 m m and 110 m m respectively [9]. T h e result ing Froude numbers range from Fr = 0.37 for small stones in open positions, that is for D50 and 6 = 0.03, to Fr — 0.69 for large stones in shielded positions, for DM and 9 = 0.09. T h i s exercise indicates that even for large, locked clasts such as those compris ing Harr is Creek's cobble patterns, the average Froude number of the flow was well below 1.0 and therefore sub-crit ical , imply ing the patterns were formed dur ing firmly sub-crit ical flows. T h i s impl icat ion is supported by observations that the cobble patterns in Harr i s Creek were destroyed dur ing peak flows a nd reformed after floods of lesser magnitude [10]. T h e calculat ion was also performed for the largest stones found in Harr is Creek, r>95 of 20 cm. T h e resulting Froude number ranges from 0.57 for stones in open positions to 0.94 for stones in shielded positions. T h i s indicates the largest stones require sub-crit ical to crit ical flows to mobilize, especially if they are shielded by neighbors. A s wi l l be discussed later, the largest stones may not need to be entrained to form the cobble patterns since they appear to be scattered about the bed whereas the smaller stones are clustered together, defining the structure. It is important to keep in m i n d that because turbulent flows cause widely fluctuating shear stresses over the bed, this simple calculation may not embody the dynamics of Harr is Creek's actual flows. T h e ridges compris ing Harr i s Creek's patterns have a sub-transverse arrangement resembling trans-verse ribs. T h a t they are distinct from transverse ribs is indicated by the oblique orientation of the cobble ridges to the mean flow direction (figures 2.3 and 2.4) and by calculations showing that flows were sub-crit ical dur ing the formation of the ridges. 12 2.4 Digitizing Harris Creek T h e Harr i s Creek study reach was digitized to determine characteristic values for parameters used in the s imulat ion and to provide data for the statistical a lgorithm, which would allow a more quantitat ive com-parison between the real phenomenon and the simulation. Dig i t iz ing also allowed the stones in the reach to be displayed according to their size, revealing dissimilar distributions of large and small stones that are incorporated into the simulation. O n an enlargement of figure 2.1 each stone's largest diameter was measured by h a n d and the stone assigned to one of five size bins on the basis of its radius. T h e stone's largest diameter corresponds roughly to the a axis, a point that can be appreciated by considering that an ellipsoidal stone wil l most likely be deposited ly ing on its side, so to speak, in the stream bed and not posit ioned end up, balanced on its smallest cross-section. T h e bins used are 0-6 .25, 6.25-12.5, 12.5-18.75, 18.75-25.0 and 25:0-31.25 c m . E a c h stone is digit ized as a circle of radius r equal to the m a x i m u m allowable radius of the b in the stone belongs to, i.e., a stone of radius 11 c m would fall into the second bin and be digitized as a circle of radius 12.5 cm. E a c h stone is writ ten to a datafile as the x and y coordinates of the center point and radius r of the covering circle. T h e lateral boundaries of the reach were digitized as a sequence of points (x,y) defining line segments in space. T h e digit ized reach is shown in figure 2.7 and separate maps showing the location of each stone in a size b in are given in figures 2.8 to 2.12. In practice, the digitizing software writes a graphic file which must be edited in order to get da ta files of (x, y, r) for each stone, one stone to a line, and of (x, y) for the points defining the lateral boundaries. Instead of using covering circles, the precise outline of each stone could have been digitized along with a point representing the center of each stone. However, the tedious task of out l ining 2500 irregular shapes was passed over in favor of the simpler method. T h e binned stones' radi i (Table 2.2) are approximately lognormally distr ibuted and correspond well wi th the size distr ibut ion measured on the riffle surface of a reach not far from the study section [9]. T h e stones in the two smallest size classes cluster together and exhibit the cobble patterns fairly well (compare figures 2.8 and 2.9 wi th figure 2.1 and 2.7). T h e stones in the three larger size classes do not delineate the cobble patterns, rather, they appear to be scattered throughout the reach (compare figures 2.10, 2.11 and 2.12 with figure 2.1 and 2.7). Separately displaying stones according to radius suggests that the larger stones are the keystones of the pattern while it is the accumulat ion of smaller stones that determines the pattern. 13 F i g u r e 2.7: Stones digitized as covering circles of five different radii. This figure is drawn to the same scale as figure 2.1. Careful comparison of the two figures will reveal a few places where the covering circle of a large stone completely overlaps a smaller adjacent stone. In these instances, computed areal coverage will be reduced slightly but stone population counts will be accurate since the smaller stone, though invisible, still remains in the digital file. 14 15 Figure 2.9: Stones of radius 6.25 cm to 12.5 cm. These stones tend to cluster together into ridges that in some places trace the cobble pattern. 16 Figure 2.10: Stones of radius 12.5 cm to 18.75 cm. These stones are less clustered and more scattered throughout the reach than smaller stones. 17 Figure 2.11: Stones of radius 18.75 to 25.0 cm. These stones appear randomly scattered throughout the reach and do not delineate the cobble patterns. 18 Figure 2.12: Stones of radius 25.0 to 31.25 cm. 19 bin size Section A Section B (cm) frequency frequency 0 - 6.25 663 713 6.25 - 12.5 438 461 12.5 - 18.75 100 74 18.75 - 25.0 16 11 25.0 - 31.25 6 0 section total 1223 1259 total clasts 2482 Table 2.2: Size distributions of identified stones in the digitized Harris Creek reach. T h e projected cross sectional area of stones in each size b in was computed from the digit izing information and is tabulated in Table 2.3. Note that these values are m a x i m a since covering circles were used in the calculations. T a k e n together, the stones cover 23% to 30% of the bed. Exper iments in sand showed that particles covering about 6% to 12% of the bed maximized the flow resistance [34]. A densely covered bed facilitates stabil ity since the flow steps above the wake interference of the rough bed, no longer exerting force close to the bed to move stones [25]. T h e form of the bed, then, can be considered a stable configuration since the flow is displaced elsewhere. Harr is Creek's dense coverage suggests the cobble patterns are an expression of the most stable configuration of the bed to shear flow, a s ituation supported by observations that the patterns are destroyed when the stresses in Harris Creek exceed about three times the threshold stress required to move the larger half of clasts [10]. Large clasts posit ioned as keystones of the patterns, the net-like or cellular arrangement of ridges, and the high populat ion of stones contr ibut ing efficient energy dissipation [25], all attest to the stability of the cobble patterns. 20 B o x 1 dimensions: 16.0 m by 5.6 m radius frequency absolute relative abs. covered (cm) area ( m 2 ) area (%) area (%) j 6.25 301 3.69 18.12 4.12 6.25 - 12.5 204 10.01 49.17 11.17 12.5 - 18.75 45 4.97 24.41 5.55 18.75 - 25.0 7 1.37 6.73 1.53 25.0 - 31.25 1 0.31 1.50 0.35 total 558 20.36 22.72 B o x 2 dimensions: 15.6 m by 4.5 m radius frequency absolute relative abs. covered (cm) area ( m 2 ) area (%) area (%) 0 - 6.25 329 4.04 18.85 5.75 6.25 - 12.5 256 12.57 58.66 17.89 12.5 - 18.75 33 3.64 16.99 5.18 18.75 - 25.0 6 1.18 5.51 1.68 25.0 - 31.25 0 0 0 0 total 624 21.43 30.50 B o x 3 dimensions: 16.5 m by 4.5 m radius frequency absolute relative abs. covered (cm) area ( m 2 ) area (%) area (%) 0 - 6.25 353 4.33 19.74 5.83 6.25 - 12.5 235 11.54 52.60 15.54 12.5 - 18.75 39 4.31 19.64 5.80 18.75 - 25.0 9 1.77 8.07 2.38 25.0 - 31.25 0 0 0 0 total 636 21.94 29.55 Table 2.3: Size distributions and area covered by stones in the digitized Harris Creek reach. Refer to figure 2.7 for location of boxes. 21 2.5 A Statistical Algorithm to Identify Ridge Structures In order to compare the real phenomena as exemplified in Harr is Creek with the computer s imulat ion results, a statistical a lgor i thm was developed that identifies ridge structures and their orientations. A close examinat ion of figure 2.2 reveals a pattern of cobble ridges which occur in and are s tructural ly associated with other structures such as pebble clusters, transverse ribs and isolated stones. W h a t exactly is a cobble ridge a n d how can it be measured? A cobble ridge, the basic constituent of the cobble patterns, is defined here as a linear or arcuate alignment of 3 or more stones. T h e interclast spacing, the (non)linearity of the ridge, and the sizes of the component stones are several factors that complicate the definition of a cobble ridge pattern. O t h e r compl icat ing factors include the proximity of other patterns (clusters, river bed boundaries , isolated stones) an d the angle of the ridge with respect to the flow. A statistical a lgor i thm was developed that would register linear accumulations of stones trending at any angle to the mean flow direction by measuring the percent of covered area in strips of different orientation (see A p p e n d i x A for a complete description of the statistical test). T h e rectangular r iverbed is rotated from — 4 5 ° to 4 5 ° about its lower left corner and at each 1 degree interval the fraction of the bed covered by stones is calculated for vertical and horizontal strips and displayed as a histogram projected onto the reference x and y axes. D u r i n g the rotat ion of the river bed, an obliquely trending ridge will show up on the his togram as a peak when it is rotated into parallel ism with the strips along x or y. A t all other angles it wil l contribute to a more uni form distr ibution. T h e areal percentage of the strip covered by stones at each angle is displayed as one pixel in a density plot of angle vs b in where a white pixel indicates zero covered area and a black pixel indicates a high percentage of covered area. T h e information in the density plot is condensed as a pseudo \ 2 (PC2) graph, so called because it uses percent covered area in the calculation instead of the absolute values needed for a true x2 graph. Peaks in the PC2 graph indicate the angle at which the d is tr ibut ion of stones is least uni form. In the Harr is Creek maps and in s imulation results, the downstream direct ion is defined as 0 ° a n d clockwise angles are positive, however, in this statistical a lgor i thm clockwise angles are negative. T h i s difference in convention necessitates that a peak occurring at some angle in the statistical graphs must be converted to find the corresponding angle of the ridge in the study reach. In the density plot, if there are alternating light and dark regions along an angle, it is likely there are alternating low an d high concentrations of stones at corresponding angles in the river bed. A candidate ridge structure is identified as a light area (low concentration of stones in the strip) changing to a dark area (high concentration of stones in the strip) and then back to a light area along a horizontal line. T h i s records the rotat ion of a ridge at approximately that b in location into and out of paral lel ism wi th the b i n strip. T h e density plot along x and the PC2 test for the cobble structures in B o x 3 are shown in figure 2.13. T h e density plot is well banded with dark and light strips. T h e PC2 graph has peaks at 1 5 ° a nd —7° corresponding to ribs oriented at — 7 5 ° and 8 3 ° to the flow direction. T h e other two boxes were also 22 analyzed with this statistical method and the results are shown in figures 2.14 and 2.15. O n l y the density plot and PC2 graphs along x are shown because the widths of the strips spanning the x axis change much less than they do along y axis therefore contr ibut ing less error to the results. Also , the test was written so that ridges with a strong transverse component wil l register in the density plot along x. u 4 a 2 -40 -20 0 20 40 t h e t a (degrees) Figure 2.13: Density plot (top) and PC2 graph (bottom) for Box 3 in Harris Creek along x. For the density plot, white pixels indicate regions with zero percentage of stones and black pixels indicate regions of high percentage of stones. The peaks in the PC2 graph correspond to angles at which the distribution of stones is least uniform. In the PC2 graph maxima appear at angles + 1 5 ° and —7° and can be seen to be coincident with the angles at which there exist regions of alternating light and dark pixels in the density plot, that is, along vertical lines in the top graph. These angles correspond to ridges aligned at —75° and 83° in the Harris Creek map. T h e rose diagrams (figure 2.3) show that there is a predominant alignment approximate ly 60 to 7 7 ° either side of the flow direction. T h e PC2 graphs for the three boxes show comparable m a x i m a corresponding to approximate ly 72 to 8 6 ° . T h i s indicates that the statistical test is detecting the predominant orientations of the ridges. 23 . , , . . • , • • • , . • • . , --45 -25 -5 15 35 theta (degrees) -40 .: II O 20 40 t h e t a (dooraea) Figure 2.14: Density plot (top) and PC2 graph (bottom) for Box 1 in Harris Creek along x. For the density plot, white pixels indicate regions with zero percentage of stones and black pixels indicate regions of high percentage of stones. The peaks in the PC2 graph at —13° and 19° correspond to ridge structures oriented at angles of 77° and —71° in the Harris Creek reach. 24 300 y 0 -45 -25 -5 15 35 theta (degrees) Figure 2.15: Density plot (top) and PC2 graph (bottom) for Box 2 in Harris Creek along x. For the density plot, white pixels indicate regions with zero percentage of stones and black pixels indicate regions of high percentage of stones. The peaks in the PC2 graph and the angles at which there are regions of alternating light and dark pixels in the density plot occur at —15° , 10° and 16° corresponding to ridge structures oriented at angles of 7 5 ° , —80° and —74° in the Harris Creek reach. 25 Chapter 3 Computer Stimulations of Cobble Patterns T h e rise in popular i ty of computer s imulation models in the physical sciences over the past ten years or so is due to the sophistication of computer model ing techniques, in part icular the development of cellular automata , lattice models and new numerical techniques. E a r t h scientists have long noted the regular a nd distinctive patterns that can develop in earth materials yet only now do they have the computat ional tools wi th which to study such phenomena [16]. Several investigators have been successful in reproducing ge-omorphic patterns such as eolian ripples [3, 20, 31], beach cusps [32], non-periglacial sorted nets [1] and stone stripes [33] using computer s imulation models. These simulations all share in c o m m o n the desire to reproduce a natural ly occurring pattern developed in noncohesive sediments of various sizes and to retain the discrete character of the sand or pebbles compris ing the pattern. In contrast to using differential equations dealing with continuous variables [20], the discrete particles are programmed to move or behave according to the investigators' understanding of the local physics of the system and usually involve randomness in the a lgor i thm. As ide from being currently in vogue, computer s imulation experiments allow certain phenomena to be studied without the logistical problems accompanying physical experimentation. In this case, cob-ble transport is s imulated without the problems of clast retrieval and equipment damage that so hamper empir ica l investigation. It is essential to realize that a computer s imulation, by virtue of its discrete nature, never completely reproduces real conditions. A successful computer simulation of a river, like the one presented here, only validates the hydrologic and sedimentologic scenario programmed into it. It is also crucial to remember that results of the R i v e r b e d model are a function of the physical scenario programmed into the code, in the form of parameter values an d rules, as well as the problem solving architecture, in the form of loops a nd logical 26 criteria. B o t h factors must be kept in m i n d as possibly contributing to s imulat ion results. Before describing the computer model in detail , reviewing the current knowledge about coarse sedi-ment transport will sharpen our intuitive understanding of the physics of the system, outline the hydrologic a n d sedimentologic scenario used in the s imulation, and uncover several principles that wil l be incorporated into the design of the model . 3.1 Facts about coarse sediment transport In a gravel bed river the size of the stones is commensurate with the depth of water flow result ing in very rough and irregular flow boundary conditions. T h e stones are often observed to be imbricated against one another and packed such that some are shielded from the shear effects of the flow by virtue of being h idden by neighbors [21]. In contrast to a stone ly ing in an open posit ion on the bed, shielded stones of equal size generally require higher shear stresses to be entrained, typical ly because they can't move unt i l the obstacle, which is usually larger, moves. Tradit ional ly , the shear stress acting parallel to the bed has been characterized as a function of the water velocity and grain size of the particle, an approach which works reasonably well for sand-sized sediment. However, on a packed gravel bed the shear stress acting on each stone must also be a function of the sheltering and stability afforded by surrounding stones (as al luded to in Section 2.3.1). A l l this is to say that the mobi l i ty of coarse clasts is a complex function of flow conditions, particle size and bed structure [21]. T h e difficulty in predict ing the onset of movement lends itself to a stochastic model ing approach and says that the movement of stones in the model should be part ia l ly dependent on whether the stone is in an open or shielded posit ion. G i v e n that a cobble has been successfully entrained into the flow, how does it travel a long the bed? Because they are so large, cobbles in transport will rol l and slide along the b o t t o m of the bed. Consequently they are transported along a b u m p y and irregular bed surface and encounter obstacles in their p a t h resulting in an irregular transport trajectory. Coarse clasts are known to engage in sporadic and discontinuous mot ion , a situation which again supports the appropriateness of a stochastic model ing approach and does not contradict the irregular trajectories that the s imulated stones travel . A stone will tend to stop if it hits an obstacle from behind, causing an elastic rebound in the upstream direction which is quickly d a m p e d by the flow downstream. If a stone hits an obstacle on the side then the combined effects of a rebound trajectory not diametrical ly opposed to the flow and the force of the water pushing the stones downstream wil l tend to keep the stone in mot ion. Large stones travel l ing downstream, by virtue of their momentum, may tend to roll over small protruberances in the stream b e d and only be stopped by impinging against a stone of comparable radius. Smal l stones travell ing downstream m a y tend to stop behind all obstacles or be drawn into the obstacle's lee by virtue of being small enough to be affected 27 by the reduced lift and drag forces caused by the obstacle's perturbat ion of the flow [7]. In the s imulat ion, the deposit ion of a stone can depend on the size of the obstacle and the angle at which it hits the obstacle. 3.2 Riverbed Computer Simulation Model T h e pr ime motivat ion in writ ing the R i v e r b e d model is to reproduce Harr is Creek's collection of obliquely and variably trending cobble ridges. T h e adopted rules reflect this goal. R i v e r b e d is writ ten in C++ and version 7 . 0 can be found in A p p e n d i x C . T h e code is heavily documented, especially at the beginning of functions and m a j o r loops. In the following description of the algorithm, P A R A M E T E R S which can be changed from one experiment to the next are indicated by small capitals and the various rules implemented are indicated by italics. Different versions of the code were written to incorporate different rules and parameters . For instance, R i v e r b e d 5 . 0 uses the rebound rule wi th clasts of equal radius while R i v e r b e d 6 . 2 uses rebound and obstacle rules with clasts of different radi i . For ease of exposition, only one result f rom any part icular set of rules and parameters is shown and is often referred to by the version of the code that produced it. In actuality, several experiments were performed for each of many combinations of rules and parameters. T h e result shown is representative of the results of these experiments. In the R i v e r b e d computer model , indiv idual cobbles in the reach are successively picked at r a n d o m , or wi th some nonuniform probability, as in the entrainment probabi l i ty rule, and entrained along a randomly chosen trajectory. U p o n impact ing an obstacle, the cobble stops or continues depending on the cobble-obstacle geometry and the rules governing deposition. If a stone hits an obstacle right f rom behind then it wil l tend to stop but if it hits at an angle, the stone could bounce off and continue downstream in a different direction. After rebounding off the obstacle the cobble continues on its new p a t h to further encounters or loses its m o m e n t u m and becomes deposited, as in the travel distance rule, otherwise it remains at the posit ion of impact . T h e rebound rule is modeled after elastic collisions (figure 3.1). G i v e n the stone's x, y coordinates, its trajectory and the coordinates of the obstacle, the rebound angle is calculated to be equal to the impact angle as measured about the normal to the tangent of the impact site (see figure 3.2 for reference d iagram). If the rebound trajectory is within 9 0 ° either side of the mean flow direction, the cobble wil l continue down the reach on this new path , in other words it bounces off the obstacle. If the rebound trajectory is greater than ± 9 0 ° from the flow, the cobble stops in the posit ion at impact . T h e rationale beh ind deposit ing those cobbles wi th a rebound trajectory greater than ± 9 0 ° is that otherwise the cobble would travel upstream. In this way the rebound rule models the general effects of the water flow, a l though hydrodynamics are not explicit ly treated. T h e rebound rule may be slightly modified from one experiment to the next. T h e simplest way to change it is to change the angular range permit t ing rebounding by specifying the parameter F U L C R U M 28 flow direction y Figure 3.1: Schematic of the rebound rule showing that an incoming stone whose rebound angle is determined to be greater than 90 degrees from the flow direction (left diagram) will stop whereas a stone whose rebound angle is computed to be less than 90 degrees from the flow direction (right diagram) will be deflected by the obstacle and continue on a new trajectory. to be other than 90. T h e rebound rule is used in all versions of R i v e r b e d . T h e neighbor rule is intended to model the shielding effects of neighboring stones: after a stone has been chosen to be randomly entrained, the neighbor rule counts the number of stones touching or within some smal l distance, typical ly 1 to 3 c m , of the chosen stone. If the number of neighbors, regardless of their spatial d is tr ibut ion, is greater than or equal to the crit ical value NN then this stone is not entrained and another candidate is chosen. If the neighbors number less than NN, the stone is entrained. T h e neighbor rule is used in R i v e r b e d 5 . 1 . T h e obstacle rule models the effects of the relative size of impinging and obstacle clasts: the obstacle a stone is deposited against must have a radius greater than or equal to half the stone's radius, otherwise the stone rolls over it. T h e rationale for this rule is that a large stone in transport will tend to be stopped by obstacles of comparable or larger radius while small stones wil l be stopped behind anyth ing they r u n into. T h e obstacle rule is used in R i v e r b e d 6 .2 . W h e n locating an appropriate obstacle, this rule takes precedence over the rebound rule. A n entrainment probability rule models the fact that larger stones are less likely to be entrained into the flow than are smaller stones at any given shear stress. For any entrained stone, the entrainment probabi l i ty is calculated as the stone's radius mult ipl ied by a uniform r a n d o m number. If the result is less t h a n some CRITERION the stone is entrained, otherwise it is not entrained. T h e entrainment probabi l i ty is p r o g r a m m e d in R i v e r b e d 7 .0 . T h e travel distance rule embodies the fact that a travelling stone wil l slow down and stop due to its inert ia and friction with the bed and water. T h i s rule is model led after the results of experiments on the travel distance of stones in open positions in gravel bed rivers [8] where a stone's travel distance was found 2 9 to vary in inverse proportion to its size as L = 1.77(1- log D)135 (2.1) where L = L i / ' i D 5 o , u r / o c e i D = Di/D50subsurfa.ee and Li and Di are respectively the travel distance and diameter of the stone in question. With this rule, each stone is assigned a maximum distance of travel derived from Equation 2.1. Once it is entrained, the stone travels along its trajectory until it stops behind an obstacle (according to the rebound rule) or until its cumulative travel distance has exceeded the maximum travel distance originally assigned to it. If the stone bounces off an obstacle, it continues on its new trajectory until the accumulated travel distance exceeds the maximum travel distance. This rule is implemented in Riverbed 7.0. The entire simulation is essentially scale free except that a representative clast size must be declared in order to supply travel distances and stones sizes. To adopt this computer model to a different prototype with a different representative stone size, the radius of the stone must be known as well as the length and width of the simulated reach. The computer simulation is written such that, all parameters and rules being equal, using the same integer to seed the random number generator will produce identical results. In this way the computer model is deterministic, determinism being used in the same sense as for deterministic chaotic systems, that is, although indivdual trajectories in a chaotic system may be impossible to predict, given exactly the same initial conditions, that trajectory will be exactly reproduced. The computer model is also stochastic in the sense that any outcome depends on the initial random number used. In fact, the model can be more stochastic if more rules and criterion incorporated random numbers. Therefore, a simulation result depends on the design of the computer program as well as on probability. Given that the simulation is fundamentally stochastic, it cannot be expected to reproduce the Harris Creek patterns exactly. At best, it can be expected to yield realistic and similar looking results for a high percentage of trials using the same rules and parameters with different random numbers. 3.2.1 Initial and Boundary Conditions To initiate a simulation, the reach is seeded with randomly placed stones. The number of seed stones is determined by calculating how many non-overlapping circles are required to cover the fraction S E E D F R A C T I O N of the reach area. The parameter S E E D F R A C T I O N is typically 0.20 to 0.30 in accord with the fractional covered area computed for the digitized boxes in the Harris Creek reach. In Riverbed 5.0a radius of 10.0 cm is used to calculate the number of seeding stones. In Riverbed 6.0 and later the average value of 10 cm is used even though the stones have different sizes. The x and y coordinates of the seed stones are determined by drawing uniform random numbers with zero mean and unit variance and scaling this number by the L E N G T H 30 Figure 3.2: Diagram showing how Riverbed calculates the impact position xnew, ynew and the rebound angle vtemp2 using knowledge of the stone's starting coordinates and trajectory, xtemp, ytemp and vtemp respectively, and the obstacles coordinates obx and oby. The equations used to calculate the impact position and rebound angle can be found in part 4 of the entrain() function in all versions of the Riverbed code. The grey filled stone is the obstacle, the black unfilled stone skims this obstacle and moves on as calculated by the black angles and vectors, and the grey unfilled stone stops behind the obstacle, as calculated with the grey angles and vectors. A l l angles and vectors have the same names as in the Riverbed code. Flow direction is along the positive x axis. The line segments labelled 10 are center-to-center lines between the obstacle and the impinging stone, both of radius 5.0. If different sized stones are used then this length must be adjusted accordingly. and W I D T H of the reach respectively. Seed stones do not touch one another. F igure 3.3 shows seeded reaches wi th fractions of 0.10, 0.20 and 0.25 of the reach area covered by randomly placed stones. T h e y are the init ial conditions for an experiment. T h e upstream and downstream boundaries of the simulated reach are open: once a stone manages to travel through the downstream boundary, a new stone is introduced on the upstream boundary. T h i s setup maintains a constant populat ion of stones in the reach at all times. T h e left and right boundaries of the reach are rigid: any stone that intersects them is rebounded off wi th a new trajectory that is equal to the incoming trajectory reflected about the normal to the boundary. 31 200 400 600 800 1000 1200 1400 x downstream (cm) 200 400 600 800 1000 1200 1400 x downstream (cm) Figu re 3.3: The init ial seeded beds for experiments with 10% (top), 20% (center), and 25% (bottom) coverage. The top bed has stones of equal radius and is the init ial condition for experiments using Riverbeds 5.0 and 5.1 . The center bed is the initial condition for experiments using Riverbeds 6.1 and 6.2. The bottom bed is the init ial condition for an experiment using Riverbed 7.0. 32 3.2.2 Experimental Setup After the initial bed has been generated with seed stones and the x and y coordinates and radius of each stone have been written to a file, one of these stones is picked at random to be entrained into the flow. This stone is given a trajectory by drawing a number from a gaussian distribution, with zero mean and unit variance, and scaling by R A N G E , allowing most stones to travel within R A N G E degrees of the mean flow direction. The chosen stone then travels along this trajectory until it encounters an obstacle in its path at which point the rebound rule determines whether the stone glances off the obstacle and continues on a new path or whether it is deposited. If it is deposited then another stone in the reach is picked at random to be entrained into the flow. Since Riverbed is a time-independent model, the length of the simulation is set by the number of random entrainements, R A N D O M C L A S T S , specified for each experiment. Typically 2000 to 3000 random entrainments are specified per simulation and once these have been done the positions of all the stones and their radii are written to a file. An experiment of length R A N D O M C L A S T S = 2000 using 200 stones corresponds to each stone in the river bed having been chosen to be entrained on average 10 times. Experiments using up to 6000 random entrainments do not reveal a pattern of stones significantly differing from the results of shorter duration. The initial and final files are read into a plotting program, in this case Mathematica, where a graph is made of the initial and final configurations of the simulated reach. To maintain a steady state population, if a stone exits the downstream margin of the reach, a new stone is introduced on the upstream margin. The new stone's position is determined by drawing a uniform random number scaled by the W I D T H , and the trajectory is randomly chosen by drawing from a gaussian distribution and scaling by R A N G E . Al l the simulations follow the same basic outline though different parameters or rules can be used. Al l experiments were performed on a Macintosh Quadra 605 personal computer, using Symantec's C compiler. Simulation run times vary from approximately one and a half hours for 300 stones and 2000 random entrainments to three hours for 500 stones and 3500 random entrainments. Omitting the commands for generous screen output in the code will reduce run times significantly. The Mathematica script takes approximately two hours to complete on the Macintosh computer. When run on Sparc 20 computers, the script could take as little as five minutes providing a minimum of 35% of one processor was available for the computing task. 3.2.3 Model Calibration and Testing The W I D T H of the simulated reach is set at 5 m, comparable to the average width of the Harris Creek study reach. The L E N G T H is set at 15 m which is shorter than the length of the study reach but long enough to 33 provide a buffer zone for upstream and downstream boundary effects. T o ameliorate the prob lem of the upstream b o u n d a r y clogging up with stones, the program is modified such that if a newly introduced stone wil l impinge on a stone which itself lies very close to the upstream boundary, the new stone is discarded a nd another one chosen. Since stones near the downstream boundary have relatively fewer potent ial obstacles in their p a t h , they tend to preferentially exit the reach. T h e result is that dur ing the course of a s imulat ion the downstream port ion of the reach tends to empty out its stones. flow direction 0 50 100 150 200 250 300 x downstream (cm) Figure 3.4: Diagram showing the river bed resulting from 20 randomly entrained stones and 20 newly injected stones in a short Riverbed test reach using the rebound rule on stones of equal size. Black circles represent the positions of stones at the initial conditions of the experiment. Entrained stones travel along vectors until they are deposited as a grey stone. Newly injected stones are transported along vectors originating on the upstream boundary (at x = 0) and are deposited in positions denoted by grey circles. Stones which travel out of the reach have vectors extending past the downstream boundary and are counted as throughput. Note that the upstream margin of the reach, at y = 350 - 400, begins to clog up with stones. Even after only 40 entrainments the stones begin to cluster. In experiments using stones of equal radius (as in R i v e r b e d 5 .0 ) , the radius is chosen to be 10 c m corresponding to the D96 (20 cm) of the riffle surface in Harr is Creek. In experiments using stones of different radi i (as in R i v e r b e d 6 .1 and later versions), the radi i are drawn from a lognormal dis tr ibut ion, truncated at a m i n i m u m size of £ > 7 5 = 10 cm, a m a x i m u m size of Dwo = 50 c m and a median diameter of D95 = 20 c m , produc ing a populat ion with a similar size distribution to Harr is Creek. A census is taken of stone size 34 in the init ial bed and the final bed and of the stones which exit the reach. C o m p a r i s o n of these three size distributions determines if stones of certain radi i are preferentially transported through the reach. Since a computer model rapidly grows into a complicated structure capable of unexpected behavior without rigorous attempts to test each subcomponent, R i v e r b e d was tested at every step to ensure that the code performed as intended. For instance, all the r a n d o m number generators were checked by graphing their distributions and calculating their mean and spread. T h e code writes to a file the values of key variables at intermediate stages of the s imulation, providing all the information necessary to follow by h a n d the newly entrained stone, draw its transport vector, determine the appropriate obstacle, the impact a n d rebound angle of this configuration, and the subsequent movement or deposition of the stone. Such an exercise was performed several times dur ing the course of programming produc ing diagrams similar to figure 3.4. T h e informat ion in the file of key values was used to manual ly draw, on the init ial bed , the vector for each entrained stone, following it along checking rebound angles and determining the appropriate obstacle and final posit ion of the stone. T h e marked-up init ial bed is then compared to the final bed to ensure a complete match . 3.2.4 Programming Details R i v e r b e d is based on l inked lists where each list member represents one stone and contains m e m o r y space for the x and y coordinates of the stone in cm, the radius in c m , the trajectory when moving , where clockwise is positive, counterclockwise is negative and downstream is zero, the m a x i m u m travel distance for a stone of that size in c m , and other characteristics as needed. L i n k e d lists are a dynamic d a t a structure meaning that the list is only as long as the number of stones in the reach and shrinks or grows if stones are lost or added; it does not need its size declared beforehand as is necessary if using arrays. L i n k e d lists have other advantages over model ing approaches using arrays, cellular automata or lattice models. For example, l inked lists use real variables in x — y space to posit ion the stones which avoids the effects of rectangular or hexagonal grids and gr id spacing. O n e peculiarity of l inked lists is that to determine the obstacle for any entrained cobble, the cobble-potent ial obstacle vector must be calculated for each member of the list. In other words, rather t h a n following the transport ing cobble along its path and seeing which stone it eventually collides with, l inked lists require that the length and direction of the vector from every stone in the river bed to the cobble are calculated. T h e obstacle chosen is that list member which has the shortest vector within grazing angle of the cobble. 35 3 . 3 S i m u l a t i o n s u s i n g e q u a l s i z e d s t o n e s In order to start simply, the first s imulation R i v e r b e d 5 . 0 uses only the rebound rule on stones of equal radius, in this case 10 cm. A typical result is shown in figure 3.5. In an effort to produce structures with a stronger transverse component and to take into account the shielding effects of neighboring clasts, R i v e r b e d 5 .1 introduces the neighbor rule, requir ing stones to have few close neighbors to be entrained. F igure 3.6 shows the result of a s imulat ion using the same parameters, r a n d o m seed and rules as in the previous experiment with the addit ion of the neighbor rule. 500 400 — 300 E ' 200 100 - v On.*** 200 400 600 800 x downstream (cm) 1000 1200 1400 Figure 3.5: F ina l bed of Riverbed 5.0 using the rebound rule after 2000 random entrainments of 238 stones of radius 10 cm with F U L C R U M = 90 and S E E D F R A C T I O N = 0.1. Flow is to the right. 1400 x downstream (cm) Figu re 3.6: F ina l bed of Riverbed 5.1 using the rebound and neighbor rules after 2000 random entrainments of the 238 stones of radius 10 cm in the bed. F U L C R U M = 90 and S E E D F R A C T I O N = 0.1. Flow is to the right. 36 3.3.1 Discussion T h e bed resulting from R i v e r b e d 5 . 0 displays linear ridges consisting of approximately 3 to 15 members which trend, an d tend to fan out, in the upstream direction. Approx imate ly 5 % of the stones remain in open positions. Similarly, the results of R i v e r b e d 5 .1 show the vast major i ty of stones compris ing upstream undulat ing , elongate clusters of approximately 3 to 20 members. H a r d l y any stones remain in open positions. C o m p a r i n g this result with the previous one, the ridges are less linear with a slightly wider spatial spread than when using only the rebound rule. Therefore, the rebound rule, wi th or without the neighbor rule, is insufficient to produce transverse bedforms with equal sized stones on a sparsely populated bed. B o t h of these simulations produce longitudinal clusters which do not resemble the Harr i s Creek patterns. T h e neighbor rule was added in the hopes of introducing lateral clustering and can be seen to produce clusters with a less pronounced linear trend upstream but it does not significantly change the orientation of the clusters. C o m p a r i n g the results of these two different sets of rules, it appears that the neighbor rule is redundant since the bi l l iard ball rules impose it: stones cannot move downstream if other stones are in front of them. T h e upstream trending clusters of stones produced with R i v e r b e d 5 . 0 resemble cluster bedforms (figure 2.5). A l t h o u g h flow perturbations are thought to be partial ly responsible for cluster bedforms, their presence in this s imulation suggests they may be a result of obstacles to clast movement on a sparsely populated bed. T h e difference in the number of stones remaining in open positions between these two experiments is a result of the design of the computer model which uses the number of r a n d o m entrainments as the m a i n counter governing the length of the s imulation. In R i v e r b e d 5 . 1 , since a proport ion of stones wil l be closely surrounded by other stones and therefore fail to meet the entrainment criterion imposed by the neighbor rule, the l inked list must be searched more frequently to find a stone that can move before the r a n d o m entrainment counter is incremented. Consequently, many more r a n d o m numbers are culled, the list is more thoroughly searched, an d any stones in open positions are likely to be entrained sooner t h a n they would have in R i v e r b e d 5 .0 . If the R i v e r b e d 5 .0 s imulation is run for a longer t ime, eventually all stones leave the open positions an d cluster together. T h i s effect demonstrates how simulat ion results m a y be affected by the design of the computer code and not be solely the expression of the physics p r o g r a m m e d into it. T h e statistical a lgor i thm was performed on the R i v e r b e d 5 .1 experiment (figure 3.7). Since the structures trend parallel to the flow direction, the density plot along x does not register t h e m very well. T h e y would show up better in the y density plot, however as stated previously, the y plots are unreliable due to excessive changes in b in widths. T h e density plot shows a low percentage of covered area for v ir tual ly all angles. T h e dark b a n d at xb in = 10 to 25 is the signature of the accumulat ion of stones on the upstream m a r g i n of the reach. T h e large white area above xb in = 270 is the signature of the empty downstream 37 marg in of the s imulated reach. T h e two m a x i m a in the PC2 graph at 2 3 ° and — 1 6 ° record the high covered area of strips trending at angles of — 6 7 ° and — 7 4 ° in the simulated bed. t h e t a ( d e g r e e s ) Figure 3.7: Density plot (top) and PC2 graph (bottom) for the Riverbed 5.1 result along x. For the density plot, white pixels indicate regions with zero percentage of stones and black pixels indicate regions of high percentage of stones. The peaks in the PC2 graph correspond to angles at which the distribution of stones is least uniform. In the PC2 graph maxima appear at angles +15° and —7° and can be seen to be coincident with the angles at which there exist regions of alternating light and dark pixels in the density plot, that is, along vertical lines in the top graph. These angles correspond to ridges aligned at —75° and 83° in the Harris Creek map. 3.4 Simulations using unequal sized stones Since Harr i s Creek's cobbles comprise a range of grain sizes and cover 20 — 30% of the reach, it is possible that a denser populat ion of unequal sized stones will introduce different structures into the s imulated reach. R i v e r b e d 6 . 1 and all later versions use stones of unequal diameter drawn from a lognormal distr ibution truncated at D75 and Dwo for Harris Creek, 10 c m and 50 c m respectively [9]. R i v e r b e d 6 .1 uses only the 38 rebound rule on stones of unequal size. A typical result is shown in figure 3.8. 500 400 • 300 100 It ~« 4* . . . . . . . . j . * . . . . 600 800 1000 x downstream (cm) Figure 3.8: Final bed of Riverbed 6.1 using 477 unequal sized stones and the rebound rule after 3000 random entrainments. SEEDFRACTION is 0.20 and FULCRUM is set to 60. Flow is to the right. R i v e r b e d 6 . 2 introduces the obstacle rule requiring an obstacle to be at least half the size of the impact ing stone before the stone will stop. T h e obstacle rule takes precedence over the rebound rule when choosing the obstacle for a stone. Consequently, if a candidate obstacle is small enough, the stone wil l pass over it unt i l s topped by a large enough stone, causing some stones to overlap smaller stones when they finally come to rest. A typical result, using the same parameters as the R i v e r b e d 6 .1 experiment, is shown in figure 3.9. 500 400 ~ 300 B o * 200 r—i—2 ' ' — I C — ' ' ' 1 ' ' ' ' 1—W—' 1 1 ' ' ' ' ' 1 1 ' ' ' 1 ' 5V*—' 1 r 200 400 500 800 x downstream (cm) 1000 1200 Figure 3.9: Final bed of Riverbed 6.2 using the obstacle and rebound rules on 477 unequal sized stones after 3000 random entrainments. SEEDFRACTION is 0.20 and FULCRUM is set to 60. Flow is to the right. Because of the primacy of the obstacle rule over the rebound rule, some overlap of stones occurs. 39 3.4.1 Discussion In b o t h R i v e r b e d 6 . 1 and 6 . 2 , the stones tend to form a variety of structures inc luding cluster bedforms and arcuate clusters trending at an angle to flow direction. G o o d examples of cluster bedforms occur in figure 3.8 at coordinates (1140, 300) and in figure 3.9 at (240, 400) and (350, 275) . A r c u a t e ridges occur in figure 3.8 wi th apexes posit ioned at (875, 130) and (1280, 340) and in figure 3.9 centered at (650, 180). N o good examples of transverse ribs occur. B o t h of these simulations produce the unrealist ic structure consisting of a large stone upstream of a string of small stones. T h e addit ion of the obstacle rule sorts the stones and tends to deposit large stones near other large stones, produc ing an unnatura l size dis tr ibut ion in the light of the Harr i s Creek size distributions shown in figures 2.8 to 2.12. Overa l l , the results of these two simulations are more realistic looking than the results using equal sized stones though unrealistic structures occur. These two simulations used twice as many stones as the R i v e r b e d 5 experiments. T h e increase in number of stones contributed to the growth of oblique clusters and ridges. U s i n g a F U L C R U M of 60 allows stones to glance off an obstacle only when the rebound angle is less than 6 0 ° from the flow direction. T h i s means that some stones stop in positions slightly more lateral to the obstacle than they would have if F U L C R U M was 90. However, this effect is not extreme: comparison of simulations in R i v e r b e d 6 . 1 using different values for the F U L C R U M parameter show that lateral clustering is not significantly enhanced by this change alone. C o m p a r i n g the histograms of size distributions (figure 3.10) for the init ial and final beds with the dis tr ibut ion of clasts exiting the reach for the R i v e r b e d 6 . 2 , no size stone is proport ional ly more likely to exit than others. T h e results of the statistical a lgor i thm on the R i v e r b e d 6 . 2 experiment can be found near the end of A p p e n d i x B , the first of two cases shown. T h e density plot is less definitely banded compared to the density plots of the three boxes in Harr is Creek (figures 2.13, 2 .14 and 2.15), yet it is significantly more banded compared to the density plot of the R i v e r b e d 5 . 1 . T h e density plot displays moderately well-developed dark and light bands indicat ing that oblique linear clusters are created better than in R i v e r b e d 5 . 1 . A g a i n , the dark b a n d at xb in = 0 to 30 and the light b a n d above xb in = 290 reflect the high density of stones at the upstream marg in and the paucity of stones at the downstream margin , respectively. T h e PC2 graph has peaks at —23° and 12° corresponding to ridges oriented at 6 7 ° and - 7 8 ° in the s imulated reach. T h e location of these ridges is h a r d to see in the in figure 3.9. T h i s is because cobble ridges in the bed extend only part way across the bed while the strips used in the statistical a lgor i thm extend the full w id th of the reach, registering any stones in that strip regardless of whether there are large gaps between adjacent stones or whether some of the stones more properly belong to other ridges trending at different angles but intersecting the strip. 40 0.3r 0.3. 0.3 0 5 10 15 20 25 0 5 10 15 20 25 0 5 10 15 20 25 radius (cm) radius (cm) radius (cm) Figure 3.10: Size distribution of stones in the initial bed (left), final bed (center) and of the stones that travel through the reach (right) for Figure 3.9. 3.5 Simulations with reluctant large and mobile small stones In gravel bed rivers, large stones are less likely to be entrained and travel shorter distances before stopping than smal l stones [8]. T o incorporate this behavior into the s imulation model Riverbed 7 . 0 uses the entrainment probabi l i ty rule, with CRITERION = 15, a n d the travel distance rule. T h e m a x i m u m travel distance for each stone is computed using equation (2.1) with D5osubsurface = 10 c m a n d I n 6 0 , „ r / „ e = 236 cm. D50 sub surface for Harr is Creek is 6 c m , however the higher value of 10 c m is used because this is the size of the smallest, a n d theoretically the most mobile, stones in the s imulation. A typical result of a s imulat ion using these rules along with the rebound rule is shown in figure 3.11. T h e bed is 25% covered to mirror the coverage of the Harr i s Creek boxes. 500 400 — 300 6 * 200 100 200 400 600 800 1000 x downstream (cm) 1200 1400 Figure 3.11: Final bed of Riverbed 7.0 using 596 unequal sized stones, the entrainment probability, travel distance and rebound rules after 3500 random entrainments. SEEDFRACTION is 0.25 and FULCRUM is set to 60. Flow is to the right. 41 3.5.1 Discussion T h i s combinat ion of rules produces a variety of structures including clusters bedforms, at (1290, 145) and (1030, 160), transverse ribs, extending from (580, 260) to (780, 0), (800, 210) to (800, 500) a n d (420, 0) to (440, 220), an d circular stone cells, centered at (440, 360) and (570, 400). Close examinat ion of figure 3.11 shows the large stones are more or less scattered throughout the bed while the smal l stones tend to cluster together in the vicinity of large stones, a distribution that mirrors Harris Creek's d is tr ibut ion (refer to figures 2.8 to 2.12). T h e stones are less tightly packed than in previous experiments because they can now stop in open positions on the bed and not only flush against an obstacle. 300 V -45 -25 -5 15 35 theta (degrees) 12 _ 10 -40 -20 0 20 40 t h e t a (degrees) Figure 3.12: Density plot (top) and PC2 graph (bottom) along x for the Riverbed 7.0 simulation result in figure 3.11. The density plot exhibits alternate regions of light and dark pixels at an angle of —19°. The P C 2 graph maxima appear at angles +15° and —35°. These angles correspond to ridges aligned at 71°, 55° and —75° in the simulated creek bed. T h e results of the statistical a lgorithm on the R i v e r b e d 7.0 bed are shown in figure 3.12. T h e density plot registers the upstream margin clogged with stones and the downstream marg in empty of stones 42 as in all previous experiments. T h e density plot and PC2 m a x i m a correspond to ridge structures oriented at angles of 7 1 ° , 5 5 ° and - 7 5 ° in the simulated reach. These orientations compare well wi th the trends of the ridges as determined in the rose diagrams (figure 2.3). T h e density plot is better banded than those of all previous simulations and best resembles the results for the Harr i s Creek boxes. In addit ion to a variety of structures, R i v e r b e d 7 . 0 also produces a bed composed almost entirely of cluster bedforms when the stone populat ion is small enough to cover only 10% of the reach area. A typica l result of an experiment using S E E D F R A C T I O N of 0.10, the same value used in R i v e r b e d 5 . 0 and 5 . 1 , is shown in figure 3.13. G o o d cluster bedforms occur at (840, 180), (1240, 80) and (320, 220). x downstream (cm) Figure 3.13: Final bed of Riverbed 7.0 using 2 3 8 unequal sized stones, the entrainment probability, travel distance and rebound rules, after 2 5 0 0 random entrainments. S E E D F R A C T I O N is 0 . 1 0 and F U L C R U M is set to 6 0 . Flow is to the right. 43 Chapter 4 Conclusions Harr i s Creek displays a net-like pattern of cobble ridges, one or two cobbles thick, wi th a variable yet strong transverse alignment. T h e patterns can be characterized as a framework of smaller stones art iculated about r a n d o m l y posit ioned large stones, as suggested by the separate maps of the stones in each size b in (figures 2.8 to 2.12). Alternatively , the patterns can be thought of as an admixture of several different structures: cluster bedforms, transverse ribs, arcuate ridges and isolated stones. C r u d e hydraul ic calculations support observations that the cobble patterns form in sub-critical flows and not in the near crit ical to crit ical flows thought to be necessary to form transverse ribs. Exper iments on a bed with 10% of its area covered by stones almost exclusively produce cluster bedforms regardless of the rules used or whether the stones are all the same size or not. T h i s f inding suggests that cluster bedforms are a result of obstacles to clast movement in a sparsely popula ted bed. In a more densely populated bed, the large obstacle clasts in cluster bedforms m a y become the nuclei f rom which oblique ridges develop out of the excess sediment. Exper iments wi th unequal sized stones introduce a greater variety of orientations to the resulting structures than those developed in equal sized stones. T h i s suggests that the magnitude of the stone populat ion and the range of sizes influences the suite of bedforms that emerge. Simulations of R i v e r b e d 7 . 0 using entrainment probabilities and m a x i m u m travel distances, bo th of which vary inversely proport ional to stone radius, produce patterns most closely resembling Harr i s Creek's patterns. These simulations produce a bed with pebble clusters, transverse and oblique ribs, arcuate ridges and isolated stones. T h e density plot is narrowly and finely banded and is similar to the density plot of the Harr i s Creek boxes. T h e density plot also registers ridges at a similar orientation as indicated by the rose diagrams. T h e rebound rule appears to be an adequate representation of how stones will impact a downstream 44 obstacle an d serves as an adequate base u p o n which to bui ld different rules. It can be substantial ly enhanced by incorporat ing several changes. M o r e randomness can be programmed into the impact to deliver the unpredictable rebound trajectories resulting from irregular shaped stones col l iding with one another. T h e m o m e n t u m exchange between obstacle and impinging stone can be calculated to determine how a smaller stone might be launched into transport as a result of gaining m o m e n t u m from the impact of a larger stone. A s programmed, the neighbor rule is unnecessary since its inclusion does not produce structures differing from using only the rebound rule. A more realistic and effective neighbor rule would discriminate between stones in front of and behind a stone and implement the decrease and increase, respectively, in the lift forces associated with such posit ioning [7]. Latera l neighbors should be discr iminated too a l though no experimental results exist on the influence of lateral neighbors on a particle's mobility. T h e obstacle rule sorts the stones according to size resulting in large stones deposit ing close to other large stones. T h e positions of the largest stones in Harr is Creek are not clustered together like the smaller stones, but scattered throughout the reach. Therefore, the obstacle rule leads to unrealistic results. T h i s implies that large stones can stop behind any size obstacle or that they stop before they collide wi th anything. T h a t stones do not necessarily need an obstacle to their continued downstream movement in order to stop is supported by the results of R i v e r b e d 7 . 0 simulations. T h e entrainment probabi l i ty rule is a realistic rule, supported as it is by empirical investigation [8] and consideration of the fact that at a given shear stress, smaller stones in open positions are more likely to move than larger stones in s imilar positions. T h i s rule uses a simple linear expression for the probabil i ty of entrainment. It m a y benefit f rom using a nonlinear probabi l i ty of entrainment. T h e rule which seems to contribute most to generating realistic cobble structures is the travel distance rule. T h e advantage of this rule is that it has an empirical basis [8] and allows stones to stop in open positions. A l l other rules deposit stones flush against other stones, achieving a strict clustering that is not seen in real river beds. A n o t h e r advantage of the travel distance rule is that it condenses several hydrologic effects into a single relation. For example, the friction between a stone and the surrounding water and between a stone and the river bed wil l act to slow a transport ing stone, eventually causing it to stop. A l s o the inert ia of a stone wil l affect its travel distance in the same manner, and the variabil ity in shear stress in a turbulent flow m a y be strong enough to entrain a stone at one instance but drop significantly the next instant, causing the stone to be deposited. T h e statistical a lgori thm does register linear ridges of stones and their corresponding orientations a l though the results are somewhat muddied by other parts of the bed that intersect the b in strips whose covered area is being calculated. T h e R i v e r b e d computer model is bo th stochastic and deterministic. Randomness enters the mode l v i a the choice of seed value used in the r a n d o m number generator. If the same value is used and all parameters 45 and rules remain the same, the outcome will be exactly duplicated. For this reason, the s imulat ion cannot be expected to exactly reproduce the Harr is Creek patterns but it can be expected to produce similar results over a large number of trials using different seed values. 4.1 Future Work A computer s imulat ion can always be improved, though usually with considerable cost in the form of in -creasing the complexity of the program and longer computat ion times. M a n y of the enhancements to the rebound an d neighbor rules and the entrainment probabil i ty as discussed in the preceeding section can be p r o g r a m m e d to improve realism. F low effects are only qualitatively modeled but could be incorporated by uti l iz ing equations governing water velocity, shear stress, and other pertinent h y d r o d y n a m i c factors. T h e size distributions of the stones that exit the downstream margin of the s imulated reach can be compared to d a t a on sediment transport collected for Harr is Creek to further calibrate and validate the model . T h e statistical a lgori thm can be improved to discriminate ridges. However the exact form of c o m p u -tat ional ly possible improvements is complicated. A t present the a lgori thm only indicates areas of possible ridges. R e p r o d u c i n g Harr is Creek's patterns was the motivat ion for the suite of experiments presented here. T h e success of R i v e r b e d 7 . 0 warrants a more systematic investigation into the effects of different parameter values. It m a y be possible to extend this version of the model to simulate coarse sediment transport in more general situations. Other changes that can be made include allowing stones to exit the reach without having to introduce new ones, no longer maintaining a constant stone populat ion. 46 Appendix A Statistical Algor i thm to Identify Ridge Structures A search for a suitable statistical test to identify transverse ridge patterns in stones uncovered comparable problems in forestry and marine biology. T h e spatial distr ibution of the clasts is s imilar to the spatial d is tr ibut ion of trees in a forest [29] and of acorn barnacles on the side of a ship [11]. T h e d a t a of these problems record the spatial coordinates (x and y) and the radius r of the d a t a points and the analyses cater to cluster formations an d their associated statistics to determine measures of spatial randomness, regularity, and clustering. T h e y do not address the spatial statistics of linear accumulations. For this reason, the following statistical a lgori thm was developed to identify cobble ridges. Suppose two cobble ridges span part of the w id th of the stream as i l lustrated in figure A . l . B y rotat ing the rectangular bounding box about the lower left corner point and calculating the fraction of area covered by the stones in vertical and horizontal strips, as the ridges rotate into parallel ism with the bins, a m a x i m u m of covered area wi l l register. B y calculat ing the covered area all dimensions are rendered unnecessary. In order to cover all possible angles in a 3 6 0 ° range that the ridges might trend, the rectangular bounding box was rotated over a span of 9 0 ° from — 4 5 ° to 4 5 ° , where 0 ° is along the positive x direct ion of the reference axes and also the downstream direction of the unrotated box. Ridges oriented transverse or sub-transverse to the flow would rotate into parallel ism with the strips spanning the reference x axis a nd produce a peak in the histogram of percent area covered along the reference x axis. T h i s is i l lustrated in figure A . l by the transverse ridge on the right. For this ridge, the bins along x register a high percentage of covered area an d produce a peak in the x h istogram but no peak in the his togram spanning the y axis. T h e more longitudinal ly oriented r ib on the left contributes to a uniform distr ibution in the x h i s togram a nd a peak in the y h istogram. 47 x histogram Figure A . l : Schematic diagram illustrating how ridges oriented at different angles register in the x and y histograms. The transverse ridge on the right causes a maximum in the x histogram since, for this angle of 9, it is parallel to the bin strips along x. The longitudinal ridge on the left causes a peak in the y histogram. Flow direction is parallel to the box and in the positive x direction when the box is rotated to 0 ° . After the x and y histograms along the reference axes have been calculated at every degree for the box rotated from — 4 5 ° to 4 5 ° , the PC2 test is calculated along x and y. T h e PC2 test is defined as PC2 = £ { 0 I ~ 6 ) 2 (o.i) where e is the fraction of the total bed covered by the entire stone populat ion , o; is the covered fraction of a single strip an d n is the number of bins. A constant number of histogram bins along the reference x and y axes was used for al l angles of rotat ion to simplify the Mathematica a lgori thm. L o o k i n g at figure A . l , it is evident that the far corners of the box wil l span a greater dstance along the reference x axis when it is rotated from 0 ° . Consequently, the w id th of the bins and the span of the histogram changes as the bed is rotated in order to mainta in a constant number of bins. U s i n g wide bins tends to smooth the histograms whereas using narrow bins tends to show finer detail . However, smaller bins require significantly more computat ion t ime. A t 6 — 0, the b in 48 1750 500 250 -34.4 -22.9 -11.5 11.5 22.9 34.4 theta 0 F i g u r e A . 2 : A graph showing how the span of the x (top curve) and y (bottom curve) histograms change depending on the angle the bounding box is rotated to. The vertical axis is the length of the histogram in cm and the horizontal axis is 0 in degrees. The span of the x histogram is more constant than that of the y histogram. w id th is set to 5 c m , a width that ensures that the smallest stone in the s imulat ion wil l register in at least two adjacent strips. It is important to note that this heuristic does not unequivocally identify cobble ridges. T h e statistical a lgor i thm is only effective at locating regions where the area covered by the stones differs from an assumed uni form distr ibut ion of stones as indicated by m a x i m a in the PC2 graph. A high percentage of covered area for a strip implies a high probabi l i ty of finding a ridge only if adjacent strips have a low percent of covered area. T h i s statistical a lgor i thm would faithfully identify the trends of cobble ridges generously spaced along a line because then the pattern is effectively one dimensional . However, Harr i s Creek's patterns range over two spatial dimensions with a fair amount of variability. Arcuate ridges which span several adjacent strips, regardless of rotat ion angle, are difficult to isolate by this method. T h e following is a description of the a lgori thm. T h e Mathematica code which performs these calcu-lations is in A p p e n d i x B . 1. R e a d in each clast as a vector of the form: x y r 2. G i v e n a value of 9, mult ip ly by the rotation matr ix A cos 9 — sin 9 0 A = sin 9 cos 9 0 0 0 1 49 while preserving the radius to get jx cos 9 — y sin 9 x sin 9 + y cos 9 r 3. Since stones might intersect the edges of the original bounding box, and in order to significantly simplify the a lgor i thm, the bounding box was increased in size in both dimensions by two times the radius of the largest contained stone. T h i s ensures every stone is completely contained in the interior of the box. T h u s for any stone [x cos 9 - y sin 9 + rmax x sin 9 + y cos 9 + rmax 4. T h e four corner points of the bounding rectangle, where xm and ym are the m a x i m u m dimensions in the x and y dimensions, respectively, are mapped as follows: a i : (0 ,0 )—> (0,0) (0.2) ai : (xm,0)—• (xm cos9,xm sinf?) (0.3) a-3 • (xm,ym)—> (xm cos9 - ymsin0,xm sin9 + ym cos9) (0.4) d4 • (0,ym)—• ( -2 / m s in0 ,2 / m cos f3 ) (0.5) 5. T h e four lines forming the bounding box have the following equations (for all non-zero 9) L i : y = {t<m9)x (0.6) L2 : y = -(cot9)x + ^ (0.7) L3 : 2/ = ( t a n 0 ) x + ^ (0.8) L 4 : y = -(cot9)x (0.9) 6. T h e project ion of the 'outermost' points of the box onto the x and y reference axes defines the span along x and y of the histograms. For 9 6 [ 0 ° , 4 5 ° ] the x h istogram has a span from the x-coordinate of the 0 4 point, denoted (a,i)x in figure A . l , to the x-coordinate of the ai point, denoted {a?)x- Expressed in absolute units, this is ymsm9 + xmcos9. T h e span of the y h is togram is from ( a i ) y to ( 0 3 ) 3 , or ym cos9 + xm sinf5. B y symmetry, for 9 6 [ - 4 5 ° , 0 ° ] the span changes in a similar manner as i l lustrated in figure A . 2 . It is important to note that the span of either his togram is a funct ion of 9 and can change dramatical ly . 50 y bin: x b i n : j k I F i g u r e A . 3 : Given a clast that spans several bins in either the x or y dimension, it is necessary to determine which bin the leftmost (bottommost) point, the centre, and the rightmost (topmost) point are contained. The area of the stone in each bin is computed by the integrals Ix or Iy . 7. T o calculate the area of each clast in any given b in , Ix: rb r + Vr2 — X2 Ix = / / dydx (0.10) rb = 2 / Vr2 -x2dx (0.11) J a = by/r2 — b2 + r2 arcsin - — a \/r2 — a2 — r2 arcsin — (0-12) r r = b \ / r 2 — b2 — a V'T2 — a2 + r2 . b . a arcsin arcsin — r r F o r horizontal strips, the area of each clast in any given b in , Iy, is the same: Iy = d \/r2 - d? - c \Jr2 - c2 + r2 . d . c arcsin arcsin -r r (0.13) (0.14) 8. T h e area of any strip can be calculated from the equations of the four perpendicular lines m a k i n g up the bounding box. Since the number of bins on either axis is specified and as the span of the histograms is known, the width and area of each b in can be calculated. 9. T h e PC2 test, as defined in equation 0.1, was redefined to minimize the contr ibut ion of bins at either end of the bounding box which, by virtue of experiencing dramat ic width changes dur ing rotat ion, can re turn anomalous values. For example, if the bounding box was rotated counterclockwise by one degree, the leftmost and rightmost bins of the x h istogram have significantly smaller areas than those in the interior of the box. It is very likely that a stone could cover a high percentage of the area of this smal l strip and return an usually high value for the fraction of covered area for that strip. T h e number of bins in b o t h histograms was set to 310 (thus for 8 = 0° the width of a b in in the x h i s togram would have an init ial width of 5) thus 80 bins on either side were ignored to calculate the PC2 test, redefined 51 230 , ,.2 p c 2 = £ (0i^£)_ i=81 (0.15) 52 Appendix B Mathematica Script for Statistical Algor i thm T h e following is an edited version of the Mathematica script that performs the statistical test for ridge structures. T h e script generates data for a range of 8 from - 4 5 ° to 45°. T h e datafile b e d l . d a t used in this example is the final result of the Riverbed 6.1 experiment shown in figure 3.8. T h e density plot a n d PC2 graph for this result is shown in the script using the data file of the final condit ion of the Riverbed 6.1 s imulat ion as shown in figure 3.8. adat=ReadList["bedl.dat".{Number.Number.Number}]; rmax=25; nx=310; ny=310; xmax=1500+2 rmax; ymax=500+2 rmax; ms=Dimensions[adat] sa=0; bdat=Table [0, {ms [ [1] ] }. {ms [ [2] ] }] ; For[i=l,i<=ms[[1]],i=i+l,bdat[[i,1]]=adat[ [i,1]]+rmax; bdat [ [ i , 2] ] =adat [ [ i , 2] ] +rmax; bdat[[i,3]]=adat[[i,3]] ; sa+=N[Pi adat [[i,3]] "2] ;] ; N[sa,10] 156653.979 A[theta_] :=N[{{Cos[N[theta Pi/180]] ,-Sin[N[theta Pi/180]] ,0}, {Sin[N[theta Pi/180]] .Cos[Nftheta Pi/180]] ,0},{0,0,1}}] ; xmat=Table[0,{91},{nx}]; ymat=Table[0,{91},{ny}] ; al={0.0}; a2={N[xmax Cos[theta Pi/180]] ,N[xmax Sinftheta Pi/180]]}; a3={N[xmax Cosftheta Pi/180]-ymax S i n f t h e t a Pi/180]], 53 N[xmax Sin[theta Pi/180]+ymax Cos[theta Pi/180]]}; a4={N[-ymax Sin[theta Pi/180]] ,N[ymax Cos[theta Pi/180]]}; Ix=Compile[{a,b,r},b Sqrt[r~2-b~2]-a Sqrt[r"2-a"2] +r~2 (ArcSin[b/r]-ArcSin[a/r])] ; Iy=Compile[{c,d,r},d Sqrt[r~2-d~2]-c Sqrt[r"2-c~2] +r~2 (ArcSin[d/r]-ArcSin[c/r])] ; For[theta=0,theta<=45,theta=theta+l, cdat=N [Transpose [A [theta] . Transpose [bdat] ] ] ; ddat=Table [0, {ms [ [1] ] }, {ms [ [2] ] }] ; al={0,0}; a2={N[xmax Cos[theta Pi/180]] ,N[xmax Sin[theta Pi/180]]}; a3={N[xmax Cos[theta Pi/180]-ymax Sin[theta Pi/180]], N[xmax Sin [theta Pi/180]+ymax Cos [theta Pi/180]]}; a4={N[-ymax Sin[theta Pi/180]] ,N[ymax Cos[theta Pi/180]]}; For [ i = l , i<=ms [ [1] ] , i=i+l,ddat [ [ i , 1] ] =cdat [ [ i , 1] ] +Abs [a4 [ [1] ] ] ; ddat [ [ i , 2] ] =cdat [ [ i , 2] ] ; ddat[[i,3]]=cdat[[i,3]] ; ] ; dx=N [ (-a4 [ [1] ] +a2 [ [1] ] ) /nx] ; dy=N[a3[[2]]/ny]; For[i=l,i<=ms[[1]],i=i+l, x=ddat [ [ i , 1] ] ; r=ddat [ [ i , 3] ] ; If[r!=0,r=r,r=0.00001]; j =Floor[(x-r)/dx]+1; k=Floor [x/dx]+l; l=Floor[(x+r)/dx]+l; del=x-(k-l) dx; eps=k dx-x; m=k-j; n=l-k; If[m>0,For[s=m,s>m-0.5,s=s-l,xmat[[theta+46,k-s]]+=Ix[-r, (k-s) dx-x,r]]; For[s=m-l,s>=l,s=s-l,xmat[[theta+46,k-s]]+=Ix[(k-s-l) dx-x, (k-s) dx-x,r]] For[s=0,s>=-0.5,s=s-l,xmat[[theta+46,k]]+=Ix[(k-l) dx-x,0,r]];, xmat [ [theta+46, k] ] +=Ix [-r, 0, r] ] ; If[n>0,For[s=n,s>n-0.5,s=s-l,xmat[[theta+46,k+s] ] +=Ix[(k+s-1) d x - x , r , r ] ] ; For[s=n-l,s>=l,s=s-l,xmat[[theta+46,k+s]]+=Ix[(k+s-l) dx-x,(k+s) dx-x,r]]; For[s=0,s>=-0.5,s=s-l,xmat[[theta+46,k]]+=Ix[0,k dx-x,r]] ; , xmat[[theta+46,k]]+=Ix[0,r,r]] ; ] ; For [i=l,i<=ms[[1]],i=i+l, y=ddat[[i,2]] ; r=ddat[[i,3]] ; If[r!=0,r=r,r=0.00001]; j=Floor[(y-r)/dy]+l; k=Floor [y/dy]+l; l=Floor[(y+r)/dy]+l; del=y-(k-l) dy; eps=k dy-y; m=k-j; n=l-k; If[m>0,For[s=m,s>m-0.5,s=s-l,ymat[[theta+46,k-s]]+=Iy[-r,(k-s) dy-y,r]] ; For[s=m-l,s>=l,s=s-l,ymat[[theta+46,k-s]]+=Iy[(k-s-l) dy-y, (k-s) dy-y,r]] For[s=0,s>=-0.5,s=s-l,ymat[[theta+46,k]]+=Iy[(k-l) dy-y,0,r]];, ymat[[theta+46,k]]+=Iy[-r,0,r]] ; If[n>0,For[s=n,s>n-0.5,s=s-1,ymat[[theta+46,k+s]]+=Iy[(k+s-1) d y - y , r , r ] ] ; For[s=n-l,s>=l,s=s-l,ymat[[theta+46,k+s]]+=Iy[(k+s-l) dy-y,(k+s) d y - y , r ] ] ; For[s=0,s>=-0.5,s=s-l,ymat[[theta+46,k]]+=ly[0,k d y - y , r ] ] ; , ymat [ [theta+46, k] ] +=Iy [0, r , r] ] ; ] ; 54 For[theta=-l,theta>=-45,theta=theta-l, cdat=N [Transpose [A [theta] . Transpose [bdat] ] ] ; ddat=Table [0, {ms [ [1] ] }, {ms [ [2] ] }] ; al={0,0}; a2={N[xmax Cos[theta Pi/180]] ,N[xmax Sin[theta Pi/180]]}; a3={N[xmax Cos[theta Pi/180]-ymax Sin[theta Pi/180]], N[xmax Sin [theta Pi/180]+ymax Cos [theta Pi/180]]}; a4={N[-ymax Sin[theta Pi/180]] ,N[ymax Cos[theta Pi/180]]}; For[i=l,i<=ms[[l]] ,i=i+l ,ddat [ [ i , 1]] =cdat [ [ i , 1] ] ; ddat [ [ i , 2] ] =cdat [ [ i , 2] ] +Abs [a2 [ [2] ] ] ; ddat[[i,3]]=cdat[[i,3]] ; ] ; dx=N[(a3[[l]])/nx] ; dy=N [ (a4 [ [2] ] -a2 [ [2] ] ) /ny] ; For[i=l,i<=ms[[1]],i=i+l, x=ddat[[i,l]] ; r=ddat[[i,3]] ; If [r!=0,r=r,r=0.00001]; j=Floor[(x-r)/dx]+1; k=Floor [x/dx]+l; l=Floor[(x+r)/dx]+l; del=x-(k-l) dx; eps=k dx-x; m=k-j; n=l-k; If[m>0,For[s=m,s>m-0.5,s=s-l,xmat[[theta+46,k-s]]+=Ix[-r,(k-s) dx-x,r]]; For[s=m-l,s>=l,s=s-l,xmat[[theta+46,k-s]]+=Ix[(k-s-l) dx-x, (k-s) dx-x,r]]; For[s=0,s>=-0.5,s=s-l,xmat[[theta+46,k]]+=Ix[(k-l) dx-x,0,r]];, xmat[[theta+46,k]]+=Ix[-r,0,r]] ; If[n>0,For[s=n,s>n-0.5,s=s-1,xmat[[theta+46,k+s] ] +=Ix[(k+s-1) d x - x , r , r ] ] ; For[s=n-l,s>=l,s=s-l,xmat[[theta+46,k+s]] +=Ix[(k+s-l) dx-x,(k+s) dx-x,r]]; For[s=0,s>=-0.5,s=s-l,xmat[[theta+46,k]]+=Ix[0,k dx-x,r]] ; , xmat[[theta+46,k]]+=Ix[0,r,r]] ; ] ; For[i=l,i<=ms[[1]],i=i+l, y=ddat[[i,2]] ; r=ddat[[i,3]] ; If[r!=0,r=r,r=0.00001]; j=Floor[(y-r)/dy]+l; k=Floor [y/dy]+l; l=Floor[(y+r)/dy]+l; del=y-(k-l) dy; eps=k dy-y; m=k-j; n=l-k; If[m>0,For[s=m,s>m-0.5,s=s-l,ymat[[theta+46,k-s]]+=Iy[-r,(k-s) dy-y, r ] ] ; For[s=m-l,s>=l,s=s-l,ymat[[theta+46,k-s]]+=Iy[(k-s-l) dy-y, (k-s) d y - y , r ] ] ; For[s=0,s>=-0.5,s=s-l,ymat[[theta+46,k]]+=Iy[(k-l) dy-y,0,r]];, ymat [ [theta+46, k] ] +=Iy [-r, 0, r] ] ; If [n>0,For[s=n,s>n-0.5,s=s-l,ymat[[theta+46,k+s] ] +=Iy[(k+s-1) d y - y , r , r ] ] ; For[s=n-l,s>=l,s=s-l,ymat[[theta+46,k+s]]+=Iy[(k+s-l) dy-y,(k+s) d y - y , r ] ] ; For[s=0,s>=-0.5,s=s-1,ymat[[theta+46,k]]+=Iy[0,k dy-y, r] ] ; , ymat[[theta+46,k]]+=ly[0,r,r]] ; ] ; 13xl4x=Compile[{a,b,theta},l/2 (b"2-a"2) (Tan[theta Pi/180]+ Cot[theta Pi/180])+(b-a) (ymax/Cos [theta Pi/180])]; 13xllx=Compile[{a,b,theta},(b-a) ymax/Cos[theta Pi/180]]; 12xllx=Compile[{a,b,theta},1/2 (a"2-b"2) (Tan[theta Pi/180]+ Cot[theta Pi/180]) +(b-a) (xmax/Sin[theta Pi/180])]; 55 normxmat=Table[0,{91},{nx}]; For[theta=0,theta<=0.5,theta=theta+l, al={0,0}; a2={N[xmax Cos[theta Pi/180]] ,N[xmax Sin[theta Pi/180]]}; a3={N[xmax Cos[theta Pi/180]-ymax Sin[theta Pi/180]], N[xmax Sin [theta Pi/180]+ymax Cos [theta Pi/180]]}; a4={N[-ymax Sin[theta Pi/180]] ,N[ymax Cos[theta Pi/180]]}; dx=N [ (-a4 [ [1] ] +a2 [ [1] ]) /nx] ; dy=N[a3[[2]]/ny] ; For [i=l,i<=nx,i=i+l, normxmat[[thet a+46,i] ] +=ymax dx; ] ; ] ; For [theta=l,theta<=45,theta=theta+l, al={0,0}; a2={N[xmax Cos[theta Pi/180]] ,N[xmax Sin[theta Pi/180]]}; a3={N[xmax Cos[theta Pi/180]-ymax Sin[theta Pi/180]], N[xmax Sin [theta Pi/180]+ymax Cos [theta Pi/180]]}; a4={N[-ymax Sin[theta Pi/180]] ,N[ymax Cos[theta Pi/180]]}; dx=N [ (-a4 [ [1] ] +a2 [ [1] ]) /nx] ; dy=N[a3[[2]]/ny] ; il=Floor[(-a4[[1]]+al[[1]])/dx]+1; i2=Floor[(-a4[[1]]+a3[ [1]])/dx]+1; F o r [ i = l , i < = i l - l , i = i + l , normxmat[[theta+46,i]]+=13xl4x[a4[[l]] + ( i - l ) dx,a4[[l]]+i dx,theta] ; normxmat [[-theta+46,i]]+=13xl4x[a4[[l]] + ( i - l ) dx,a4[[l]]+i dx,theta] ; ] ; For [i=il,i<=il+0.5,i=i+l, normxmat[[theta+46,i]]+=13xl4x[a4[[l]] + ( i - l ) dx,0,theta] ; normxmat [ [theta+46, i ] ] +=13xllx [0, a4 [ [1] ] +i dx,theta] ; normxmat [[-theta+46,i]]+=13xl4x[a4[[l]] + ( i - l ) dx,0,theta] ; normxmat [[-theta+46,i]]+=13xllx[0,a4[[l]]+i dx,theta] ; ] ; For [ i = i l + l , i < = i 2 - l , i = i + l , normxmat [[theta+46,i]]+=13xllx[a4[[l]] + ( i - l ) dx,a4[[l]]+i dx.theta] ; normxmat[[-theta+46,i]]+=13xllx[a4[[l]] + ( i - l ) dx,a4[[l]]+i dx,theta] ; ] ; For [i=i2,i<=i2+0.5,i=i+l, normxmat[[theta+46,i]]+=13xllx[a4[[l]] + ( i - l ) dx,a3[[l]] ,theta] ; normxmat[[theta+46,i]]+=12xllx[a3[[l]] ,a 4 [ [ l ] ] + i dx.theta] ; normxmat [ [-theta+46, i]]+=13xllx[a4[[l]] + ( i - l ) dx,a3[[l]] ,theta] ; normxmat [[-theta+46,i]]+=12xllx[a3[[l]] ,a 4 [ [ l ] ] + i dx.theta] ; ] ; For [i=i2+l,i<=nx,i=i+l, normxmat [[theta+46,i]]+=12xllx[a4[[l]] + ( i - l ) dx,a4[[l]]+i dx.theta] ; normxmat[[-theta+46,i]]+=12xllx[a4[[l]] + ( i - l ) dx,a4[[l]]+i dx.theta] ; ] ; ] ; 14ylly=Compile[{c,d,theta},1/2 (d~2-c"2) (Tan[theta Pi/180]+Cot[theta Pi/180])]; 14yl2y=Compile[{c,d,theta},xmax/Cos[theta Pi/180] (d-c)]; 13yl2y=Compile[{c,d,theta},(xmax/Cos[theta Pi/180]+ymax/Sin[theta Pi/180]) (d-c) -1/2 (d"2-c*2) (Tan[theta Pi/180]+Cot[theta Pi/180])]; 13ylly=Compile[{c,d,theta},ymax/Sin[theta Pi/180] (d-c)]; normymat=Table[0,{91},{ny}]; For[theta=0,theta<=0.5,theta=theta+1, 56 al={0,0}; a2={N[xmax Cos[theta Pi/180]] ,N[xmax Sin[theta Pi/180]]}; a3={N[xmax Cos[theta Pi/180]-ymax Sin[theta Pi/180]], N[xmax Sin [theta Pi/180]+ymax Cos [theta Pi/180]]}; a4={N[-ymax Sin[theta Pi/180]] ,N[ymax Cos[theta Pi/180]]}; dy=N[a3[[2]]/ny] ; For[i=l,i<=ny,i=i+l, normymat[[thet a+46,i]]+=xmax dy; ] ; For[theta=l,theta<=19,theta=theta+l, al={0,0}; a2={N[xmax Cos[theta Pi/180]] ,N[xmax Sin[theta Pi/180]]}; a3={N[xmax Cos[theta Pi/180]-ymax Sin[theta Pi/180]], N[xmax Sin [theta Pi/180]+ymax Cos [theta Pi/180]]}; a4={N[-ymax Sin[theta Pi/180]] ,N[ymax Cos[theta Pi/180]]}; dy=N[a3[[2]]/ny]; il=Floor[(a2[[2]])/dy]+1; i2=Floor[(a4[[2]])/dy]+1; F o r [ i = l , i < = i l - l , i = i + l , normymat[[theta+46,i]]+=14ylly[(i-l) dy.i dy,theta]; normymat[[-theta+46,i]]+=14ylly[(i-l) d y , i dy,theta]; ] ; For[i=il,i<=il+0.5,i=i+l, normymat[[theta+46,i]]+=14ylly[(i-l) dy,a2[[2]] ,theta] ; normymat [[theta+46,i]]+=14yl2y[a2[[2]] , i dy,theta] ; normymat[[-theta+46,i]]+=14ylly[(i-l) dy,a2[[2]] ,theta] ; normymat[[-theta+46,i]]+=14yl2y[a2[[2]],i dy,theta]; ] ; F o r [ i = i l + l , i < = i 2 - l , i = i + l , normymat [ [theta+46, i]]+=14yl2y[(i-l) d y , i dy, theta] ; normymat[[-theta+46,i]]+=14yl2y[(i-l) dy.i dy,theta] ; ] ; For [i=i2,i<=i2+0.5,i=i+l, normymat[[theta+46,i]]+=14yl2y[(i-l) dy,a4[[2]] ,theta] ; normymat[[theta+46,i]]+=13yl2y[a4[[2]] , i dy,theta] ; normymat[[-theta+46,i]]+=14yl2y[(i-l) dy,a4[[2]] ,theta] ; normymat[[-theta+46,i]]+=13yl2y[a4[[2]] , i dy,theta] ; ] ; For [i=i2+l,i<=ny-l,i=i+l, normymat[[theta+46,i]]+=13yl2y[(i-l) dy.i dy,theta]; normymat[[-theta+46,i]]+=13yl2y[(i-l) d y , i dy,theta]; ] ; For [i=ny,i<=ny+0.5,i=i+l, normymat[[theta+46,i]]+=13yl2y[(i-l) dy,a3[[2]] ,theta] ; normymat[[-theta+46,i]]+=13yl2y[(i-1) dy,a3[[2]] ,theta] ; ] ; ] ; For[theta=20,theta<=45,theta=theta+l, al={0,0}; a2={N[xmax Cos[theta Pi/180]] ,N[xmax Sin[theta Pi/180]]}; a3={N[xmax Cos[theta Pi/180]-ymax Sin[theta Pi/180]], N[xmax Sin[theta Pi/180]+ymax Cos[theta Pi/180]]}; a4={N[-ymax Sin[theta Pi/180]] ,N[ymax Cos[theta Pi/180]]}; dy=N[a3[[2]]/ny] ; i l = F l o o r [ ( a 4 [ [2]])/dy]+1; i2=Floor [ (a2[ [2]])/dy]+1; F o r [ i = l , i < = i l - l , i = i + l , normymat[[theta+46,i]]+=14ylly[(i-l) d y . i dy,theta]; 57 normymat[[-theta+46,i]]+=14ylly[(i-l) d y , i dy,theta]; ] ; For[i=il,i<=il+0.5,i=i+l, normymat[[theta+46,i]]+=14ylly[(i-l) dy,a4[[2]] ,theta] ; normymat[[theta+46,i]]+=13ylly[a4[[2]] , i dy,theta] ; normymat[[-theta+46,i]]+=14ylly[(i-l) dy,a4[[2]] ,theta] ; normymat[[-theta+46,i]]+=13y1ly[a4[[2] ] , i dy,theta]; ] ; F o r [ i = i l + l , i < = i 2 - l , i = i + l , normymat[[theta+46,i]]+=13ylly[(i-l) dy.i dy,theta]; normymat[[-theta+46,i]]+=13ylly[(i-l) d y , i dy,theta]; ] ; For[i=i2,i<=i2+0.5,i=i+l, normymat[[theta+46,i]]+=13ylly[(i-l) dy,a2[[2]] ,theta] ; normymat[[theta+46,i] ] +=13yl2y[a2[[2] ] , i dy,theta]; normymat [ [-theta+46, i ] ] + = 1 3 y l l y [ ( i - l ) dy,a2[[2]] ,theta] ; normymat [[-theta+46,i]]+=13yl2y[a2[[2]] , i dy,theta] ; ] ; For [i=i2+l,i<=ny-l,i=i+l, normymat[[theta+46,i]]+=13yl2y[(i-l) d y , i dy,theta] ; normymat[[-theta+46,i]]+=13yl2y[(i-l) dy.i dy,theta]; ]; For [i=ny,i<=ny+0.5,i=i+l, normymat[[theta+46,i]]+=13yl2y[(i-l) dy,a3[[2]] ,theta] ; normymat[[-theta+46,i]]+=13yl2y[(i-l) dy,a3[[2]] ,theta] ; ] ; ] ; xresmat=Table[0,{91},{nx}] ; For[i=l,i<=91,i=i+l, For[j=l,j<=nx,j=j+l.xresmat[[i,j]]=N[xmat[[i,j]]/normxmat[[i,j]]]] ] ; xresmax=Max[xresmat] 0.634253 f [x_] :=Which[ 0 == x , 1.0, 0 0 < x <= 0.05, 0 95, 0 05 < x <= 0.1, 0 9, 0 1 < x <= 0.15, 0 85, 0 15 < x <= 0.2, 0 8, 0 2 < x <= 0.25, 0 75, 0 25 < x <= 0.3, 0 7, 0 3 < x <= 0.35, 0 65, 0 35 < x <= 0.4, 0 6, 0 4 < x <= 0.45, 0 55, 0 45 < x <= 0.5, 0 5, 0 5 < x <= 0.55, 0 45, 0 55 < x <= 0.6, 0 4, 0 6 < x <= 0.65, 0 35, 0 65 < x <= 0.7, 0 3, 0 7 < x <= 0.75, 0 25, 0 75 < x <= 0.8, 0 2, 0 8 < x <= 0.85, 0 15, 0 85 < x <= 0.9, 0 1, 0 9 < x <= 0.95, 0 05, 0 95 < x <= 1.0, 0 0]; xtest=Table[0,{91},{nx}]; For[i=l,i<=91,i=i+l, For[j=l,j<=nx,j=j+l,xtest[[i,j]]=N[xresmat[[i,j]]/xresmax]]; 58 testx=Map[f,xtest, {2}]; Show[DensityGraphics[Transpose[testx]], Mesh->False, AspectRatio->0.618,Frame->True, FrameLabel->{"theta (degrees) ", "x bin"}, PlotRange->{{-l ,92},-[-5,nx+5} ,{0.0,1.0}}] -45 -25 -5 15 35 theta (degrees) xmidpoints=Table[-0.5+j-45,{j,91}]; xchitest=Table[0,{91}]; xoffset=10; For[i=l,i<=91,i=i+l, For [j = l+xoffset,j<=nx-xoffset,j=j+l, x c h i t e s t [ [i]]+=(xresmat[[i,j]]-sa/(xmax ymax))~2/(sa/(xmax ymax))]; ] ; ListPlot[Transpose[{xmidpoints.xchitest}], PlotRange->{5,20}, Frame->True, FrameLabel->{"theta (degrees)","pc squared ( i n x)"}, Axes0rigin->{-50,0}, PlotJoined -> True, PlotStyle->AbsoluteThickness[0.5],AspectRatio->0.4] 12 10 8 -40 -20 0 20 40 t h e t a (degrees) 60 yresmat=Table[0,{91},{ny}]; For[i=l,i<=91,i=i+l, For[j=l,j<=ny,j=j+l, yresmat [ [ i , j ] ] =N[ymat [ [ i , j ] ] /normymat [ [ i , j ] ] ] ] ; ] ; yresmax=Max[yresmat] 0.671906 ytest=Table[0,{91},{ny}]; For[i=l,i<=91,i=i+l, For[j=l,j<=ny,j=j+1,ytest[[i,j]]=N[yresmat[[i,j]]/yresmax]]; ] testy=Map [f , yt e s t , {2}] ; Show[DensityGraphics[Transpose[testy]],Mesh->False, AspectRatio->0.618, Frame->True, FrameLabel->{"theta (degrees)","y bin"}, PlotRange->{{-l,92},{-5,ny+5},{0.0,1.0}}] -45 -25 -5 15 35 theta (degrees) ymidpoints=Table[-0.5+j-45,{j,91}]; ychitest=Table [0,{91}]; yoffset=10; For[i=l,i<=91,i=i+l, For[j=l+yoffset,j<=ny-yoffset,j=j+l, ychitest[[i]]+=(yresmat[[i,j]]-sa/(xmax ymax))~2/(sa/(xmax ymax))]; ] ; ListPlot[Transpose[{ymidpoints.ychitest}],PlotRange->{0,20}, Frame->True, FrameLabel->{"theta (degrees)","pc squared ( i n y)"}, Axes0rigin->{-50,0}, PlotJoined -> True, PlotStyle->AbsoluteThickness[0.5],AspectRatio->0.4] 61 10 -40 -20 0 20 40 t h e t a (degrees) 62 Appendix C Computer Code for Riverbed #include <iostream.h> #include <fstream.hpp> #include <math.h> #include <stdio.h> #include <stdlib.h> #define Width 500.0 #define Length 1500.0 #define SeedFraction 0.20 fldefine RandomClasts 3001 #define Range 10.0 #define Lambda 4.0 #define Epsilon 0.0 #define D_75 10.0 #define D_95 20.0 #define D.100 50.0 #define D_50sub 10 #define L_D50surf 236 #define FULCRUM 60.0 #define Criterion 15.0 //for f i l e class //for exit //width of reach in cm //length of reach in cm //fraction of riverbed covered by seed stones //belonging to D_95. //number of random entrainments, main counter //for program //for angular range of trajectories //for lognormal distribution of radii //for lognormal distribution of radii //Harris Creek minimum r i f f l e surface grain size, cm. //grain size of the 95th coarsest percentile //in Harris Creek in cm. //maximum grain size in cm allowed in Harris Creek, //mean diameter of subsurface clasts in cm. //mean distance of movement in cm for D_50 surface //material. //skim/stop criterion in entrainO in degrees //either side of downstream, //criterion for clast entrainment. // RIVERBED 7.0 // Riverbed 7.0 models coarse clast transport in a stream bed with the rebound rule along with // an entrainment probability and travel distance for each stone that varies inversely // proportionally to the radius. A stone is introduced into the simulated reach from a random // position on the upstream margin or is randomly entrained from the seeded stone population // produced in the i n i t i a l conditions. It travels downstream with a randomly chosen trajectory // and may encounter obstacles. Upon impact with an obstacle, the stone's rebound angle is // calculated by the rebound rule to be equal to the incident angle measured about the normal // to tangent of the impact site. It stops i f rebound angle is greater than FULCRUM degrees; // otherwise the stone moves along the new trajectory to further encounters. If a stone exits // the reach or travels through the bed without encountering obstacles, a new stone is // introduced upstream to take i t s place. Riverbed 7.0 also uses different radii drawn from // a truncated lognormal distribution, and an anti-clogging rule (Vanish) for the upstream // margin. The river bed width and length are in cm. Detailed notes on the algorithm appear at // the beginning of crucial loops and for each function. // Written August 5-13, 1996 by Selina Tribe. 63 struct NodeType { float x; float y; float r; float v; float 1; NodeType* link; }; typedef NodeType* NodePtr; NodePtr head; NodePtr endPtr; long seed=-34581; float SeedNumber; long *idum; float tr,g; int totalClasts; int throughput; int adjust; int input; int rebound; int vanishedRock; int reversal; int deleteClast; int stop; int vanish; int Number; int exitClast; float xtemprc, ytemprc, rtemprc; float vtemprc, ltemprc; float D; float L; streampos savedPos3=0; float rani(long *idum); float gasdev(long *idum); void seedbedO ; void entrainQ ; void f ileSeedBedO ; void filebed (); //dynamic memory allocation to make a linked l i s t //x position //y position //clast radius //clast trajectory //maximum travel length //pointer to next struct //head of l i s t //points to end of l i s t / / i n i t i a l seed for rani //number of stones in the bed at i n i t i a l conditions. //for random number generator //results of ranlO and gasdevO . //counter: total number of stones in the reach, //= SeedNumber + input - throughput. //counter: number of stones passing through the //downstream boundary of reach. //counter: how many times new positions were adjusted //because they f e l l out of the reach, //counter: number of new stones introduced into the //active bed. //counter: total number of rebounds on reach margins //counter for number of repicks of new clast to avoid //clogging the upstream margin. //counter for number of trajectory reversals //performed on bed margins. //flag signifying a re-entrained clast. //flag signifying a clast has been intercepted and //stopped, to get out of entrainQ . //flag signalling that obstacle picked is too close //to upstream margin so repick new stone position, //counter: number of nodes in l i s t . //flag:inject new stone since one has passed through, //coordinates, radius of randomly entrained stone, //angle, maximum distance of randomly entrained stone. //D=2*radius/D_50sub. //L= l/L_D50surf. //for f i l e "values.dat" void main () { int r c , i ; NodePtr NodePtr int Position; int moveable; delPtr; prevPtr; //counters //pointer to randomly entrained node to be deleted //from l i s t //pointer to node before entrainPtr allowing for //easy deletion //position in l i s t chosen to be deleted=stone chosen //to be re-entrained. //flag indicating whether stone has f u l f i l l e d //entrainment criterion. // Initialize riverbed by seeding i t with clast. idum=&seed; Set counters and flags to zero. 64 seedbed(); cout«"bed seeded. "«"\n"; f ileSeedBedO ; cout«"bed filed"«"\n"; throughput = 0; input = 0; adjust = 0; rebound = 0; r e v e r s a l = 0; stop = 0; deleteClast = 0; vanishedRock = 0; vanish = 0; moveable = 0; // random entrain loop // 1. This procedure w i l l randomly chose a stone located at some p o s i t i o n i n the ri v e r b e d // to be entrained into the flow. Each stone picked randomly w i l l have a p r o b a b i l i t y of // entrainment inversely proportional to s i z e : a large stone i s l e s s l i k e l y to be picked. // The stone i s chosen from the e x i s t i n g l i s t by p i c k i n g a uniform random number ranging // from 1 to Number, the length of the l i s t . A t r a j e c t o r y i s picked from a normal // d i s t r i b u t i o n with the same variance as used f o r the new stones. Any node i n the l i s t // can be entrained: i f the f i r s t node happens to be chosen then the head pointer and a l l // other pointers, which Eire always repositioned to the l i s t head, are s h i f t e d one down. // The pointer delPtr i s positioned at the node i n the l i s t representing the chosen // stone then x, y, r , t r a v e l distance, v are written to temporary var i a b l e s (xtemprc, // ytemprc, rtemprc, ltemprc, vtemprc). The f l a g deleteClast i s set to 1. The pointer // prevPtr t r a v e l s down the l i s t to the node immediately p r i o r to de l P t r . The pointers // point to new po s i t i o n s and delPtr i s deleted r e s u l t i n g i n the l i s t being one node // shorter. Xtemprc, ytemprc, rtemprc, vtemprc are fed into the e n t r a i n O fu n c t i o n to // compute where the stone ends up. for(rc=l;rc<RandomClasts;rc++) { //begin r c loop fstream f i l e C v a l u e s . d a t " ) ; i f ( ! f i l e ) { cerr«"failed to open f i l e 'values-nc'\n"; exit(EXIT.FAILURE); } while (moveable == 0) { file.seekp(savedPos3); t r = ranl(idum); P o s i t i o n = int(tr*Number) ; //chosen c l a s t i s a random number scaled //by l i s t s i z e cout«"rc= " « r c « " tr= " « t r « " Position= "«Position«"\n"; file«"rc= " « r c « " tr= " « t r « " Position= "«Position«"\n"; delPtr = head; i f ( P osition > 1) { f o r (i=2; i<(Position+1); i++) { delPtr = delPtr->link; i f (delPtr == NULL) { cout«"ERR0R : delPtr POINTS TO NULL"«"\n"; file«"ERROR : delPtr POINTS TO NULL"«"\n"; } } } tr=ranl(idum); i f ((delPtr->r*2*tr) < C r i t e r i o n ) { moveable = 1; xtemprc = delPtr->x; ytemprc = delPtr->y; / / P r o b a b i l i t y of entrainment: l a r g e r //stones l e s s l i k e l y to be entrained. //coordinates are copied from l i s t / / i n t o the bedload. 65 rtemprc = delPtr->r; ltemprc = delPtr->1; g = gasdev(idum); vtemprc = g*Range; //get t r a j e c t o r y from gaussian d i s t . i f (ytemprc <= 0.0 kk vtemprc > 0.0) { vtemprc = - vtemprc; } i f (ytemprc >= Width && vtemprc < 0.0) { vtemprc = - vtemprc; } deleteClast = 1; cout«"xtemprc= "«xtemprc«" ytemprc= "«ytemprc<<" ltemprc= "«ltemprc « " rtemprc= "«rtemprc«" vtemprc= "«vtemprc«"\n"; file«"xtemprc= "«xtemprc«" ytemprc= " « y t e m p r c « " ltemprc= "«ltemprc « " rtemprc= "«rtemprc«" vtemprc= "«vtemprc«"\n"; i f (delPtr == head) { head = head->link; / / r e p o s i t i o n head and a l l other //pointers one node down. } else { prevPtr = head; while (prevPtr->link != delPtr) { prevPtr = prevPtr->link; • } prevPtr->link = delPtr->link; } i f (prevPtr->link == NULL) { endPtr = prevPtr; } delete delPtr; Number — ; cout«"Randomly entrain # "«Position«" " « N u m b e r « " nodes l e f t i n l i s t " « " i= " « i « " \ n " ; file«"Randomly entrain # "«Position«" " « N u m b e r « " nodes l e f t i n l i s t " « " i= " « i « " \ n " ; savedPos3=f i l e . tellgO ; f i l e . c l o s e ( ) ; e n t r a i n 0 ; } } moveable = 0; stop = 0; i f ( e x i t C l a s t == 1) { e x i t C l a s t = 0; e n t r a i n ( ) ; stop = 0; } } // end rc loop t o t a l C l a s t s = SeedNumber + input - throughput; . filebedO; fstream f i leC'values.dat"); i f ( ! f i l e ) { cerr<<"failed to open f i l e ' i n i t ' \ n " ; exit(EXIT JFAILURE); > file.seekp(savedPos3); c o u t « " rc= " « r c « " t o t a l clasts= "«totalClasts«"\n"; f i l e « " rc= "«rc«"total clasts= "«totalClasts«"\n"; c o u t « " input = " « i n p u t « " throughput= "«throughput«"\n"; f i l e « " input = " « i n p u t « " throughput= "«throughput«"\n"; c o u t « " adjusts= "<<adjust<<" rebounds= " « r e b o u n d « " vanished= "<<vanishedRock « " reversals= "<<reversal«"\n"; 66 file«" adjusts= "«adjust<<" rebounds1 savedPos3=f i l e . t e l l g O ; file.close(); "«rebound«"\n"; // end main } // void entrain() // This function projects a stone along a trajectory and determines which obstacle i t w i l l // impact, how i t will rebound and where i t will eventually rest. Stop must be zero to enter // this function and stop = 1 to exit the function. Detailed notes are at the beginning of each // numbered section. The stone's i n i t i a l position is xtemp and ytemp, i t s radius rtemp, and the // trajectory is vtemp in degrees where negative vtemp sends a stone counterclockwise down the // river and positive clockwise. A stone does not leave this function until i t is stopped at // which point its coordinates are determined and added to the end of the l i s t . Values are // written to f i l e values.dat. An anti-clogging rule has been written to avoid clogging up of // the upstream margin. When an obstacle is found, i f i t has an x position less than the radius // of the incoming clast, then the clast vanishes and vanishedRock counter is incremented. NodePtr NodePtr float xtemp, tempPtr; obstaclePtr; ytemp, rtemp, vtemp, ltemp, mileage; float oldDist, touchDistance, ratio; float a, b, hy, k, obx, oby; float A, B, C, G, M, Ob, Q, Z; float hnew, xnew, ynew, alpha, xlag, ylag; float xfree, yfree, hfree, xN, yN, IN, xS, yS, IS float xtemp2, ytemp2, vtemp2, ltemp2; int reboundedClastN, reboundedClastS, int reboundsN, reboundsS; int obstacle, skim, Case; int checkedNode ; shortmargin; //pointer to find obstacle clast //pointer to obstacle clast stone hits //temporary coordinates, radius, //angle, travel length, distance //travelled for entrained clast. //vector distance between new clast //and potential obstacle //vectors in the encounter triangles //angles in the encounter triangles //vectors for position of new clast //at encounter point //for boundary conditions //temporary coordinates and angle for //skimming clast. //flags. //counters: number of rebounds along //north and south margin, //flag: obstacle found; skimming stone //orientation of encounter triangle. //counts number of l i s t nodes checked. reboundsN = 0; reboundsS = 0; reboundedClastN = 0; reboundedClastS = 0; shortmargin = 0; skim = 0; while(stop == 0) { //stop condition fstream fileC'values.dat") ; i f (Ifile) { cerr<<"failed to open f i l e 'values-entrainO .'\n"; exit(EXIT.FAILURE); } file.seekp(savedPos3); // 1. The entrained stone can be a randomly entrained clast (when deleteClast = 1) // a clast rebounding off the northern and southern rigid boundaries, a clast which // has just skimmed an obstacle and continues to move downstream, a clast repicked // since the previous one would have clogged the upstream border, or a new clast // introduced on upstream margin in which case xtemp = 0, ytemp is drawn from a // uniform distribution and vtemp is drawn from a gaussian distribution, i f (reboundedClastN == 1) { //begin margin&skim loop ytemp = yN; xtemp = xN; vtemp = alpha; ltemp = IN; 67 reboundedClastN = 0; } else i f (reboundedClastS == 1) { ytemp = yS; xtemp = xS; vtemp = alpha; ltemp = IS; reboundedClastS = 0; } else i f (skim ==1) { ytemp = ytemp2; xtemp = xtemp2; vtemp = vtemp2; ltemp = ltemp2; skim = 0; i f (vtemp > 0 && ytemp == 0) { //to avoid errors vtemp = -vtemp; re v e r s a l ++; } else i f (vtemp < 0 && ytemp == Width) { vtemp = -vtemp; re v e r s a l ++; } } else i f (deleteClast ==1) { ytemp = ytemprc; xtemp = xtemprc; rtemp = rtemprc; vtemp = vtemprc; ltemp = ltemprc; deleteClast = 0; } else i f (vanish == 1){ //repick to avoid clogging upstream border. tr=ranl(idum); ytemp = tr*Width; xtemp = 0.0; g = gasdev(idum); rtemp = Lambda*exp(g) + Epsilon; / / r a d i i drawn from lognormal d i s t r i b u t i o n , while (rtemp < (D_75*0.5) II rtemp > (D_100*0.5)) { g = gasdev(idum); rtemp = Lambda*exp(g) + Epsilon; } g = gasdev(idum); vtemp = g*Range; D = 2*rtemp/10; L =1.77*(pow(l-loglO(D),1.35)); ltemp = L*L_D50surf; f i l e « " \ t " « t r « " \ t " « g « " \ n " ; c p u t « " \ t " « t r « " \ t " « g « " \ n " ; vanish =0; vanishedRock ++; y else { / / i n j e c t a new c l a s t at upstream margin. tr=ranl(idum); ytemp = tr*Width; xtemp = 0.0; g = gasdev(idum); rtemp = Lambda*exp(g) + Epsilon; while (rtemp < (D_75*0.5) II rtemp > (D_100*0.5)) { g = gasdev(idum); rtemp = Lambda*exp(g) + Epsilon; } D = 2*rtemp/10; L =1.77*(pow(l-loglO(D),1.35)); ltemp = L*L_D50surf; g = gasdev(idum); vtemp = g*Range; f i l e « " \ t " « t r « " \ t " « g « " \ n " ; 68 cout«"\t"«tr«"\t"«g«"\n"; input ++; } //end margin&skim loop tempPtr = head; obstaclePtr = head; obstacle =0; vanish = 0; checkedNode = 0; oldDist = 2000; // 2. Starting at the head of the l i s t , compare the stone's i n i t i a l position and // trajectory with each member of the l i s t positioned downstream from the obstacle // using tempPtr. If the stone encounters an obstacle and i f the distance to this // obstacle is less than what is was previously then the pointer obstaclePtr is // positioned at that point in the l i s t and the obstacle flag is set to one. while(tempPtr != NULL && vanish == 0) { //begin while l i s t i f ((tempPtr->x) > xtemp) { //begin forward loop, //only looks at stones in front of i t oby' = sqrt ((tempPtr->y-ytemp) * (tempPtr->y-ytemp)) ; hy = sqrt((tempPtr->x-xtemp)*(tempPtr->x-xtemp) + (tempPtr->y-ytemp)*(tempPtr->y-ytemp)); 0b = asin(oby/hy)*57.3; i f (tempPtr->y > ytemp) { 0b = - 0b; } touchDistance = rtemp + tempPtr->r; Q = (atan(touchDistance/hy))*57.3; i f ((fabs(vtemp - 0b) < Q) && hy < oldDist) {//if new clast w i l l touch the //obstacle and i f distance //gets closer obstaclePtr = tempPtr; oldDist = hy; obstacle = 1; tempPtr = tempPtr->link; i f (obstaclePtr->x < rtemp) { obstacle = 0; vanish = 1; //set flag to redraw to avoid upstream clog cout«"0bstacle too close to upstream margin so vanish. "«"\n"; f ile«"0bstacle too close to upstream margin so vanish. "«"\n"; } } else { tempPtr = tempPtr->link; } y else { tempPtr = tempPtr->link; } // end forward loop checkedNode ++; } // end while l i s t tempPtr = head; // reposition to head of l i s t to avoid deletion d i f f i c u l t i e s . // 3. The entire l i s t has been checked against the trajectory of the entrained stone // and no obstacle has been found (obstacle=0). Therefore we must check to see i f the // stone wi l l travel i t s maximum distance before hitting the north or south boundaries // or pass through the downstream margin of the bed. If the stone passes through, i t is // counted as throughput, stop is set to 1, and the entrain function is exited. If i t // hits the north or south boundaries the stone is given marginal coordinates xN, yN or // xS, yS and rebound flags are set to 1. These coordinates are sent back to part 1 where // the l i s t is checked for potential obstacles to the rebounder. i f (obstacle == 0 ft& vanish == 0) { //begin obstacle loop cout<<"No obstacle, boundary conditions"«"\n"; file«"No obstacle, boundary conditions"«"\n"; cout«"xtemp= "«xtemp«" ytemp= "«ytemp«" ltemp="«ltemp«" rtemp= " 69 «rtemp«" vtemp= "«vtemp«"\n"; file«"xtemp= "«xtemp«" ytemp= "«ytemp«" ltemp="«ltemp«" rtemp= " «rtemp«" vtemp= "«vtemp«"\n"; i f (vtemp <= 0) { //begin rebound loop yfree = Width - ytemp; xfree = xtemp + yfree/(tan(-vtemp*0.01745)); hfree = yfree/(sin(-vtemp*0.01745)); i f (ltemp <= hfree) { shortmargin = 1; stop = 1; xlag = ltemp*cos(-vtemp*0.01745); ylag = ltemp*sin(-vtemp*0.01745); } else i f ((xfree > Length) && (ltemp > hfree)) { throughput ++; stop = 1; //to get out of stop condition exitClast = 1; //flag to get a new stone injected into bed. cout«"throughput"«"\n"; f ile«"throughput"«"\n"«"\n"; } else { yN = Width; xN = xtemp + xfree; IN = ltemp - hfree; reboundedClastN = 1; reboundsN ++; rebound ++; alpha = 0 - vtemp; cout«"N rebounder"«"\n"; file«"N rebounder"«"\n"; } } else i f (vtemp > 0) { yfree = ytemp; xfree = xtemp + yfree/(tan(vtemp*0.01745)); hfree = yfree/(sin(vtemp*0.01745)); i f (ltemp < hfree) { shortmargin = 1; stop = 1; xlag = ltemp*cos(vtemp*0.01745); ylag = ltemp*sin(vtemp*0.01745); } else i f ((xfree > Length) && (ltemp > hfree)) { throughput ++; stop = 1; exitClast = 1; cout«"throughput"«"\n"; f ile«"throughput"«M\n"; } else { yS = 0.0; xS = xtemp + xfree; IS = ltemp - hfree; reboundedClastS = 1; reboundsS ++; rebound ++; alpha = 0 - vtemp; cout«"S rebounder"«"\n"; f ile«"S rebounder"«"\n"; } y else { cout«"ERR0R IN REBOUND L00P"«"\n"; file«"ERR0R IN REBOUND L00P"«"\n"; } //end rebound loop } else i f (obstacle ==1 && vanish == 0) { //continue obstacle loop // 4. The entire l i s t has been checked and an obstacle has been found. 70 // Trigonometry on the t r i a n g l e s formed by the stone and the obstacle give // parameters which w i l l be used i n part 5 to get new coordinates at point of // contact. C a p i t a l l e t t e r s are angles always converted to degrees and small // l e t t e r s are lengths of t r i a n g l e sides. Case keeps track of which set of // t r i a n g l e s i s being calculated. Reference diagram given i n f i g u r e 15 of // thesis text, obx = sqrt((obstaclePtr->x-xtemp)*(obstaclePtr->x-xtemp)); oby = sqrt((obstaclePtr->y-ytemp)*(obstaclePtr->y-ytemp)); cout«"xtemp= " « x t e m p « " ytemp= " « y t e m p « " ltemp= " « l t e m p « " rtemp= " « r t e m p « " vtemp= "«vtemp«"\n"; cout«"obstx= "«obstaclePtr->x«" obsty= "«obstaclePtr->y«" obrad= " «obstaclePtr->r«" obx= " « o b x « " oby= " « o b y « " \ n " ; file«"xtemp= " « x t e m p « " ytemp= " « y t e m p « " ltemp= " « l t e m p « " rtemp= " « r t e m p « " vtemp= "«vtemp«"\n"; file«"obstx= "«obstaclePtr->x«" obsty= "«obstaclePtr->y«" obrad= " «obstaclePtr->r«" obx= " « o b x « " oby= " « o b y « " \ n " ; hy = sqrt((obstaclePtr->x-xtemp)*(obstaclePtr->x-xtemp) + (obstaclePtr->y-ytemp)*(obstaclePtr->y-ytemp)); Ob = asin(oby/hy)*57.3; touchDistance = rtemp + obstaclePtr->r; i f (obstaclePtr->y > ytemp) { Ob = - Ob; } Q = (atan(touchDistance/hy))*57.3; cout«"hy= " « h y « " 0b= " « 0 b « " Q= " « Q « " touchDist= "«touchDistance«"\n" file«"hy= " « h y « " 0b= "<<0b«" Q= " « Q « " touchDist= "«touchDistance«"\n" Case = 0; i f ((vtemp<0.0 &ft 0b<0.0) && fabs(vtemp) < fabs(Ob)) {//start case loop, case C = 90.0 + Ob; Z = - Ob + vtemp; b = oby - obx*tan(-vtemp*0.017452); Case = 1; } else i f ((vtemp<0.0 && 0b<0.0) && fabs(vtemp) > fabs(Ob)) { //case 2 C = 90.0 - 0b; Z = - vtemp + 0b; b = obx*tan(-vtemp*0.017452) - oby; Case = 2; } else i f (vtemp>0.0 && 0b>0.0 ft& vtemp<0b) { //case 3 C = 90.0 - 0b; Z = 0b - vtemp; b = oby - obx*tan(vtemp*0.017452); Case = 3; } else i f (vtemp>0.0 && 0b>0.0 && vtemp>0b) { //case 4 C = 90.0 + 0b; Z = vtemp - Ob; b = obx*tan(vtemp*0.017452) - oby; Case = 4; } else i f (vtemp>0.0 && 0b<= 0.0) { //cases 5 C = 90.0 + Ob; Z = - Ob + vtemp; b = oby + obx*tan(vtemp*0.017452); Case = 5; } else i f (vtemp<0.0 && 0b>= 0.0) { //case 6 C = 90.0 - 0b; Z = 0b - vtemp; b = oby + obx*tan(-vtemp*0.017452); Case = 6; } else { cout«"ERR0R IN CASE L00P"«"\n"; f ile«"ERR0R IN CASE L00P"«"\n"; } //end case loop cout«"C= " « C « " Z= " « Z « " b= " « b « " case= "«Case«"\n"; 71 file«"C= " « C « " Z= " « Z « " b= " « b « " case= "«Case«"\n"; G = 180.0 - Z - C; B = (asin((sin(G*0.017452))*b/touchDistance))*57.3; A = 180.0 - B - G; a = touchDistance*sin(A*0.017452)/sin(G*0.017452); k = hy*sin(C*0.017452)/sin(G*0.017452); cout«"G= " « G « " B= " « B « " A= " « A « " a= " « a « " k= "«k«"\n"; file«"G= " « G « " B= " « B « " A= " « A « " a= " « a « " k= "«k«"\n"; hnew = k - a; obstaclePtr = head; //reposition to head of l i s t to avoid deletion problems // 5. Using the triangles solved above, the coordinates of the center of the // stone are determined at its point of contact with the obstacle. If the stone // happens to end up slightly out of the digital reach then i t is nudged back to // the zero border, a message signals this adjustment and adjust is incremented. // At the point of contact, the stone will rebound at an angle equal to the // incident angle as measured from the normal to the encounter surface. If the // rebound angle is greater than FULCRUM degrees the stone stop's and is // deposited at the impact position. If the angle is less than FULCRUM then the // stone skims the obstacle, stop flag is zero so entrainO cannot be exited. // The new coordinates, trajectory, and remaining travel distance, xtemp2, ytemp2, // vtemp2 and ltemp are sent back to part 1. i f ((Case ==1) II (Case ==2)) { //begin skim/stop loop ynew = hnew*sin(-vtemp*0.017452); xnew = hnew*cos(-vtemp*0.017452); xtemp2 = xtemp + xnew; ytemp2 = ytemp + ynew; i f (Case ==1) { vtemp2 = 180 - 2*B + vtemp; } else { vtemp2 = vtemp -180 + 2*B; } } else i f ((Case ==3) II (Case ==4)) { ynew = hnew*sin(vtemp*0.017452); xnew = hnew*cos(vtemp*0.017452); xtemp2 = xtemp + xnew; ytemp2 = ytemp - ynew; i f (Case == 3) { vtemp2 = -180 + 2*B + vtemp; } else { vtemp2 = 180 -2*B + vtemp; } } else i f ((Case ==5) II (Case ==6)) { i f (Case == 5) { ynew = hnew*sin(vtemp*0.017452); xnew = hnew*cos(vtemp*0.017452) ; xtemp2 = xtemp + xnew; ytemp2 = ytemp - ynew; vtemp2 = 180 - 2*B + vtemp; } else { ynew = hnew*sin(-vtemp*0.017452); xnew = hnew*cos(-vtemp*0.017452); xtemp2 = xtemp + xnew; ytemp2 = ytemp + ynew; vtemp2 = vtemp -180 + 2*B; } } else { cout«"ERR0R IN 0UTERSKIM/ST0P L00P"«"\n"; file«"ERR0R IN OUTER SKIM/STOP L00P"«"\n"; y //end skim/stop loop i f (vtemp2 < FULCRUM && vtemp2 > -FULCRUM) { //skim/stop assignment skim = 1; 72 } else i f (vtemp2 > FULCRUM I I vtemp2 < -FULCRUM) { stop = 1; } else { cout«"ERRQR IN SKIM/STOP ASSIGNMENT"«"\n"; file«"ERROR IN SKIM/STOP ASSIGNMENT"«"\n"; } i f (xtemp2 < 0.0) { //adjustments i f stone i s s l i g h t l y out of bounds xtemp2 = 0.0; adjust ++; cout«"x adjustment made to 0."«"\n"; file«"x adjustment made to 0."«"\n"; } i f (ytemp2 < 0.0) { ytemp2 = 0.0; adjust ++; cout«"y adjustment made to 0."«"\n"; file«"y adjustment made to 0."«"\n"; } else i f (ytemp2 > Width) { ytemp2 = Width; adjust ++; cout<<"y adjustment made to width. "«"\n"; file«"y adjustment made to width. "«"\n"; } ltemp2 = ltemp - (hy - touchDistance); //re-evaluate allowable t r a v e l distance cout«"hnew= " « h n e w « " ynew= " « y n e w « " xnew= " « x n e w « " xtemp2= "«xtemp2 « " ytemp2= " « y t e m p 2 « " ltemp2= " « l t e m p 2 « " vtemp2= "«vtemp2«"\n"; i f (skim ==1) { cout«" skim"«"\n"; } else {cout«" stop"«"\n"; } file«"hnew= " « h n e w « " ynew= " « y n e w « " xnew= " « x n e w « " xtemp2= "«xtemp2 « " ytemp2= " « y t e m p 2 « " ltemp2= " « l t e m p 2 « " vtemp2= "«vtemp2«"\n"; i f (skim ==1) { f i l e « " skim"«"\n"; > else {file«" stop"«"\n"; } // 6. The stone ends up stopping therefore a node i s added to the end of the // l i s t i n to which the stone's coordinates are written. Number, representing // the length of the l i s t , i s incremented. The entrain function i s completed // and exited. i f (stop ==1|| ltemp2 < 0.0) { //begin stop portion of obstacle loop endPtr->link = new NodeType; endPtr->link->link = NULL; endPtr = endPtr->link; i f (ltemp2 < 0.0 && vtemp < 0.0) { endPtr->x = xtemp + ltemp*cos(-vtemp*0.01745); endPtr->y = ytemp +ltemp*sin(-vtemp*0.01745); stop = 0; } else i f (ltemp2 < 0.0 && vtemp >= 0.0) { endPtr->x = xtemp + ltemp*cos(vtemp*0.01745); endPtr->y = ytemp - ltemp*sin(vtemp*0.01745); stop = 0; } else i f (stop ==1) { endPtr->x = xtemp2; endPtr->y = ytemp2; } stop = 1; endPtr->r = rtemp; D = 2*rtemp/10; L =1.77*(pow(l-loglO(D),1.35)); endPtr->l = L*L_D50surf; //reset ltemp endPtr->v = 1.0; Number ++; cout«"x= "«endPtr->x«" y="«endPtr->y«" ltemp="«endPtr->l«" r=" «endPtr->r«" v="«endPtr->v«"\n"; 73 file«"x= "«endPtr->x«" y="«endPtr->y«" ltemp="«endPtr->l«" r=" «endPtr->r«" v="«endPtr->v«"\n"; cout«"\n"; file«"\n"; } //end stop portion of obstacle loop } //end obstacle loop i f (shortmargin == 1) { i f ((xtemp + xlag)<= LengthH //stopped before hitting boundary endPtr->link = new NodeType; endPtr->link->link = NULL; endPtr = endPtr->link; endPtr->x = xtemp + xlag; i f (vtemp<=0){ ylag = -ylag; } endPtr->y = ylag; endPtr->r = rtemp; D = 2*rtemp/10; L =1.77*(pow(l-loglO(D),1.35)); endPtr->l = L*L_D50surf; //reset ltemp endPtr->v = 1.0; stop = 1; //to exit entrainO cout«"x= "«endPtr->x«" y="«endPtr->y«" ltemp="«endPtr->l«" r=" «endPtr->r«" v="«endPtr->v«"\n"; file«"x= "«endPtr->x«" y="«endPtr->y«" ltemp="«endPtr->l«" r=" «endPtr->r«" v="«endPtr->v«"\n"; cout«"\n"; file«"\n"; y else i f ((xtemp + xlag)>Length){ //passed through downstream boundary throughput**; exitClast=l; cout«"throughput"«"\n"; f ile«"throughput"«"\n"; } else { cout«"error at downstream margin"<<"\n"; file<<"error at downstream margin"«"\n"; } } savedPos3=f i l e . t e l l g O ; f ile.closeO ; } //end stop condition void seedbed() // This function- generates a river bed covered with SeedFraction non-touching stones whose // positions are drawn from a uniform distribution in x and y, whose radii are drawn from a // lognormal distribution and whose maximum travel distance, 1, is calculated. { NodePtr checkPtr; int overlap; float xtemp, ytemp, rtemp, xdiff, ydiff; float root, Distance, touchDistance; SeedNumber = int((Width*Length*SeedFraction)/(D_95*0.5*D_95*0.5*3.14159)); Number = 0; head = new NodeType; head->link = NULL; endPtr = head; //positions pointer to top of l i s t tr=ranl(idum); / / f i l l in f i r s t item in li s t = f i r s t clast endPtr->x = tr*Length; tr=ranl(idum); 74 endPtr->y = tr*Width; g = gasdev(idum); endPtr->r = Lambda*exp(g) + Epsilon; while ( endPtr->r < (D_75*0.5) I IendPtr->r > (D_100*0.5)) { g = gasdev(idum); endPtr->r = Lambda*exp(g) + Epsilon; } endPtr->v = 0.0; //seeded clast signifier D = 2*endPtr->r/10; L =1.77*(pow(l-loglO(D),1.35)); endPtr->l = L*L_D50surf; endPtr->link = new NodeType; endPtr->link->link=NULL; endPtr = endPtr->link; Number ++; overlap = 0; while (Number < SeedNumber) { tr=ranl(idum); / / f i l l in f i r s t item in l i s t = f i r s t clast xtemp = tr*Length; tr=ranl(idum); ytemp = tr*Width; g = gasdev(idum); rtemp = Lambda*exp(g) + Epsilon; while (rtemp < (D_75*0.5) II rtemp > (D_100*0.5)) { g = gasdev(idum); rtemp = Lambda*exp(g) + Epsilon; } checkPtr = head; overlap = 0; while (checkPtr->link !=NULL) { xdiff = checkPtr->x - xtemp; ydiff = checkPtr->y - ytemp; root = (xdiff*xdiff) + (ydiff*ydiff); Distance = sqrt (root); touchDistance = checkPtr->r + rtemp; i f (Distance < touchDistance) { checkPtr = endPtr; overlap =1; y else i f (Distance >= touchDistance) { checkPtr=checkPtr->link; } else { cout«"ERR0R in DISTANCE L00P"«"\n"; } } i f (overlap == 0) { endPtr->x = xtemp; endPtr->y = ytemp; endPtr->r = rtemp; endPtr->v = 0.0; D = 2*endPtr->r/10; L =1.77*(pow(l-loglO(D),1.35)); endPtr->l = L*L_D50surf; Number ++; } i f (overlap == 0 && Number < SeedNumber) { endPtr->link = new NodeType; endPtr->link->link=NULL; endPtr = endPtr->link; } } cout«"SeedNumber= "«SeedNumber«" Number= "«Number«"\n"; } 75 // void fileSeedBed () // Prints the i n i t i a l conditions to a f i l e resulting in a table of x,y coordinates, radii and // travel length. It opens f i l e (init.dat) and uses the pointer visitPtr to write the contents / / o f the l i s t to f i l e . Includes error handling and returns visitPtr to head after use to // avoid problems when deleting nodes. { NodePtr visitPtr; streampos savedPosl=0; //for f i l e "init.dat" fstream f i l e ( " i n i t . d a t " ) ; i f (!file) { cerr<<"failed to open f i l e 'init'\n"; exit(EXIT FAILURE); } file.seekp(savedPosl); visitPtr = head; while (visitPtr != NULL) { f ile«visitPtr->x«" "«visitPtr->y«" "«visitPtr->r«"\n" ; visitPtr = visitPtr->link; } savedPosl= f i l e . t e l l g O ; visitPtr = head; f i l e . closeO ; cout«SeedNumber«"\n"; ' } // void filebed () // Prints the f i n a l conditions to f i l e resulting in a table of x,y coordinates, radii and travel // length. It opens f i l e bed.dat and uses visitPtr to v i s i t each node in the l i s t and write // contents to f i l e . Includes error handling and returns visitPtr to head after use. { NodePtr visitPtr; streampos savedPos2=0; //for f i l e "bed.dat" fstream fileC'bed.dat") ; i f (Ifile) { cerr«"failed to open f i l e 'bed"'«"\n"; exit(EXIT.FAILURE); } file.seekp(savedPos2); visitPtr = head; while (visitPtr != NULL) { f ile«visitPtr->x«" "«visitPtr->y«" "«visitPtr->r«"\n"; visitPtr = visitPtr->link; } savedPos2=f i l e . t e l l g O ; visitPtr = head; file. c l o s e O ; cout<<totalClasts; } // #define IA 16807 #define IM 2147483647 #define AM (1.0/IM) #define IQ 127773 #define IR 2836 #define NTAB 32 #define NDIV (1+(IM-1)/NTAB) #define EPS 1.2e-7 #define RNMX (1.0-EPS) float rani(long *idum) // Random number generator RANI from Numerical Recipes in C. Returns a uniform random deviate // between (0.0,1.0). Call with idum a negative integer to i n i t i a l i z e and do not change 76 // idum between successive deviates in a sequence. RNMX should approximate the largest // floating value that is less than 1. i int j ; long k; static long iy=0; static long iv[NTAB]; float temp; i f (*idum<=OII!iy) { i f (-(*idum)<l) *idum=l; else *idum=-(*idum); for (j=NTAB+7;j>=0;j—) { k=(*idum)/IQ; *idum=IA*(*idum-k*IQ)-IR*k; i f (*idum,0) *idum+=IM; i f (j<NTAB) iv[j]=*idum; } iy=iv[0] ; } k=(*idum)/IQ; *idum=IA*(*idum-k*IQ)-IR*k; i f (*idum<0) *idum+=IM; j=iy/NDIV; iy=iv[j] ; iv[j]=*idum; i f ((temp=AM*iy)>RNMX) return RNMX; else return temp; } #undef IA #undef IM #undef AM #undef IQ #undef IR #undef NTAB #undef NDIV #undef EPS #undef RNMX // float gasdev(long *idum) // Returns normally distributed deviate with zero mean and unit variance using rani(idum) // as the source of uniform deviates. From Numerical Recipes in C. { float rani(long *idum); static int iset = 0; static float gset ; float fac, rsq, v l , v2; i f (iset == 0) { do { vl=2.0*ranl(idum)-1.0; v2=2.0*ranl(idum)-1.0; rsq=vl*vl+v2*v2; } while (rsq >= 1.01 I rsq == 0.0); fac=sqrt(-2.0*log(rsq)/rsq); gset=vl*fac; iset=l; return v2*fac; } else { iset=0; return gset; } 77 Bibliography [1] A h n e r t , F . , 1994. M o d e l l i n g the development of non-periglacial sorted nets. Catena 23, 43-63. [2] A l l e n , J . R . L . , 1983. A simplified cascade model for transverse stone-ribs in gravelly streams. Proceedings, Royal Society of London A,385, 253-266. [3] A n d e r s o n , R . S . , 1990. E o l i a n ripples as examples of self-organization in geomorphological systems. Earth-Science Reviews 29, 77-96. [4] B luck , B . J . , 1987. Bedforms and clast size changes in gravel-bed rivers, in R ichards , K . J . ( E d . ) , R iver channels; environment and process. Special Publication, Institute of British Geographers 18, 159-178. [5] Brayshaw, A . C . , 1984. Characterist ic and origin of cluster bedforms in coarse-grained al luvia l channels, in Kos ter , E . H . , and Steel, R . J . (Eds . ) , Sedimentology of Gravels and Conglomerates Canadian Society of Petroleum Geologists Memoir 10, 77-85. [6] Brayshaw, A . C . , 1985. B e d microtopography and entrainment thresholds in gravel-bed rivers. Geological Society of America Bulletin 96, 218-223. [7] Brayshaw, A . C . , Frostick, L . E . , and I. Re id , 1983. T h e hydrodynamics of particle clusters and sediment entrainment in coarse al luvial channels. Sedimentology 30, 137-143. [8] C h u r c h , M . and M . A . Hassan, 1992. Size and distance of travel of unconstrained clasts on a streambed. Water Resources Research 28 (1), 299-303. [9] C h u r c h , M . , Wolcot t , J . F . , and W . K . Fletcher, 1991. A test of equal mobi l i ty in f luvial sediment transport: behavior of the sand fraction. Water Resources Research 27, (11), 2941-2951. [10] C h u r c h , M . 1996. Personal communicat ion. [11] Cressie, N . A . C . , 1993. Statistics for Spatial Data, Revised edition, J o h n W i l e y & Sons, Inc. New Y o r k . [12] D a l C i n , R . , 1968. "Pebble Clusters": their origin and uti l ization in the study of paleocurrents. Sedi-mentary Geology 2, 233-241. [13] D e Jong , C , 1991. A reappraisal of the significance of obstacle clasts in cluster bedform dispersal . Earth Surface Processes and Landforms 16, 737-744. 78 [14] D e Jong , C . and P. Ergenzinger, 1995. T h e interrelations between mounta in valley form and r iver-bed arrangement, in H i c k i n , E . J . (Ed . ) , River Geomorphology, J o h n W i l e y & Sons. [15] Gus tavson , T . C . , 1974. Sedimentation on gravel outwash fans, M a l a s p i n a Glac ier Fore land , Alaska . Journal of Sedimentary Petrology 44, 374-389. [16] Hal let , B . , 1990. Spat ia l self-organization in geomorphology: from periodic bedforms a n d patterned ground to scale-invariant topography. Earth-Science Reviews 29, 57-75. [17] Hassan , M . A . , 1993. B e d material and bedload movement in two ephemeral streams. Special Publications of the International Association of Sedimentologists 17, 37-49. [18] Hassan , M . A . and M . C h u r c h , 1992. T h e movement of indiv idual grains on the s treambed in B i l l i , P . , Hey. R . D . , T h o m e , C . R . , and P. Taccon i (Eds. ) , Dynamics of Gravel-bed Rivers, J o h n W i l e y & Sons L t d . [19] Koster , E . H . , 1978. Transverse Ribs: T h e i r characteristic origin and paleohydraul ic significance. Cana-dian Society of Petroleum Geologists Memoirs 5, 161-186. [20] L a n d r y , W . and B . T . Werner, 1994. C o m p u t e r simulations of self-organized wind ripple patterns. Physica D, 238-260. [21] L a r o n n e , J . B . and M . A . C a r s o n , 1978. Interrelationships between bed morphology and bed-material transport for a small , gravel-bed channel. Sedimentology 23, 67-85. [22] M a r t i n i , L P . , 1977. Grave l ly flood deposits of Irvine Creek, Ontar io , C a n a d a . Sedimentology 24, 603^622. [23] M c D o n a l d , B . C . and I. Banerjee, 1971. Sediments and bedforms on a bra ided outwash p la in . Canadian Journal of Earth Sciences 8, 1282-1301. [24] M c D o n a l d , B . C . and T . J . Day, 1978. A n experimental flume study on the formation of transverse ribs. Geological Survey of Canada Paper 78-1A, 441-451. [25] Nowell , A . R . and M . C h u r c h , 1979. Turbulent flow in a depth-l imited boundary layer. Journal of Geo-physical Research 84, 4816-4824. [26] Press, W . H . , Teukolsky, S . A . , Vetterl ing, W . T . and B . P . Flannery , 1994. Numerical Recipes in C. Second E d i t i o n . Cambr idge Univers i ty Press. [27] R e i d , I and L . E . Frostick, 1984. Part ic le interaction and its effect on the thresholds of in i t ia l and final bedload mot ion in coarse al luvial channels, in Koster , E . H . , and Steel, R . J . (Eds . ) , Sedimentology of Gravels and Conglomerates Canadian Society of Petroleum Geologists Memoir 10, 61-68. [28] R e i d , I. and M . A . Hassan, 1992. T h e influence of microform bed roughness elements on flow and sediment transport in gravel bed rivers. Earth Surface Processes and Landforms 15, 739-750. 79 [29] Ross i , R . E . , M u l l a , D . J . , Journel , A . G . , a n d E . H . Franz , 1992. Geostatist ical Tools for M o d e l i n g a n d Interpreting Ecologica l Spat ial Dependence, Ecological Monographs, 62,(2), 217-314. [30] Weatherby, H . , 1996. Unpubl i shed M . S c . thesis, Department of E a r t h a n d O c e a n Science, Univers i ty of B r i t i s h C o l u m b i a . [31] Werner , B . T . , 1995. E o l i a n dunes: computer simulation and attractor interpretation. Geology 23 (12), 1107-1110. [32] Werner , B . T . and T . M . F i n k , 1993. Beach cusps as self-organized patterns. Science 260, 968-971. [33] Werner , B . T . and B . Hallet , 1993. Numer ica l s imulation of self-organized stone stripes. Nature 361, 142-145. [34] W h i t e , C M . , 1940. T h e equi l ibr ium of grains on the bed of a stream. Proceedings of the Royal Society Series A, 174, 322-338. 80
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Computer simulations of cobble structure on a gravel...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Computer simulations of cobble structure on a gravel stream bed Tribe, Selina 1996
pdf
Page Metadata
Item Metadata
Title | Computer simulations of cobble structure on a gravel stream bed |
Creator |
Tribe, Selina |
Date Issued | 1996 |
Description | In Harris Creek, a gravel bed river in southcentral British Columbia, the large stones protruding from the surface armour of the bed are arranged as curvilinear ridges oriented at varying angles to the flow, defining a net-like surface structure. A computer simulation, using various rules governing the entrainment and deposition of stones, was written to reproduce these structures. The working hypothesis for the simulation is that the patterns are the result of the mutual interaction of stones as they are mobilized into the flow and collide with downstream obstacles. An opposing hypothesis would be that the pattern being imposed on the stones by the flow. A statistical algorithm was developed to identify cobble ridges and their orientations and to allow a comparison between the simulation results and the real example. Simulations using different sized stones in a bed with 10% of its area covered by stones produce flow-parallel, linear clusters resembling cluster bedforms. Simulations using a more densely populated bed, up to 25% covered area, produce a variety of gravel bedforms including cluster bedforms, transverse ribs, some trending at oblique orientations, and arcuate ridges. These results indicate that the magnitude of the stone population influences the suite of bedforms that emerge. Cluster bedforms may be the result of obstacles to the downstream movement of stones in a sparsely populated bed. Furthermore, cluster bedforms may be nuclei from which oblique ridges develop out of excess sediment. Finally, simulations in which large stones are rarely entrained and travel shorter distances than small stones, produce structures very similar to those found in Harris Creek. |
Extent | 7653961 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-02-17 |
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.0099040 |
URI | http://hdl.handle.net/2429/4663 |
Degree |
Master of Science - MSc |
Program |
Geography |
Affiliation |
Arts, Faculty of Geography, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1996-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1996-0544.pdf [ 7.3MB ]
- Metadata
- JSON: 831-1.0099040.json
- JSON-LD: 831-1.0099040-ld.json
- RDF/XML (Pretty): 831-1.0099040-rdf.xml
- RDF/JSON: 831-1.0099040-rdf.json
- Turtle: 831-1.0099040-turtle.txt
- N-Triples: 831-1.0099040-rdf-ntriples.txt
- Original Record: 831-1.0099040-source.json
- Full Text
- 831-1.0099040-fulltext.txt
- Citation
- 831-1.0099040.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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0099040/manifest