UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Episodic processes in fluvial landscape evolution 2000

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
ubc_2000-0377.pdf [ 9.23MB ]
Metadata
JSON: 1.0089851.json
JSON-LD: 1.0089851+ld.json
RDF/XML (Pretty): 1.0089851.xml
RDF/JSON: 1.0089851+rdf.json
Turtle: 1.0089851+rdf-turtle.txt
N-Triples: 1.0089851+rdf-ntriples.txt
Citation
1.0089851.ris

Full Text

EPISODIC PROCESSES IN FLUVIAL L A N D S C A P E EVOLUTION by Simon Dadson B.A. (Hons), University of Oxford, 1998 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF Master of Science in T H E F A C U L T Y OF G R A D U A T E STUDIES (Department of Geography) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y OF BRITISH C O L U M B I A September 2000 © Simon Dadson, 2000 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of kEOfrWH 7 The University of British Columbia Vancouver, Canada D a t e (o S&PTCMOK XQOQ DE-6 (2/88) Abst rac t The fractal dimension provides a succinct summary of topographic intricacy and a measure of landscape roughness. The work presented in this thesis is an attempt to examine the spatial variation of fractal properties of topography. It is hypothesized that distinct zones, or domains, may be identified in the landscape, within which the fractal dimension is constant, owing to the dominance of a characteristic geomorphic process. The variogram method (structure-function analysis) is used to measure the fractal properties of digital elevation models of the Queen Charlotte Islands and small tributaries of the Capilano River (in the Southern Coast Mountains). Process do- mains are identified using upslope-area criteria, and are shown to have distinct fractal dimensions. In the Capilano basin, the fractal dimension is found to be higher on upper slopes; whereas on the Queen Charlotte Islands, it is higher on lower slopes. This difference is explained with reference to the details of specific surficial processes operating in each of the areas, and to their different glacial histories. A multiple-process landscape evolution model is presented to explore the link between geomorphic process and fractal form. The model is carefully calibrated using field data from several recent surveys within the study areas. Non-linear diffusion is used to represent surficial slope failure; fluvial processes are represented using a simplified Bagnold equation; and deep-seated rock slides are introduced at random on the steepest, highest slopes. The model results confirm the smoothing effect of diffusion, which lowers the fractal dimension of the resulting topography. Structures introduced by rock-slope failures introduce visible short-term complexity on hillslopes; however, this is quickly incorporated into the channel network so does not appear to increase the fractal dimension of areas defined as hillslopes using the threshold slope- area criterion. Further work is suggested to improve modelling of fluvial processes in steep, headwater channels, and to investigate the details of hillslope-channel coupling at higher spatial resolution. ii Contents Abstract ii List of Tables vii List of Figures viii Acknowledgements xi Chapter 1 Introduction 1 1.1 Aims 1 1.1.1 Spatial Organization of Landscape Processes 2 1.2 Fractal Models of Fluvial Landscapes 5 1.2.1 Fractal Properties of Topography 5 1.3 Objectives 11 Chapter 2 Surface Processes and Tectonic History in the Queen Char- lotte Islands and Southern Coast Mountains of British Columbia 12 2.1 Location 12 2.2 Tectonic background and Geological History 16 2.2.1 Tectonics and Topography 16 2.3 Queen Charlotte Islands 17 2.3.1 Bedrock Geology 18 ii i 2.3.2 Glacial Legacy 18 2.3.3 Contemporary Processes 19 2.4 G V R D Capilano Watershed . . . 20 2.4.1 Bedrock Geology 20 2.4.2 Tectonic Uplift 21 2.4.3 Glacial Legacy 21 2.4.4 Contemporary Processes 22 2.5 Summary 23 Chapter 3 Sediment Flux due to Episodic Events 24 3.1 Process rates 24 3.2 Sources of Data 24 3.3 Surficial Slides: magnitude-frequency 25 3.3.1 Results 25 3.3.2 Power-law Scaling and Self-organized Criticality 26 3.3.3 Power-law Scaling and Sampling Bias 28 3.3.4 Finite Scaling Effects 28 3.3.5 G V R D Slides 31 3.3.6 Implications 32 3.4 Slope and Volumetric Transport Rate 33 3.4.1 Approach 33 3.4.2 Parameterization 33 3.4.3 Discussion 36 3.5 Rock Slides 36 3.5.1 Rock Slide Inventory 37 3.5.2 Discussion 37 iv 3.6 Summary 39 Chapter 4 Fractal properties of topography 40 4.1 Geomorphic Processes and Fractal Landscapes 40 4.2 Fractal Properties of Landscape 41 4.3 Measuring Fractal Properties of Landscape 43 4.4 Process-domains and the Fractal Dimension 44 4.5 Results 46 4.5.1 G V R D Capilano Watershed 46 4.5.2 Queen Charlotte Islands 60 4.6 Discussion 69 Chapter 5 Synthesis: Process-Based Modelling with Episodic Events 72 5.1 Introduction 72 5.2 Constructing Testable Landscape Models 73 5.3 Framework for Landscape Modelling 74 5.3.1 Continuum Model of Landscape Evolution 75 5.3.2 Parameters 78 5.3.3 Canonical Scaling 80 5.3.4 Boundary Conditions 81 5.3.5 Numerical Schemes 81 5.3.6 Stability 82 5.4 Results: Preliminary Runs 83 5.5 Results: Detailed Model Runs 88 5.5.1 Lembke Creek 89 5.5.2 Gregory Creek 94 5.6 Discussion 95 v Chapter 6 Conclusion 98 6.1 Practical Implications 99 6.2 Theoretical Implications 100 Bibliography 102 Appendix A Computation of Fractal Dimension 112 A . l Slope-Area Plot Constructions 112 A. 2 Construction of Hillslope Variograms 116 Appendix B Derivation of Landscape Modelling Equations 132 B. l Continuity Equation 132 B.2 Fluvial Transport Rule 133 B.3 erode2d.cc 135 vi List of Tables 2.1 Sediment Flux due to Landsliding on Queen Charlotte Islands (data from Rood 1984) 20 3.1 Slope-Volume Classes for G V R D Landslides 34 3.2 Parameters for Non-Linear Slope-Volume Relation 34 4.1 G V R D Capilano Channel Network Fractal Dimensions 47 4.2 G V R D Capilano Hillslope Fractal Statistics 59 4.3 Drainage Basin Characteristics and Channel Network Fractal Dimen- sions: Queen Charlotte Islands 63 4.4 Queen Charlotte Islands Hillslope Fractal Statistics 68 5.1 Parameter Estimates for Landscape Evolution Model 79 5.2 Parameters Used in Model Runs 88 5.3 Fractal Dimensions of Model Results 90 vii List of Figures 1.1 Continuum of Geomorphic Processes (after Montgomery and Foufoula- Georgiou, 1993). The thick solid line is typical slope-area plot for a steepland landscape; the thick broken line is a typical plot for a lowland area 3 1.2 Simulated Fractal Landscapes: fractional Brownian surfaces (Klinken- berg and Goodchild, 1993) 6 1.3 Moment statistics of a land-surface (part of Daniels Creek in the Capi- lano Watershed; see Chapter Four). The lower of the two plots shows a linear moment-order relation, indicating that this particular landscape is a mono-fractal 9 2.1 Location of the two Study Areas (also showing principal tectonic fea- tures, described in this chapter) 13 2.2 Queen Charlotte Islands (study basins are described here and in Chap- ter Four) 14 2.3 G V R D Capilano Watershed (study basins are described here and in Chapter Four) 15 3.1 Queen Charlotte Islands: Landslide Magnitude Frequency 26 3.2 Extreme Event Distribution 29 viii 3.3 Resolution Function 30 3.4 Landslide Magnitude Frequency modified by Resolution Function . . . 30 3.5 G V R D Slope-Volume Relation 35 3.6 Rock Slide Magnitude Frequency (data from Savigny 1991) 38 4.1 Zones for Computation of D 48 4.2 Slope-Area Plots: G V R D Daniels Creek 49 4.3 Slope-Area Plots: G V R D Sisters Creek 50 4.4 Slope-Area Plots: G V R D Enchantment Creek 51 4.5 Hillslope Variograms: G V R D Daniels Creek (D l - D2) 53 4.6 Hillslope Variograms: G V R D Daniels Creek (D3 - D4) 54 4.7 Hillslope Variograms: G V R D Sisters Creek (SI - S2) 55 4.8 Hillslope Variograms: G V R D Sisters Creek (S3 - S4) 56 4.9 Hillslope Variograms: G V R D Enchantment Creek ( E l - E2) 57 4.10 Hillslope Variograms: G V R D Enchantment Creek (E2 - E3) 58 4.11 Slope-Area Plots: Queen Charlotte Islands (Rennel Sound Area) . . . . 61 4.12 Slope-Area Plots: Queen Charlotte Islands (Skidegate Channel Area) . 62 4.13 Hillslope Variograms: Queen Charlotte Islands (Rennel Sound: Riley and Gregory Creeks) 64 4.14 Hillslope Variograms: Queen Charlotte Islands (Rennel Sound: Moun- tain and Hangover Creeks) 65 4.15 Hillslope Variograms: Queen Charlotte Islands (Skidegate Channel: Ar- mentieres and Government Creeks) 66 4.16 Hillslope Variograms: Queen Charlotte Islands (Skidegate Channel: Ja- son and Security Creeks) 67 5.1 Artificial Model Surface after lOOka 86 ix 5.2 Model Rock Slides: Magnitude-Frequency (100 year slide interval; N = 1000) 87 5.3 Model Rock Slides: Magnitude-Frequency (400 year slide interval; N = 250) 87 5.4 Lembke Creek: 10 ka no rock slides 91 5.5 Lembke Creek: 10 ka with rock slides 92 5.6 Gregory Creek: 10 ka with rock slides 96 x Acknowledgements I am grateful for the assistance and encouragement given to me by Michael Church, who supervised this work throughout. Olav Slaymaker made several suggestions during the for- mulation of ideas for this work and gave friendly encouragement. Brian Klinkenberg provided help with GIS, and invaluable discussions of fractal concepts. Garry Clarke provided excellent help with the formulation and implementation of the landscape evolution model described in Chapter Five; this numerical model has benefitted from insights obtained by reference to the recent work of Yvonne Martin (acknowledged in the text), which have greatly speeded its implementation. I am grateful to Russell White for sharing some of the benefit of his field experience in the Capilano watershed. Darren Ham helped with A R C / I N F O and ArcView GIS software. Catherine Griffiths gave several computer graphics tutorials. Vincent Kujala generously pro- vided access to computing facilities and gave assistance with the I ^ T E X typesetting software. The analyses presented here would have been impossible had not the following people and organizations provided access to data: Tony Cheong (Queen Charlotte Islands D E M data); Derek Bonin, Dave Dunkley, Lome Gilmour (Greater Vancouver Regional District D E M and landslide data); Ken Rood (Queen Charlotte Islands landslide data); Wayne Savigny (Fraser Valley landslide data). This work was completed while the author was recipient of the Wendy Fan Memorial Scholarship and a University Graduate Fellowship at the University of British Columbia. Finally, I am grateful to Elaine Cho and the office staff, who made the University easier to deal with, and to my family and Emma McKenzie for their support. S I M O N D A D S O N The University of British Columbia August 2000 xi Chapter 1 Introduct ion 1.1 A i m s Studies of fluvial landscape evolution have demonstrated that fluvial networks create fractal landscapes1 (see Rodriguez-Iturbe and Rinaldo, 1997, for a recent, comprehen- sive review). The intuitive justification for this result is that branching river networks dissect landscapes over a wide range of scales. The effects of hillslope processes such as rock slides, shallow landslides and debris flows on the fractal dimension of the land- scape have been less well studied. The aim of the present work is to determine the extent to which hillslope processes modify the effects of dissection by fractal fluvial networks. In order to achieve this goal, the rates and spatio-temporal variability of episodic, structure-generating hillslope processes are examined in conjunction with an analysis of the spatial variation of the fractal dimension in different process-domains within the landscape. Finally, an empirically calibrated model of fluvial and hillslope processes is used to make explicit the links between episodic processes and the fractal dimension of 1A preliminary, non-technical definition of a fractal is a shape that is composed of parts of itself, with each part in some way similar to the whole. A fractal landscape is therefore a landscape that shows similar variability at all spatial scales. This definition is elaborated in the following sections. 1 topography. With these objectives in mind, there follows some theoretical justification for the link between geomorphic processes and fractal landscapes. 1.1.1 Spatial Organization of Landscape Processes The properties of topography are determined by geomorphic processes; it is therefore reasonable to expect that spatial variations in the fractal dimension of topography will be related to spatial variations in landscape-forming processes. Several dynamical considerations may be used to explain the spatial organization of landscape processes. Fluvial landscapes comprise water and differentially erodible sediment in spatially and temporally variable quantities. The relative amounts of water and sediment lead to dif- ferent geomorphic processes, which cause changes in the configuration of the landscape. Geomorphic processes have been organized according to relative amounts of water and sediment by Takahashi (1981), Savage (1993), and Montgomery and Foufoula-Georgiou (1993) inter alia. The scheme of Montgomery and Foufoula-Georgiou is the most re- cent and most pertinent to the present study area, so it will form a prototype for the discussion presented here. Montgomery and Foufoula-Georgiou's scheme (illustrated in Figure 1.1) puts various fluvial and hillslope processes into categories controlled by upslope area (serv- ing as a proxy for both the amount and persistence of hydrological inputs), and slope. They suggest that it is instructive to consider hillslope and fluvial transport as end-members of a continuum of landscape processes ranging from turbulent fluvial washload transport, through bedload transport and gravity-driven processes in low- order streams, to shallow surficial hillslope failures and soil creep. Upon this range of processes it is necessary to impose other processes, such as deep-seated rockslides, which are contingent on the local configuration of bedrock features, and pervasive processes such as chemical denudation (by solution). 2 - 5 - 4 - 3 - 2 - 1 0 1 2 3 10 10 10 10 10 10 10 10 10 drainage area (km2) Figure 1.1: Continuum of Geomorphic Processes (after Montgomery and Foufoula- Georgiou, 1993). The thick solid line is typical slope-area plot for a steepland land- scape; the thick broken line is a typical plot for a lowland area. These modifications notwithstanding, one may expect to find testable evidence in the landscape for this interpretation of the fundamental controls on sediment transporting processes, and therefore landscape form. Several authors have outlined empirically-justified schemes for ordering landscape process and their associated mor- phological changes in this way (Ryder, 1981; Culling, 1988; Hewitt, 1989). The relation between altitudinally-zoned processes has often stood aside from altitudinally-zoned form (as observed by Outcalt and Melton, 1992; Evans and Mc- Clean, 1995); their re-amalgamation is central to the attempt within this thesis to relate fractal properties of topography to the processes that produce them. Among writers who do discuss the relation between the spatio-temporal structure of process and form at the landscape scale, few discuss the fractal properties of form. However Veneziano and Iacobellis (1999) consider fractal properties in some detail: they quickly point out that the scaling properties of hillslope and fluvial processes are not the same. They suggest that channel networks are fractal because channel incision varies as a 3 fractional power of upslope area, which takes values over a range of scales between zero and the size of the area of interest, leading to changes in the landscape with no characteristic scale. Furthermore, Veneziano and Iacobellis (1999) claim that fractal behaviour should not be expected on hillslopes because it is not inherent in the gov- erning equation for hillslope evolution, which is usually supposed to be gradient-driven linear diffusion. Culling (1960), Kirkby (1971), Tucker and Slingerland (1994), Martin and Church (1997) and Martin (1998) have shown that diffusion is a good representa- tion of hillslope processes over long time-scales, however it may not be adequate at reproducing the properties of landscape form at all spatial scales. In particular, diffu- sion does not correctly reproduce the fractal properties of real landscapes (as discussed explicitly by Turcotte, 1998). Diffusion can be used to represent processes that can be averaged over timescales of interest (Martin and Church, 1997); it therefore correctly reproduces the average changes of a landscape. However, observations of form (e.g., digital elevation data) are instantaneous samples that reveal the effects of events that occur below the timescales over which diffusion is valid. If these events (e.g., landslides) are not smooth, they might be expected to create landscape zones with fractal dimensions greater than two (i.e., non-planar hillslopes). Indeed, Hovius et al. (1997), Pelletier et al. (1997), Dens- more et al. (1997) and Turcotte (1999) have shown that individual mass-movement processes responsible for hillslope evolution are power-law distributed (over a limited range of scales), which suggests that fractal scaling ought to be an inherent feature of hillslope morphology. In order to present these arguments in greater detail, it is necessary to define the concept of a fractal landscape in a more rigorous way. 4 1.2 Fractal Models of Fluvial Landscapes 1.2.1 Fractal Properties of Topography The fractal dimension of a topographic surface is a measure of the degree to which elevation changes over different spatial scales. It is an indicator of landscape intricacy and an indirect measure of landscape roughness (Outcalt and Melton, 1992). A fractal is a shape without a characteristic length scale. It has a non-integer Hausdorff (fractal2) dimension, D, indicating that it is a space-filling curve (if 1 < D < 2), or a volume-filling surface (if 2 < D < 3). As geometric objects, fractals are characterized by the fact that they are invariant under scale transformations (enlarge- ment or reduction). Examples of fractional Brownian surfaces with fractal dimensions between two and three are shown in Figure 1.2. The history of the use of fractal concepts to characterize geomorphic features begins with L. F. Richardson3, who showed that coastlines vary over many length scales. Later, writing about branching river networks, Horton (1945) presented laws showing relations between statistics of streams of various sizes. None of these authors' works were interpreted in the context of fractal geometry until Mandelbrot (1967) pointed out the utility of this erstwhile theoretical branch of pure mathematics for describing intricate features of the natural world, which he elaborated considerably in later work4. His Fractal Geometry of Nature (1983) presents several geomorphological examples of fractal geometry including branching river networks, intricate island coasts and topographic transects. 2The Hausdorff and fractal dimensions are not always equivalent. Exceptions to the equivalence assumed here are uncommon, so will not be discussed (see Barnsley, 1988, for more rigorous definitions of each of these measures). described in Mandelbrot (1983). 4While Mandelbrot was the first to use the term fractal (fractus (Lat.) having been broken), various mathematicians (e.g., von Koch, Peano, Sierpinski) had studied the theoretical properties of curves that would now be called fractals (see Mandelbrot, 1983, for historical perspectives on this). 5 Figure 1.2: Simulated Fractal Landscapes: fractional Brownian surfaces (Klinkenberg and Goodchild, 1993) 6 Burrough (1981) was among the first geographers to apply fractal concepts to landscapes. His interim assessment was optimistic, but he urged practitioners to refrain from using the fractal concept when it was valid over only a small range of length scales and to seek physical reasons for the appearance of fractal forms. This concurs with the view of Leo Kadanoff (1986), a physicist whose initial reaction was cautiously optimistic. Both Burrough and Kadanoff recognized that fractal geometry is well suited to describing geographical and physical patterns, but can be poor at describing the processes by which these patterns are generated. A similar view was put by Church and Mark (1980), although fractal geometry was not the primary focus of their discussion. A contemporary reading of Church and Mark provokes the same view of fractals as of scale relations: they are a useful start but, through paucity of process, a poor conclusion. Throughout the 1980s, attempts were made to measure the fractal dimension of geomorphic features. Naturally, most success was found in the prototypes set out by Mandelbrot in 1983: • Branching river networks (Tarboton et al, 1991; Rodriguez-Iturbe and Rinaldo, 1997). • River meanders (Snow, 1989; Nikora, 1991; St0lum, 1996). • Topographic transects and fields (Klinkenberg and Goodchild, 1992). In general, mixed results were obtained: some were encouraging to supporters of the fractal geometry of nature, others showed results that discounted fractal hypotheses. If the 1980s were characterized by the measurement of fractal objects, the 1990s were dominated by attempts to link fractal form and physical or geomorphic processes. In the context of river basins, much of this work is summarized by Rodriguez-Iturbe 7 and Rinaldo (1997), two of the major participants in researching the fractal-process link. The current state of fractal research in geomorphology can be summarized in the following three views: • The fractal concept is useful in process-response models in so far as simple mono- fractal (or unifractal5) scaling results from the underlying physics of some land- scape processes; • Fractals are useful, as above, but landscapes are not simple mono-fractals; • The fractal dimension is a useful geomorphometric statistic that describes surface roughness, but it cannot be related to the physics of process-geomorphology. The first view, that landscapes are simple mono-fractals, was prevalent following Mandelbrot's thesis (1983). Turcotte (1992) demonstrated that topographic transects are examples of self-affine fractals (a self-affine fractal is one that scales in different ways in different directions). Further work has shown that in most cases landscape variability is not this straightforward. Lavallee et al. (1993) claim that multifractal concepts are useful in characterizing land-surfaces. A multifractal is a set on a geometric support; it may be formed by a multiplicative cascade, or by sampling a set with a spatially variable fractal dimension (Rodriguez-Iturbe and Rinaldo, 1997). Multifractal statistics are a natural extension of the fractal dimension, wherein higher order moments of the distribution under investigation are studied. A set is multifractal if the relation between moment order and moment variogram-slope is non-linear. Figure 1.3 shows a linear moment- slope relation, indicating mono-fractal scaling; for a multifractal, the moment-slope plot would be non-linear. 5Evans and McClean (1995) point out that "unifractal" is more consistent etymologically, how- ever "monofractal" is more common in the fractal literature: "mono-fractal" will be used here to accommodate both views. 8 1e+08 1e+07 h M1 (H = 0.858356) M M 2 (H = 1.69294) x separat ion (h) M 3 (H = 2.49445) • M 4 ( H = 3.25169) + 1000 M 5 (H = 3.96022) o Fractal R e g r e s s i o n — - 3.5 h Q. O <" 2.5 h E o E moment order Figure 1.3: Moment statistics of a land-surface (part of Daniels Creek in the Capilano Watershed; see Chapter Four). The lower of the two plots shows a linear moment-order relation, indicating that this particular landscape is a mono-fractal 9 In parallel to the positive reaction to Mandelbrot's (1983) Fractal Geometry of Nature came a more sceptical reaction (e.g., Kadanoff, 1986; Shenker, 1994). The central tenets of this objection are summarized neatly by Avnir et al. (1998) in a re- view of a sample of empirical fractal claims published in the physics journal Physical Review Letters. These objections are: (i) empirical fractals cannot continue indefi- nitely without violating other physical laws; and (ii) most studies of empirical fractals demonstrate scaling over a range of fewer than three orders of magnitude. Avnir et al. (1998) suggest that an understanding of these finite-scaling structures may be impaired by their treatment as fractals: they appear to be fractal, but are probably not. In landscape-scale geomorphology, Burrough (1981, 1993) has introduced the term pseudo-fractal to describe phenomena that appear fractal, and for which there is practical utility in measuring the fractal dimension, D, but the genetic processes of which are not fractal. In the present context this is appealing: erosive processes remove landscape material within a wide, but inherently limited, range of length scales. It seems impossible for a process that does have a characteristic length scale (or limited range of characteristic length scales) to create a land-surface with no characteristic length scale (or much wider range of spatial variability). However the cumulative effect of many scale-bound processes operating on the same surface may be to create a pseudo-fractal landscape. One goal of the present study is to demonstrate the effect on the fractal dimension of a surface of perturbations caused by large slope failures, and therefore to determine whether or not landscape processes with a wide but finite range of length scales can create the pseudo-fractal landscape that Burrough (1993) suggests. Veneziano and Iacobellis propose a compelling merger of fractal and pseudo- fractal theories of landscape that avoids the complexities and dissatisfactions of a 10 multifractal approach. They show that the combination of fluvial and hillslope pro- cesses can create complex fractal landscapes across which the fractal dimension is not stationary. They suggest that fractal analyses of topography should segregate between fluvial and hillslope process-domains. Their method will be used here in order to ex- amine the link between the nature and rate of hillslope mass-wasting processes and the fractal dimension of topography in the hillslope process-domains. 1.3 Object ives The first objective of this thesis is to test whether distinct fractal domains can be identified in the landscape. It is hypothesized that if landscape processes can be organized into process-domains, then their impact on landscape morphology ought to be to create areas of homogeneous fractal dimension. In order to test this hypothesis, several digital elevation models are analysed, after having been partitioned into the process domains suggested by Veneziano and Iacobellis (1999). The second objective of the present work is to examine whether process do- mains with similar fractal dimensions may be reproduced using a field-calibrated model of landscape evolution incorporating fluvial erosion, diffusive hillslope processes and episodic deep-seated bedrock landslides. These constitute some of the most important geomorphic processes in the areas chosen for empirical study, which are described in Chapter Two. 11 Chapter 2 Surface Processes and Tectonic History in the Queen Charlotte Islands and Southern Coast Mountains of British Columbia 2.1 Location Two areas are used to test the work reported here: parts of the Queen Charlotte Islands (Haida Gwaii), British Columbia, and the Greater Vancouver Regional Dis- trict (GVRD) Capilano Watershed, within the Southern Coast Mountains of British Columbia. Figure 2.1 gives the regional situation for these two sites, which form part of the Canadian Western Cordillera. Figures 2.2 and 2.3 show details of the location of the study areas. The purpose of this chapter is to provide tectonic, geological and climatic con- text in order that modelling and analysis decisions may be justified, and results in- terpreted. Beginning with the geological background and progressing to a review of surficial processes, the differences and similarities between the two areas are dis- cussed. The importance of tectonic structure and earth history in the development of the drainage network is also considered. 12 Figure 2.1: Locat ion of the two Study Areas (also showing principal tectonic features, described in this chapter) 13 Figure 2.2: Queen Charlotte Islands (study basins are described here and in Chapter Four) 14 Figure 2.3: G V R D Capilano Watershed (study basins are described here and in Chap- ter Four) 15 2.2 Tectonic background and Geological History The two areas lie within the same general tectonic region on the western edge of the Canadian Cordillera, close to the triple junction between the Pacific, Juan de Fuca/Explorer, and North America plates, which lies approximately 100 km west of Vancouver Island (Parrish, 1983). The tectonic regime is influenced by the relative rates of motion of these plates, given in Figure 2.1. South of the triple junction, subduction occurs offshore, where the Juan de Fuca/Explorer plate passes under the North America plate (Figure 2.1). North of the triple junction, is the Queen Charlotte Fault, which separates the Pacific and North America plates. This transform fault is the most seismically active in Canada, a fact that has profound consequences for slope stability on the Queen Charlotte Islands. The most recent period of mountain building in the Coast Range began 10 Ma BP, with a maximum average rate on the order 0.5 mm a - 1 Parrish (1983). More recent uplift rates of up to 2 mm a - 1 have been inferred from tidal records. In some heavily glaciated areas this uplift may be isostatic in origin, whereas in others the movements are mainly tectonic (Clague, 1989). In the following sections, these general rates will be supplemented with information about the nuances of each site. 2.2.1 Tectonics and Topography The long-term interactions between tectonics and topography have been the subject of much recent work (Molnar and England, 1990; Beaumont et ai, 2000). Uplift increases the amount of potential energy available for erosion by surface processes. This is important because, over long timescales, lowering of the land-surface (resulting from denudation by surface processes) is accompanied by a considerable removal of landscape mass. For normal crustal rock densities, this leads to further (isostatic) uplift amounting to approximately four-fifths of the original lowering. 16 Over short timescales, the relation between tectonic processes and landscape processes is also relevant to studies of the observed topography. The influence of lin- ear tectonic structures (lineaments, including faults and thrusts) on large landslides has been the subject of several studies in the Coast Mountains. Leir (1995) and En- glish (1996) have shown that the locations of structure-forming landslides are strongly influenced by the location of lineaments. In the present landscape, deep-seated land- slides occur with a lower rate than is likely either after glaciation (Ryder, 1981) or in the very early history of the Cordilleran orogen (cf. Hovius et a/., 1998). In the current study area, the influence of early slope failures on the present drainage network may have been obfuscated by the later occurrence of several widespread glaciations, which provide an immediate explanation of many important features of the contemporary landscape. Nevertheless, the pre-glacial topography in regions such as this may have been strongly influenced by widespread landsliding associated with orogenic activity. If this is true, landslides would have provided significant depressions in which the first cirque glaciers of subsequent glaciations may have formed. Early controls on the devel- opment of the drainage network are beyond the temporal scope of this study, however the influence of large contemporary slope failures on the drainage network is evident in the current landscape. 2.3 Queen Charlotte Islands The Queen Charlotte Islands (Haida Gwaii), lie approximately 100 km off the west coast of British Columbia. They are part of the insular tectonic belt of the Canadian Cordillera and one of the most seismically active regions in Canada, being close to the Queen Charlotte Fault. The islands have been uplifted by a maximum of 1 km in the last 10 Ma, according to the fission-track observations of Parrish (1983). 17 2.3.1 Bedrock Geology The Queen Charlotte Islands comprise three principal physiographic provinces: the Queen Charlotte Ranges, the Lowlands, and the Skidegate Plateau. Details of the bedrock geology are given in Sutherland-Brown (1968): the majority of the study area is underlain by early Tertiary lavas, which are gently bedded or flat. Sutherland- Brown (1968) divides the volcanics into several formations. The "soft" volcanics of the Masset and Yakoun formations are common: they have highly developed joints (especially in the columnar basalts), which render the rocks highly permeable and therefore susceptible to extensive sub-surface weathering. As a consequence, the area is prone both to deep-seated and superficial slides; indeed, the majority of recent rock slides have occurred in this formation (Alley and Thompson, 1978). Another important formation for slope failures is the "hard" volcanic Karmutsen group, found along the Skidegate Channel between Graham Island in the north and Moresby Island in the south. A number of recently intruded granitic plutons are of local significance, being responsible for the relatively high incidence of slope failure near Pivot Mountain and in Riley and Bonanza Creeks. Most notably, there is an area of sedimentaries in the north (Kunga and Longarm formations); it is not well-weathered, but is deeply gullied owing to a high degree of jointing. Weathering rates are higher for the Longarm sedimentaries, so some surficial slides have been observed (Alley and Thompson, 1978). 2.3.2 Glacial Legacy The Fraser Glaciation had an important effect on the landscapes of the area. Alley and Thompson (1978) note that at the peak of this most recent Cordilleran glacial episode there was an ice cap approximately 900 m thick over the area, which coalesced with the main Cordilleran ice sheet. The retreat of this ice cap, which occurred in line 18 with the retreat of the mainland ice sheet (between 14 and 12 ka BP) , left a veneer of t i l l masking the surface and has resulted in at least 18 m of isostatic rebound since 5 ka B P (Sutherland-Brown, 1968). During deglaciation, a period of low sea level (from 15 ka to 10 ka BP) was followed by a rise until 8 ka B P leaving some shorelines 15 m higher than at present Clague (1989). Clague (1989) suggests that ice loads on the Queen Charlotte Islands were insufficient to cause isostatic depression. In addition, the eastward migration of mantle material following the decay of the mainland icesheet provides a further explanation for the general rise in sea-level on the Queen Charlotte Islands during the Holocene. 2.3.3 Contemporary Processes The contemporary climate is one of cool summers and mild winters. Annual precip- itation at lowland settlements is generally 1,000-1,250 mm, but some western slopes receive up to 7,600 mm per year. This falls mainly as rain, with snow in winter at altitudes above 600 m. A further meteorological feature that is worthy of note is the prevalence of violent winds, gusting in excess of 190 km/h, which frequently lead to tree-throw. Seismic activity has direct consequences for slope stability. Alley and Thompson (1978) comment on 1,268 recorded earthquakes between the years 1899 and 1974, with 31 between Richter magnitudes 5 and 7; 7 between magnitudes 7 and 8; and 1 of magnitude greater than 8 (in August 1949). The combination of steep slopes, highly jointed rock, a failure-prone glacial veneer, high levels of intense winter storm precipitation, violent gusts of wind, and frequent seismic shocks results in a dynamic landscape in which slope failures are 19 frequent and sediment yields are high. Rates of sediment flux due to landsliding were measured by Rood (1984, 1990) and Gimbarzevsky (1988). A summary of the rates in each formation is given in Table 2.1. Formation Type Average Yield m3 ha-1 a - 1 Masset soft volcanics 1.5 Karmutsen hard volcanics 2.0 Haida sedimentaries 2.0 P T P granites 1.3 Table 2.1: Sediment Flux due to Landsliding on Queen Charlotte Islands (data from Rood 1984) 2.4 G V R D Capi lano Watershed The Capilano Watershed, in North Vancouver, is part of the Pacific Ranges of the Coast Mountains. The area is characterized by steep slopes (commonly greater than 35 degrees) and high relief (between 900 and 1,300 m over 2 or 3 km), which lead to active hillslope processes and high sediment yields. Resulting from the needs of watershed planners, there is a well-presented record of these processes available for use in the present study, including gully processes, shallow landsliding and (to a lesser degree) larger rock slides. The following sections are based on this inventory (referred to as G V R D , 1999), and on the earlier synthesis of (Ryder, 1981). 2.4.1 Bedrock Geology Bedrock in the area comprises intrusive igneous granodiorite, quartz diorite, diorite and lesser amounts of gabbro and migmatite, with great intact strength and widely- spaced joints. There are small areas of volcanic rocks, which are similar but more 20 highly jointed than the intrusives. Bedrock is exposed at mid and upper elevations; at lower elevations, there is a mantle of glacial drift and colluvial deposits. 2.4.2 Tectonic Uplift Uplift history in the Capilano watershed is similar to that of the Queen Charlotte Islands, and was the subject of the same regional study by Parrish (1983). However, in the last 10 Ma, the G V R D Watershed has seen between 2 and 3 km of uplift (cf. 1 km on the Queen Charlotte Islands). Ryder (1981) highlights the importance of a pause in uplift late in the Tertiary period, which allowed the development of an erosion surface that was further dissected upon subsequent uplift (but is currently visible in some upland parts). 2.4.3 Glacial Legacy The main Cordilleran ice sheet covered the G V R D Water Supply Area, which has been subject to a larger number of Pleistocene glacial repetitions than documented on the Queen Charlotte Islands (where the ice cover was separate from the mainland). Evidence of three major glaciations has been found in the Lower Mainland, although the most recent (the Fraser Glaciation) dominates the landscape. The Fraser Glaciation in the study area is described by Lian and Hickin (1993). It began roughly 30 ka B P with slow development of cirque glaciers in high, north- facing depressions, expanding until about 20 ka BP, when large valley glaciers had developed (Coquitlam Stade). Between 18 and 16 ka BP, the Port Moody Interstade saw a recession of the ice and reforestation of lowlands. Subsequently, the ice re- advanced until the glacial maximum at 14.5 ka BP, when a 1.5 km thick ice-tongue extended well into Washington State, USA. Ice retreated rapidly and it is likely that 21 it had melted from even the higher parts of the Capilano watershed before 10 ka BP. Whereas sea levels have risen since deglaciation on the Queen Charlotte Islands, on the mainland a period of isostatic rebound immediately followed the retreat of ice, leading to a marine regression in spite of global eustacy. In the first 2-4 ka since the retreat of ice, isostatic rebound of formerly glaciated areas was larger than the coeval eustatic rise, leading to a reduction of shoreline elevation by 100-200 m. However, in the last 8-10 ka shoreline elevation has risen again by about 20 m to its present level, probably as a result of global eustatic rise in sea level and forebulge collapse (Clague, 1989). Erosion rates in the watershed prior to the revegetation of slopes were consid- erably higher than at present, with a much higher incidence of both deep-seated rock slides and shallow slides in the veneer of glacial material (Ryder, 1981; G V R D , 1999). 2.4.4 Contemporary Processes The present climate is characterized by large frontal depressions that track succes- sively across the area between October and March, and summer high pressure that is occasionally interrupted by small intense storms that may produce extreme flooding. These conditions, combined with contemporary tectonic uplift, steep slopes and a residual veneer of glacial deposits result in high rates of hillslope failure, which are documented by G V R D (1999), Savigny (1991) and English (1996). Campbell (1999) recently determined a sediment yield due to landsliding of 168 m 3 k m - 2 a - 1 (= 1.7 m 3 h a - 1 a - 1 , similar to the values reported for the Queen Charlotte Islands) for the nearby Lynn Valley. 22 2.5 Summary The purpose of this chapter is to present the geological and geographical characteristics of the study areas, giving some published estimates of background rates and summary process statistics. In particular, it highlights the findings of Ryder (1981), Rood (1984, 1990), English (1996) and G V R D (1999), that contemporary processes and forms are strongly influenced by the tectonic structures in the area, and the glacial overprint. In Chapter Three, the rates and magnitude-frequency characteristics will be examined in more detail, so that their influence on landscape form may be assessed. 23 Chapter 3 Sediment F l u x due to Episod ic Events 3.1 Process rates Estimates of sediment transport by episodic processes presented in the previous chap- ter indicate that landslides and debris flows play an important role in the long-term denudation of the G V R D Water Supply area and the Queen Charlotte Islands. Quan- tification of sediment flux due to episodic processes is an important step in the con- struction of testable models of long-term landform evolution (Martin, 1998, 2000), and has important consequences for the resulting landscape form. In this chapter, the record of slope-failures in the two study areas will be examined in order to determine whether they follow a power-law magnitude-frequency distribution, and to discuss how they can be incorporated into a landscape evolution model. 3.2 Sources of Da ta The statistical analyses presented in this chapter are based on three previously pub- lished landslide datasets: (i) a set of 1,337 failures mapped by Rood (1984, 1990) for parts of the Queen Charlotte Islands; (ii) a set of 865 landslides collected as part of 24 the G V R D Ecological Inventory Program (GVRD, 1999); and (iii) a set of 49 major rock slides collected by Savigny (1991). 3.3 Surficial Slides: magnitude-frequency 3.3.1 Results The results presented here are from the G V R D and Queen Charlotte Islands study areas. Figure 3.1 shows a magnitude-frequency distribution for 1,337 failures on the Queen Charlotte Islands (the dashed line). A similar distribution for the G V R D land- slide inventory is also shown in Figure 3.1 (the dot-dashed line), although in this latter case there is little preserved evidence of a power law. In each case the slide volumes extend over a range of two to three orders of magnitude. For the Queen Charlotte Islands landslides, the power-law in Equation 3.1 (the solid line on Figure 3.1) was fitted to the most robust central part of the curve (between 400 and 7,000 m 3 ) . Frequency = 0.8 * Magn i tude - 0 6 3 (3.1) The exponent here is negative, indicating that frequency of events declines with their magnitude; it is also less than one, which shows that large events are more important than the cumulative effect of small events (Hovius et al., 1997). This is of particular interest because it implies that small events may be modelled with diffusion provided that large events are considered explicitly. There are two interesting features at each end of the distribution: (i) a relative lack of small events (< 400 m 3), probably a result of the limited sampling resolution and rapid revegetation of landslide remains, but possibly a dynamical feature of the landsliding process; and (ii) a break of slope for large events (> 7,000 m 3), which may indicate a change in controlling conditions or another sampling limitation. 25 0.1 0.01 k 0.001 0.0001 1 ' ' ' 1 ' 1 1 — ' — ' *. QCI (Frequency) Power Law (QCI) • . . . - O £H GVRD (Frequency) : E / / / i / / i / / \ \ \ \ / / ,° x \ ^ \ i . . . i . 1 . . . . . . . . 100 1000 10000 Slide Volume (m3) 100000 Figure 3.1: Queen Charlotte Islands: Landslide Magnitude Frequency 3.3.2 Power-law Scaling and Self-organized Criticality Several hypotheses have been offered to account for power-law scaling in natural phe- nomena. These include the dynamical state of self-organized criticality, and the pos- sibility that a power-law may arise from imperfect sampling of rare events. These competing explanations are elaborated below. Power-law scaling may be a natural feature of hillslope instability (Hovius et al, 1997; Pelletier et al, 1997; Turcotte, 1999). Failures resulting from the exceedence of critical thresholds commonly have power-law distributed magnitude distributions; self- organized criticality has been offered as an explanation for this by Bak et al (1987, 1988). 26 Self-organized criticality (SOC) is a feature of systems maintained far from equi- librium that organize themselves through internal feedbacks to a state in which a sys- tem parameter fluctuates with power-law distributed adjustments around a marginally- stable critical value (Bak et al, 1987; Turcotte, 1992). In the present context, slope- failures may be the power-law distributed adjustments to which Bak et al. refer. This analogy follows naturally from their definitive example of a growing pile of dry sand. In this demonstration, successive grains of sand are added to a circular table. Initially the sand grains form a pile that is stable because its angle is below the angle of repose for sand. As grains continually fall the material in the pile reaches the angle for fail- ure. The claim of Bak et al. is that further grains will not simply roll down the stable surface of the pile. Instead, owing to the configuration of the grains, the pile will grow beyond the critical angle by some amount until an avalanche (sand-slide) occurs. Bak et al. present a variety of experimental evidence to show that the magnitude-frequency distribution of the resulting avalanches follows a power-law. The direct analogy of the sandpile to a layer of surficial hillslope material is initially appealing, although Mehta (1992) shows that there are significant differences. Recently, a more physically-based justification of the power-law distribution of slide magnitudes has been presented by Hergarten and Neugebauer (1999), who show that power-law distributed landslides result from a simple model of threshold instability on an inclined plane with one Dirichlet boundary (i.e., constant fluvial downcutting along one side). These theoretical conclusions concur with the empirical findings and explanations of Hovius et al. (1997) and Pelletier et al. (1997). An alternative explanation of the proposed power-law is obtained through a consideration of the natural timescales over which contingent hillslope stability criteria apply. In most cases, slope failure is a result of interactions between slope loading by soils and vegetation, strengthening through the development of a root network and 27 strength-reducing factors (especially pore water pressure in soils). Since these each have different natural timescales (pore water pressure typically changes more rapidly than vegetation changes can occur), it is likely that observations over a discrete time- interval will show a large number of small events and a small number of large events; possibly resulting in a finite-scaling power-law magnitude frequency distribution of slide volumes. Nevertheless, there are some differences between a theoretical power- law and the distributions obtained from field and air photo data, which relate to sampling biases discussed below. 3.3.3 Power-law Scaling and Sampling Bias The nature of rare and episodic events is that the largest events are usually the most rare, and the smallest events are the most common. This introduces a natural bias to any data collected on such events: all observations must be collected over a limited time period and a restricted area, therefore a large number of small events and a small number of large events will be recorded. Further temporal or spatial sampling biases may result if the temporal range of observations is too short, in which case a survey may underestimate the number of large events and therefore overestimate the importance of small events. Similarly, if the sample area is too small, the survey may miss large, rare slides (Figure 3.2). 3.3.4 Finite Scaling Effects The kink in the Queen Charlottes slide distribution, at 400 m 3 , may be another sam- pling effect. This represents a relative paucity of events with small volumes (roughly a cube of 5-7 m on each side). One explanation of this feature is that a proportion of these events is missed in the air photo analysis either because the events are too 28 Size Size (A) High Resolution, (B) Low Resolution, Small Extent Large Extent Figure 3.2: Extreme Event Distribution small to be observed, or because small events are more rapidly and easily obscured by revegetation. Smith et al. (1986) have shown that the persistence of slide evidence is a log- arithmic function of time, suggesting that the probability of observing a slide can be written as a function of slide size: Pr(obs) = exp[(logP(re/) * cutoff)/size] (3.2) where cutoff is the size (in m 3) of a reference slide visible with probability P(ref). An example is plotted in Fig. 3.3, where cutoff = 2,500 and P(ref) = 0.5. If this resolution function is convolved with an artificially generated set of slide magnitudes chosen to obey the expected power-law the observed landslide distribution also shows a kink at low magnitudes similar to that found in the real data set (compare Fig. 3.4 with Fig.3.1). Departures from power-law distribution may be due to sampling biases. How- ever the observed kink at low magnitudes may be a genuine feature of the mechanics of slope failure: it may result from the need to accumulate a critical amount of sedi- ment in order that its weight may overcome frictional resistance, consequently fewer 29 ^ ' ' CUTOFF = 2500 j . . . . . . . - J . . •—I 1 10 100 1000 10000 100000 Volume (m3) Figure 3.3: Resolution Function 1000 10000 100000 Volume (m3) gure 3.4: Landslide Magnitude Frequency modified by Resolution Function o 1e-50 1e-100 1e-150 1e-200 30 slides occur with volume lower than a critical amount which may be represented by the observed kink. To decide between these two hypotheses is beyond the scope of the present work, however it is anticipated that a critical test between the two could be performed in an experimental study of failure of a naturally occurring surficial material: if fail- ures occurred on all scales down to the dimensions of individual grains then the kink observed in the distribution of real landslides is likely to be an artifact; if there is a measurable lower limit to the size of failure that can occur in the experimental case, then confidence increases that finite scaling is an inherent feature of slope failure in natural materials. 3.3.5 G V R D Slides The magnitude-frequency distribution of the G V R D dataset (the dot-dashed line in Figure 3.1) has a strikingly different shape to the Queen Charlotte Islands dataset. The curve has no linear part, so these data do not follow a power-law distribution. Instead they appear to have a characteristic scale. Two explanations for this result are possible. The first explanation is that a different resolution function would be capable of creating such an extreme lower-end kink to reproduce the curve shown from an original power-law. The combination of an enhanced lower-end kink and the relative lack of large events (exacerbated by the finite sampling space scale and time interval) leads to the hump observed for the G V R D data. This hypothesis implies that these landslides occur with a power law distribution over a wide range of scales, but it is impracticable to make observations with sufficient resolution, over a wide enough or long enough time, to see the underlying power-law. 31 The second explanation is that the underlying data are in fact exponentially distributed. This implies that there is a characteristic scale at which landslides occur. Distinction between the two theories is important in the study of the dynamics of mass- wasting. It is suggested that the most fruitful approach would combine experiment (as above) and theory, to determine a suitable distribution for episodic events, which is consistent with their known mechanisms. 3 . 3 . 6 Implications The practical implications of surveys of the magnitude-frequency distribution of slope failure are readily apparent in the quantification of risk to humans, their possessions and the structures placed close to landslide hazards. It is worth note that the approach taken here provides probabilistic predictions of landslide statistics, rather than deter- ministic predictions of the spatial and temporal incidence of the events in question. It is therefore of value in the evaluation of generalized risk, and has applications in actu- arial risk analysis, and in the optimal design of risk-mitigating engineering structures (for a recent example in the same area, see Hungr et al., 1999). The demonstration of a power-law landslide magnitude-frequency distribution also has several theoretical implications: (i) the statistics can be transformed for incor- poration into long-term landscape evolution models (see Chapter 5), and (ii) power-law perturbations to hillslope areas may lead to an increase in the fractal dimension of their topography. The second of these points forms the basis of the subsequent chap- ters of this thesis. In what remains of the current chapter, the landslide statistics are transformed into a non-linear slope stability relation and the influence of very large bedrock landslides is assessed. 32 3.4 Slope and Volumetr ic Transport Ra te 3.4.1 Approach The relation between slope and sediment flux is important because it reveals slope- stability thresholds governed by the dynamics of slope failure. In addition, a constitu- tive relation (Eq. 3.3) between slope (dh/dxi) and sediment flux (q) forms an integral part of the approach to long-term landform evolution taken in Chapter 5. 3.4.2 Parameterization In order to predict sediment flux from slope, a suitable function, / , must be found, and parameters must be fit using the observed landslide data. Martin (1998, 2000) has suggested that an hyperbolic tangent function is appropriate: Martin (1998) parameterized Equation 3.4 for the same Queen Charlotte Islands land- slides described here, with Kx = 0.43, hx0 = 1.06, hscaie = 0.19. The procedure adopted by Martin (2000) for the development of a slope-transport rate for shallow landslides and debris flows on the Queen Charlotte Islands was re- peated for the G V R D landslide inventory. Landslide width, length, thickness and initiation slope were reported in the G V R D Ecological Inventory (1999), with measure- ments made from air photos and calibrated using field data. Volumes were computed using a relation presented in the G V R D analysis, and the slide population was divided into ten slope-classes as shown in Table 3.1. The area falling into each slope class was q = f(dh/dxi) (3.3) (3.4) 33 estimated from a digital elevation model of the area (also supplied by G V R D ) . The temporal range of the observed data was estimated as 30 years, using dendrochronol- ogy of woody debris found within a sample of the deposits ( G V R D , 1999). Slope Class Yield m/m 3 —2 —1 0.00-0.25 0.004 0.25-0.50 0.011 0.50-0.75 0.035 0.75-1.00 0.068 1.00-1.25 0.113 1.25-1.50 0.112 1.50-1.75 0.060 1.75-2.00 0.028 2.00-2.25 0.011 2.25-2.50 0.000 Table 3.1: Slope-Volume Classes for G V R D Landslides A function of the form in Equation 3.4 was fit by eye to the data in Table 3.1 to obtain initial "seed" values for a non-linear least-squares curve fitting routine nlm in the statistical package R, which resulted in the relation shown in Figure 3.5, with r2 = 0.86 (The fitted parameters are shown in Table 3.2). This relation is valid for slopes with gradients between zero and approximately 1.5 (i.e., 60 degrees). On steeper slopes, surficial failures of the type analysed here give way to bedrock-controlled failure mechanisms, such as deep seated rock sliding. Parameter Queen Charlotte Islands G V R D Martin (1998) this study 0.43 0.12 1.06 0.81 hscale 0.19 0.32 Table 3.2: Parameters for Non-Linear Slope-Volume Relation 34 35 3.4.3 Discussion The results reported in the previous section provide strong support for the general applicability of the non-linear transport function proposed by Martin (1998, 2000). They form the basis of the constitutive relation between slope and sediment flux used in the landscape model of Chapter 5. Other non-linear functions have subsequently been proposed (Tucker and Bras, 1998; Roering et al, 1999), however these are a step function and exponential function respectively. The advantage of the hyperbolic tan- gent, proposed by Martin (1998, 2000) and subsequently used here, is that it correctly identifies the maximum slope on which (shallow) landsliding is important. Beyond this value of slope, stability is governed by factors other than the internal friction of surficial material. 3.5 Rock Slides Large rock slides, which remove material on scales comparable to local relief, are important features in the evolution of landscapes in this area (Savigny, 1991; English, 1996; Hungr et al, 1999). The morphological effect of a rock slide may persist if it is large compared to the rate of its subsequent obliteration (by weathering and erosion). If rock slides occur over a wide range of scales their effect may be to increase the fractal dimension of hillslopes, because they will impart to a hillslope topographic variation over a range of scales. It is appropriate to consider large rock-slope failures as an extraneous shot noise in the present work, since none of their controlling factors besides slope is codified within the landscape evolution model. In this section, the statistics of a regionally-significant rock slide inventory in the region are examined, in order to provide a basis upon which to compare model-generated landslide statistics with those observed in nature. 36 3.5.1 Rock Slide Inventory One of the largest inventories of rock slides in the region was made by Savigny (1991), who catalogued 49 rock slides in the Fraser Canyon and vicinity. This is the best record of very large slope failures close to the areas in which the present work is based. Savigny (1991) produced a map of slide areas based on the interpretation of 1:50,000 black and white aerial photography, which was digitized to find the slides areas. The magnitude-frequency statistics of this dataset are shown in Figure 3.6. The record of large landslides extends over more than two orders of magnitude (0.2 - 35 Mm 3 ) . Their cumulative distribution takes the general form of a power-law, which is similar to that found by Hungr et al. (1999) for a dataset from a nearby area. However, the power-law does not appear to be a good fit. There are several distinct breaks of slope at 1, 4 and 20 M m 3 , which are also present in the ordinary magnitude-frequency distribution. 3.5.2 Discussion Large slope failures occur over a range of scales, however within this range, there ap- pear to be characteristic scales, showing that slope failures tend to have characteristic volume ranges. The multimodal magnitude-frequency plot may owe its explanation to the particular range of geological configurations in the area, or to different rock-slope stability controls (e.g., seismicity, hydrometeorological factors). The implications for landscape form and the fractal dimension cannot be determined from process observa- tions alone. The landslides may increase the fractal dimension by imposing a "noise" on the diffusive processes; however, this noise may reduce the fractal dimension be- cause it occurs within a finite range of scales (2.25 orders of magnitude). 37 Cumulative Frequency Frequency - - *- •-JK-l- 100000 1e+06 Magnitude (m 1e+07 1e+08 Figure 3.6: Rock Slide Magnitude Frequency (data from Savigny 1991) 38 3.6 Summary Magnitude-frequency statistics of geomorphic processes indicate the relative impor- tance of different sized events. In this chapter, it has been shown that on the Queen Charlotte Islands larger events are more important than the cumulative effect of small events. This has important implications for the way that these processes are incor- porated into models of long-term landscape evolution: it means that small events can be treated diffusively, leaving only large mass-wasting processes to be represented explicitly. A second use to which magnitude-frequency statistics can be put is to deter- mine whether particular processes obey a power-law. These observations may be used to justify conjectures about the system-wide dynamics of the processes in question, and their effects on landscape morphology. Recent claims of power-law distributed landslides (Hovius et ai, 1997; Pelletier et ai, 1997) are supported by the findings of the present study from the Queen Charlotte Islands, although the power-law found here is valid over just one order of magnitude. It is premature to conclude that the slope-failures studied here are scale-invariant, although they do occur over three orders of magnitude. This does not rule out the possibility that landsliding is a power-law process that may result from the self-organization of mountain slopes to a marginally stable state (Turcotte, 1999), however it does require the restriction that this be con- sidered a valid representation only within a limited range of topographic scales (cf. Burrough, 1993; Avnir et al, 1998). The influence of landsliding on the fractal prop- erties of hillslope topography will be the subject of the next chapter. 39 Chapter 4 Fracta l propert ies of topography 4.1 Geomorph ic Processes and Fracta l Landscapes Fluvial processes are frequently thought to increase the fractal dimension of topogra- phy, whereas it is often supposed that hillslope processes reduce the fractal dimension of landscapes (see Rodriguez-Iturbe and Rinaldo, 1997, for examples of this view). In the previous chapter, hillslope processes in the study-area were shown to act over a wide range of length scales. In this chapter, digital elevation data are used to examine the relative roles of fluvial and hillslope processes in creating fractal landscapes. The fractal dimension of a land surface is a measure of the degree to which elevation varies over different spatial scales. It is an indicator of topographic auto- correlation, and an indirect measure of topographic roughness. Fluvial features are expected to scale with upslope area because, in general, a smaller upslope area leads to a smaller river. At some point, this argument must cease to hold: i.e., for some value of upslope area the landscape no longer contains a channel. In areas of the landscape that are beyond this threshold for channel initiation, any fractal property of topography must be due to hillslope processes. 40 4.2 Fractal Properties of Landscape Three types of fractal behaviour are widely acknowledged (see Lam and de Cola, 1993, for geographical examples): mono-fractal scaling: self-similar scaling with a single fractal dimension, or self- affine scaling after an anisotropic geometric transformation (Mandelbrot, 1983); multifractal scaling: characterized by an infinite number of fractal dimensions on a geometric support (first described by Frisch and Parisi, 1985, cited by Rodriguez- Iturbe and Rinaldo, 1997); pseudofractal scaling: having measurable fractal characteristics but without scale- invariant cause, or scaling over only a limited range of scales (often caused by sampling limitations; Burrough, 1993). Since the introduction of the term "fractal", each of these concepts has been applied to the study of landscapes. River networks are examples of self-similar mono- fractals (Mandelbrot, 1983; Tarboton et al., 1991); some topographic transects are examples of self-affine mono-fractals, although Veneziano and Iacobellis (1999) have suggested that scaling in the transect is not this simple. Lavallee et al. (1993) contend that landscapes (along with "most geographical and geophysical fields", p. 158) are multifractal because they exhibit the mixed effects of several different processes. Re- sponses to the claim that topography is multifractal can be divided into three groups: • topography is a true multifractal because it is produced by a multifractal process, such as a multiplicative random cascade (Lavallee et al, 1993; Rodriguez-Iturbe and Rinaldo, 1997). 41 • topography appears to be a multifractal because a typical landscape is produced by several processes, each imparting certain scaling characteristics to the sur- face. Thus non-stationarity of the fractal dimension leads to a spurious (and misleading) claim of multifractality (see Veneziano and Iacobellis, 1999). • fractal and multifractal concepts are not useful descriptions of landscapes, be- cause scaling occurs over only a limited range of scales and is not always related to landscape-forming processes (Burrough, 1993; Evans and McClean, 1995; Avnir et al, 1998). Veneziano and Iacobellis (1999) suggest that some parts of the landscape close to river channels do have fractal characteristics, while hillslopes do not. They show that if these two zones are mixed in fractal analysis then multifractality is observed; whereas if the two zones are analysed separately, simple mono-fractal scaling is found. They conclude that an individual watershed may be partitioned into spatial zones of homogeneous fractal dimension, based on area upslope of each point. Once these definitions have been made, elevation in the domains so created scales as a mono- fractal (i.e., with a single fractal dimension), because a single geomorphic process characteristically dominates the domain. In the remainder of this chapter, the methods with which this result may be established are outlined, and the results of such analyses are shown for several drainage basins within the study areas described in Chapter Two. 42 4.3 Measur ing Fracta l Propert ies of Landscape Several methods are available for estimating the fractal dimension of a topographic sur- face (Klinkenberg and Goodchild, 1992; Rodriguez-Iturbe and Rinaldo, 1997). These are: • dividers method for a topographic cross-section. • box-counting method for the boundary of a contour polygon. • structure function analysis (variogram method). Klinkenberg and Goodchild (1992) have performed a comparison of these meth- ods for the estimation of fractal dimensions from digital elevation data. They conclude in favour of the variogram method because it is most stable. The variogram method will be used in the present work owing to its stability, the ease with which irregularly- spaced elevation sample points may be analysed, and the straightforward modification required to examine multifractal characteristics of the surface. The following descrip- tion of this approach is based on Veneziano and Iacobellis (1999). The fractal dimension of a sampled field of elevations may be obtained by supposing that the distribution of elevations between points at a given separation- distance is related to the separation-distance according to Equation 4.1. Ah(L) = r~H Ah(rL) (4.1) where Ah(L) is the increment in elevation between two points at distance L, = denotes equality in distribution, r is some positive number, and H is a constant called the index of self-similarity. If this is the case, then the 5th general moment of the distribution will behave according to Equation 4.2, for all moments that exist. 43 E[\Ah(L)\q] = r~qHE[\Ah(rL)\q] (4.2) The fractal dimension, D = 2 — H, which lies between 2 and 3 for a topographic surface, may be obtained from a log-log plot of the qth moment of Ah(L) against separation distance L. The intercept of this plot is also useful in the present context: it is the lacunarity of the surface (70). It is independent of D and expresses the patchiness of the surface variability. A surface with a high lacunarity has many patches of different surface variability and so appears "moth-eaten". In addition to obtaining D and 70 from the variogram plot, it is possible to examine the scaling behaviour of higher order moments of the dataset. A plot of H(q) against q will be linear if the surface is a mono-fractal, or non-linear if it is multifractal. Data from the Queen Charlotte Islands and Greater Vancouver Regional Dis- trict (GVRD) Capilano watershed were analysed using these tools, after first parti- tioning the surface into possible process-domains. 4.4 Process-domains and the Fracta l D imens ion Several authors have proposed the existence of landscape areas dominated by differ- ent processes (e.g., landsliding, soil creep, solution, fluvial erosion). Ryder (1981) and Hewitt (1989) propose classifications based on altitudinal zonation of landscape processes. These observations are interesting in the present context because they are supported by claims that the fractal dimension changes with altitude (Burrough, 1993; Evans and McClean, 1995). Altitudinal zonation of landscape processes is intuitively appropriate because processes at higher altitudes have available lower concentrations of water than at lower altitudes. Reduced temperatures and different biotic processes result in a slower rate of weathering, which leads to a reduced amount of easily mo- 44 bilized material. Consequently, sediment transport is less frequent, or more episodic, and takes the form of surficial landslides and larger rockslides. This spectrum of processes is summarized in the state-space diagrams of Taka- hashi (1981), Savage (1993), and Montgomery and Foufoula-Georgiou (1993), each of whom demonstrates that the relative concentration of water in a landscape unit is re- sponsible for defining the dominant processes that occur within the unit (see Chapter One). While altitude and concentration of water are correlated, a more appropriate measure of (time-averaged) water concentration that may easily be approximated from digital elevation data is the area upslope of a given point. Following Veneziano and Iacobellis (1999), landscapes may be divided into zones that are influenced by the channel network, and those that are unchanneled hillslopes. The boundary between the two is defined by examining plots of average slope against upslope area (see Fig 1.1). Mark and Aronson (1984), Tarboton et al. (1991) and Montgomery and Foufoula-Georgiou (1993) have each observed a break of slope at upslope areas of the order ~10 2 - 103 m 2 . It is claimed that this corresponds to the hillslope-fluvial transition; thus it is similar to the constant of channel maintenance defined by Schumm (1956). In the present work, slope was determined from a 10 m D E M of the study area. Upslope area was found using A R C / I N F O ' s flow accumulation algorithm, from the same elevation data. Each slope-area plot shown in this chapter demonstrates an inflection at its lower end. This was used to identify the constant of channel maintenance, allowing the hillslope and channel process domains to be separated1. Channel network fractal dimensions were measured from the slope-area plots (Veneziano and Iacobellis, 1999), using the Octave/Matlab script mfslarea.m given 1For the GVRD locations, slope-area plots created using this procedure were found to be similar to plots from field data (collected by White, pers comm., 2000); the two methods also produced similar estimates of the constant of channel maintenance. 45 in Appendix A . The fractal and multifractal characteristics of the hillslope areas were then studied using mf s, a C++ program also given in Appendix A . Preliminary inves- tigation indicated that the area above the channel network was slightly multifractal when treated as a homogeneous unit. To correct for this, the hillslope region was fur- ther subdivided into an upper region (dominated by failures in bedrock) and a lower region (dominated by shallow failures in surficial material). A n example of the typical zonation pattern is shown in Figure 4.1. Results are presented below for twelve sub-basins of the Capilano River and eight drainage basins on the Queen Charlotte Islands (see Chapter Two). 4.5 Resul ts In total, twenty drainage basins were investigated. Twelve sub-basins of the Capilano River were chosen; they are underlain primarily by diorite and granodiorite with a surficial cover of glacial drift and colluvial deposits (see Chapter Two). Table 4.1 shows the area of each of these; small sub-basins were chosen in an attempt to analyse areas with consistent controlling conditions, and to preserve computational tractability in the analysis of this high-resolution data-set. On the Queen Charlotte Islands, eight watersheds were chosen in the central portion of the islands; their geologies and other data are given in Table 4.3. Of these, only in Riley, Mountain and Armentieres creeks is there a significant amount of logging (10-20% by area, according to Rood, 1990). Maps showing the location of each watershed are given in Figures 2.3 and 2.2. 4.5.1 G V R D Capilano Watershed Slope-area plots for a 10 m D E M of the Capilano Watershed (supplied by G V R D Watershed Planning Department) are shown in Figures 4.2, 4.3 and 4.4. The most 46 common hillslope-channel threshold is 2,000 m 2 . This indicates that that only 2,000 m 2 (=0.002 km 2) of contributing area is required for channel initiation in the G V R D watersheds. Fractal dimensions of the stream networks in each of these watersheds are given in Table 4.1. The average of these is 1.33, suggesting that channel networks in this area are not as space-filling as usual (typical channel network fractal dimension is between 1.6 and 2.0; Veneziano and Iacobellis, 1999). A l l the DEMs studied had mono-fractal channel networks, as shown by a linear moment-slope relation (see Figures 4.2, 4.3, 4.4). The low channel fractal dimension combined with small upslope area required for channel initiation is consistent with the area's Quaternary history, the legacy of which is a mantle of highly erodible glacial deposits on which a steep, linear channel network has recently re-formed. Watershed Area Slope-Channel Transition Fractal Dimension (Network) km? m? Daniels (1) 4.5 2000 1.39 Daniels (2) 1.1 sample too small Daniels (3) 3.4 2000 1.52 Daniels (4) 1.0 2000 1.30 Sisters (1) 4.2 2000 1.26 Sisters (2) 6.5 2000 1.33 Sisters (3) 4.1 2000 1.23 Sisters (4) 2.1 2000 1.30 Enchantment (1) 4.3 2000 1.43 Enchantment (2) 2.0 sample too small Enchantment (3) 1.8 2000 1.40 Enchantment (4) 2.3 2000 1.18 Table 4.1: G V R D Capilano Channel Network Fractal Dimensions 47  uj/iu ado|s |EOO| LU/UI edo|s |eoo| Ul/W 3dO|S |EOO| o o o Figure 4.2: Slope-Area Plots: G V R D Daniels Creek 49   Points on the channel network (i.e., with upslope area greater than 2,000 m 2 (20 pixels)) were removed from the elevation data in all subsequent analyses in order to avoid the problems of non-stationarity documented by Veneziano and Iacobellis (1999). To separate the upper bedrock slopes from the lower soil mantled slopes, the remaining elevation points were divided into two sets: the set U, with upslope area < 600 m 2 (< 6 pixels), and the set L, with 600 m 2 - 2,000 m 2 (6-20 pixels). Fractal dimensions were measured using mf s (detailed in Appendix A) which measures fractal dimension using the variogram or structure function method. It is capable of estimating D with precision ±0.05 for a sample of 2,000 data points. It also plots moment-slopes, which are used to determine whether or not the elevation field is multifractal (see Section 4.3). Where available, 2,000 elevation points were used; in all cases the fractal dimension reported is the average of three independent measurements of a randomly drawn subsample of the dataset in question. Variograms are plotted for separations between 10 m and 100 m (the 10 m lower limit is given by the resolution of the D E M ; the upper limit is imposed to ensure that sample points remain close to the same hillslope, following the advice of Veneziano and Iacobellis 1999). The resulting variograms are shown in Figures 4.5 through 4.10 and the results are summarized in Table 4.2. The results show that in all cases D is significantly greater than 2.0, indicating that there is scaling (over a limited range of one order of magnitude) for the hillslope regions. In all cases, D for upper hillslopes is greater than or equal to that for lower hillslopes (although it remains unclear whether the difference is within the resolution of the algorithm for estimating D). In all cases but one, the scaling was mono-fractal (as shown by the linear increase of higher order moment-slopes in Figures 4.5 through 4.10. The one case where a slight multifractality was noticed is for Sisters Creek (3). This creek has a particularly low channel network fractal dimension (cf. drainage density), 52 1U8LU0LU )U91U0UI Figure 4.5: Hillslope Variograms: G V R D Daniels Creek ( D l - D2) 53 I" " • I' • I I' • • I" 1 X It C L C L Q "_ 5 & CM to r~ w co cy CT ll ll co^r 9 s 11 tr X — "—C\J ^ 2 o « s i ™ co lire S.-S 8" Z.E 5^ IS era r-~ co ,_: Figure 4.6: Hillslope Variograms: G V R D Daniels Creek (D3 - D4) 54 h x * o H- Q x k CL 8 8 o * - " 1 Sis II rr co OJ CO as O CM II £ 2,2 Si s i CO —' CM CO SI Si CT> © © O g O i - Figure 4.7: Hillslope Variograms: G V R D Sisters Creek (Si - S2) 55 luaiuow luauuow Figure 4.8: Hillslope Variograms: G V R D Sisters Creek (S3 - S4) 56 . ii t r C\J i n cn u i RS CO to " — O J ^ 2 8 4 CNJ £ r - cn co $ csi <ri II II S CO o Ol CO '—CM to if M • x k C L C L CM LU • - •- - - •- - g 2 X X * -I ° CD 1 CM LU s s s 8 o o g I o Figure 4.9: Hillslope Variograms: G V R D Enchantment Creek ( E l - E2) 57 C K M + CJ X k -I ° _o co LU l C O w CM CO co r- O CO CM <*? CO — ' 5 5 CM e II rr § 8 3 S O CM S" <° ̂ <5 C M CO II II C O S CO w co oir-. co • 8 £ Figure 4.10: Hillslope Variograms: G V R D Enchantment Creek (E2 - E3) 58 Watershed Lower Slopes Upper Slopes D To D 7o Daniels (1) 2.10 -0.26 2.15 -0.27 Daniels (2) 2.10 -0.27 2.11 -0.28 Daniels (3) 2.10 -0.17 2.12 -0.20 Daniels (4) 2.10 -0.27 2.14 -0.27 Sisters (1) 2.10 -0.19 2.13 -0.15 Sisters (2) 2.07 -0.27 2.11 -0.22 Sisters (3) 2.09 -0.22 2.14 -0.10 Sisters (4) 2.10 -0.17 2.12 -0.16 Enchantment (1) 2.12 -0.16 2.17 -0.14 Enchantment (2) 2.08 -0.12 2.13 -0.07 Enchantment (3) 2.11 -0.14 2.15 -0.06 Enchantment (4) 2.10 -0.14 2.13 -0.11 Average 2.10 -0.20 2.13 -0.17 Table 4.2: G V R D Capilano Hillslope Fractal Statistics indicating that hillslope and channel processes may be less distinct than in areas with a more extensively developed channel network. The appearance of multifractality may therefore have resulted from inhomogeneity of the processes occurring within this sub-basin. The Capilano sub-basins have non-planar hillslopes. This is contrary to the common assumption that hillslopes are dominated by diffusive, homogenizing pro- cesses that destroy fractal properties introduced by fluvial processes. Furthermore, in some cases there is a distinction between the "rougher" upland hillslopes that may be dominated by rarer events and the lower hillslopes, which are dominated by more frequent shallow failures leading to a lower fractal dimension. The results presented here suggest that there is some characteristic of hillslope mass-wasting that modifies D in areas not influenced by stream channels. This may be a sediment transport process that operates over a wide range of length scales (such as the hillslope processes described in Chapter Three); alternatively, it may be 59 a scale-bound process that leads to a pseudo-fractal effect captured in this analysis (see Burrough, 1993). Since this question is relevant to both the G V R D Capilano watersheds and the Queen Charlotte Islands, it will be discussed after the presentation of results from the latter. 4.5.2 Queen Charlotte Islands Slope-area plots for eight watersheds on the Queen Charlotte Islands are shown in Figures 4.11 and 4.12. (These were created from a 30 m resolution grid interpolated from B C T R I M data using the A R C / I N F O T0P0GRID routine.) The most common hillslope-channel threshold is 6,000 m 2 (= 0.006 km 2). This value is three times higher than the hillslope-fluvial transition in the Capilano sub-basins, suggesting that there is a difference in the flow-properties and erodibility of the substrate2. The fractal dimensions of the channel networks in the Queen Charlotte Islands watersheds are given in Table 4.3. These are considerably higher than those for the G V R D , indicating that the channel networks are better developed on the Queen Char- lotte Islands, since they fill the plane more completely. They are all mono-fractals, as shown by the linear relation between moment-slopes for higher order moments. The combination of the greater upslope area required for channelized flow and higher fractal dimension of streams on the Queen Charlotte Islands indicates that these land- scapes are more permeable and more highly erodible than those of the G V R D Capilano watershed, which is consistent with the geology and glacial history of the two areas (the Queen Charlotte Islands have had longer ice-free periods during which channel network development has occurred; see Chapter Two). 2 The hillslope-channel threshold increased three-fold when the DEM resolution decreased by the same proportion. In the present case, this is considered to be genuine, rather than an artifact of the flow accumulation routine, since that algorithm does not use information on cell sizes. 60 i i v E U l / U I 8dO|S |B00 | ui/ui edo|s |EOO| UI/UJ SdO|S | E 0 0 | U I / U J ado | s |Bao | Figure 4.11: Slope-Area Plots: Queen Charlotte Islands (Rennel Sound Area) 61 Ui/lU 8dO|S |BQO| UI/IU 8dO|S |B00| ui/LU sdo|S |eoo| UI/UJ edo|s |BOO| Figure 4.12: Slope-Area Plots: Queen Charlotte Islands (Skidegate Channel Area) 62 Watershed Area Physiography Major Geological Slope-Channel Fractal km2 Formation Transition 9 vnr Dimension (Network) Riley 28.7 SP Yakoun 6000 1.57 Gregory 36.7 SP Masset 6000 1.56 Hangover 21.2 SP Masset 6000 1.62 Mountain 12.8 Q C R P T P 3000 1.61 Government 16.1 Q C R Karmutsen 6000 1.67 Armentiere 4.0 Q C R Karmutsen 6000 1.70 Jason 11.9 Q C R Karmutsen 6000 1.71 Security 33.9 Q C R Karmutsen 6000 1.54 Table 4.3: Drainage Basin Characteristics and Channel Network Fractal Dimensions: Queen Charlotte Islands (compiled from Rood, 1984) Key: SP=Skidegate Plateau; Q C R = Queen Charlotte Ranges; P T P = Post-tectonic Plutons (see Chapter Two) The procedure for zoning elevations by upslope area was carried out: fluvial points with upslope area > 6,000 m 2 (> 6 pixels) were removed; those pixels remaining were divided into two zones: upper hillslopes with area < 1,800 m 2 (< 2 pixels) and a lower hillslope zone with 1,800 m 2 < upslope area < 6,000 m 2 . These zones formed clip regions for the original B C T R I M spot heights, whose fractal dimensions were then analysed using mf s, with a separation range 60 m - 600 m (again, the lower limit set by the data resolution and the upper limit imposed to ensure that hillslope boundaries were not crossed, following Veneziano and Iacobellis 1999). The resulting variograms and moment-slope plots are shown in Figures 4.13 through 4.16 and summarized in Table 4.4. In all cases the fractal dimension was significantly greater than two, which is the base value that would be expected for planar slopes resulting from diffusive processes. In general, hillslope fractal dimensions on the Queen Charlotte Islands are higher than in the Capilano watershed, which is consistent with the wider range of mass-movement events and more rapid rates of sediment flux presented in Chapter Three. 63 I jueujouj Figure 4.13: Hillslope Variograms: Queen Charlotte Islands (Rennel Sound: Riley and Gregory Creeks) 64 II cr £ (OC4 CM C O II II X X cn cn co 3 ; U 8 L U O L U J U D L U O L U Figure 4.14: Hillslope Variograms: Queen Charlotte Islands (Rennel Sound: Mountain and Hangover Creeks) 65 iicc in o O CO c\i CM ll ll X X cofC CM CM CO * - * - LO "—'CM x in TJ ll II I I C O ^ J - Figure 4.15: Hillslope Variograms: Queen Charlotte Islands (Skidegate Channel: Ar- mentieres and Government Creeks) 66 E E 2 CNJ CO II II CO LT) CO O CO en co Figure 4.16: Hillslope Variograms: Queen Charlotte Islands (Skidegate Channel: Ja- son and Security Creeks) 67 Watershed Lower Slopes Upper Slopes N D 7o N D 7o Riley 1314 2.35 -0.20 2000 2.22 -0.17 Gregory 1818 2.30 -0.07 2000 2.23 -0.15 Mountain 805 2.25 -0.16 1368 2.18 0.00 Hangover 1027 2.24 0.01 1711 2.15 -0.27 Government 866 2.48 0.49 1449 2.28 -0.04 Armentiere 250 2.47 0.66 309 2.32 0.28 Jason 723 2.31 0.28 1288 2.17 -0.12 Security 1604 2.35 0.26 2000 2.23 -0.10 Average 2.34 0.16 2.22 -0.07 Table 4.4: Queen Charlotte Islands Hillslope Fractal Statistics There is an important difference between the characteristics of the Queen Char- lotte Islands slopes and those in the Capilano watershed: without exception, the upper hillslopes on the Queen Charlotte Islands have a lower fractal dimension than the lower hillslopes. For some watersheds the difference is less than the precision of the mea- surement, but for many it is clear. There may be several explanations for this result, which is counter-intuitive and opposite to the G V R D result: 1. there is a lower incidence of large structure-generating landslides on the Queen Charlotte Islands than in the G V R D mountains. This is supported by Alley and Thompson (1978) and Ryder (pers. comm.), who have found only localized areas of deep-seated bedrock landsliding on the Queen Charlottes. 2. weathering rates are higher on the Queen Charlotte Islands, so slide evidence is obliterated more rapidly. 3. there is more hillslope activity on the lower slopes of the Queen Charlotte Islands. However, the hillslope processes expected to act on these slopes are predomi- nantly shallow landslides, which have been assumed to smooth the surface over 68 long timescales. Their cumulative effect is considered to be diffusive, so they ought to reduce the fractal dimension rather than increase it. The increased fractal dimension observed in this process-domain may be a result of the "noise" generated by individual landslides. The smoothing effect of diffusion is revealed only after a long time; meanwhile, individual landslides (which are power-law distributed over a finite range) increase the fractal dimension. 4. human activity (logging) has increased the tendency for midslope regions to fail, leading to a higher D in those areas. Of these explanations, (1), (2) and (3) seem the most realistic, since (4) implies that its effect will be seen only in areas with high areal proportion of logging, which is not the case. Point (3) is significant because it implies that, while it may correctly reproduce process rates, diffusion may be less successful at reproducing the resulting landscape form. A further explanation which may be important in the present context relates to the different glacial history of the Queen Charlotte Islands (see Chapter Two). This region was not fully glaciated, leaving higher elevations ice-free. Conversely, the G V R D area was fully glaciated. The absence of ice from upper parts of the Queen Charlotte Islands landscape may therefore be responsible for the lower fractal dimension seen there, compared with the glaciated lowlands. 4.6 Discussion The two areas studied in this chapter share several fractal characteristics: the most elementary is that hillslopes in the study area have a fractal dimension consistently higher than 2.0, which is the base dimension that would be expected from planar 69 slopes. This is contrary to the widespread view that hillslope processes destroy the fractal properties of topography; The present work has found that hillslope processes have distinct fractal characteristics when studied separately from channel networks. This evidence supports the view of Veneziano and Iacobellis (1999) that confounding hillslope and channel behaviour leads to multifractal scaling similar to that observed by Lavallee et al. (1993) and Rodriguez-Iturbe and Rinaldo (1997). It is at present unclear whether the hillslope part of the landscape is truly fractal. It does not appear to have multifractal characteristics in this study, as shown by the linear relation between moment order and moment slope for higher general moments of the variogram. However, the possibility remains that the hillslope region is a pseudo-fractal, owing to the small range over which the observations have been shown to scale (one order of magnitude). If this is the case then the fractal appearance of the landscape should be interpreted as a finite scaling effect. A distinction can be made between the findings of Veneziano and Iacobellis (1999) and those of the present study. Veneziano and Iacobellis found mono-fractal channel networks and hillslope regions with D close to 2.0. In the areas studied here, a mono-fractal channel network was found, but D for the upper and lower hillslope regions is significantly different from 2.0. This implies that in this environment (a recently glaciated, collisional orogen), there is a separate process responsible for in- creasing the hillslope fractal dimension. It is suggested here that power-law-distributed landslides may be responsible for the observed fractal properties of hillslopes. Over sufficiently long timescales the effect of shallow landslides may be generalized, yielding a slope-dependent diffusive transport relation (e.g., Culling, 1960; Martin and Church, 1997). However, for in- stantaneous digital elevation samples, the effects of the small-scale process that make up the diffusive behaviour may introduce detectable noise (as in the case of the Queen 70 Charlotte Islands). The effect of large events such as deep-seated rockslides is more permanent: they cannot be treated diffusively unless the timescales of the study are very long (~106 a - 109 a). In the final chapter of this thesis, these conjectures will be tested using a numer- ical model of landform evolution that incorporates hillslope diffusion, fluvial erosion and stochastic deep-seated bedrock landsliding. 71 Chapter 5 Synthesis: Process-Based Mode l l i ng w i th Episod ic Events 5.1 Int roduct ion The formulation of integrated, multiple-process models of the long-term evolution of landscapes has been a recurrent goal in geomorphology (Davis, 1899; Penck, 1924; Culling, 1960; Merrits and Ellis, 1994; Martin, 1998; Beaumont et ai, 2000). A recent trend has been to develop models that incorporate generalized physical representations of measurable small-scale processes, seeking to demonstrate their cumulative effects over much longer timescales. This is an attractive approach because it allows models to be calibrated with respect to time at the outset, a feature that was deficient in earlier models (especially those of Davis). The use of a landscape evolution model in the present work is to simulate the projected effects of measured rates of episodic processes over timescales longer than the available record. In addition, it is possible to control the processes incorporated in the model in order to decide more conclusively which are responsible for the resulting structure of the landscape. In this chapter the use of such a model is described. Land- scape evolution equations are formulated and solved numerically using a computer, in 72 order to determine the effects of various landscape processes (episodic and otherwise) on the fractal dimension of topography. This is a strictly uniformitarian approach, since it assumes that process rates calibrated from contemporary field observations will persist over the timescales of interest. 5.2 Const ruc t ing Testable Landscape Mode ls A numerical model of landform evolution that includes several processes is a type of complex hypothesis (Oreskes et al., 1994); the scientific use of this type of model is therefore worthy of discussion. The following example serves to illustrate the problems involved in testing a multiple process landform evolution model: Suppose that a geomorphologist constructs a model (complex hypothesis) of the form: "there are processes P, Q, R, and S, which are linked in the way T, and which produce the outcome X, where X is an observable property of the resulting landscape". (For clarity, one may suppose that P is soil creep, Q is shallow landsliding, R is fluvial erosion, S is rock sliding and T is a means of linking the processes such as a continuity equation or sediment budget. X is a measurable geomorphometric property such as the average slope, hypsometry or fractal dimension that arises from the modelled processes.) Verification of a model such as this conventionally involves a demonstration of its truth (i.e., that it is a good representation of reality), most obviously through a correct prediction of the property X. It is widely acknowledged that a model or hypothesis is not conclusively verified by a correct prediction; its truth is merely made more likely. Falsification of a theory, model or hypothesis is generally taken to have a more conclusive effect on the device's perceived veracity (Haines-Young and Petch, 1986, give a clear illustration of these claims in the context of physical geography). 73 However, as Oreskes et al. (1994) explain, falsification is also problematic because even simple single-hypothesis models are never closed systems; they rely on theoretical terms that are laden with inferences and assumptions. As a consequence, any failure of the model to agree with empirical observations may be explained and eliminated by modification of an auxiliary hypothesis or initial condition. This is an example of the problem of the underdetermination of theories by evidence (i.e., the Duhem-Quine thesis, for discussion of which see Duhem 1905, Quine 1953). In short, it is not possible to devise a critical experiment. At most, one can enumerate all the recognized processes that might be responsible for creating fractal topography and construct sufficiently critical tests between them, ultimately to arrive at a tentative claim that one process combination is more (or less) successful at pro- ducing the relevant observed feature of topography. This is a more rigorous scientific method than testing a fixed suite of processes against a fixed set of observables in order to decide wherein lies the ultimate truth. In the present work, the following are compared: diffusive hillslope processes, fluvial erosion and deep-seated bedrock land- slides. In the remainder of this chapter, the methods used to represent these processes will be explained. 5.3 Framework for Landscape Mode l l i ng Two common frameworks for long-term landform evolution modelling are discrete and continuous models1. A n example of a discrete model is a cellular automaton, in which a set of rules is applied to cells in a (usually regular) grid (Bak et al., 1987; Chase, 1994; Rodriguez-Iturbe and Rinaldo, 1997). In an atomist world discrete models ought to rule, however continuous models are often applied in situations where particles or 1 Naturally, all numerical models must be solved in discrete space and time if a (conventional) computer is used; the distinction at issue here is the framework in which the model is formulated. 74 quanta are very small compared to the spatial scales of interest. Owing to the large size difference between landscape particles and river basins, a continuum approach will be taken here; it is described below2. 5.3.1 Continuum Model of Landscape Evolution The most axiomatic step in the construction of a continuum model of landscape evo- lution is to consider elevation as a scalar field z(x,y), where z is elevation and x,y are mutually-orthogonal horizontal space variables. Following this initial step, it is possible to derive a transport theorem for elevation (as if it were energy, mass or momentum), based on the assumption that volume is conserved3. The local form of this (Eq. 5.1) balances the change in elevation over time at a point, (dz/dt), with the introduction to (or removal from) the model domain of landscape material, $, and the flux of sediment from that point, (q = —K(z, zXi,t)dz/dXi): + $ u + $F + ri(z,Zxi) (5-1) where t = time; X{ = space variables; K = diffusivity function; and zXi = dz/dxi. The production function, $, is divided into three parts, each to represent a different landscape process: §u represents uplift, §p represents fluvial erosion, and rj(z,zXi) represents stochastic deep-seated bedrock landslides. 2The distinction between discrete and continuous models is therefore one of convenience: Ger- shenfeld (1999) demonstrates analytically (using statistical mechanics) that a gas-lattice cellular automaton (discrete) exhibits exactly the behaviour of a fluid obeying the Navier-Stokes equations (continuous). A statistical approach is interesting in the context of modelling sediment transport in geomorphology because, although this is an area where particles are often much smaller than the space-scales of ultimate interest, the results of non-trivial interactions between these particles has been difficult to capture using traditional continuum mechanics. 3See Appendix B for a full derivation of this and subsequent results stated in this section. The assumption of conservation of volume is necessary to recover landscape form. Over longer timescales, it is more appropriate to enforce conservation of mass, however landscape form would not easily be modelled using this approach. Since morphology is central to the present work, a volume-conservative approach is more suitable here. dz dt _d_ dxi 75 The individual components of Equation 5.1 are discussed below. Diffusion The first term represents the flux of sediment away from the point: it is a diffusive term and can be used to account for processes that can be averaged over the time-scales of interest. These include creep, shallow landsliding and debris flows (see Chapter Three for justification of diffusion in the present study area). Physical intuition suggests that the rate of sediment transport by these processes should be related to slope, since they are driven primarily by gravity. The diffusivity function can be a constant, or a function of space and time vari- ables and their derivatives. Empirical work by Martin and Church (1997) and Martin (1998, 2000) has shown that a suitable relation between slope and sedi- ment flux is a non-linear sigmoidal curve (see Chapter Three). This is intuitively appealing because it suggests that, while sediment flux is low at low gradients, it increases rapidly until slopes reach approximately 45 degrees. Above this angle, failure is through different mechanisms that cannot be approximated by diffu- sion over geomorphic timescales. Martin (2000) has shown that an hyperbolic tangent function gives a suitable fit, which she has parameterized using data from the Queen Charlotte Islands. Details of a parameterization of the same function for the G V R D Capilano watershed are given in Chapter Three. Upl i f t Tectonic and isostatic uplift may be included as a production term. In the present study uplift is included only in long model runs; uplift rates are poorly constrained compared with the detailed denudation observations used in this work, therefore in shorter runs (~ 10 ka) the effects of uplift are not studied. F l u v i a l Eros ion Fluvial erosion and deposition are included in the model based on a functional relation between slope, S, and stream discharge, Q, obtained by 76 simplifying the Bagnold equation presented by Martin and Church (2000). The simplified version, derived in Appendix B, is: 6.44Q 0 3 5 S1-5 if S < 0.01 (5.2) $F = < 40.32 Q 0 3 4 S2-22 if S > 0.01 Following Martin (1998), a correction is made to account for configurational stability in steep headwater channels (i.e., those with S > 0.01); ordinary fluvial processes are assumed not to operate where S > 0.1. At regular intervals the drainage network is rerouted using the D8 flow-routing algorithm (Jenson and Domingue, 1988). Discharge is computed from upslope area, A, using the results of Church (1997): where kfi is a parameter to represent hydroclimate. Deep-Seated. Landslides In order to incorporate the largest landslides, the effects of which cannot be averaged over the time-scales of interest, the stochastic term rj(z, zXi) is added. The use of a stochastic term to represent large landslides is justified on the basis that the precise location of these events cannot be determined from variables resolved in the model. Large, deep-seated bedrock landslides occur as a result of a combination of weathering, cumulative hydrometeorologic conditions and seismic vibration; their location is determined by the distribution of these forcings in space and time. It it not practical to resolve these forcings over the time-scales of the present study. In addition, the exact location of rock slides is contingent Q = kflA{ 0.68 (5.3) 77 on local features of the bedrock configuration, such as bedding planes and the stratigraphic sequence, which are not codified or analysed. Despite these claims, it is possible to constrain the areas most likely to fail on the basis of local slope, relief and elevation. Consequently, in the present scheme each node is assigned a probability of failure, Ps\, based on local slope, relief and elevation criteria. At regular intervals, -IB\, model failure sites are chosen at random from a subset of model nodes where Ps\ exceeds a critical value, PCnt- The failed material is then routed downslope in the direction of steepest descent so that its angle of repose is (j)t. The stochastic term makes Equation 5.1 a non-linear Langevin equation similar to the type that Sornette and Zhang (1993) discuss. This type of equation is well known to create surfaces with interesting fractal properties. Rodriguez-Iturbe and Rinaldo (1997) suggest a similar formulation, but derive the stochastic forc- ing from hydrological parameters rather than from hillslope-process measure- ments. 5.3.2 Parameters As stated in Equation 5.1, the model involves several parameters. These are listed in Table 5.1, along with methods for their estimation from empirical data, and order of magnitude estimates (exact values used in model runs are given in Section 5.4). 78 62 CD" i - l P 3 CD c t - CD i-i H co >—1 • 3 P c t - CD CO CL CO o P as o o P o CL CD > 2 O o P [inim u i-i _ cr [inim u p CL ST o [inim u o' £L nne rocl CD 3 w CO 3 CL CD p c t - CD o CL CO CD T) O CD i-i CO Inte oba c t > CL CD rva cr •—1 CD P de c t - cr >-!_ co' 1 CD i - l P TO CD CD P i - l P 3 CD c t - CD i - l to CO o h-1 O h-1 o o o o CL CD o ra 3 t—' TO p CL to Cn sr ^ CL 3 CL rjq sr |-* OS ps r- p - CL CL CD 3' < CD p p i ffi H 2 5' ^ h-» P " S" t—i 3 ^ co CD CD o o o CD co 4^ 8 c f O P O p co ra. „ co O o p- i = tr tr Pr " o >—1 sr CD 5 co 93 ^ P5 o o to CO O p ?r co a? 3- 2. TO CO. o & 3' o o ^CL pj o to p CD 00 CO , CO 3 3f co 00 p- ^ O c t - h-i o CO co CO 03 i ? 5 CL S p i P CD 1-1 sr c t - 5' o - cr I—L CD CD £L CD ^ 00 p ' ~_ TO O < O co CD co 3 o CO p 1-1 o CD o p c + P 5.3.3 Canonical Scaling Using the order-of-magnitude estimates given in Table 5.1, and the structure of the model in Equation 5.1 it is possible to arrive at preliminary conclusions regarding the relative importance of the components of the model by making the individual terms of Equation 5.1 dimensionless. This is achieved with the following scalings: z = [z]z*, t = [t]t\ xt = [Xi]xl K = [K]K*, $u = [^u]^u, $ F = [$F]$*f, V = [v]v*- Applying these to Equation 5.1 gives: from which the natural time constant for diffusion can be obtained. The time (in years) for elevation perturbations in a 1 km 2 area to become completely flat, assuming that only diffusive processes are operating, is: Further results can be derived from Equation 5.4, including an expression that com- pares the relative roles of uplift and diffusion using the quotient of their respective characteristic scales: I r , f t p ,. [$r/]fo]2 (0-002)(1000)2 Uplift Ratio = = ( 1 0 Q 0 ) ( 0 2 ) = 10 (5-6) In this case, the relative importance of uplift is greater than 1.0, indicating that landscapes in the present study-area are dominated by uplift, which is consistent with their tectonic setting and current relief. 80 5.3.4 Boundary Conditions Conditions must be specified for the geometric boundaries of the model. Typical choices include constant-value boundaries, constant-derivative boundaries and combi- nations of these. Since neither fixed elevation nor fixed flux is a physically correct description of the situation a more complex set of boundary conditions is required. The model will be run for individual watersheds, so the physical feature at the bound- ary will be a ridge. The most realistic model boundary may therefore be a "reflection" of the topography within the model domain. Using the method of images to produce a reflected ridge is an appropriate approach; it is also necessary to include a larger area of D E M than is of interest (i.e., to include a buffer zone within which results are computed but then discarded). In this way, the natural structure of the landscape isolates the area of interest from spurious effects at the model edge. 5.3.5 Numerical Schemes For simple cases with only linear hillslope diffusion, the problem can be solved ana- lytically, yielding equilibrium solutions (e.g., Kirkby, 1971). Owing to the complex topography of a real drainage basin, a numerical approach will be taken here. The choice of numerical schemes is first between finite difference and finite el- ement models. The former discretize the governing equation on a regular grid; the latter use non-rectangular facets (commonly but not necessarily triangles). The ad- vantages and disadvantages of each method are similar to those between lattice and triangular irregular network (TIN) representations of digital elevation data (see Kum- ler, 1994, for a recent comparison). In general, lattices (finite difference schemes) are simpler and easier to manipulate, however TINs (finite element schemes) are able to model complex surfaces more economically. The TIN structure is also appropriate for 81 modelling with different resolutions on hillslopes and channels. However, algorithms for flow-routing on TINs are in their infancy, which make them less suitable for simple models at present (McAllister and Snoeyink 1999; McAllister pers. comm., 1999). Consequently, a finite difference scheme was employed for this work. 5.3.6 Stability A further choice is between explicit and implicit solutions on the finite difference grid. The conditionally-stable explicit solution takes the form of a Taylor series expansion of the differential equation resulting in an equation for the height of each node in the domain at time t + 1, which may be written solely in terms of the values of its neighbouring nodes at the earlier time, t. Conversely, the implicit method produces an equation at each node for time t + 1, which is in terms of the value of that node at time t, and at several neighbouring nodes at time t + 1. Since, in this case, the future value of each node depends on the future value of every other node, the implicit method requires the solution of a set of simultaneous algebraic equations. Imaginably, this method is more computationally demanding especially for the non-linear case, but it is unconditionally stable. Owing to its simplicity, an explicit solution is used here. In order to achieve this, it is necessary to use a forward-difference approximation for the time-derivative at the same time as a central-difference approximation for the spatial derivatives. This introduces errors of different orders of magnitude, and a problem of instability arises if the model is advanced for large values of At (although the solution is accurate for small values of At). The largest time-step for which the model can be run is therefore inversely proportional to the square of the chosen space step. Typically this time is much 82 smaller than the time constant for the problem (which is the time for any initial structure to be made uniform). Happily, in the present context, one is content to stop the model before the topography has disappeared, so one can tolerate the small time-step that an explicit scheme requires, particularly since the approach is very fast and a simple modification can take care of the non-linear case. The final C++ code is in Appendix B. 5.4 Resul ts : Pre l im inary Runs Preliminary runs were conducted to assess the initial plausibility of the model. Results from the most comprehensive run are shown in Figure 5.1, which is the final surface resulting from the evolution of a triangular prism (to represent a mountain ridge) with diffusion, fluvial processes and deep-seated bedrock landsliding over 100 ka. The initial surface was planar with no perturbations: all structure was derived from the erosion rules4. The results of this run show the emergence of linear channels, characteristic of the present study area (and others with high uplift ratios; e.g., parts of New Zealand, Taiwan and other Pacific Islands). In addition, 1,000 deep-seated bedrock landslides occurred over the 100 ka period. These have altered the drainage network across the model domain, but particularly in the central upland portion where significant tributaries have been created. An ensemble of 30 model runs was made, yielding an average magnitude- frequency distribution of rock slides, and allowing an estimate of its variability to be obtained. This distribution, which is shown in Figure 5.2, spans approximately 2 4This example is not intended to reproduce the exact topography in any natural location, since few landscapes begin as plane-faced triangular prisms, and it is not reasonable to expect that climate parameters or uplift rates will have remained the same throughout the 100 ka time period; they will almost certainly have changed. 83 orders of magnitude, however it does not follow a power-law. Figure 5.2 is similar to the distribution of the G V R D surficial slides (Fig. 3.1), showing a low-magnitude kink and a rounded central portion. Two suggestions follow from this: that the lower-end kink is an emergent feature of threshold exceedence failures, rather than an obser- vational inadequacy; and that the distribution of modelled rock slides is exponential rather than power-law. The observations presented for this model-run differ from the findings reported by Hergarten and Neugebauer (1999) in their simpler (but similar) model, and the experimental findings of Densmore et al. (1997). In these cases, the observation of power-law distributed rock slides has been used to suggest that landslides are a results of self-organized critical behaviour. A n alternative explanation is suggested here. In the present case, the rock slide interval was seen to have the most important effect on the resulting rockslide distribution. A further model ensemble was gath- ered with a higher rock-slide interval of 400 years. This group (Figure 5.3) has a magnitude-frequency distribution that more closely resembles a power-law (though it is not perfect). The effect of an increased slide interval is to reduce the total number of slides that occur, and to allow slopes to become steeper before they are likely to fail. The result of these differences is a reduction in the proportion of average-sized events, leading to a more extreme distribution. It is, at present, unclear whether the power-law-like distribution results from the limitation imposed on the number of rock slides in the latter run (Figure 5.3, and cf. the sampling limitation imposed on a human observer), or whether it is a result of an increase in the degree to which slopes may be steepened between slides. The latter situation is important in the present regional context since it implies that over- steepening introduced by rapid tectonic uplift and as part of the legacy of glaciation is a necessary condition for the observation of power-law (or perhaps exponentially) 84 distributed slope-failures. In the search for answers to the questions raised so far, it is anticipated that a modelling approach may prove useful. 85 s . s • % •i % •i % • % • S11 % •i • . " i i f." '• .•• j"« d"• ̂  • rf-• <• ^. W H S • Iii j> i i BVi'ijiVi'ii'i nftSiftft?1" i J J W J -,S • S • % • S • • [S-IS'S'SJ" Figure 5.1: Artificial Model Surface after lOOka 86 100 0.01 Frequency (Ensemble Average) 1 Standard Deviation 10000 100000 , 1e+06 Volume (m3) Figure 5.2: Model Rock Slides: Magnitude-Frequency (100 year slide interval; N 1000) 100 0.1 0.01 10000 Frequency (Ensemble Average) 1 Standard Deviation 100000 Volume (m3) 1e+06 1e+07 Figure 5.3: Model Rock Slides: Magnitude-Frequency (400 year slide interval; N 250) 87 5.5 Resul ts : Detai led M o d e l Runs A more realistic use of the landform evolution model presented here is to examine the predicted changes over a shorter time period for a real landscape. In this section, the results of such an application are discussed, and the fractal dimensions of various hypothetical process-domains are presented in order to approach an answer the pre- viously posed questions regarding the effect of episodic processes such as deep-seated bedrock landslides on the fractal dimension of topography. One sub-basin was chosen within each of the study areas: Lembke Creek, a 3.6 km 2 sub-basin in the western side of the Capilano watershed (see Fig. 2.3); and a 4.4 km 2 tributary of Gregory Creek, which is in the Rennel Sound area of the Queen Charlotte Islands (Fig. 2.2). Small basins were chosen in order to reduce topographic complexity and to preserve computational tractability. Table 5.2 gives the parameters used in each model run, the details of which are discussed below. A l l runs were for 10 ka, using the appropriate non-linear diffusion rule (see Chapter Three) on a regular square grid of 30 m on a side. The effects of tectonics were not included, since they obscure the subtleties of other processes. Watershed Hydroclimate Hillslope-Fluvial Rock Slide Threshold Church (1997) Threshold (m 2) (slope * elevation) m (see Chapter Four) (see Chapter Four) Lembke 0.15 2000 1000 Gregory 0.17 6000 450 Table 5.2: Parameters Used in Model Runs 88 5.5.1 Lembke Creek Results for Lembke Creek without deep-seated rock slides are shown in Figure 5.4. This map shows contours of the final surface. Changes between the start and finish of the model run are shown in reds (erosion) and greens (deposition)5. Areas of net fluvial erosion are shown in blue; net fluvial deposition is marked in grey (N.B.: values for fluvial changes have been divided over the whole cell in which the river lies (900 m 2)). The steep slopes to the north-west have yielded considerable amounts of sediment that has been transported onto the valley floor by headwater streams. Most of this sediment is from areas where slopes have been over-steepened by glaciation or uplift. The effect of fluvial processes is to homogenize the valley flat, by redistributing sediment that is supplied from the hillsides. The area to the south-east shows a lower amount of erosion by diffusive processes: it is a steep rock slope and therefore reasonably stable in this model run. Figure 5.5 shows the same situation as the previous run, but with the addition of stochastic deep-seated bedrock landslides. Slide heads are shown in yellow, and tracks are clearly visible as areas of high erosion. Deep-seated slides are common in the south-eastern rock slopes, where they have altered the drainage network and the rock face. Areas of deposition on lower slopes beneath large landslides (shown in green) may be interpreted as talus cones (or alluvial fans in the isolated cases along the central and eastern slopes of the southern valley wall). The interaction of channel and hillslope processes is especially noteworthy: channel processes are able to respond to changes in the topography induced by hill- slope processes (episodic and otherwise) because the channel network is recomputed from the topography after each rock slide (every 100 years). The effect of this is to 5Shading and symbols are used in place of colour on black-and-white archival versions; in each case the appropriate key is given on the figure. 89 re-route flow into the crenulations created by rock slides, and to allow the passage of colluvial deposits into river channels. Examples of the former can be seen along the southern valley wall; of the latter on the most southerly talus cone. The primary reason for having created this model is to assess the effect of the erosion rules on the fractal dimension of the topography. Table 5.3 shows the fractal dimensions of the channel networks at 2,500 year intervals under the influence of diffusion and fluvial processes and with the addition of rock slides. It also shows the division of the non-fluvial part of the landscape into hypothetical process-domains (upland areas with upslope area 0 - 1,800 m 2 ; and lowland areas with upslope area 1,800 - 5,400 m 2 ) . The procedures used for the estimation of these parameters were identical to those used in Chapter 4 for the original D E M data. Time Dcfiannei ^upland Dlowland a With Rock Slides 0 1.62 2.17 2.21 2500 1.64 2.14 2.19 5000 1.46 2.13 2.19 7500 1.46 2.13 2.19 10000 1.50 2.14 2.19 No Rock Slides 0 1.62 2.17 2.21 2500 1.74 2.15 2.20 5000 1.69 2.15 2.20 7500 1.71 2.15 2.20 10000 1.58 2.15 2.20 Table 5.3: Fractal Dimensions of Model Results 90 91 8 CO O 9 2 The fractal dimension of the channel network responds as expected: in the case without rock slides, the channel network fractal dimension fluctuates about an average value of 1.67. The addition of rock slides liberates large quantities of material from rock slopes, filling the channel network and restricting the areas where evidence of fluvial processes may be found. This leads to a systematic decrease of the fractal dimension of the channel network over the course of 10 ka. The effect of hillslope processes is less obvious. Of particular note is that the upper hillslopes in each case have a lower fractal dimension than the lower hillslopes. This is the opposite of the result found in Chapter Four. This discrepancy may have resulted from the necessity to resample elevation data onto a coarser grid for modelling work in order to relieve computational constraints (this was achieved using A R C / I N F O s cubic interpolation algorithm). The effects of D E M resolution on the determination of fractal properties of landscape have not, to my knowledge, been examined in the literature. The effect of diffusion, which can be seen most clearly in Figure 5.4, is to smooth the surface; this is captured in the fractal dimension, which declines from 2.17 to 2.15 (upland) and from 2.21 to 2.20 (lowland) over the 10 ka period. The change may be smaller than the assumed precision of the technique for estimating D; nevertheless, the decline in D is consistent. The addition of bedrock landslides has a roughening effect visible in Figure 5.5, however they cause a decrease in the fractal dimension, indicating a relative smoothing effect. There may be several explanations for this result: (i) bedrock landsliding occurs over only a limited range of scales and therefore destroys part of the fractal scaling introduced by fluvial processes; or (ii) the traces of landslides quickly become part of the channel network and therefore are excluded from the fractal analysis presented here (using the method of Veneziano and Iacobellis, 1999). The first explanation is 93 consistent with the idea of pseudo-fractal topography suggested by Burrough (1993); it implies that, while individual landslides decrease the fractal dimension, their cumu- lative interference may lead to an increase in the fractal dimension over much longer time-scales; 10 ka is therefore not sufficient to create pseudo-fractal landscapes. The outright acceptance of the second explanation is difficult because, if it were true, one would expect deep-seated landslides to contribute to an increased channel network fractal dimension, which is not the case. A third explanation might therefore be that the deep-seated landslides create a different type of sediment flow-path that is missed by the analysis presented here. The landslide flow paths have a higher slope than would be expected if they were fluvial channels, but still have sufficiently high upslope areas that they are influenced by the flow of water. These configurationally- stable channels on steep slopes do not fit well into the distinction between hillslope and fluvial processes, and further work is necessary to enable these hypotheses to be tested. The first steps are likely to involve, better representation of channel processes on slopes greater than 0.1, since they were not considered in the present study. 5.5.2 Gregory Creek Results from Gregory Creek are shown in Figure 5.6. The most notable feature here is that although the landscape was tested for areas likely to slide with the same regularity as the model of Lembke Creek (above), the occurrence of deep-seated landslides was confined to the most southerly area where two small rock slides occurred; they were not large enough to create a contour crenulation. This was found to be the case for all physically realistic values of the rock slide stability threshold, suggesting that there is an underlying topographic reason for the lack of large rock slides in Gregory Creek. This conjecture is supported by the evidence of Ryder (pers. comm., 2000), whose 94 recent aerial-photograph inventories of slope processes on the Queen Charlotte Islands show little evidence of large deep-seated landslides. Most slope failures in Gregory Creek are shallow landslides in surficial material: for long-term landform modelling in this area non-linear diffusion alone can account for all sediment flux from hillslopes. The fluvial network shown in the results from Gregory Creek contains a greater amount of deposition than that for Lembke Creek. In particular there is a large area of deposition at the confluence of several streams as they emerge from the upland area in the eastern part of the map. The reproduction of this feature, which represents an alluvial fan, gives confidence in the way that fluvial processes have been incorporated into the model. The fractal properties of this Gregory Creek tributary were not studied because the rarity of rock slides means that there is no basis for comparison between models with and without them. 5.6 Discussion The results presented in this chapter show that diffusive hillslope processes lead to a reduction in the fractal dimension of topography. In cases where rock slides do not occur, the fluvial network removes material in such a way as to preserve the fractal dimension of the channel network. The addition of episodic events has an effect that is more complex than simply roughening the surface (although roughening is clearly observed in Figure 5.5). Indeed, the fractal dimensions of hillslopes in model runs with structure generating landslides were lower than those of hillslopes without these perturbations. Two possible reasons for this result are suggested. The first is that while slope failures often appear to be power-law distributed, this holds over only a limited range 95  of scales; therefore rockslides do have a characteristic range of length scales and their occurrence leads to a reduction in the fractal dimension of parts of the landscape in which they occur. This may, over long time periods, create a pseudo-fractal landscape of the type described by Burrough (1993). A n alternative explanation is that rock slides quickly modify the drainage net- work by creating steep headwater channels or gullies, which are not identified on slope-area plots (even using the procedure of Veneziano and Iacobellis) because they are too steep to be channels, yet have too large a contributing area to be hillslopes. The solution of this problem is non-trivial and would require a study of the morphom- etry of headward channels. Further work on hillslope-fluvial process coupling is also necessary in order to decide between the these two explanations. This has also been suggested by Martin (1998), Hovius et al. (2000) and Beaumont et al. (2000), owing to the realization of the importance of headwater channels in long-term landscape evolution, which is summarized in Church (in press). 97 Chapter 6 Conclus ion The work in this thesis has attempted to test two hypotheses concerning geomorphic processes and the fractal properties of topography: (i) that there exist areas of the landscape with distinct fractal dimensions; and (ii) that these fractal domains are created by the operation of characteristic geomorphic processes. Several conclusions can be made. In Chapter Three, landslides on the Queen Charlotte Islands were shown to have a power-law magnitude-frequency distribution similar to those observed in other regions by Hovius et al. (1997), Pelletier et al. (1997) and Turcotte (1999). The slope-volume relation for a landslide inventory in the G V R D Capilano River basin was shown to be non-linear, in agreement with the hyperbolic tangent relation found by Martin (1998, 2000) for the Queen Charlotte Islands landslides. In Chapter Four, it was shown that evidence of power-law landslide statistics may be reflected in the non-integer fractal dimension observed for hillslopes in both the G V R D and Queen Charlotte Islands watersheds. Using a modified version of Veneziano and Iacobellis' (1999) procedure for separating areas dominated by hillslope and channel processes, it was shown that, in several cases, hillslope process-domains have a higher fractal dimension. In all cases, there is a consistent difference between the fractal dimension of topography in the hillslope and fluvial domains, supporting 98 the first hypothesis that process domains have distinct fractal dimensions. It is further suggested that details of the spatial variability of the fractal dimension depend on the glacial history of individual areas. The second hypothesis, that spatial differences in fractal dimension can be re- lated to geomorphic processes, was examined in Chapter Five. Results from a land- scape evolution model show that smoothing by diffusion is observed in the fractal dimension of simulated topography. Large episodic landslides are shown to have a short-term impact on the perceived "roughness" of the surface, however the pertur- bations that they create are quickly incorporated into the drainage network, where they appear as steep headwater channels (which are a common feature of landscapes in the area of study). The influence of deep-seated bedrock landslides on the fractal dimension of topography remains unclear: landslides reduce the fractal dimension in the short term, but their cumulative effect may be to increase it over longer periods of time, creating a pseudo-fractal landscape. In conclusion, the second hypothesis requires further work before it can be fully accepted. It is suggested that a numeri- cal modelling approach is appropriate, but that further attention should be given to geomorphic processes in mountain streams. 6.1 Prac t i ca l Implicat ions The practical utility of the findings reported here is in the area of long-term planning for land-use management and risk mitigation. The magnitude-frequency distributions presented for slope failure may be used to assess the likelihood of failure of a given volume over a representative time period. The landform modelling approach used in Chapter Five is also of relevance to land-use planning: its results are not deterministic predictions of landscape change; rather, they indicate areas of potential long-term 99 instability and allow the identification of important fluvial and hillslope sediment pathways. At time-scales longer than 100 years deterministic predictions are of limited use: generalized probabilistic predictions are more informative, and of greater utility in applications in which a measure of the uncertainty of predicted events is required. A probabilistic approach is especially useful in view of non-steady-state variability in erosion rates over long timescales (cf. Benda, 1994). In this respect, the ensemble approach taken in Chapter Five is useful: the model results shown there are generalized probabilistic predictions, with an attached measure of uncertainty. A n important line of further work in the use of the current approach is to assess the ways in which modification by human activity can be included. Initially, this might include modification of slope-failure rates resulting from timber harvesting or other changes from natural land use. Modification of fluvial processes in response to climate changes is also possible, although the correct calibration of such scenarios may present difficulty. 6.2 Theoret ica l Impl icat ions The theoretical implications of the work presented here concern the use of fractal prop- erties of topography to understand geomorphic processes. It has been confirmed that splitting landscapes into fluvial and hillslope process-domains can lead to the creation of areas with homogeneous fractal dimension, however some landscape features do not fit well into the hillslope-fluvial distinction. These include head ward channels and the remains of deep-seated bedrock landslides. It is also possible that landscape processes with a limited range of length scales can create non-fractal areas of landscape in the short term, yet over a long time their influence may lead to pseudo-fractal effects, as suggested by Burrough (1993). 100 A fundamental need in contemporary geomorphology is to relate fractal land- forms to geomorphic processes. Several recent attempts have produced wide-reaching theories of drainage-basin-scale fluvial geomorphology (these include Rodriguez-Iturbe and Rinaldo, 1997; Veneziano and Iacobellis, 1999). Theories of long-term landscape evolution are essential to understanding the ways in which information is transferred across scales in the landscape. It is important that they include processes of regional significance, such as deep-seated bedrock landsliding, the controls of which are contin- gent on factors such as bedrock geology and landscape history. 101 Bib l iography Alley, N . F. and Thompson, B . (1978). Aspects of Environmental Geology, Parts of Graham Island, Queen Charlotte Islands. British Columbia Ministry of Environ- ment, Resource Analysis Branch (Bulletin No. 2), Victoria, B C , Canada. Avnir, D., Biham, 0., Lidar, D., and Malcai, 0 . (1998). Is the geometry of nature fractal? Science, 279(5347), 39-40. Bak, P., Tang, O , and Weisenfeld, K . (1987). Self-organized criticality: an explanation of 1/f noise. Physical Review Letters, 59, 381-384. Bak, P., Tang, C , and Weisenfeld, K . (1988). Self-organized criticality. Physical Review A, 38(1), 364-374. Barnsley, M . F. (1988). Fractals Everywhere. Academic Press, San Diego, C A . Beaumont, C , Kooi, H. , and Willet, S. (2000). Coupled tectonic-surface process mod- els with applications to rifted margins and collisional orogens. In M . A . Summerfield, editor, Geomorphology and Global Tectonics, pages 29-55. John Wiley & Sons Ltd., Chichester. Benda, L. E. (1994). Stochastic Geomorphology in a Humid Mountain Landscape. Ph.D. Thesis, Department of Geological Sciences, The University of Washington. 102 Burrough, P. A . (1981). The fractal dimension of landscapes and other environmental data. Nature, 294, 240-242. Burrough, P. A . (1993). Fractal and geostatistical methods in landscape studies. In N . S. N . Lam and L. de Cola, editors, Fractals in Geography, pages 87-121. Prentice Hall, Englewood Cliffs, New Jersey. Campbell, D. (1999). Sediment Sources and Transfers in Lynn Valley, B.C.. B. Sc. Honours Thesis, Department of Geography, The University of British Columbia. Chase, C. G. (1994). Fluvial landsculpting and the fractal dimension of topography. Geomorphology, 5, 39-57. Church, M . (1997). Regionalised Hydrological Estimates for British Columbia: First Approximation of Scale Effects. Resource Inventory and Data Management Branch, British Columbia Ministry of Environment, Lands and Parks, Victoria, Canada. Church, M . (in press). Mountains and montane channels. In R. J. Allison and T. P. Burt, editors, Sediment Cascades: An Integrated Approach. Church, M . and Mark, D. M . (1980). On size and scale in geomorphology. Progress in Physical Geography, 4, 342-390. Clague, J. J. (1989). Quaternary geology of the Canadian Cordillera. In R. J. Ful- ton, editor, Quaternary Geology of Canada and Greenland, pages 15-96. Geological Survey of Canada, Ottawa. Cruden, D. M . (1985). Rock slope movements in the Canadian Cordillera. Canadian Geotechnical Journal, 22, 528-540. Culling, W. E. H. (1960). Analytical theory of erosion. Journal of Geology, 68, 336-344. 103 Culling, W. E. H . (1988). A unified theory of particulate flows in geomorphic settings. Earth Surface Processes and Landforms, 13, 431-585. Davis, W. (1899). The geographical cycle. Geographical Journal, 14, 481-504. Day, T. J. (1969). The Channel Geometry of Mountain Streams. M.Sc. thesis, The University of British Columbia. Densmore, A . L. , Anderson, R. S., McAdoo, B. G., and Ellis, M . A . (1997). Hillslope evolution by bedrock landslides. Science, 275(5298), 369-372. Duhem, P. (1905). Physical Theory and Experiment. Princeton University Press, Princeton, (trans. P. P. Weiner, 1954). English, R. R. (1996). Lineament Control on Drainage Basin Development, Large Rock Landslides and Mountain Slope Deformation in the SW Coast Mountains, BC, Canada. M.Sc. thesis, Department of Geological Sciences, The University of British Columbia. Evans, I. S. and McClean, C. J. (1995). The land surface is not a multifractal. Zeitschrift fur Geomorphologie N.F. Suppl.-Bd., 101, 127-147. Gershenfeld, N . A . (1999). The Nature of Mathematical Modelling. Cambridge Uni- versity Press, Cambridge. Gimbarzevsky, P. (1988). Mass Wasting on the Queen Charlotte Islands (Land Man- agement Report No. 29). British Columbia Ministry of Forests and Lands, Victoria, Canada. G V R D (1999). Greater Vancouver Regional District Watershed Management Report No. 5. Ecological Inventory Programme, Vancouver, Canada. 104 Haines-Young, R. H. and Petch, J. R. (1986). Physical Geography: its nature and methods. Harper and Row, London. Hergarten, S. and Neugebauer, H. J. (1999). Self-organized criticality in landsliding processes. In S. Hergarten and H. J. Neugebauer, editors, Process Modelling and Landform Evolution, pages 231-251. Springer-Verlag, Berlin. Hewitt, K . (1989). The altitudinal organization of Karakoram geomorphic processes and depositional environments. Zeitschrift fur Geomorphologie N.F. Suppl.-Bd., 76, 9-32. Horton, R. E. (1945). Erosional development of streams and their drainage basins: Hydrophysical approach to quantitative morphology. Bulletin, Geological Society of America, 56, 275-370. Hovius, N . , Stark, C. P., and Allen, P. A . (1997). Sediment flux from a mountain belt derived by landslide mapping. Geology, 25(3), 231-234. Hovius, N . , Stark, C. P., Tutton, M . A. , and Abbot, L. D. (1998). Landslide-driven drainage network evolution in a pre-steady-state mountain belt: Finisterre Moun- tains, Papua New Guinea. Geology, 26(12), 1071-1074. Hovius, N . , Stark, C. P., Chu, H.-T., and Lin, J.-C. (2000). Supply and removal of sediment in a landslide-dominated mountain belt: Central Range, Taiwan. Journal of Geology, 108, 73-89. Hungr, O., Evans, S. G., and Hazzard, J . (1999). Magnitude and frequency of rock falls and rock slides along the main transportation corridors of southwestern British Columbia. Canadian Geotechnical Journal, 36(2), 224-238. 105 Jenson, S. K . and Domingue, J. O. (1988). Extracting topographic structure from digital elevation data for geographic information system analysis. Photogrammetric Engineering and Remote Sensing, 54(11), 1593-1600. Kadanoff, L . P. (1986). Fractals: Where's the physics? Physics Today, February, 6-7. Kirkby, M . J. (1971). Hillslope process-response models based on the continuity equa- tion. Transactions, Institute of British Geographers, Special Publication No. 3. Klinkenberg, B. and Goodchild, M . F. (1992). The fractal properties of topography: A comparison of methods. Earth Surface Processes and Landforms, 17, 217-234. Klinkenberg, B. and Goodchild, M . F. (1993). Statistics of channel networks on frac- tional Brownian surfaces. In N . S. N . Lam and L. de Cola, editors, Fractals in Geography. Prentice Hall, Englewood Cliffs, New Jersey. Knighton, D. (1998). Fluvial Forms and Processes: A New Perspective. Arnold, London. Kumler, M . P. (1994). An intensive comparison of triangular irregular networks (TINs) and digital elevation models (DEMs). Cartographica, 31(2), 1-99. Lam, N . and de Cola, L. (1993). Fractal in Geography. Prentice Hall, Englewood Cliffs, New Jersey. Lavallee, D., Lovejoy, S., Schertzer, D., and Ladoy, P. (1993). Nonlinear variability of landscape topography: Multifractal analysis and simulation. In N . S. N . Lam and L. de Cola, editors, Fractals in Geography, pages 158-192. Prentice Hall, Englewood Cliffs, New Jersey. 106 Leir, M . C. (1995). Airborne Synthetic Aperture Radar, Digital Terrain Models and Geographic Information Systems: Tools for Mapping and Managing Large Land- slide Hazards in South-Western British Columbia. M.A.Sc. thesis, Department of Geological Sciences, University of British Columbia. Lian, O. B. and Hickin, E. J. (1993). Late Pleistocene stratigraphy and chronology of lower Seymour Valley, southwestern British Columbia. Canadian Journal of Earth Sciences, 30, 841-850. Mandelbrot, B. B. (1967). How long is the coast of Britain? statistical self-similarity and fractional dimension. Science, 156, 636-638. Mandelbrot, B. B. (1983). The Fractal Dimension of Topography. Freeman, New York. Mark, D. M . and Aronson, P. B. (1984). Scale-dependent fractal dimensions of to- pographic surfaces: A n empirical investigation with applications in geomorphology and computer mapping. Mathematical Geology, 16(7), 671-683. Martin, Y . (1998). Modelling Geomorphology in Long Term Landform Evolution. Ph.D. thesis, Department of Geography, University of British Columbia. Martin, Y . (2000). Modelling hillslope evolution: Linear and non-linear transport relations. Geomorphology, in press. Martin, Y . and Church, M . (1997). Diffusion in landscape development models: On the nature of basic transport relations. Earth Surface Processes and Landforms, 22, 273-279. Martin, Y . and Church, M . (2000). Re-examination of Bagnold's empirical bedload formulae. Earth Surface Processes and Landforms, 25. 107 McAllister, M . and Snoeyink, J. (1999). The geometry of watershed on TINs. Paper Presented at GIS'99, Vancouver, Canada. Mehta, A . (1992). Real sandpiles: Dilatancy, hysteresis and cooperative dynamics. Physica A, 186, 121-153. Merrits, D. and Ellis, M . (1994). Introduction to special section of tectonics and topography. Journal of Geophysical Research, 99(B6), 12,135-12,141. Middleton, G. V . and Wilcock, P. R. (1994). Mechanics in the Earth and Environ- mental Sciences. Cambridge University Press, Cambridge. Molnar, P. and England, P. (1990). Late Cenozoic uplift of mountain ranges and global climate change:chicken or egg. Nature, 346, 29-34. Montgomery, D. R. and Foufoula-Georgiou, E. (1993). Channel network source repre- sentation using digital elevation models. Water Resources Research, 29(12), 3925- 3934. Nikora, V . (1991). Fractal structure of river planforms. Water Resources Research, 27, 1327-1333. Oreskes, N . , Shrader-Frechette, K. , and Belitz, K . (1994). Verification, validation and confirmation of numerical models in the earth sciences. Science, 263 , 641-646. Outcalt, S. I. and Melton, M . A. (1992). Geomorphic applications of the Hausdorff- Besicovitch dimension. Earth Surface Processes and Landforms, 17, 775-787. Parrish, R. R. (1983). Cenozoic thermal evolution and tectonics of the Coast Moun- tains of British Columbia: 1. fission track dating, apparent uplift rates and patterns of uplift. Tectonics, 2(6), 601-631. 108 Pelletier, J. D., Malamud, B. D., Blodgett, T., and Turcotte, D. L. (1997). Scale- invariance of soil moisture variability and its implications for the frequency-size distribution of landslides. Engineering Geology, 48, 255-268. Penck, W. (1924). Morphological Analysis of Landforms. MacMillan, London, trans. H. Geek and K . C. Boswell, 1953. Quine, W. v. 0 . (1953). From a Logical Point of View. Harvard University Press, Cambridge, Mass. Rodriguez-Iturbe, I. and Rinaldo, A . (1997). Fractal River Basins. Cambridge Uni- versity Press, Cambridge. Roering, J. J. , Kirchner, J. W., and Dietrich, W. E. (1999). Evidence for non-linear, diffusive sediment transport on hillslopes and implications for landscape morphology. Water Resources Research, 35(3), 853-870. Rood, K . M . (1984). An Aerial Photograph Inventory of the Frequency and Yield of Mass-Wasting on the Queen Charlotte Islands, British Columbia (LandManagement Report No. 34)- British Columbia Ministry of Forests, Victoria, Canada. Rood, K . M . (1990). Site Characteristics and Landsliding in Forested and Clearcut Terrain, Queen Charlotte Islands, British Columbia (Land Management Report No. 64)- British Columbia Ministry of Forests, Victoria, Canada. Ryder, J. M . (1981). Geomorphology of the southern part of the Coast Mountains of British Columbia. Zeitschrift fur Geomorphologie N.F. Suppl.-Bd., 37, 120-147. Savage, S. B. (1993). Mechanics of granular flows. In K . Hutter, editor, Continuum Mechanics in Environmental Sciences and Geophysics, pages 467-522. Springer- Verlag, Wien. 109 Savigny, K . W. (1991). Engineering geology of large landslides in the Lower Fraser River valley transportation corridor southwestern Canadian Cordillera. In S. G. Evans, editor, Landslide Hazards in the Canadian Cordillera. Geological Association of Canada. Schumm, S. A . (1956). Evolution of drainage systems and slopes in badlands at Perth Amboy, New Jersey. Bulletin, Geological Society of America, 67, 597-646. Shenker, O. R. (1994). Fractal geometry is not the geometry of nature. Studies in the History and Philosophy of Science, 25(6), 967-981. Smith, R. B. , Commandeur, P. R., and Ryan, M . W. (1986). Soils, Vegetation and Forest Growth on Landslides and Surrounding Logged and Old-Growth Areas on the Queen Charlotte Islands (Land Management Report No. 41)- British Columbia Ministry of Forests, Victoria, Canada. Snow, R. S. (1989). Fractal sinuosity of stream channels. Pure and Applied Geophysics, 131, 99-109. Sornette, D. and Zhang, Y . - C . (1993). Non-linear Langevin model of geomorphic erosion processes. Geophysical Journal International, 113, 382-386. St0lum, H.-H. (1996). River meandering as a self-organization process. Science, 271, 1710-1713. Sutherland-Brown, A. (1968). Geology of the Queen Charlotte Islands, British Columbia (Bulletin 54)- British Columbia Department of Mines and Petroleum Resources, Victoria. Takahashi, T. (1981). Debris flow. Annual Review of Fluid Mechanics, 13, 57-77. 110 Tarboton, D. G., Bras, R. L. , and Rodriguez-Iturbe, I. (1991). On the extraction of channel networks from digital elevation data. Hydrological Processes, 5, 81-100. Tucker, G; E. and Bras, R. L. (1998). Hillslope processes, drainage density and land- scape morphology. Water Resources Research, 34(10), 2751-2764. Tucker, G. E. and Slingerland, R. L. (1994). Erosional dynamics, flexural isostasy and long-lived escarpments: A numerical modelling study. Journal of Geophysical Research, 99(B6), 12,229 - 12,243. Turcotte, D. L. (1992). Fractals and Chaos in Geology and Geophysics. Cambridge University Press, Cambridge. Turcotte, D. L. (1998). Fractals and Chaos in Geology and Geophysics. Cambridge University Press, Cambridge (Second Edition). Turcotte, D. L. (1999). Self-organized criticality. Reports on Progress in Physics, 62(10), 1377-1429. Veneziano, D. and Iacobellis, V . (1999). Self-similarity and multifractality of topo- graphic surfaces at basin and subbasin scales. Journal of Geophysical Research, 104(B6), 12,797-12,812. I l l Append i x A Computa t ion of Fracta l D imens ion A . l S lope-Area P lo t Construct ions mfslarea.m was written to construct slope-area plots from grid DEMs. It runs in Octave/Matlab, and takes as input the output from ARCINFO's GRIDASCII command (with header removed). ## SLOPEAREA.M Simon Dadson, 16th May 2000 ## Octave script to produce a slope-area plot from ARC/INFO grids of ## slope and cumulative area (named slope.<fname> and area.<fname>) ## user i d e n t i f i e s hillslope-channel scale cutoff and ## channel network f r a c t a l dimension i s computed. ## Based on algorithms i n Montgomery and FoufoulaTGeorgiou (1993) and ## Veneziano and Iacobellis (1999) clear; closeplot; disp("Slope-Area S t a t i s t i c s . Simon Dadson (2000)."); dispC'Reading ARC/INFO data"); fflush(stdout); name = inputC'Enter a name for t h i s analysis: ", " s " ) ; l s l = strcatC'load area.", name); ls2 = strcatC'load slope.", name); 112 e v a l ( l s l ) ; eval(ls2); slope = slope ./ 100; ## because ARC/INFO gives gradient as a percentage area = area .* 100; ## because GVRD DEMs have 100 sq. m. pixels dimensions = size(area); area_max = max(max(area)); area_min = 100; ## i . e . , gri d c e l l area divs =20; ## number of divisions step_len =3; ## step length results = zeros(divs +1,8); results(1:divs,2) = rot90(logspace(logl0(area_min),logl0(area_max),divs),3); results(2:divs + 1, 1) = rot90(logspace(logl0(area_min),logl0(area_max),divs),3) r e s u l t s ( l . l ) = 0; results(:,3) = (results(:,1) + results(:,2)) ./ 2; ## do average for each area class disp("Computing averaged slope-area s t a t i s t i c s " ) ; fflush(stdout); for ypos = 1:step_len:dimensions(2) for xpos = 1:step_len:dimensions(1) tmp_area = area(xpos, ypos); ## to use as an index tmp_slope = slope(xpos, ypos); i f tmp_area != -999900 & tmp_slope != -99.99 ## skip N0DATA values for i = l:divs i f tmp_area >= r e s u l t s ( i , l ) & tmp_area < results(i,2) results(i,4) = results(i,4) + 1; results(i,5) = results(i,5) + tmp_slope; results(i,6) = results(i,6) + tmp_slope"2; results(i,7) = results(i,7) + tmp_slope~3; resu l t s ( i , 8 ) = results(i,8) + tmp_slope"4; endif endfor endif endfor endfor results(:,5) = results(:,5)./results(:,4); results(:,6) = results(:,6)./results(:,4); 113 results(:,7) = results(:,7)./results(:,4); results(:,8) = results(:,8)./results(:,4); ## produce slope area plot sl_min = min(results(2:divs,8)); sl_max = max(max(results(2:divs,5:8))); a x i s ( [ l , area_max, sl_min, sl_max]); hold on; title(strcat("Slope-Area Plot: ", name)); xlabel("drainage area (square metres)"); y l a b e l ( " l o c a l slope m/m"); gset key l e f t bottom t i t l e "Select Cutoff..." box; loglog(results(2:divs,3), results(2:divs,5), "-@x;Ml;"); loglog(results(2:divs,3), results(2:divs,6), "-@+;M2;"); loglog(results(2:divs,3), results(2:divs,7), "-@*;M3;"); loglog(results(2:divs,3), results(2:divs,8), "-<3o;M4;"); disp("(N.B.: a message about NaN indicates that the sample size taken from the DEM i s too small for the number of class divisions chosen. From t h i s point on, empty class-bins w i l l be ignored.");fflush(stdout); cutoff = input("Enter point for f l u v i a l i n f l e x i o n : " ) ; for count = l:divs i f cutoff >= results(count,1) & cutoff < results(count,2) f i r s t _ p o i n t = count; endif endfor ## compute f r a c t a l dimension disp("Computing f r a c t a l dimension of r i v e r network"); fflush(stdout); n = zeros(l,4); sumx = zeros(1,4); sumy = zeros(1,4); sumxy = zeros(1,4); sumxsq = zeros(1,4); for m = 1:4 for i = f i r s t _ p o i n t : d i v s 114 tmpx = logl0(results(i,3)); tmpy = logl0(results(i,4+m)); i f isnan(tmpy) == 0 sumx(m) = sumx(m) + tmpx; sumy(m) = sumy(m) + tmpy; sumxy(m) = sumxy(m) + (tmpx * tmpy); sumxsq(m) = sumxsq(m) + tmpx"2; n(m) = n(m) + 1; endif endfor endfor xmean = zeros(1,4); ymean = zeros(1,4); reg_slope = zeros(1,4); reg_intcpt = zeros(1,4); xmean = sumx ./ n; ymean = sumy ./ n; reg_slope = ( (n .* sumxy) - (sumx .* sumy) ) ./ ( (n . * sumxsq) - (sumx.~2) ); reg_intcpt = ymean - (reg_slope .* xmean); xfit = linspace(cutoff, area_max, 10); loglog(xfit, 10."(reg_intcpt(l) + loglO(xfit) * reg_slope(1)), \ ';Fractal Regression;'); loglog(xfit, 10."(reg_intcpt(2) + loglO(xfit) * reg_slope(2)), \ >..,\. loglog(xfit, 10.~(reg_intcpt(3) + loglO(xfit) * reg_slope(3)), \ )..>•). I I / i loglog(xfit, 10."(reg_intcpt(4) + loglO(xfit) * reg_slope(4)), \ i i J > theta = -1 * reg_slope(l); H_cn = 1 - (2 * theta) D_cn = 2 - H_cn tstr = ['gset key left bottom t i t l e "H_cn = ' num2str(H_cn) ' D_cn = ' num2str(D_cn) ""] ; eval(tstr); hold off; replot; disp("Saving plot to ./slopearea.ps"); fflush(stdout); 115 gset term postscript landscape; outstr = s t r c a t C g s e t output "slopearea' .name, '.ps'"); eval(outstr); replot; gset term x l l ; disp("The next plot w i l l be a straight l i n e i f the r i v e r network i s a se l f - s i m i l a r f r a c t a l ; i f the r i v e r network i s a m u l t i f r a c t a l , a non-linear curve w i l l be seen."); fflush(stdout); dispC'Press a key to diagnose m u l t i f r a c t a l i t y . " ) ; fflush(stdout); foo = kbhitO ; figure(2); axis([1,4,min(reg_slope), max(reg_slope)]); gset x t i c s 1; gset nokey; t i t l e ( s t r c a t ( " S l o p e Area Moments: ", name)); xlabeK"moment order"); ylabel("moment slope"); plot(linspace(l,4,4), reg.slope, "-@+;;"); disp("Saving plot to ./mmtslp.ps"); fflush(stdout); gset term postscript landscape; outstr2 = s t r c a t C g s e t output "mmtslp'.name, '.ps'"); eval(outstr2); replot; gset term x l l ; disp("Done."); A . 2 Const ruc t ion of Hi l ls lope Var iograms Hillslope variograms were constructed using mfs .cc , given below. It produces variograms of the first through fifth order moments of an elevation field defined in a file of x, y, z coordinates in the format output by A R C / I N F O . /*mfs.cc Simon Dadson 2000 */ 116 /•Compute f r a c t a l dimension from */ /•ARC/INFO point coverage using */ /•variogram method */ /***********************************/ /* Algorithm set up parameters from user input or command-line read input f i l e into l i n k e d - l i s t of Point2d (see objgeom) select random subsample by purging points from l i s t determine boundaries for log-spaced bins compute variogram for each point pair compute distance (and angle) f i t pair into a bin and include i t i n the moment calculation complete moment calculation and return plot the variogram with octave ask user to i d e n t i f y range over which to compute D (unless running i n batch mode, i n which case choose the largest range) do linear regression & convert results to f r a c t a l s t a t i s t i c s write output f i l e s */ #include <cmath> #include <cstring> #include <iostream.h> #include <fstream.h> #include <string> #include <strstream.h> #include <list> #include "objgeom" #include "stats" #include "random" using namespace std; // (std) Basic Maths // (std) i/o (suns needs .h on some includes) // (std) F i l e i/o // (std) C++ string class // (std) support for streaming operations on strings // (std) support for linked l i s t s // geometry algorithms and Point2D class // least-squares linear regression algorithms // random number generator from Knuth 1981 typedef list<Point2D> p l i s t ; // Defines Point2D l i s t type typedef list<Point2D>::iterator p l i s t _ i t ; // See "objgeom" include f i l e typedef vector<double> dvec; // define a double-precision vector class string runref; int h_start, h_stop, h_num, h_t o l , anisodirn, adirn_tol = 0; fl o a t sample = -1; 117 void UsageMsgO ; void ReadParameters(); void Summary(plistfe inputdata); void OctavePlot(dvec h, dvec gamma); void OctaveFDPlot(dvec h, dvec gamma, dvec regstats, dvec f_bounds); void WriteLog(const stringfe runref); void WriteOut(const stringfe runref, const dvecfe h, const dvecfe gamma); void WriteFrac(string runref, dvec l h , dvec lgamma, dvec regstats, dvec fr a c s t a t s ) ; p l i s t readinputfile(const string& runref); dvec logspace(int s t a r t , int stop, int nsteps); dvec variogram(const dvec hvalues, plistfc inputdata); dvec askbounds(const bool&); dvec regress(const dvecfe y, const dvecfe x, const dvecfe range); dvec fconvert(dvec regstats, dvec range); int str2int(char* s); p l i s t purge(plistfe, f l o a t ) ; int f i t ( c o n s t doublefe d _ i j , const dvecfe hvalues); dvec loglO (dvec v); int main(int argc, char* argv[]) { bool iac t = true; i f (argc == 1) UsageMsgO ; i f (string(argv[l]) == s t r i n g ( " - i " ) ) ReadParameters(); else i f (string(argv[l]) == ("-b")) { // interpret command l i n e parameters i f (argc != 9) UsageMsgO; iact = fa l s e ; runref = argv[2]; h_start = str2int(argv[3]); h_stop = str2int(argv[4]); h_num = str2int(argv[5]); h_tol = str2int(argv[6]); anisodirn = str2int(argv[7]); adirn_tol = str2int(argv[8]); } else UsageMsgO ; p l i s t b i g l i s t = readinputfile(runref); p l i s t inputdata = purge(biglist, sample); dvec h = logspace(h_start, h_stop, h_num); dvec gamma = variogram(h, inputdata); 118 i f (iact == true) OctavePlot(h, gamma); dvec f_bounds = askbounds(iact); dvec regstats = regress(loglO(gamma), loglO(h), f_bounds); dvec fracstats = fconvert(regstats,f_bounds); i f (iact == true) OctaveFDPlot(h, gamma, regstats, f_bounds); WriteLog(runref); WriteOut(runref, h, gamma); WriteFrac(runref, loglO(h), loglO(gamma), regstats, f r a c s t a t s ) ; > void ReadParametersO { cout « "mfs v.2 SD 21 January 2000" « endl; cout « "Interactive Session" « endl; cout << "Enter runref (input f i l e w i l l be <runref>.in) "; cin >> runref; cout « "Enter sample size (fraction to keep) "; cin » sample; cout « "Enter i n i t i a l value for h (h_start): "; cin » h_start; cout << "Enter f i n a l value for h (h_stop): "; cin >> h_stop; cout « "Enter number of steps (h_num): "; cin » h_num; cout << "Enter h-tolerance (h_tol): "; cin » h_t o l ; cout << "Enter anisotropy dirn (360-degree bearing from N; -1 = is o t r o p i c ) : "; cin » anisodirn; i f (anisodirn != -1) { cout << "Enter angle tolerance (adirn_tol): "; cin » adirn_tol; } } void UsageMsgO { cout « "mfs v2" « endl; cout « "Usage: f d - i (interactive mode)" << endl; cout « " f d -b runref h_start h_stop h_num h_tol anisodirn adirn_tol (batch mode)" « endl; exit(0); 119 void OctavePlot(dvec h, dvec gamma){ // Generate graphs using Octave string dotm = runref + ".m"; ofstream octfile(dotm.c_str()); o c t f i l e « "load 'mmt." « runref « "';" « endl; o c t f i l e « "h = mmt(:,l);" « endl; o c t f i l e << "gamma = mmt(:,2);" « endl; o c t f i l e « "second = mmt(:,3);" « endl; o c t f i l e « "t h i r d = mmt(:J4);" « endl; o c t f i l e « "fourth = mmt(:,5);" « endl; o c t f i l e « " f i f t h = mmt(:,6);" « endl; o c t f i l e « "title('Variogram: " « runref « ' " ) ; " « endl; o c t f i l e « "xlabelC separation (h)');" « endl; o c t f i l e « "ylabel('moment');" « endl; o c t f i l e « "gset key below;" « endl; o c t f i l e « "loglog(h,gamma,'@33;Ml;')" « endl; o c t f i l e « "hold on;" « endl; o c t f i l e « "loglog(h,second,'@33;M2;')" « endl; o c t f i l e « "loglog(h,third,'@33;M3;')" « endl; o c t f i l e « "loglog(h,fourth,'033;M4;')" « endl; o c t f i l e « "loglogQi,fifth,'@33;M5;')" « endl; o c t f i l e << "pause" « endl; string callstr="octave -q " + dotm; system(callstr.c_str()); // c a l l s octave to create graphs void OctaveFDPlot(dvec h, dvec gamma, dvec regstats, dvec f_bounds){ // plot f r a c t a l regression on variogram using Octave string f r a c l i n e = runref + "fin.out"; ofstream f l n ( f r a c l i n e . c _ s t r ( ) ) ; for (int count = (int) f.bounds[0]; count <= (int) f_bounds[l]; count++) { fl o a t x, y; x = h[count]; y = pow(10,(regstats[0]*loglO(x))+regstats[l]); f i n « x « " " « y « endl; } string dotm = runref + "f.m"; ofstream octfile(dotm.c_str()); o c t f i l e « "load 'mmt." « runref « "';" « endl; 120 'load "' « runref « "fln.out';" « endl; 'h = mmt(:,1);" « endl; 'gamma = mmt(:,2);" « endl; 'second = mmt(:,3);" « endl; 'third = mmt(:,4);" « endl; 'fourth = mmt(:,5);" « endl; " f i f t h = mmt(:,6);" « endl; "x =" « runref « " f l n ( : , l ) ; " « endl; "y =" « runref « " f l n ( : , 2 ) ; " « endl; 'title('Variogram: " « runref « ' " ) ; " « endl; "xlabelC separation (h)');" « endl; 'ylabel('moment');" « endl; "gset key below;" << endl; "hold on" « endl; 'loglog(h,gamma,'033;Ml;')" « endl; 'loglog(x,y,';Fractal Regression;')" << endl; 'loglog(h.second,'033;M2;')" « endl; 'loglog(h,third,'033;M3;')" « endl; 'loglog(h,fourth,'033;M4;')" « endl; 'loglog(h,fifth,'033;M5;')" « endl; "hold o f f " « endl; 'pause" « endl; // o c t f i l e << "gset term postscript landscape" « endl; // o c t f i l e « "gset output '" « runref « ".ps"' « endl; // o c t f i l e « "replot" « endl; // o c t f i l e « "gset term x l l " « endl; string callstr="octave -q " + dotm; system(callstr.c_str()); // c a l l s octave to create graphs o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e << o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « o c t f i l e « void WriteLog(const stringfe runref){ // write parameters to log f i l e s t r i n g dotlog = runref + ".log"; ofstream l o g f i l e ( d o t l o g . c _ s t r ( ) ) ; l o g f i l e « "fd.cc produced t h i s data using the following parameters:" « endl; l o g f i l e << "h_start, h_stop, h_num, h_t o l , anisodirn, a d i r n _ t o l , sample" « endl; l o g f i l e « h.start « "," « h_stop « "," « h_num « "," « h_tol « "," « anisodirn « "," « adirn.tol « "," « sample « "," « endl; 121 void WriteOut(const stringfc runref, const dvec& h, const dvecfe gamma){ // write variogram to f i l e s t r i n g dotout = runref + ".out"; ofstream out f i l e ( d o t o u t . c _ s t r ( ) ) ; for (int count = 0 ; count <= (h_num-l); count++) { o u t f i l e << h[count] « " " « gamma[count] « endl; } void WriteFrac(string runref, dvec l h , dvec lgamma, dvec regstats, dvec fracstats){ // write f r a c t a l stats to f i l e s t r i n g dotfrc = runref + " . f r c " ; ofstream f r c f i l e ( d o t f r c . c _ s t r ( ) ) ; for (int count = 0 ; count <= (h_num-l); count++) { f r c f i l e « lh[count] « " " << lgamma[count] « endl; } f r c f i l e << "Regression S t a t i s t i c s : slope=" « regstats [ 0 ] « " intercept=" « regstats [ 1 ] « endl; f r c f i l e « "Fractal S t a t i s t i c s : D=" « fracstats [ 0 ] « " Gamma(0)=" « fracstats [ 1 ] « " Bounds " « fracstats [2] « " <= h <=" « fracstats [ 3 ] « endl; cout « "Regression S t a t i s t i c s : slope=" << regstats [ 0 ] « " intercept=" << regstats [ 1 ] « endl; cout « "Fractal S t a t i s t i c s : D=" « fracstats [ 0 ] « " Gamma(0)=" « fracstats [ 1 ] « " Bounds " « fracstats [ 2 ] « " <= h <=" « fracstats [ 3 ] << endl; p l i s t readinputfile(const stringfe runrefM /*** open data f i l e for input ***/ string dotin = runref + " . i n " ; ifstream i n f i l e ( d o t i n . c _ s t r ( ) ) ; i f ( U n f i l e ) { cerr « "Cannot open: " « dotin « " for input. Aborting." << endl; e x i t ( l ) ; } else cout « "Reading Datafile: "; 122 /*** create l i s t of points and read data f i l e into t h i s l i s t ***/ str i n g dataline; p l i s t r e a d l i s t ; p l i s t _ i t read_it = readlist.begin(); // starts traverse of r e a d l i s t on // read_.it int size =0; // number of lin e s read while ( g e t l i n e ( i n f i l e , dataline)) // read one dataline { size++; i f (size == 1) /*do nothing (to eat header line) */; else { f l o a t x, y, z, rubbish; //temp vars to decode stream const char* buffer = dataline.c_str(); //dataline -> c-string for //streaming istrstream iss(buffer,80); i s s » rubbish >> z >> x >> y; // rubbish catches ARC rec num Point2D temp_point(x,y,z); r e a d l i s t . i n s e r t ( r e a d _ i t , temp_point); // cout << temp_point « endl; // useful for testing... } } /*** Perform niceties and return linked l i s t of points ***/ cout « size-1 « " data lines read." « endl; return r e a d l i s t ; > dvec logspace(int start, int stop, int nsteps){ // returns a dvec of nsteps values, logarithmically spaced between // start and stop (inclusive) dvec logbins; cout « "Computing logarithmically spaced h-values: "; double l s t a r t = (double) loglO(start); double lstop = (double) loglO(stop); double incr = (lstop - l s t a r t ) / ((double)nsteps - 1); for (double lgh = l s t a r t ; lgh < lstop; lgh += incr) { double temp = pow(10,lgh); logbins.push_back(temp); } logbins.push_back((double)stop); //cheat for la s t value to cover 123 / / any roundoff e r ro r cout « "Done." « end l ; r e tu rn l o g b i n s ; } dvec askbounds(const bool& i a c t ) { / / ask user fo bounds of f r a c t a l regress ion / / choose e n t i r e range i f running i n command l i n e mode dvec bounds; i f ( i a c t == t rue) { double f i r s t , l a s t = 0; cout « "Enter f i r s t point fo r f r a c t a l est imate: " ; c i n » f i r s t ; cout « "Enter l a s t point fo r f r a c t a l est imate: " ; c i n » l a s t ; bounds .push_back(f i rs t ) ; bounds.push_back(last) ; } else { bounds.push_back(0); bounds.push_back(h_num-l); } r e tu rn bounds; > dvec regress(const dvec& y , const dvecfe x , const dvec& range){ / / regress y on x over range, r e tu rn s lope , in te rcep t i n a dvec dvec x f ; dvec y f ; copy(xf , x , ( i n t ) range[0] , ( i n t ) r a n g e [ l ] ) ; copy(yf, y , ( i n t ) range[0] , ( i n t ) r a n g e [ l ] ) ; dvec x f sq = squares(xf) ; dvec y f sq = squares(yf) ; dvec xyfprod = products (xf , y f ) ; i n t nx = xf . s i z e O ; i n t ny = y f . s i z e O ; i f (nx!=ny) e x i t ( 0 ) ; double sumxf = sum(xf); / / w e ' l l use t h i s a l o t double slope = (ny*sum(xyfprod)-sumxf*sum(yf))/ (nx*sum(xfsq)-sumxf*sumxf); double in te rcep t = mean(yf) - (slope * mean(xf)); dvec output; 124 output.push_back(slope); output.push_back(intercept); return output; dvec fconvert(dvec regstats, dvec range){ // convert regression results into f r a c t a l s t a t i s t i c s dvec fracstats; fracstats.push_back(3 - regstats[0]); fracstats.push_back(regstats[1] ); fracstats.push_back(range[0]); fracstats.push_back(range[1]); return fracstats; } int str2int(char* s) { // return an integer representation of a numeric string int sign = 1; int value = 0; int place = 1; int length = s t r l e n ( s ) ; i f (length > 10) { cout « "str2int f a i l e d : number too large" << endl; e x i t ( l ) ; } for (int counter = 0; counter < length; counter++) { i f (s[counter]=='-') { sign = -1; place++; } else { int number = int(s[counter])-48; i f ((number >= 0) && (number <= 9)) { int power = (int)pow(10,(length - place)); value += number * power; place++; } else { cout « "str2int: non-integer argument." « endl; e x i t ( l ) ; } 125 } value *= sign; return value; } p l i s t purge(plistfe input, f l o a t proportion) { // create subsample of points based on user specified proportion // i f running i n batch mode, always choose 2000 points ( i f available) p l i s t purged = input; int l i s t s i z e = input.size(); i f (proportion == -1) proportion = ((float)2000 / ( f l o a t ) l i s t s i z e ) ; i f (proportion > 1) proportion = 1; cout << "Original d a t a f i l e has: " << l i s t s i z e << " points." « endl; cout « "Random sampling..." « endl; Random choose; for (int count = 1; count <= listsize*(1-proportion); count++) { int purgelocn = choose.integer(purged.size()-1,0); p l i s t _ i t d e l _ i t = purged.begin(); for (int count2 = 1; count2 <= purgelocn; count2++) { ++del_it; } purged.erase(del_it); } cout << "Subsample of: " « purged.sizeO « " points taken." << endl return purged; dvec variogram(const dvec hvalues, plistfe inputdata){ // return variogram of inputdata over bins hvalues. // see thesis for d e t a i l s of this algorithm // (write out a l l moments for p l o t t i n g (added 25/05/00)) int hsize = hvalues.size(); dvec N; dvec f i r s t ; dvec second; dvec t h i r d ; dvec fourth; dvec f i f t h ; dvec gamma; 126 / / f i l l N, sum and gamma with zeros for ( int count = 0; count < hsize; count++) { N.push_back(0); f irst .push_back(0); second.push_back(0); third.push_back(0); fourth.push_back(0); f i f th.push_back(0); gamma.push_back(0); } / / attempt to f i t each point into an h-bin for ( p l i s t _ i t t l l= inputdata .begin() ; t i l != inputdata.end(); ++tll) { p l i s t _ i t temp = inputdata .beginO; / / allow advance by one next l ine temp = t i l ; temp ++; for ( p l i s t _ i t tl2=temp; t l 2 != inputdata.end(); ++tl2) { / / Finds distance and angle between point p a i r s . double d_ij = d i s t ( * t l l , * t l 2 ) ; double theta_i j = a n g l e ( * t l l , * t l 2 ) ; double z d i f f _ i j = z d i f f ( * t l l , * t l 2 ) ; i f (anisodirn != -1) { cout « "I have not yet learned to make anisotropic variograms!" « endl; e x i t ( l ) ; > else { in t l oca t ion = f i t ( d _ i j , h v a l u e s ) ; i f ( locat ion != 9999) { f i r s t [ l o c a t i o n ] += z d i f f _ i j ; second[location] += pow(zdi f f_ i j , 2); th ird[ locat ion] += pow(zdi f f_ i j , 3 ) ; fourth[ locat ion] += pow(zdi f f_ i j , 4 ) ; f i f th [ l oca t ion ] += pow(zdi f f_ i j , 5 ) ; N[location]++; } } } } for ( int count = 0; count < hsize; count++) { i f (N[count] > 100) { 127 gamma[count] = f i r s t [ c o u n t ] / ( N [ c o u n t ] ) ; second[count] /= N[count] * 2 ; th i rd [coun t ] /= N[count] * 3 ; fourth[count] /= N[count] * 4 ; f i f t h [ c o u n t ] /= N[count] * 5 ; } else { cout « "warning: only " « N[count] « " input po in t s i n h_bin « count « " . " « end l ; } } s t r i n g o u t f i l e = "mmt." + runref ; ofstream o u t ( o u t f i l e . c _ s t r ( ) ) ; fo r ( i n t count = 0; count < h s i z e ; count++) { out « hvalues [count] << " " « gamma[count] « " " « second[count] « " " « th i rd[count ] « " " « fourth[count] « " " « f i f t h [ coun t ] « end l ; } r e tu rn gamma; > i n t f i t ( c o n s t doublefe d _ i j , const dvecft hvalues) { / / r e tu rn the array index i n hvalues i n t o which a po in t w i th / / separa t ion dis tance d_ i j should be p laced , i n t l o c a t i o n ; bool housed = f a l s e ; i n t count = 0; whi le (housed == f a l s e ) { double h = hvalues[count] ; i f ( (d_ i j <= (h + h_ to l ) ) && (d_i j >= (h - h _ t o l ) ) ) { l o c a t i o n = count; housed = t rue ; } i f (count == h v a l u e s . s i z e ( ) ) { l o c a t i o n = 9999 ; housed = t rue ; } count++; } r e tu rn l o c a t i o n ; 128 } dvec loglO (dvec v) { / / return loglO of a dvec ( in one go) / / same as octave/matlab dvec logs; int v_size = v . s i z e O ; for (int i = 0; i < v_size; i++) { logs.push_back(loglO(v[i])); / / take loglO of each element > return logs; } / l i e * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * / /* class defn for Point2D * / /* generic functions on p l i s t s */ class Point2D { friend double dist(const Point2D&, const Point2D&); friend double angle(const Point2D&, const Point2D&); friend double zsq(const Point2D&, const Point2D&); friend double zdiff(const Point2D&, const Point2D&); friend istreamfe opera tor»( i s t ream& i s t r , Point2D& p) ; friend ostreamfe operator«(ostream& ostr , const Point2D& p); public: Point2D(double x, double y, double z) { _xcoord = x; _ycoord = y _zvalue = z; } Point2D() { _xcoord = -9999; _ycoord = -9999; _zvalue = -9999; } double xposO { return _xcoord; > double yposO { return _ycoord; } double zval() { return _zvalue; } void p r i n t O { cout « "(" « _xcoord « " , " « _ycoord « ") = « _zvalue;} private: double _xcoord; double _ycoord; double _zvalue; >; / / streaming operators (to write to screen and read from keyboard) 129 ostreamfe operator<<(ostream& o s t r , const Point2D& p) { i f ( (p ._xcoord == -9999) I I (p._ycoord — -9999)) o s t r << " P o s i t i o n Undefined"; e l se i f (p._zvalue == -9999) os t r « " (" « p ._xcoord « " , " << p._ycoord « ") = no va lue" ; e lse o s t r « " ( " « p._xcoord « " , " « p ._ycoord « ") = " « p ._zva lue ; r e tu rn o s t r ; } istreamfe o p e r a t o r » ( i s t r e a m & i s t r , Point2D& p) { / / format fo r input i s x , y , z char b i n ; / / used to eat ( and , i s t r » p ._xcoord; i f ( i s t r . p e e k ( ) = = ' , ' ) i s t r >> b i n » p ._ycoord; i f ( i s t r . p e e k ( ) = = ' , ' ) i s t r >> b i n » p ._zva lue ; r e tu rn i s t r ; } double d i s t ( c o n s t Point2D& p i , const Point2D& p2){ / / Finds s t r a i g h t l i n e dis tance from p i to p2 i f ( ( p l . _xcoord == -9999) I I (p l ._ycoord == -9999) I I (p2._xcoord == -9999) | | (p2._ycoord == -9999)) r e tu rn -9999; e lse { double xdtemp = p2._xcoord - p l . _ x c o o r d ; double ydtemp = p2._ycoord - p l . _ y c o o r d ; xdtemp *= xdtemp; ydtemp *= ydtemp; double dtemp = sqrt(xdtemp+ydtemp); r e tu rn dtemp; } } double angle(const Point2D& p i , const Point2D& p2){ i f ( (p l . _xcoord == -9999) I I (p l ._ycoord == -9999) I I (p2._xcoord == -9999) | | (p2._ycoord == -9999)) r e tu rn -9999; e l se { / / Finds 360 degree bear ing from North of pt 2 from pt 1 double xdtemp = p2 ._xcoord-p l ._xcoord ; double ydtemp = p2 ._ycoord-p l ._ycoord; double theta=0; 130 / / Correct for quadrant to get radian bearing from north i f (xdtemp > 0) theta = 1.57079632679 - atan(ydtemp/xdtemp); / / = p i / 2 else i f (xdtemp < 0) theta = 4.71238898038 - atan(ydtemp/xdtemp); / / = 3 p i by 2 / / Convert to degrees theta *= 57.2957795131; / / 180 / p i (rad -> deg) return theta; } } double zsq(const Point2D& p i , const Point2D& p2){ i f ( (pl ._xcoord == -9999) I I (pl ._ycoord == -9999) I I (p2._xcoord == -9999) I I (p2._ycoord == -9999)) return -9999; else { / / f i n d square of distance between two points double z d i f f = (pi ._zvalue)-(p2._zvalue); double zsq = z d i f f * z d i f f ; return zsq; } } double zdi f f (const Point2D& p i , const Point2D& p2){ / / return absolute value of z difference i f ( (pl ._xcoord == -9999) I I (pi ._ycoord == -9999) I I (p2._xcoord == -9999) I I (p2._ycoord == -9999)) return -9999; else { double z d i f f = (pi ._zvalue)-(p2._zvalue); i f ( zd i f f < 0) z d i f f *= -1; return z d i f f ; } } 131 Appendix B Derivation of Landscape Modelling Equations In this appendix, the rules used to represent hillslope and fluvial landscape processes are derived from assumptions about their physical mechanisms. B . l Continuity Equation Assuming that Z , the elevation at a point in the landscape, comprises the sum of landscape elements, z: Z The time variation of Z is: / zd?r (B.l) Jv dZ_ _ d_ dt dt JV j zd3r (B.2) v = / v { l * + A - ( B ' 3 ) Where Vj is a velocity vector of the boundary of the landscape element. Suppose that: § = f * * r . - f p - # r (B.4) at Jv Jv OXJ i.e., that changes in elevation are due to (i) "introduction" or "removal" of landscape mate- rial, by tectonic or fluvial processes, and (ii) flux of sediment, qj, away from the point. 132 Combining B3 and B4, and applying conservation of landscape volume, the local form of the balance equation becomes: dt dxj f dxj J ^ ( 5 ) B.2 Fluvial Transport Rule The rationalized Bagnold equation written by Martin and Church (2000) is: ib = [u- uof^/V^l/^y V/4)] (B.6) where if, = bedload (dry) mass flux per unit width, u — stream power per unit width, u>o = critical stream power, D = characteristic grain size diameter, d = channel depth, pr = relative density of sediment in water, and g = acceleration due to gravity. Assuming UJ = (pwgQS)/w, UJQ ~ 0, and D « 1, then this expression simplifies to: I l = 0 . 0 1 3 9i<^H^ (B.7) a which is multiplied by channel width, w, to obtain the mass flux of sediment from the model cell. This is then divided by the density of sediment, ps, and the area of the channel bed (length of the channel represented in the cell ( « As, the cell dimension) multiplied by channel width) to obtain the change in river-bed elevation for the model cell, Azb, in metres per second: 0.0139 w dps Ax Sediment transporting discharge, Q, is assumed to flow for one day each year (=86,400 s), and may be approximated from upslope area using the relation obtained by Church (1997): Q = kfiA0-68, where kji is a parameter to represent hydroclimate. The terms w and d 133 are correlates of discharge, so may be approximated by the hydraulic geometry relations: w = 3.26 Q05, and d = 0.42 Q0A (Kellerhals, in Knighton, 1998) for slopes less than 0.01: Azb (per year) = 6.44 Q 0 3 5 S 1 ' 5 (B.9) Note that channel change has been "smeared" over the grid cell in order to account for the fact that the river is a linear feature in an areal model unit. For slopes greater than 0.01, the hydraulic geometry relations of Day (1969) were used: w = 0.95 Q 0 7 6 , and d « 0.92 Q05. A correction to the slope exponent is also made, following Mart in (1998), to account for configurational stability in steep channels, giving: Azb (per year) = 40.32 Q 0 3 4 5 2 2 2 (B.10) 134 B.3 erode2d.cc / * Landscape Evolut ion Model * / / * Simon Dadson 25.02.00 * / / * PROGRAM HEADERS * / #include <iostream.h> #include <cmath> #include <fstream.h> #include <strstream.h> #include <string> #include <vector> #include "random" #include "landnode" #define SIDELENGTH 30 #define TSTOP 100000 #define NSIDE 50 #define NT 10000 #define HMAX 1500 using namespace std; / / (std) support for streaming operations on str ings / / length of dem g r i d s ize / / length of model run (years) / / nodes on each side of square model domain / / number of timesteps / / maximum height (used i n generating / / a r t i f i c i a l topography) / / Support for custom data types typedef double Matrix[NSIDE][NSIDE]; typedef vector<double> dvec; ofstream s l i d e f i l e C ' l s c a p e . s l i " ) ; ofstream s l idedataC' l scape . s ld") ; ofstream sedtranspC'lscape.sed"); Random b i l l y ; Matrix f i l lnumber; / / keeps record of sinks i n f low-routing algorithm. Matrix chaimelch; / / keeps record of net aggradation/degradation i n channels / * * * * * * * * * * * * * * * * / / * DECLARATIONS * / 135 / / Global Constants const int flowback[8] = {16,32,64,128,1,2,4,8}; / / look-up arrays const int xlookup[8] = {1 ,1 ,1 ,0 , -1 , -1 , -1 ,0}; / / for flow- const int ylookup[8] = {-1,0,1,1,1,0,-1,-1}; / / rout ing const double rt2 = pow(2, 0.5); const double P i = 3.1415927; / / Global Parameters double K, hxO, hscale; double r ; double cr i t_prob; double alpha = 35*(Pi/180); double k _ f l u v i a l ; double hftrans; double Umax; int s l ide_tes t_ int = 1; s t r ing i n p u t f i l e ; s t r i n g upl i f tmodel; double rho_rock = 2.65; double rho_sed = 1.60; bool diagnose = fa l se ; / / Functions void t i t l e ( ) ; void getparamsO; void expnonlinear(Matrix i n i t i a l , Matrix f i n a l , Matrix production, Matrix u p l i f t , f loat s t a r t , f loat stop, f loat dt , s t r i n g kmodel); void renewuplift(Matrix u p l i f t , f loat t , s t r ing upl i f tmodel , double Umax); void dnet(Matrix landscape, Matrix slope, Matrix r i v e r n e t , Matrix flow, double k _ f l u v i a l , f loat t ) ; void f i l l s i n k s ( M a t r i x landscape, Matrix flow); void flowroute(Matrix i n a l t g r i d , Matrix slope, Matrix o u t d i r g r i d ) ; void accumulate(Matrix i n d i r g r i d , Matrix outareagrid); void Bagnold(Matrix acc, Matrix slope, Matrix transport , double k _ f l u v i a l ) ; void fluxsep(Matrix riverchange, Matrix r i v e r f l u x , Matrix flow, f loat t ) ; void dsbl(Matrix landscape, Matrix s l idescape, Matrix s l iderecord , Matrix slope, Matrix flow, f loat t ) ; void export(Matrix output, bool d i sp lay , int i d ) ; void matrix_copy(Matrix o l d , Matrix newmatrix); void matrix_add(Matrix f i r s t , Matrix second); void threshold(Matrix i n , Matrix out, double cutoff ) ; void f i l K M a t r i x input, s t r ing f i l e ) ; / / non-l inear d i f fus ion parameters / / s t a b i l i t y parameter for so lut ion / / c r i t i c a l p r o b a b i l i t y for dsbl s l ides / / repose angle for talus (radian) / / proport iona l i ty const, for f l u v i a l transport / / upslope area where channel i n i t i a t i o n occurs / / maximum u p l i f t rate / / minimum i n t e r v a l between dsbl s l ides / / source of i n i t i a l surface / / u p l i f t regime / / density of rock / / density of sediment / / turns on verbose output for diagnostic work 136 v o i d w r i t e l o g O ; bool i s_a_zero(Mat r ix g r i d ) ; bool i s_even( in t t e s t i n t ) ; bool is_mult(double a, double b ) ; i n t d r a in (Ma t r ix i n g r i d , i n t x i , i n t y i ) ; i n t dsxpos( int o l d x , double f l o w d i r ) ; i n t dsypos( int o ldy , double f l o w d i r ) ; double volume(Matrix landscape); /****************/ / * MAIN PROGRAM * / i n t main( int argc, char* argv[]) { t i t l e O ; getparams () ; Mat r ix landscape, r iverchange, s lope , f low, s l i d e r e c o r d ; f i l l ( l a n d s c a p e , i n p u t f i l e ) ; f i l K s l i d e r e c o r d , "nodata"); f i l l ( f i l l n u m b e r , "zeros") ; f i l l ( c h a n n e l c h , "ze ros" ) ; export( landscape, diagnose, 0 ) ; f l o a t dt = TSTOP/NT; r = (K * dt) / SIDELENGTH; cout « "r = " « r « end l ; / / cout « " S t a b i l i t y parameter = " « (D * dt) / (NSIDE*NSIDE) « endl i n t i d = 1; w r i t e l o g O ; ofstream v f i l e ( " l s c a p e . v o l " ) ; fo r ( f l o a t t = 0; t <= TSTOP-dt; t+=dt) { i f (t==0) dnet(landscape, s lope , r iverchange, f low, k _ f l u v i a l , t ) ; i f ( i s _ m u l t ( t , s l i d e _ t e s t _ i n t ) ) { Mat r ix s l idescape ; dnet(landscape, s lope , r iverchange, f low, k _ f l u v i a l , -1); / / before s l i d i n g (to refresh) dsbl ( landscape, s l idescape , s l i d e r e c o r d , s lope , f l ow, t ) ; mat r ix_copy(s l idescape , landscape); dnet(landscape, s lope , r iverchange, f low, k _ f l u v i a l , t ) ; / / a f te r s l i d i n g } 137 Matrix output; v f i l e << t << " " « volume(landscape) « endl; Matrix u p l i f t ; renewupl i f t (up l i f t , t , upl i f tmodel , Umax); expnonlinear(landscape, output, riverchange, u p l i f t , t , t+dt, 1, "constant"); i f ( is_mult(t , TSTOP/40)) { cout « "time = " « t+(TSTOP/40) « " yrs . " « "Writing lscape" << i d « endl; export(output, diagnose, i d ) ; export(s l iderecord, diagnose, id+50); id++; } matrix_copy(output, landscape); / / feed th i s i t e r a t i o n ' s output back as input for next > / / c a l l Octave scr ip t to display animated 3D projec t ion of landscape export(s l iderecord, diagnose, 99); systemC 'mv lscape99.out l scape.s lh"); export(f i l lnumber, diagnose, 98); systemC'mv lscape98.out lscape.snk"); export(channelch, diagnose, 97); systemC 'mv lscape97.out l scape .cnl"); system("./vscape 2> /dev /nu l l" ) ; } / / TECTONICS void renewuplift(Matrix u p l i f t , f loat t , s t r ing upl i f tmodel , double Umax) { i f (upliftmodel == "none") { / / NO UPLIFT for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { u p l i f t [xpos] [ypos] = 0; > } } i f (upliftmodel == "const") { / / UPLIFT CONSTANT IN SPACE AND TIME for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { uplift[xpos][ypos] = Umax; } 138 } } i f (upliftmodel == "thrust") { / / THRUST for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { i f (xpos < (NSIDE / 2)) upl i f t [xpos] [ypos] = Umax; else uplift[xpos][ypos] = 0 ; > } } i f (upliftmodel == "pluton") { / / PLUTON double step_diff = (Umax / NSIDE) * 2; for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { i f (xpos < (NSIDE / 2)) u p l i f t [xpos][ypos] = xpos * s tep_dif f ; else upl i f t [xpos] [ypos] = (2 * Umax) - (xpos * s tep_di f f ) ; } } } i f (upliftmodel == "rebound") { / / ISOSTATIC REBOUND WITH (LINEAR) DECREASING RATE for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { upl i f t [xpos] [ypos] = Umax * ((TST0P - t) / TSTOP); } } } / / FLUVIAL PROCESSES void dnet(Matrix landscape, Matrix slope, Matrix riverchange, Matrix flow, double k _ f l u v i a l , f loat t) { * II Compute the drainage network from elevations using D8 algorithm / / each c e l l i n Q i s the annual average discharge at that po int , / / based on the area, integrated upslope, above the po int . / / Each c e l l i n r ivernet i s the amount removed by the r i v e r ( i . e . , / / negative production of sediment) 139 Matrix acc, r i v e r f l u x ; flowroute(landscape, slope, flow); bool tampered = fa l se ; double before_tamper = volume(landscape); while (is_a_zero(flow)) { tampered = true; f i l l s i n k s ( l a n d s c a p e , flow); / / f i l l s small sinks only flowroute(landscape, slope, flow); } i f (tampered) { / / pr in t a warning message i f any closed p i t s have been f i l l e d / / as a resu l t s of enforcing the f low-routing rules / / never happens with usual erosion rules because no / / process i s capable of producing a p i t . double after_tamper = volume(landscape); i f ((after.tamper - before_tamper) > (after_tamper * 0.01)) { cout « "lscape (dnet): WARNING! e levat ion f i e l d modified by > 17. to route flow" « endl; } } accumulate(flow, acc); Bagnold(acc, slope, r i v e r f l u x , k _ f l u v i a l ) ; f luxsep(riverchange, r i v e r f l u x , flow, t ) ; i f (t != - 1 ) matrix_add(channelch, riverchange); i f (t = TSTOP) { Matrix network; f i l K n e t w o r k , "nodata"); threshold(acc, network, k_fluvial / (NSIDE * NSIDE)); export(network, diagnose, 96); system("mv lscape96.out lscape.net"); } } void f i l l s i n k s ( M a t r i x landscape, Matrix flow) { . / * for every c e l l , determine i t s lowest neighbour i f the lowest neighbour i s higher than the o r i g i n a l c e l l then increase the height to equal the height of the lowest neighbour 140 (cannot happen with usual erosion rules) */ int sinkcount = 0; bool sinks = true; while (sinks) { sinks = fa l se ; for ( int ypos = 2; ypos < NSIDE-2; ypos++) { / / skip the edges / / where there i s always a sink for ( int xpos = 2; xpos < NSIDE-2; xpos++) { double low_height = 1000000; / / a very large number bool i s_s ink = true; for ( int neighbour = 0; neighbour < 8; neighbour++) { int xj = xpos + xlookup[neighbour]; int yj = ypos + ylookup[neighbour]; i f (landscape[xpos] [ypos] < landscape[xj][yj]) { i f ( landscape[xj][yj] < low_height) low_height = landscape[xj][yj]; } else i s_s ink = fa l se ; } i f ( is_sink) { landscape[xpos] [ypos] = low_height; sinks = true; > } } } / / For a l l c e l l s with undefined flow d i r e c t i o n / / f i n d the global pour point double global_min = 10000; / / a r b i t r a r i l y large value for ( int ypos = 2; ypos < NSIDE - 2; ypos++) { for ( int xpos = 2; xpos < NSIDE - 2; xpos++) { i f (flow[xpos][ypos] == 0) { bool has_higher = fa l se ; for ( int neighbour = 0; neighbour < 8; neighbour++) { int xj = xpos + xlookup[neighbour]; int yj = ypos + ylookup[neighbour]; i f ((landscape[xpos][ypos] != landscape[xj][yj]) kk ( landscape[xj][yj] < global_min)) 141 global_min = landscape[xj][yj]; > } } } / / assign global_min to a l l c e l l s with undefined flow d i r e c t i o n for ( int ypos = 2; ypos < NSIDE - 2; ypos++) { for ( int xpos = 2; xpos < NSIDE - 2; xpos++) { i f (flow[xpos][ypos] == 0) { landscape[xpos] [ypos] = global_min; fillnumber[xpos][ypos]++; } } } } void flowroute(Matrix i n a l t g r i d , Matrix slope, Matrix outdirgr id) { / / kernel of D8 flow routing algorithm. / / Input must be a depressionless DEM; i f i t i s not, depressions / / are returned with zero drainage d i r e c t i o n / / Make edges flow outward; assumes that we are well ins ide area of interest / / a two-ce l l s t r i p i s used because the image boundary conditions create / / spurious features at the edge, which are of no interest for ( int xpos = 0; xpos < NSIDE; xpos++) { outdirgrid[xpos] [0] = 128; outdirgrid[xpos] [1] = 128; outdirgrid[xpos][NSIDE-1] = 8; outdirgrid[xpos][NSIDE-2] = 8; > for ( int ypos = 0; ypos < NSIDE; ypos++) { outdirgrid[0][ypos] = 32; outdirgrid[1][ypos] = 32; outdirgrid[NSIDE-1][ypos] = 2; outdirgrid[NSIDE-2][ypos] = 2; } / / Correct corners outdirgrid[0][0] = 64; outdirgrid[1][1] = 64; 142 outdirgrid[0][NSIDE-1] =16; outdirgrid[1][NSIDE-1] = 16; outdirgrid[NSIDE-1][0] = 1; outdirgrid[NSIDE-1] [1] = 1; outdirgrid[NSIDE-1] [NSIDE-1] = 4; outdirgrid[NSIDE-2] [NSIDE-2] = 4; / / for each i n t e r i o r node compare i t s e levat ion with those of i t s neighbours double diag_div = pow(2 * pow(SIDELENGTH, 2), 0.5); / / hor izonta l / / length for points on diagonal double diag_corr = SIDELENGTH / d iag .d iv ; / / correct ion for diagonals for ( int ypos = 2; ypos < NSIDE-2; ypos++) { for ( int xpos = 2; xpos < NSIDE-2; xpos++) { double ce l lhe ight = inaltgrid[xpos][ypos]; double local_max = 0; int local_flow = 0; for ( int neighbour = 0; neighbour < 8; neighbour++) { int xj = xpos + xlookup[neighbour]; in t yj = ypos + ylookup[neighbour]; double n_slope = (cel lheight - i n a l t g r i d [ x j ] [ y j ] ) / SIDELENGTH; i f (is_even(neighbour)) { / / i f on a diagonal , divide by sqrt(2) n_slope = n_slope * diag_corr ; / / correct for diagonals } i f (n_slope > local_max) { local_flow = ( int) pow(2, neighbour); local_max = n_slope; } else i f (n_slope == local_max) { local_flow += ( int) pow(2, neighbour); } } i f (local_max == 0) local_flow = 0; slope[xpos][ypos] = local_max; outdirgrid[xpos] [ypos] = (double) local_flow; } > / / END OF FIRST PASS / / Choose d i r e c t i o n randomly for c e l l s that could have several for ( int ypos = 2; ypos < NSIDE-2; ypos++) { for ( int xpos = 2; xpos < NSIDE-2; xpos++) { 143 in t local_flow = ( int) outdirgrid[xpos][ypos]; i f ( local_flow == 0) ; / / do nothing: wait for next sweep / / assign central c e l l i f a whole side i s equally lower else i f ( local_flow == 7) local_flow = 2; else i f ( local_flow == 28) local_flow = 8; else i f ( local_flow == 112) local_flow = 32; else i f ( local_flow == 193) local_flow = 128; else { / / unpack flow direct ions and make a random choice i f / / several are equally v a l i d dvec choices; int pips = local_flow; i f (pips >= 128) { choices.push_back(128); pips -= 128; } i f (pips >= 64) { choices.push_back(64); pips -= 64; > i f (pips >= 32) { choices.push_back(32); pips -= 32; } i f (pips >= 16) { choices.push_back(16); pips -= 16; } i f (pips >= 8) { choices.push_back(8); pips -= 8; } i f (pips >= 4) { choices.push_back(4); pips -= 4; } i f (pips >= 2) { choices.push_back(2); pips -= 2; } i f (pips >= 1) { choices.push_back(l); 144 pips -= 1; } i f (choices .s ize() == 1) local_flow = ( int) choices[0]; else { int index = b i l l y . i n t e g e r ( c h o i c e s . s i z e ( ) - l , 0); local_flow = ( int) choices[index]; } > outdirgrid[xpos][ypos] = local_flow; } } / / END OF SECOND PASS / / Resolve flow for f l a t nodes by i t e r a t i o n / / flow may be to any node that has an assigned d i r e c t i o n provided / / that node does not flow back to the o r i g i n a l . / / the while loop may get stuck i f there are extensive areas of / / perfect f latness i f i t does, ta i l_chaser b a i l s i t out and dnet / / w i l l c a l l f i l l _ s i n k s . Few landscapes are perfec t ly f l a t / / so th i s i s not usual ly a problem. int ta i l_chaser = 0; while ( ( is_a_zero(outdirgrid)) kk ( tai l_chaser < NSIDE*NSIDE)) { tail_chaser++; for ( int ypos = 2; ypos < NSIDE-2; ypos++) { for ( int xpos = 2; xpos < NSIDE-2; xpos++) { i f (outdirgrid[xpos][ypos] == 0) { in t f l a t . f l o w = 0; dvec f lat_choices; for ( int neighbour = 0; neighbour .< 8; neighbour++) { int xj = xpos + xlookup[neighbour]; int yj = ypos + ylookup[neighbour]; i f ( ( inaltgrid[xpos][ypos] - i n a l t g r i d [ x j ] [ y j ] == 0) kk (outd irgr id[xj ] [y j ] != 0) kk (outd irgr id[xj ] [y j ] != flowback[neighbour])) { double d i m = pow(2, neighbour); f lat_choices .push_back(dirn); } } i f ( f la t_choices . s ize ( ) == 0) f lat_f low = 0 ; / / n o change else i f ( f lat_choices . s i z e O == 1) f lat_f low = ( int) f lat_choices[0] ; 145 else { in t index = b i l l y . i n t e g e r ( f l a t _ c h o i c e s . s i z e ( ) - l , 0); f lat_f low = ( int) f lat_choices[ index]; > outdirgrid[xpos][ypos] = f lat_f low; } } } / / END OF THIRD PASS } } in t drain(Matrix i n g r i d , int x i , int y i ) { / / f i n d upslope area of point x i , y i on the / / f lowdirect ion g r i d i n g r i d / / (recursive) int cel lcount = 0; for ( int neighbour = 0; neighbour < 8; neighbour++) { int xj = x i + xlookup[neighbour]; int yj = y i + ylookup[neighbour] ; i f ((xj > 0) && (xj < NSIDE) && (yj > 0) && (yj < NSIDE)) { i f ( ingr id[x j ] [y j ] == flowback[neighbour]) { int upstream = d r a i n ( i n g r i d , x j , y j ) ; cel lcount += 1 + upstream; } > } return ce l lcount; } void accumulate(Matrix i n d i r g r i d , Matrix outareagrid) { / / given a g r i d of f low-direct ions , return upslope area / / of each c e l l for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { outareagrid[xpos][ypos] = d r a i n ( i n d i r g r i d , xpos, ypos); } } > void Bagnold(Matrix acc, Matrix slope, Matrix r i v e r f l u x , double k_f luv ia l ) { / / convert upslope area into (bedload) sediment f lux 146 / / define parameters that give regional context double k _ f l = k _ f l u v i a l ; double c_chan_maint = hftrans; / / m"2 -> km"2 for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { double upslparea = acc[xpos][ypos] * (SIDELENGTH * SIDELENGTH); / / -> m"2 double q = k _ f l * pow(upslparea / 1000000, 0.68); i f (upslparea >= c_chan_maint) { i f (slope[xpos][ypos] < 0.01) { / / compute sediment transport rate using / / normal Bagnold r e l a t i o n / / (Martin and Church, 1999; Knighton, 1998) double p o i n t d i f f = 6.44 * pow(q, 0.35) * pow(slope[xpos] [ypos], 1.5); riverflux[xpos][ypos] = p o i n t d i f f ; } else i f (slope[xpos] [ypos] < 0.15) { / / compute sediment transport rate using / / headward Bagnold r e l a t i o n (Martin, 1998; Day, 1969) double p o i n t d i f f = 40.32 * pow(q, 0.34) * pow(slope [xpos] [ypos], 2.22); riverflux[xpos][ypos] = p o i n t d i f f ; } else { riverflux[xpos][ypos] = 0; } } else r i v e r f l u x [xpos] [ypos] = 0; / / no f l u v i a l processes above channel head } } > void fluxsep(Matrix riverchange, Matrix r i v e r f l u x , Matrix flow, f loat t) { 147 / / use predicted t o t a l f lux at each point and input f lux from upstream / / neighbours to compute amount of erosion or deposit ion at each point / / change = input - output for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { double inflow = 0; for ( int neighbour = 0; neighbour < 8; neighbour++) { int xj = xpos + xlookup[neighbour]; int yj = ypos + ylookup[neighbour]; i f ((xj > 0) kk (xj < NSIDE) kk (yj > 0) kk (yj < NSIDE)) { i f (flow[xj] [yj] == flowback[neighbour]) { inflow += r i v e r f l u x [ x j ] [ y j ] ; } } } double netchange = inflow - r iverf lux[xpos][ypos]; riverchange[xpos] [ypos] = netchange / (SIDELENGTH * SIDELENGTH); } } / / write t o t a l bedload sediment f lux to sedtransp, as a row of / / NSIDE*NSIDE+1 columns / / ( f i r s t column i s time) / / . . . records one row every time dnet i s run from main. . one row / / each 100 years i f (t != -1) { sedtransp « t « " "; for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { sedtransp << riverflux[xpos][ypos] << " " ; } } sedtransp « endl; } > 148 / / DIFFUSION void expnonlinear(Matrix i n i t i a l , Matrix f i n a l , Matrix production, Matrix u p l i f t , f loat s t a r t , f loat stop, f loat dt , s t r i n g kmodel) { / / E x p l i c i t 2-D solver for non- l iear landscape d i f f u s i o n equation / / with i n i t i a l surface, from start to stop with i n t e r v a l dt / / using kmodel spec i f ied and including terms for / / production of material ( tectonics , r i v e r s etc . ) ( s p a t i a l l y dependent) / / with resu l t s written to f i n a l (which must exist when fn i s cal led) Matrix U, V; matr ix_copy( in i t i a l ,U) ; f i l l ( V , "zeros"); int toggle = 1; double Kl = K / 2; for ( f loat time = s tar t ; time <= stop; time+=dt) { i f (toggle == 1) { / / set boundary using method of images for ( int ypos = 1; ypos < NSIDE-1; ypos++) { / / l e f t and r ight bnds V[0] [ypos] = U[2] [ypos] ; V[NSIDE-1] [ypos] = U[NSIDE-3] [ypos]; } for ( int xpos = 1; xpos < NSIDE-1; xpos++) { / / top and bottom bnds V[xpos] [0] = ULxpos] [2] ; V[xpos][NSIDE-1] = U[xpos][NSIDE-3]; } / / set diagonal image at corners V[0] [0] = U[2] [2] ; V[0] [NSIDE-1] = U[2] [NSIDE-3] ; V[NSIDE-1] [0] = U[NSIDE-3][2]; V[NSIDE-1][NSIDE-1] = U[NSIDE-3][NSIDE-3]; for ( int ypos = 1; ypos < NSIDE-1; ypos++) { for ( int xpos = 1; xpos < NSIDE - 1; xpos++) { / / Compute some space derivat ives 149 double XA = (U[xpos+l][ypos] - U[xpos][ypos]) / SIDELENGTH; double YA = (U[xpos+l][ypos+1] + U[xpos] [ypos+1] - U[xpos+l] [ypos] - U [xpos] [ypos]) /(2*SIDELENGTH); double XB = (U[xpos][ypos] - U[xpos-l][ypos]) / SIDELENGTH; double YB = (U[xpos-l] [ypos+1] + U[xpos][ypos+1] - U[xpos-l][ypos] - U[xpos][ypos]) /(2*SIDELENGTH); double YC = (U[xpos] [ypos+1] - U[xpos][ypos]) / SIDELENGTH; double XC = (U[xpos][ypos] + U[xpos][ypos+1] - U[xpos-l] [ypos+1] - U[xpos-l][ypos]) /(2*SIDELENGTH); double YD = (U[xpos][ypos] - U[xpos][ypos-1]) / SIDELENGTH; double XD = (U[xpos] [ypos-1] + U[xpos] [ypos] - U[xpos-l] [ypos-1] - U[xpos-l][ypos]) /(2*SIDELENGTH); / / compute sediment f lux using non-l inear function double slpA = XA + YA double slpB = XB + YB double slpC = XC + YC double slpD = XD + YD double QA = K l * (1 + tanh( (abs( slpA ) - hxO) / hscale)) * (slpA) double QB = K l * (1 + tanh( (abs( slpB ) - hxO) / hscale)) * (slpB) double QC = K l * (1 + tanh( (abs( slpC ) - hxO) / hscale)) * (slpC) double QD = K l * (1 + tanh( (abs( slpD ) - hxO) / hscale)) * (slpD) / / compute e levat ion change double dif fuse = (dt/SIDELENGTH)*(QA-QB+QC-QD); i f ((slpA > 1.6) | | (slpB > 1.6) II 150 (slpC > 1.6) I I (slpD > 1.6)) { / / turn off d i f fus ion on slopes > " 60 degrees diffuse = 0; } double change = diffuse + (production[xpos][ypos] * dt) + (uplift[xpos] [ypos] * dt ) ; V[xpos] [ypos] = U[xpos] [ypos] + change; } } } i f (toggle == -1) { / / set boundary using method of images for ( int ypos = 1; ypos < NSIDE-1; ypos++) { / / l e f t and r ight bnds U[0] [ypos] = V[2] [ypos] ; U[NSIDE-1][ypos] = V[NSIDE-3][ypos]; } for ( int xpos = 1; xpos < NSIDE-1; xpos++) { / / top and bottom bnds U[xpos] [0] = V[xpos] [2] ; U[xpos] [NSIDE-1] = V[xpos] [NSIDE-3]; } / / set diagonal image at corners U[0] [0] = V[2] [2] ; U[0][NSIDE-1] = V[2][NSIDE-3]; U[NSIDE-1] [0] = V[NSIDE-3] [2]; U[NSIDE-1][NSIDE-1] = V[NSIDE-3][NSIDE-3]; for ( int ypos = 1; ypos < NSIDE-1; ypos++) { for ( int xpos = 1; xpos < NSIDE - 1; xpos++) { / / Compute some space derivat ives double XA = (V[xpos+l][ypos] - V[xpos][ypos]) / SIDELENGTH; double YA = (V[xpos+l][ypos+1] + V [xpos] [ypos+1] - V[xpos+l] [ypos] - V[xpos][ypos]) /(2 * SIDELENGTH); double XB = (V[xpos] [ypos] - V[xpos-1] [ypos]) / SIDELENGTH; 151 double YB = (V[xpos-l][ypos+1] + V[xpos][ypos+1] - V[xpos-l][ypos] - V[xpos] [ypos] ) /(2*SIDELENGTH); double YC = (V[xpos] [ypos+1] - V[xpos][ypos]) / SIDELENGTH; double XC = (V[xpos][ypos] + V[xpos][ypos+1] - V[xpos-l][ypos+1] - V[xpos-l] [ypos]) /(2*SIDELENGTH); double YD = (V[xpos][ypos] - V[xpos] [ypos-1]) / SIDELENGTH; double XD = (V[xpos] [ypos-1] + V[xpos] [ypos] - V[xpos-l][ypos-1] - V[xpos-l] [ypos]) /(2+SIDELENGTH); / / compute sediment f lux using non-l inear function double slpA double slpB double slpC double slpD = XA + YA = XB + YB = XC + YC = XD + YD double QA = K l * (1 + tanh( (abs( slpA ) - hxO) / hscale)) * (slpA) double QB = K l * (1 + tanh( (abs( slpB ) - hxO) / hscale)) * (slpB) double QC = K l * (1 + tanh( (abs( slpC ) - hxO) / hscale)) * (slpC) double QD = K l * (1 + tanh( (abs( slpD ) - hxO) / hscale)) * (slpD) / / compute e levat ion change double diffuse = (dt/SIDELENGTH)*(QA-QB+QC-QD); i f ((slpA > 1.6) I I (slpB > 1.6) | | (slpC > 1.6) I I (slpD > 1.6)) { / / turn off d i f fus ion on slopes > " 60 degrees diffuse = 0; } double change = diffuse + (production[xpos] [ypos] * dt) + (uplift[xpos][ypos] * dt ) ; 152 U[xpos][ypos] = V[xpos][ypos] + change; } } } toggle *= -1; } i f (toggle == 1) matrix_copy(U, f i n a l ) ; else i f (toggle == -1) matrix_copy(V, f i n a l ) ; else { cout « "lscape (expsolv): in terna l error" << endl; e x i t ( l ) ; } > / / DEEP-SEATED BEDROCK LANDSLIDES void dsbl(Matrix landscape, Matrix s l idescape, Matrix s l i d e r e c o r d , Matrix slope, Matrix flow, f loat t) { / / modify landscape to create slidescape based on the occurrence of / / deep-seated bedrock landsl ides; s l ides occur at random. / / most l i k e l y s i tes are i d e n t i f i e d using model nodes where gradient / / e levat ion and length of slope are highest / / put time i n s l i d e f i l e s l i d e f i l e « t « endl; / / Generate s l i d i n g p r o b a b i l i t i e s Matrix s l ideprob; matrix_copy(landscape, s l ideprob); matrix_copy(landscape, s l idescape); for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { slideprob[xpos] [ypos] *= slope [xpos] [ypos]; } } / / Make a l i s t of p r o b a b i l i t i e s Landnode nodelist[NSIDE*NSIDE]; for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { 153 nodelist[xpos*NSIDE+ypos].x(xpos); nodelist[xpos*NSIDE+ypos].y(ypos); i f ((xpos <= 5) | | (xpos >= NSIDE-5) II (ypos <= 5) | | (ypos >= NSIDE-5)) { / / zero p r o b a b i l i t y of s l i d i n g near edge nodelist[xpos*NSIDE+ypos].elev(O); } else { nodelist[xpos*NSIDE+ypos].elev(slideprob[xpos][ypos]); } } } / / Sort the l i s t of p r o b a b i l i t i e s sor t (node l i s t , (NSIDE*NSIDE)); / / Choose those p r o b a b i l i t i e s greater than a c r i t i c a l threshold in t l i m i t = -1; bool steep_end = fa l se ; for ( int probcount = 0; probcount < (NSIDE*NSIDE); probcount++) i f (nodelist[probcount] <= crit_prob) limit++; } / / cout « l i m i t « " ," « (NSIDE*NSIDE)-1 « endl; i f ( l imi t == (NSIDE*NSIDE)-1) { / / b a i l out here i f we arent having a s l ide s l i d e f i l e << "no s l ide occurred" « endl; } else { / /otherwise choose a s l ide locat ion at random int xpos, ypos, cougar; bool selected = fa l se ; cougar = bi l ly. integer((NSIDE*NSIDE)-1, l i m i t ) ; xpos = ( int) nodel i s t [cougar] .xcoordO; ypos = ( int) nodel i s t [cougar] .ycoordO; s l i d e f i l e << nodelist[cougar] « endl; s l idedata « nodelist[cougar] « " "; sliderecord[xpos] [ypos] = 1; / / Calculate the maximum f a i l u r e size at the chosen s i t e double current_height = slidescape[xpos][ypos]; 154 s l i d e f i l e << "current height = " « current_height << endl; int nextx = dsxpos(xpos, flow[xpos][ypos]); int nexty = dsypos(ypos, flow[xpos][ypos]); s l i d e f i l e « "Next c e l l : " « nextx « " ," « nexty « endl; double next_cell_height = slidescape[nextx][nexty]; s l i d e f i l e « "Next c e l l height = " « next_cell_height << endl; double biggest_sl ide = current.height - next_cell_height - (SIDELENGTH * tan(alpha)); i f (biggest_sl ide < 0) biggest_sl ide = 0 ; / / i f slope not steep enough s l i d e f i l e << "biggest s l ide = " « biggest_sl ide << endl; double ndmag = 1 ; / / a l l avai lable r e l i e f f a i l s s l i d e f i l e « "weibull proportion = " « ndmag « endl; double s l ide_s ize = ndmag*biggest_slide; s l i d e f i l e << "actual s l ide size = " << s l ide_s ize « endl; / / Begin volume check during s l ide double pre_vol = volume(landscape); / / reduce e levat ion at f a i l u r e s i t e slidescape[xpos][ypos] -= s l ide_s ize ; / / decide how much talus we have generated double stash = sl ide_size*(rho_rock / rho_sed); / / allow for density change s l i d e f i l e « "available r e g o l i t h = " « stash << endl; double t t l_dep_vol = stash; / / track t t l amount moved double ds_difference = SIDELENGTH*tan(alpha); / / difference to maintain / / angle of repose / / move the material downslope u n t i l i t runs out int dscount = 0 ; / / how many times have we moved downslope with debris int current_xpos = nextx; / / s tart at node downslope of or ig ina t ing node int current_ypos = nexty; while (stash > 0) { s l i d e f i l e << "Moved downstream" << endl; s l i d e f i l e << "at locat ion: " « current_xpos << ", " « current_ypos « " "; dscount++; double current.height = slidescape[current_xpos][current_ypos]; int new_xpos = dsxpos(current_xpos, flow[current_xpos][current_ypos]); int new_ypos = dsypos(current_ypos, flow[current_xpos][current_ypos]); 155 s l i d e f i l e << "done backflow" « endl; i f ((new_xpos <= 0) I I (new_xpos >= NSIDE-1) I I (new_ypos <= 0) I I (new_ypos >= NSIDE-1)) { / / i n case we go off the edge s l i d e f i l e « "lscape (dsbl): a freak s l ide has occurred: " << stash << " units of sediment have l e f t the model domain."; stash = 0; > else { double ds_height = slidescape[new_xpos][new_ypos]; double ideal_height = ds_height + ds_difference; double required_height = ideal_height - current_height; s l i d e f i l e << "downslope height = " << ds_height « " ;ideal_height = << ideal_height « " ; reqdheight = " « required_height « endl; i f (required_height >= 0) { / / i f we are i n a deposi t ional node s l i d e f i l e << "Depositing: "; i f (stash >= required_height) { / / i f there i s enough material deposit to ideal_height slidescape[current_xpos][current_ypos] = ideal_height; stash -= required_height; s l i d e f i l e << required_height << " units" << endl; } else { / / otherwise just leave the entire amount that we have slidescape[current_xpos][current_ypos] += stash; s l i d e f i l e « stash « " units" « endl; stash = 0; } } else { / / otherwise the s l ide scours the downslope area s l i d e f i l e « "Scouring "; slidescape[current_xpos][current.ypos] = ideal_height; double surplus = slidescape[current_xpos][current_ypos] - ideal_height; 156 s l i d e f i l e « surplus << " units" << endl; surplus *= (rho_rock / rho_sed); stash += surplus; t t l_dep_vol += surplus; / / add to t o t a l s l ide volume } current_xpos = new_xpos; current_ypos = new_ypos; > } double post_vol = volume(slidescape); / / to check conservation of volume s l i d e f i l e « "Total Volume: " « t t l_dep_vol « endl; s l i d e f i l e « "Runout length: " << dscount « endl; s l i d e f i l e « t t l_dep_vol « " " « dscount « endl; s l i d e f i l e « "Pre-vol: " « pre_vol « " . Post -vo l : " « post_vol « endl s l idedata « t t l_dep_vol « " " « dscount « endl; cout « "DSBL @ (" « xpos « " ," « ypos « "); Vol = " « t t l_dep_vol * (SIDELENGTH * SIDELENGTH) « " m~3 ; Runout = " « dscount * SIDELENGTH « " m." « endl; } } / / INPUT-OUTPUT FUNCTIONS void export(Matrix output, bool d i sp lay , int id) { / / write a f i l e containing output, named lscape<id>.out. / / i f the f l a g display i s t rue , the output i s also writ ten to the screen. s t r i n g pref ix = "lscape"; int tens = ( in t ) ( id /10) ; int units = i d - (10*tens); char tnum = tens + 48; char unum = units + 48; s t r ing filename = pref ix + tnum + unum + ".out"; ofstream out f i l e ( f i l ename.c_s tr ( ) ) ; for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { o u t f i l e « output[xpos][ypos] « " "; 157 } o u t f i l e « endl; } i f (display == true) { for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { cout « output[xpos][ypos] << " "; } cout « endl; } } } void f i l l ( M a t r i x input, s t r ing f i l e ) { / / F i l l a matrix, input, with values taken from the f i l e / / i f f i l e i s "random" or "zeros" then no f i l e i s read / / and the matrix i s f i l l e d with uniform random values or / / zeros accordingly i f ( f i l e == "random") { for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { input[xpos] [ypos] = b i l ly . integer(HMAX,0) ; > } } else i f ( f i l e == "zeros") { for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { input[xpos][ypos] = 0; } } } else i f ( f i l e == "nodata") { for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { input[xpos] [ypos] = -9999; } } } else i f ( f i l e == "sine") { 158 for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { i f ((xpos ==0) II (ypos ==0) I I (xpos == NSIDE - 1) I I (ypos == NSIDE - 1)) input[xpos][ypos] = 0; else input[xpos][ypos] = (2*(NSIDE)) - ((xpos) * sin((4*Pi*ypos)/(NSIDE))); } } } else i f ( f i l e == "comb") { for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { i f ((xpos ==0) II (ypos ==0) | | (xpos == NSIDE - 1) | | (ypos == NSIDE - 1)) input[xpos][ypos] = 0; else i i f (is_even(xpos)) input[xpos] [ypos] = (xpos*NSIDE + ypos); else input[xpos] [NSIDE - ypos - 1] = (xpos*NSIDE + ypos); } } } } else i f ( f i l e == "ramp") { double step_diff = HMAX / (NSIDE / 2); for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { i f (xpos < (NSIDE / 2)) input[xpos][ypos] = xpos * s tep_dif f ; else input[xpos] [ypos] = (2 * HMAX) - (xpos * s tep_di f f ) ; > } } else i f ( f i l e == "test") { / / test data from Jenson and Domingue (1988) double testinput[NSIDE][NSIDE] {{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, { 0,778,765,750,740,747,759,765,766,769,776,786,795, 0}, { 0,770,758,745,737,741,751,753,761,777,789,802,814, 0}, { 0,777,763,747,736,735,743,750,767,787,806,820,832, 0}, { 0,786,767,750,737,733,739,752,769,785,797,808,822, 0}, { 0,794,773,756,741,733,733,744,759,772,779,789,806, 0}, { 0,799,782,763,750,737,733,733,745,757,767,782,801, 0>, { 0,802,788,771,761,751,736,733,738,751,764,779,798, 0}, { 0,799,790,780,772,762,746,733,737,754,770,784,794, 0}, 159 { 0,811,799,787,771,757,741,728,730,745,765,779,783, 0}, { 0,823,807,790,774,762,748,733,725,733,750,764,763, 0}, { 0,830,814,801,787,776,761,743,728,725,737,748,751, 0}, { 0,822,818,811,801,791,776,757,739,726,725,735,751, 0}, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}; for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { / / rotate 90 degrees because i t easier to enter l i k e that above! input[xpos][ypos] = testinput[ypos][xpos]; / / / 832; > } } else { i fstream i n f i l e ( f i l e . c _ s t r ( ) ) ; i f ( U n f i l e ) { cerr << "lscape ( f i l l ) : cannot open input f i l e : " << f i l e << endl; ex i t (1) ; } cout << "lscape ( f i l l ) : reading i n i t i a l values from: " « f i l e « endl; s t r ing datal ine; for ( int ypos = 0; ypos < NSIDE; ypos++) { { get1ine(inf i l e , dat a l i n e ) ; double temp; const char* buffer = d a t a l i n e . c _ s t r ( ) ; / / d a t a l i n e -> c - s t r i n g for streaming istrstream i s s (buf fer , 10000); for ( int xpos = 0; xpos < NSIDE; xpos++) { i s s » temp; input[xpos][ypos] = temp; } } > } } / / AUXILLIARY FUNCTIONS double volume(Matrix landscape) { / / compute the approximate volume of landscape double landvolume = 0; 160 double base_area = SIDELENGTH * SIDELENGTH; for ( int ypos = 0; ypos < NSIDE - 2; ypos+=2) { for ( int xpos = 0; xpos < NSIDE - 2; xpos+=2) { double t t l _ e l e v = landscape[xpos][ypos] + landscape[xpos+1] [ypos] + landscape[xpos] [ypos+1] + landscape[xpos+1][ypos+1]; landvolume += (0.25*ttl_elev*base_area); } } return (landvolume); > void t i t l e O { cout « "erode2d: non-l inear d i f fus ion with threshold lands l id ing" « endl; cout « "Diffusive-advective landform evolution model" << endl; cout « "Written by Simon Dadson, A p r i l 2000" « endl; } void writelogO { // write log f i l e ofstream l o g f i l e C l s c a p e . l o g " ) ; l o g f i l e « "ER0DE2D Log F i l e l scape . log: " « endl; l o g f i l e « " I n i t i a l surface = " << i n p u t f i l e « endl; l o g f i l e « "SIDELENGTH = " « SIDELENGTH « endl; l o g f i l e « "TST0P = " « TST0P « endl; l o g f i l e « "NSIDE = " « NSIDE « endl; l o g f i l e « "NT = " « NT « endl; l o g f i l e « "HMAX = " « HMAX « endl; l o g f i l e « "D = " « D « endl; l o g f i l e « "K, hxO, hscale" « endl; l o g f i l e « K « " " « hxO « " " « hscale « endl; l o g f i l e « " s t a b i l i t y parameter = " « r « endl; l o g f i l e « "Pcrit = " « cr i t_prob « endl; l o g f i l e « "alpha = " « alpha « endl; l o g f i l e « "K_f luv ia l = " « k _ f l u v i a l « endl; l o g f i l e « "hftrans = " « hftrans « endl; l o g f i l e « "upl i f t model = " « upliftmodel « endl; l o g f i l e « "max u p l i f t = " « Umax << endl; l o g f i l e << "slide test i n t e r v a l = " « s l ide_tes t_ int « endl; l o g f i l e << "rho_rock = " « rho_rock « endl; l o g f i l e << "rho_sed = " « rho_sed « endl; l o g f i l e << "diagnose = " << diagnose << endl; 161 int dsxpos(int oldx, double f lowdir) { / / f i n d x coordinate of point downstream of oldx, given f lowdir int newx = -1; for ( int nc = 0; nc < 9; nc++) { i f (f lowdir == pow(2, nc)) { newx = oldx + xlookup[nc]; / / back to front ! > } i f (newx == -1) { cout « "lscape (dsxpos): error rout ing sediment" « endl; ex i t (1) ; } return (newx); > in t dsypos(int o ldy, double f lowdir) { / / f i n d y coordinate of point downstream of oldy, given f lowdir in t newy = -1 ; for ( int nc = 0; nc < 9; nc++) { i f (flowdir == pow(2, nc)) { newy = oldy + ylookup[nc]; > } i f (newy == -1) {. cout « "lscape (dsypos): error routing sediment" « endl; e x i t ( l ) ; } return (newy); > void matrix_copy(Matrix o l d , Matrix newmatrix) { / / copy o ld into new / / new must exist p r i o r to the c a l l for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { newmatrix[xpos][ypos]=old[xpos][ypos]; } 162 } } void matrix_add(Matrix f i r s t , Matrix second) { // add second to f i r s t , result i s i n f i r s t , // second i s unchanged // both must exist prior to the c a l l for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { first[xpos][ypos]+=second[xpos][ypos] ; } } } bool is_even(int t e s t i n t ) { // return true i f t e s t i n t i s even; false i f odd fl o a t teststat = 0; teststat = testint°/,2; i f (teststat != 0) return f a l s e ; else i f (teststat == 0) return true; else exit(1); bool is_mult(double a, double b) { // return true i f a i s a multiple of b double teststat = 0; teststat = int(a)'/,int(b) ; i f (teststat != 0) return f a l s e ; else i f (teststat == 0) return true; else exit(1); void getparamsO { cout << "Enter name of input f i l e : " ; cin » i n p u t f i l e ; int location; cout « "Use nonlinear d i f f u s i o n rule for GVRD (1) or QCI (2) cin >> location; i f (location == 1) { K = 0.12; hxO = 0.81; 163 hscale = 0.32; } else i f ( locat ion == 2) { K = 0.43; hxO = 1.06; hscale = 0.19; } else { cout « "Invalid location" « endl; ex i t (1) ; } cout << "Enter u p l i f t s tyle (none, const, thrust , p luton, rebound): "; c in >> upl i f tmodel; cout « "Enter maximum u p l i f t rate (metres per year): "; c i n » Umax; cout « "Enter coef f ic ient for f l u v i a l transport (k): "; c i n » k _ f l u v i a l ; cout « "Enter h i l l s l o p e - f l u v i a l t r a n s i t i o n ( in upslope area units (m"2)) "; c i n » hftrans; cout << "Enter minimum dsbl recurrence i n t e r v a l (years): "; c in » s l ide_tes t_ int ; cout << "Enter c r i t i c a l s l ide probab i l i ty (angle * e levat ion): "; c i n » cr i t_prob; } bool is_a_zero(Matrix grid) { / / return true i f g r i d contains any zero, fa l se otherwise, bool found_one = fa l se ; for ( int ypos = 0; ypos < NSIDE; ypos++) { for ( int xpos = 0; xpos < NSIDE; xpos++) { i f (grid[xpos][ypos] == 0) found_one = true; } } return found_one; } void threshold(Matrix i n , Matrix out, double cutoff) { / / produce out from i n , replacing every value greater than cutoff with 1 / / and every value lower or equal to with -9999 (arc N0DATA) for ( int ypos = 0; ypos < NSIDE; ypos++) { 164 or ( int xpos = 0; xpos < NSIDE; xpos++) { i f (in[xpos] [ypos] > cutoff) out[xpos][yp< else out[xpos][ypos] = -9999; 165

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
United States 48 2
Japan 5 0
France 4 0
Netherlands 3 2
Romania 3 0
China 2 2
City Views Downloads
Kansas City 39 0
Unknown 10 1
Tokyo 5 0
Washington 4 0
Mountain View 3 0
Hangzhou 2 0
Redmond 1 0
Ashburn 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}

Share

Share to:

Comment

Related Items