UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Episodic processes in fluvial landscape evolution Dadson, Simon 2000

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

Item Metadata


ubc_2000-0377.pdf [ 9.23MB ]
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
Original Record: 1.0089851 +original-record.json
Full Text

Full Text

EPISODIC PROCESSES IN FLUVIAL LANDSCAPE EVOLUTION by Simon Dadson B.A. (Hons), University of Oxford, 1998 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in THE FACULTY OF GRADUATE STUDIES (Department of Geography) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA 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 Date (o S&PTCMOK XQOQ DE-6 (2/88) Abstract 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 viList of Figures viii Acknowledgements xChapter 1 Introduction 1 1.1 Aims1.1.1 Spatial Organization of Landscape Processes 2 1.2 Fractal Models of Fluvial Landscapes 5 1.2.1 Fractal Properties of Topography1.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 Tectonic background and Geological History 16 2.2.1 Tectonics and Topography 12.3 Queen Charlotte Islands 17 2.3.1 Bedrock Geology 8 iii 2.3.2 Glacial Legacy 18 2.3.3 Contemporary Processes 19 2.4 GVRD Capilano Watershed . . . 20 2.4.1 Bedrock Geology 22.4.2 Tectonic Uplift 1 2.4.3 Glacial Legacy 22.4.4 Contemporary Processes 22 2.5 Summary 23 Chapter 3 Sediment Flux due to Episodic Events 24 3.1 Process rates 23.2 Sources of Data3.3 Surficial Slides: magnitude-frequency 25 3.3.1 Results 23.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 23.3.5 GVRD Slides 31 3.3.6 Implications 2 3.4 Slope and Volumetric Transport Rate 33 3.4.1 Approach 33.4.2 Parameterization 3 3.4.3 Discussion 36 3.5 Rock Slides3.5.1 Rock Slide Inventory 37 3.5.2 Discussion 3iv 3.6 Summary 39 Chapter 4 Fractal properties of topography 40 4.1 Geomorphic Processes and Fractal Landscapes4.2 Fractal Properties of Landscape 41 4.3 Measuring Fractal Properties of Landscape 43 4.4 Process-domains and the Fractal Dimension 4 4.5 Results 46 4.5.1 GVRD Capilano Watershed 44.5.2 Queen Charlotte Islands 60 4.6 Discussion 69 Chapter 5 Synthesis: Process-Based Modelling with Episodic Events 72 5.1 Introduction 75.2 Constructing Testable Landscape Models 73 5.3 Framework for Landscape Modelling 4 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 1 5.3.5 Numerical Schemes 85.3.6 Stability 82 5.4 Results: Preliminary Runs 83 5.5 Results: Detailed Model Runs 8 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 9 6.2 Theoretical Implications 100 Bibliography 102 Appendix A Computation of Fractal Dimension 112 A.l Slope-Area Plot Constructions 11A. 2 Construction of Hillslope Variograms 116 Appendix B Derivation of Landscape Modelling Equations 132 B. l Continuity Equation 13B.2 Fluvial Transport Rule 3 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 GVRD Landslides 34 3.2 Parameters for Non-Linear Slope-Volume Relation 34 4.1 GVRD Capilano Channel Network Fractal Dimensions 47 4.2 GVRD 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 GVRD Capilano Watershed (study basins are described here and in Chapter Four) 5 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 GVRD Slope-Volume Relation 35 3.6 Rock Slide Magnitude Frequency (data from Savigny 1991) 38 4.1 Zones for Computation of D 44.2 Slope-Area Plots: GVRD Daniels Creek 49 4.3 Slope-Area Plots: GVRD Sisters Creek 50 4.4 Slope-Area Plots: GVRD Enchantment Creek 51 4.5 Hillslope Variograms: GVRD Daniels Creek (Dl - D2) 53 4.6 Hillslope Variograms: GVRD Daniels Creek (D3 - D4) 54 4.7 Hillslope Variograms: GVRD Sisters Creek (SI - S2) 55 4.8 Hillslope Variograms: GVRD Sisters Creek (S3 - S4) 56 4.9 Hillslope Variograms: GVRD Enchantment Creek (El - E2) 57 4.10 Hillslope Variograms: GVRD 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)5.4 Lembke Creek: 10 ka no rock slides 91 5.5 Lembke Creek: 10 ka with rock slides 2 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 ARC/INFO 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^TEX 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 DEM data); Derek Bonin, Dave Dunkley, Lome Gilmour (Greater Vancouver Regional District DEM 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. SIMON DADSON The University of British Columbia August 2000 xi Chapter 1 Introduction 1.1 Aims 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-10 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 M2 (H = 1.69294) x separation (h) M3 (H = 2.49445) • M4(H = 3.25169) + 1000 M5 (H = 3.96022) o Fractal Regression — -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 Objectives 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: Location 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: GVRD 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 till masking the surface and has resulted in at least 18 m of isostatic rebound since 5 ka BP (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 BP 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 PTP granites 1.3 Table 2.1: Sediment Flux due to Landsliding on Queen Charlotte Islands (data from Rood 1984) 2.4 GVRD Capilano 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 GVRD, 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 GVRD 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 GVRD 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 BP 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; GVRD, 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 GVRD (1999), Savigny (1991) and English (1996). Campbell (1999) recently determined a sediment yield due to landsliding of 168 m3 km-2 a-1 (= 1.7 m3ha-1a-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 GVRD (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 Flux due to Episodic 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 GVRD 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 Data 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 GVRD 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 GVRD 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 GVRD 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 m3). Frequency = 0.8 * Magnitude-063 (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 m3), 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 m3), 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 m3, 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 m3) 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 GVRD Slides The magnitude-frequency distribution of the GVRD 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 GVRD 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 Volumetric Transport Rate 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 GVRD landslide inventory. Landslide width, length, thickness and initiation slope were reported in the GVRD Ecological Inventory (1999), with measure ments made from air photos and calibrated using field data. Volumes were computed using a relation presented in the GVRD 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 GVRD). 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 (GVRD, 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 GVRD 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 GVRD 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 Mm3). 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 Mm3, 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 Fractal properties of topography 4.1 Geomorphic Processes and Fractal 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 Measuring Fractal Properties 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 Fractal Dimension 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 ~102 - 103 m2. 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 DEM of the study area. Upslope area was found using ARC/INFO'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). An 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 Results 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 GVRD Capilano Watershed Slope-area plots for a 10 m DEM of the Capilano Watershed (supplied by GVRD Watershed Planning Department) are shown in Figures 4.2, 4.3 and 4.4. The most 46 common hillslope-channel threshold is 2,000 m2. This indicates that that only 2,000 m2 (=0.002 km2) of contributing area is required for channel initiation in the GVRD 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). All 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: GVRD 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: GVRD Daniels Creek 49 Points on the channel network (i.e., with upslope area greater than 2,000 m2 (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 m2 (< 6 pixels), and the set L, with 600 m2 - 2,000 m2 (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 DEM; 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: GVRD Daniels Creek (Dl - D2) 53 I" " • I' • I I' • • I" 1 X It CL CL Q "_ 5 & CM to r~ w co cy CT ll ll co^r 9 s 11 tr X — "—C\J ^2 o « si ™ co lire S.-S 8" Z.E 5^ IS era r-~ co ,_: Figure 4.6: Hillslope Variograms: GVRD 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 si CO —' CM CO SI Si CT> © © O g O i-Figure 4.7: Hillslope Variograms: GVRD Sisters Creek (Si - S2) 55 luaiuow luauuow Figure 4.8: Hillslope Variograms: GVRD Sisters Creek (S3 - S4) 56 . ii tr C\J in cn ui RS CO to "—OJ ^2 8 4 CNJ £ r- cn co $ csi <ri II II S CO o Ol CO '—CM to if M • x k CL CL 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: GVRD Enchantment Creek (El - E2) 57 C KM + CJ X k -I ° _o co LU l CO w CM CO co r-O CO CM <*? CO —' 55 CM e II rr § 83 S O CM S" <° ^ <5 CM CO II II CO SCO w co oir-. co • 8 £ Figure 4.10: Hillslope Variograms: GVRD 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: GVRD 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 GVRD 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 BC TRIM data using the ARC/INFO T0P0GRID routine.) The most common hillslope-channel threshold is 6,000 m2 (= 0.006 km2). 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 GVRD, 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 GVRD 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 ii v E Ul/UI 8dO|S |B00| ui/ui edo|s |EOO| UI/UJ SdO|S |E00| UI/UJ 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 QCR PTP 3000 1.61 Government 16.1 QCR Karmutsen 6000 1.67 Armentiere 4.0 QCR Karmutsen 6000 1.70 Jason 11.9 QCR Karmutsen 6000 1.71 Security 33.9 QCR 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; QCR = Queen Charlotte Ranges; PTP = Post-tectonic Plutons (see Chapter Two) The procedure for zoning elevations by upslope area was carried out: fluvial points with upslope area > 6,000 m2 (> 6 pixels) were removed; those pixels remaining were divided into two zones: upper hillslopes with area < 1,800 m2 (< 2 pixels) and a lower hillslope zone with 1,800 m2 < upslope area < 6,000 m2. These zones formed clip regions for the original BC TRIM 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 CO II II XX cn cn co 3 ;U8LUOLU JUDLUOLU 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 XX cofC CM CM CO *-*- LO "—'CM x in TJ ll II II CO^J-Figure 4.15: Hillslope Variograms: Queen Charlotte Islands (Skidegate Channel: Ar-mentieres and Government Creeks) 66 EE 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 GVRD result: 1. there is a lower incidence of large structure-generating landslides on the Queen Charlotte Islands than in the GVRD 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 GVRD 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 Modelling with Episodic Events 5.1 Introduction 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 Constructing Testable Landscape Models 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 Modelling Two common frameworks for long-term landform evolution modelling are discrete and continuous models1. An 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 GVRD Capilano watershed are given in Chapter Three. Uplift 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. Fluvial Erosion 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.44Q035 S1-5 if S < 0.01 (5.2) $F = < 40.32 Q0 34 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 ct-CD i-i H co >—1 • 3 P ct-CD CO CL CO o P as o o P o CL CD > 2 O o P [inimu i-i _ cr [inimu p CL ST o [inimu o' £L nne rocl CD 3 w CO 3 CL CD p ct-CD o CL CO CD T) O CD i-i CO Inte oba ct > CL CD rva cr •—1 CD P de ct-cr >-!_ co' 1 CD i-l P TO CD CD P i-l P 3 CD ct-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 pi ffi H 2 5' ^ h-» P " S" t—i 3^ co CD CD o o o CD co 4^ 8 cf 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 ct-h-i o CO co CO 03 i? 5 CL S pi P CD 1-1 sr ct-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 km2 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: Ir ,ftp ,. [$r/]fo]2 (0-002)(1000)2 Uplift Ratio =  (10Q0)(02) = 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 DEM 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 Results: Preliminary 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 GVRD 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. An 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 • ."ii f." '• .•• j"« d"• ^ • rf-• <• ^. WHS • 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 Results: Detailed Model 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 km2 sub-basin in the western side of the Capilano watershed (see Fig. 2.3); and a 4.4 km2 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. All 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 (m2) (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 m2)). 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 m2; and lowland areas with upslope area 1,800 - 5,400 m2). The procedures used for the estimation of these parameters were identical to those used in Chapter 4 for the original DEM 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 92 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 ARC/INFOs cubic interpolation algorithm). The effects of DEM 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). An 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 Conclusion 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 GVRD 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 GVRD 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 Practical Implications 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. An 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 Theoretical Implications 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 Bibliography 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, BC, 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, CA. 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. GVRD (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: An 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. Ill Appendix A Computation of Fractal Dimension A.l Slope-Area Plot Constructions 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 identifies hillslope-channel scale cutoff and ## channel network fractal dimension is computed. ## Based on algorithms in Montgomery and FoufoulaTGeorgiou (1993) and ## Veneziano and Iacobellis (1999) clear; closeplot; disp("Slope-Area Statistics. Simon Dadson (2000)."); dispC'Reading ARC/INFO data"); fflush(stdout); name = inputC'Enter a name for this analysis: ", "s"); lsl = strcatC'load area.", name); ls2 = strcatC'load slope.", name); 112 eval(lsl); 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., grid cell 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) results(l.l) = 0; results(:,3) = (results(:,1) + results(:,2)) ./ 2; ## do average for each area class disp("Computing averaged slope-area statistics"); 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); if tmp_area != -999900 & tmp_slope != -99.99 ## skip N0DATA values for i = l:divs if tmp_area >= results(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; results(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))); axis([l, area_max, sl_min, sl_max]); hold on; title(strcat("Slope-Area Plot: ", name)); xlabel("drainage area (square metres)"); ylabel("local slope m/m"); gset key left bottom title "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 is too small for the number of class divisions chosen. From this point on, empty class-bins will be ignored.");fflush(stdout); cutoff = input("Enter point for fluvial inflexion: "); for count = l:divs if cutoff >= results(count,1) & cutoff < results(count,2) first_point = count; endif endfor ## compute fractal dimension disp("Computing fractal dimension of river 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 = first_point:divs 114 tmpx = logl0(results(i,3)); tmpy = logl0(results(i,4+m)); if 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)), \ )..>•). II /i loglog(xfit, 10."(reg_intcpt(4) + loglO(xfit) * reg_slope(4)), \ ii J > theta = -1 * reg_slope(l); H_cn = 1 - (2 * theta) D_cn = 2 - H_cn tstr = ['gset key left bottom title "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 = strcatCgset output "slopearea' .name, '.ps'"); eval(outstr); replot; gset term xll; disp("The next plot will be a straight line if the river network is a self-similar fractal; if the river network is a multifractal, a non-linear curve will be seen."); fflush(stdout); dispC'Press a key to diagnose multifractality."); fflush(stdout); foo = kbhitO ; figure(2); axis([1,4,min(reg_slope), max(reg_slope)]); gset xtics 1; gset nokey; title(strcat("Slope 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 = strcatCgset output "mmtslp'.name, '.ps'"); eval(outstr2); replot; gset term xll; disp("Done."); A.2 Construction of Hillslope Variograms 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 ARC/INFO. /*mfs.cc Simon Dadson 2000 */ 116 /•Compute fractal dimension from */ /•ARC/INFO point coverage using */ /•variogram method */***********************************/ /* Algorithm set up parameters from user input or command-line read input file into linked-list of Point2d (see objgeom) select random subsample by purging points from list determine boundaries for log-spaced bins compute variogram for each point pair compute distance (and angle) fit pair into a bin and include it in the moment calculation complete moment calculation and return plot the variogram with octave ask user to identify range over which to compute D (unless running in batch mode, in which case choose the largest range) do linear regression & convert results to fractal statistics write output files */ #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) File i/o // (std) C++ string class // (std) support for streaming operations on strings // (std) support for linked lists // geometry algorithms and Point2D class // least-squares linear regression algorithms // random number generator from Knuth 1981 typedef list<Point2D> plist; // Defines Point2D list type typedef list<Point2D>::iterator plist_it; // See "objgeom" include file typedef vector<double> dvec; // define a double-precision vector class string runref; int h_start, h_stop, h_num, h_tol, anisodirn, adirn_tol = 0; float 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 lh, dvec lgamma, dvec regstats, dvec fracstats); plist readinputfile(const string& runref); dvec logspace(int start, 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); plist purge(plistfe, float); int fit(const doublefe d_ij, const dvecfe hvalues); dvec loglO (dvec v); int main(int argc, char* argv[]) { bool iact = true; if (argc == 1) UsageMsgO ; if (string(argv[l]) == string("-i")) ReadParameters(); else if (string(argv[l]) == ("-b")) { // interpret command line parameters if (argc != 9) UsageMsgO; iact = false; 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 ; plist biglist = readinputfile(runref); plist inputdata = purge(biglist, sample); dvec h = logspace(h_start, h_stop, h_num); dvec gamma = variogram(h, inputdata); 118 if (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); if (iact == true) OctaveFDPlot(h, gamma, regstats, f_bounds); WriteLog(runref); WriteOut(runref, h, gamma); WriteFrac(runref, loglO(h), loglO(gamma), regstats, fracstats); > void ReadParametersO { cout « "mfs v.2 SD 21 January 2000" « endl; cout « "Interactive Session" « endl; cout << "Enter runref (input file will be <runref>.in) "; cin >> runref; cout « "Enter sample size (fraction to keep) "; cin » sample; cout « "Enter initial value for h (h_start): "; cin » h_start; cout << "Enter final 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_tol; cout << "Enter anisotropy dirn (360-degree bearing from N; -1 = isotropic): "; cin » anisodirn; if (anisodirn != -1) { cout << "Enter angle tolerance (adirn_tol): "; cin » adirn_tol; } } void UsageMsgO { cout « "mfs v2" « endl; cout « "Usage: fd -i (interactive mode)" << endl; cout « " fd -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()); octfile « "load 'mmt." « runref « "';" « endl; octfile « "h = mmt(:,l);" « endl; octfile << "gamma = mmt(:,2);" « endl; octfile « "second = mmt(:,3);" « endl; octfile « "third = mmt(:J4);" « endl; octfile « "fourth = mmt(:,5);" « endl; octfile « "fifth = mmt(:,6);" « endl; octfile « "title('Variogram: " « runref « '");" « endl; octfile « "xlabelC separation (h)');" « endl; octfile « "ylabel('moment');" « endl; octfile « "gset key below;" « endl; octfile « "loglog(h,gamma,'@33;Ml;')" « endl; octfile « "hold on;" « endl; octfile « "loglog(h,second,'@33;M2;')" « endl; octfile « "loglog(h,third,'@33;M3;')" « endl; octfile « "loglog(h,fourth,'033;M4;')" « endl; octfile « "loglogQi,fifth,'@33;M5;')" « endl; octfile << "pause" « endl; string callstr="octave -q " + dotm; system(callstr.c_str()); // calls octave to create graphs void OctaveFDPlot(dvec h, dvec gamma, dvec regstats, dvec f_bounds){ // plot fractal regression on variogram using Octave string fracline = runref + "fin.out"; ofstream fln(fracline.c_str()); for (int count = (int) f.bounds[0]; count <= (int) f_bounds[l]; count++) { float x, y; x = h[count]; y = pow(10,(regstats[0]*loglO(x))+regstats[l]); fin « x « " " « y « endl; } string dotm = runref + "f.m"; ofstream octfile(dotm.c_str()); octfile « "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; "fifth = mmt(:,6);" « endl; "x =" « runref « "fln(:,l);" « endl; "y =" « runref « "fln(:,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 off" « endl; 'pause" « endl; // octfile << "gset term postscript landscape" « endl; // octfile « "gset output '" « runref « ".ps"' « endl; // octfile « "replot" « endl; // octfile « "gset term xll" « endl; string callstr="octave -q " + dotm; system(callstr.c_str()); // calls octave to create graphs octfile « octfile « octfile « octfile « octfile « octfile « octfile « octfile « octfile « octfile « octfile « octfile « octfile « octfile « octfile << octfile « octfile « octfile « octfile « octfile « octfile « octfile « void WriteLog(const stringfe runref){ // write parameters to log file string dotlog = runref + ".log"; ofstream logfile(dotlog.c_str()); logfile « "fd.cc produced this data using the following parameters:" « endl; logfile << "h_start, h_stop, h_num, h_tol, anisodirn, adirn_tol, sample" « endl; logfile « 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 file string dotout = runref + ".out"; ofstream outfile(dotout.c_str()); for (int count = 0; count <= (h_num-l); count++) { outfile << h[count] « " " « gamma[count] « endl; } void WriteFrac(string runref, dvec lh, dvec lgamma, dvec regstats, dvec fracstats){ // write fractal stats to file string dotfrc = runref + ".frc"; ofstream frcfile(dotfrc.c_str()); for (int count = 0; count <= (h_num-l); count++) { frcfile « lh[count] « " " << lgamma[count] « endl; } frcfile << "Regression Statistics: slope=" « regstats[0] « " intercept=" « regstats[1] « endl; frcfile « "Fractal Statistics: D=" « fracstats[0] « " Gamma(0)=" « fracstats[1] « " Bounds " « fracstats [2] « " <= h <=" « fracstats[3] « endl; cout « "Regression Statistics: slope=" << regstats[0] « " intercept=" << regstats[1] « endl; cout « "Fractal Statistics: D=" « fracstats[0] « " Gamma(0)=" « fracstats[1] « " Bounds " « fracstats[2] « " <= h <=" « fracstats[3] << endl; plist readinputfile(const stringfe runrefM /*** open data file for input ***/ string dotin = runref + ".in"; ifstream infile(dotin.c_str()); if (Unfile) { cerr « "Cannot open: " « dotin « " for input. Aborting." << endl; exit(l); } else cout « "Reading Datafile: "; 122 /*** create list of points and read data file into this list ***/ string dataline; plist readlist; plist_it read_it = readlist.begin(); // starts traverse of readlist on // read_.it int size =0; // number of lines read while (getline(infile, dataline)) // read one dataline { size++; if (size == 1) /*do nothing (to eat header line) */; else { float x, y, z, rubbish; //temp vars to decode stream const char* buffer = dataline.c_str(); //dataline -> c-string for //streaming istrstream iss(buffer,80); iss » rubbish >> z >> x >> y; // rubbish catches ARC rec num Point2D temp_point(x,y,z); readlist.insert(read_it, temp_point); // cout << temp_point « endl; // useful for testing... } } /*** Perform niceties and return linked list of points ***/ cout « size-1 « " data lines read." « endl; return readlist; > 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 lstart = (double) loglO(start); double lstop = (double) loglO(stop); double incr = (lstop - lstart) / ((double)nsteps - 1); for (double lgh = lstart; lgh < lstop; lgh += incr) { double temp = pow(10,lgh); logbins.push_back(temp); } logbins.push_back((double)stop); //cheat for last value to cover 123 // any roundoff error cout « "Done." « endl; return logbins; } dvec askbounds(const bool& iact){ // ask user fo bounds of fractal regression // choose entire range if running in command line mode dvec bounds; if (iact == true) { double first, last = 0; cout « "Enter first point for fractal estimate: "; cin » first; cout « "Enter last point for fractal estimate: "; cin » last; bounds.push_back(first); bounds.push_back(last); } else { bounds.push_back(0); bounds.push_back(h_num-l); } return bounds; > dvec regress(const dvec& y, const dvecfe x, const dvec& range){ // regress y on x over range, return slope, intercept in a dvec dvec xf; dvec yf; copy(xf, x, (int) range[0], (int) range[l]); copy(yf, y, (int) range[0], (int) range[l]); dvec xfsq = squares(xf); dvec yfsq = squares(yf); dvec xyfprod = products(xf, yf); int nx = xf .sizeO ; int ny = yf.sizeO; if (nx!=ny) exit(0); double sumxf = sum(xf); // we'll use this a lot double slope = (ny*sum(xyfprod)-sumxf*sum(yf))/ (nx*sum(xfsq)-sumxf*sumxf); double intercept = 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 fractal statistics 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 = strlen(s); if (length > 10) { cout « "str2int failed: number too large" << endl; exit(l); } for (int counter = 0; counter < length; counter++) { if (s[counter]=='-') { sign = -1; place++; } else { int number = int(s[counter])-48; if ((number >= 0) && (number <= 9)) { int power = (int)pow(10,(length - place)); value += number * power; place++; } else { cout « "str2int: non-integer argument." « endl; exit(l); } 125 } value *= sign; return value; } plist purge(plistfe input, float proportion) { // create subsample of points based on user specified proportion // if running in batch mode, always choose 2000 points (if available) plist purged = input; int listsize = input.size(); if (proportion == -1) proportion = ((float)2000 / (float)listsize); if (proportion > 1) proportion = 1; cout << "Original datafile has: " << listsize << " 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); plist_it del_it = 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 details of this algorithm // (write out all moments for plotting (added 25/05/00)) int hsize = hvalues.size(); dvec N; dvec first; dvec second; dvec third; dvec fourth; dvec fifth; dvec gamma; 126 // fill N, sum and gamma with zeros for (int count = 0; count < hsize; count++) { N.push_back(0); first.push_back(0); second.push_back(0); third.push_back(0); fourth.push_back(0); fifth.push_back(0); gamma.push_back(0); } // attempt to fit each point into an h-bin for (plist_it tll=inputdata.begin(); til != inputdata.end(); ++tll) { plist_it temp = inputdata.beginO; // allow advance by one next line temp = til; temp ++; for (plist_it tl2=temp; tl2 != inputdata.end(); ++tl2) { // Finds distance and angle between point pairs. double d_ij = dist(*tll,*tl2); double theta_ij = angle(*tll,*tl2); double zdiff_ij = zdiff(*tll,*tl2); if (anisodirn != -1) { cout « "I have not yet learned to make anisotropic variograms!" « endl; exit(l); > else { int location = fit(d_ij,hvalues); if (location != 9999) { first[location] += zdiff_ij; second[location] += pow(zdiff_ij, 2); third[location] += pow(zdiff_ij, 3); fourth[location] += pow(zdiff_ij, 4); fifth[location] += pow(zdiff_ij, 5); N[location]++; } } } } for (int count = 0; count < hsize; count++) { if (N[count] > 100) { 127 gamma[count] = first[count]/(N[count]); second[count] /= N[count] * 2; third[count] /= N[count] * 3; fourth[count] /= N[count] * 4; fifth[count] /= N[count] * 5; } else { cout « "warning: only " « N[count] « " input points in h_bin « count « "." « endl; } } string outfile = "mmt." + runref; ofstream out(outfile.c_str()); for (int count = 0; count < hsize; count++) { out « hvalues [count] << " " « gamma[count] « " " « second[count] « " " « third[count] « " " « fourth[count] « " " « fifth[count] « endl; } return gamma; > int fit(const doublefe d_ij, const dvecft hvalues) { // return the array index in hvalues into which a point with // separation distance d_ij should be placed, int location; bool housed = false; int count = 0; while (housed == false) { double h = hvalues[count]; if ((d_ij <= (h + h_tol)) && (d_ij >= (h - h_tol))) { location = count; housed = true; } if (count == hvalues.size()) { location = 9999; housed = true; } count++; } return location; 128 } dvec loglO (dvec v) { // return loglO of a dvec (in one go) // same as octave/matlab dvec logs; int v_size = v.sizeO; for (int i = 0; i < v_size; i++) { logs.push_back(loglO(v[i])); // take loglO of each element > return logs; } /lie**************************************/ /* class defn for Point2D */ /* generic functions on plists */ 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 operator»(istream& istr, 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 printO { 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& ostr, const Point2D& p) { if ((p._xcoord == -9999) II (p._ycoord — -9999)) ostr << "Position Undefined"; else if (p._zvalue == -9999) ostr « "(" « p._xcoord « "," << p._ycoord « ") = no value"; else ostr « "(" « p._xcoord « "," « p._ycoord « ") = " « p._zvalue; return ostr; } istreamfe operator»(istream& istr, Point2D& p) { // format for input is x,y,z char bin; // used to eat ( and , istr » p._xcoord; if (istr.peek()==',') istr >> bin » p._ycoord; if (istr.peek()==',') istr >> bin » p._zvalue; return istr; } double dist(const Point2D& pi, const Point2D& p2){ // Finds straight line distance from pi to p2 if ((pl._xcoord == -9999) II (pl._ycoord == -9999) II (p2._xcoord == -9999) || (p2._ycoord == -9999)) return -9999; else { double xdtemp = p2._xcoord - pl._xcoord; double ydtemp = p2._ycoord - pl._ycoord; xdtemp *= xdtemp; ydtemp *= ydtemp; double dtemp = sqrt(xdtemp+ydtemp); return dtemp; } } double angle(const Point2D& pi, const Point2D& p2){ if ((pl._xcoord == -9999) II (pl._ycoord == -9999) I I (p2._xcoord == -9999) || (p2._ycoord == -9999)) return -9999; else { // Finds 360 degree bearing from North of pt 2 from pt 1 double xdtemp = p2._xcoord-pl._xcoord; double ydtemp = p2._ycoord-pl._ycoord; double theta=0; 130 // Correct for quadrant to get radian bearing from north if (xdtemp > 0) theta = 1.57079632679 - atan(ydtemp/xdtemp); // = pi/2 else if (xdtemp < 0) theta = 4.71238898038 -atan(ydtemp/xdtemp); // = 3 pi by 2 // Convert to degrees theta *= 57.2957795131; // 180 / pi (rad -> deg) return theta; } } double zsq(const Point2D& pi, const Point2D& p2){ if ((pl._xcoord == -9999) II (pl._ycoord == -9999) I I (p2._xcoord == -9999) II (p2._ycoord == -9999)) return -9999; else { // find square of distance between two points double zdiff = (pi._zvalue)-(p2._zvalue); double zsq = zdiff*zdiff; return zsq; } } double zdiff(const Point2D& pi, const Point2D& p2){ // return absolute value of z difference if ((pl._xcoord == -9999) II (pi._ycoord == -9999) I I (p2._xcoord == -9999) II (p2._ycoord == -9999)) return -9999; else { double zdiff = (pi._zvalue)-(p2._zvalue); if (zdiff < 0) zdiff *= -1; return zdiff; } } 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) Jv = /v{l*+A- (B'3) Where Vj is a velocity vector of the boundary of the landscape element. Suppose that: § =f**r.-fp-#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: Il = 0.0139i<^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 Q035 S1'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 Q076, and d « 0.92 Q05. A correction to the slope exponent is also made, following Martin (1998), to account for configurational stability in steep channels, giving: Azb (per year) = 40.32 Q034 5 2 2 (B.10) 134 B.3 erode2d.cc /* Landscape Evolution 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 strings // length of dem grid size // length of model run (years) // nodes on each side of square model domain // number of timesteps // maximum height (used in generating // artificial topography) // Support for custom data types typedef double Matrix[NSIDE][NSIDE]; typedef vector<double> dvec; ofstream slidefileC'lscape.sli"); ofstream slidedataC'lscape.sld"); ofstream sedtranspC'lscape.sed"); Random billy; Matrix fillnumber; // keeps record of sinks in flow-routing algorithm. Matrix chaimelch; // keeps record of net aggradation/degradation in 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}; // routing const double rt2 = pow(2, 0.5); const double Pi = 3.1415927; // Global Parameters double K, hxO, hscale; double r; double crit_prob; double alpha = 35*(Pi/180); double k_fluvial; double hftrans; double Umax; int slide_test_int = 1; string inputfile; string upliftmodel; double rho_rock = 2.65; double rho_sed = 1.60; bool diagnose = false; // Functions void title(); void getparamsO; void expnonlinear(Matrix initial, Matrix final, Matrix production, Matrix uplift, float start, float stop, float dt, string kmodel); void renewuplift(Matrix uplift, float t, string upliftmodel, double Umax); void dnet(Matrix landscape, Matrix slope, Matrix rivernet, Matrix flow, double k_fluvial, float t); void fillsinks(Matrix landscape, Matrix flow); void flowroute(Matrix inaltgrid, Matrix slope, Matrix outdirgrid); void accumulate(Matrix indirgrid, Matrix outareagrid); void Bagnold(Matrix acc, Matrix slope, Matrix transport, double k_fluvial); void fluxsep(Matrix riverchange, Matrix riverflux, Matrix flow, float t); void dsbl(Matrix landscape, Matrix slidescape, Matrix sliderecord, Matrix slope, Matrix flow, float t); void export(Matrix output, bool display, int id); void matrix_copy(Matrix old, Matrix newmatrix); void matrix_add(Matrix first, Matrix second); void threshold(Matrix in, Matrix out, double cutoff); void filKMatrix input, string file); // non-linear diffusion parameters // stability parameter for solution // critical probability for dsbl slides // repose angle for talus (radian) // proportionality const, for fluvial transport // upslope area where channel initiation occurs // maximum uplift rate // minimum interval between dsbl slides // source of initial surface // uplift regime // density of rock // density of sediment // turns on verbose output for diagnostic work 136 void writelogO; bool is_a_zero(Matrix grid); bool is_even(int testint); bool is_mult(double a, double b); int drain(Matrix ingrid, int xi, int yi); int dsxpos(int oldx, double flowdir); int dsypos(int oldy, double flowdir); double volume(Matrix landscape); /****************/ /* MAIN PROGRAM */ int main(int argc, char* argv[]) { titleO ; getparams () ; Matrix landscape, riverchange, slope, flow, sliderecord; fill(landscape, inputfile); filKsliderecord, "nodata"); fill(fillnumber, "zeros"); fill(channelch, "zeros"); export(landscape, diagnose, 0); float dt = TSTOP/NT; r = (K * dt) / SIDELENGTH; cout « "r = " « r « endl; // cout « "Stability parameter = " « (D * dt) / (NSIDE*NSIDE) « endl int id = 1; writelogO ; ofstream vfile("lscape.vol"); for (float t = 0; t <= TSTOP-dt; t+=dt) { if (t==0) dnet(landscape, slope, riverchange, flow, k_fluvial, t); if (is_mult(t, slide_test_int)) { Matrix slidescape; dnet(landscape, slope, riverchange, flow, k_fluvial, -1); // before sliding (to refresh) dsbl(landscape, slidescape, sliderecord, slope, flow, t); matrix_copy(slidescape, landscape); dnet(landscape, slope, riverchange, flow, k_fluvial, t); // after sliding } 137 Matrix output; vfile << t << " " « volume(landscape) « endl; Matrix uplift; renewuplift(uplift, t, upliftmodel, Umax); expnonlinear(landscape, output, riverchange, uplift, t, t+dt, 1, "constant"); if (is_mult(t, TSTOP/40)) { cout « "time = " « t+(TSTOP/40) « " yrs. " « "Writing lscape" << id « endl; export(output, diagnose, id); export(sliderecord, diagnose, id+50); id++; } matrix_copy(output, landscape); // feed this iteration's output back as input for next > // call Octave script to display animated 3D projection of landscape export(sliderecord, diagnose, 99); systemC'mv lscape99.out lscape.slh"); export(fillnumber, diagnose, 98); systemC'mv lscape98.out lscape.snk"); export(channelch, diagnose, 97); systemC'mv lscape97.out lscape.cnl"); system("./vscape 2> /dev/null"); } // TECTONICS void renewuplift(Matrix uplift, float t, string upliftmodel, double Umax) { if (upliftmodel == "none") { // NO UPLIFT for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { uplift [xpos] [ypos] = 0; > } } if (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 } } if (upliftmodel == "thrust") { // THRUST for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { if (xpos < (NSIDE / 2)) uplift[xpos] [ypos] = Umax; else uplift[xpos][ypos] =0; > } } if (upliftmodel == "pluton") { // PLUTON double step_diff = (Umax / NSIDE) * 2; for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { if (xpos < (NSIDE / 2)) uplift [xpos][ypos] = xpos * step_diff; else uplift[xpos] [ypos] = (2 * Umax) - (xpos * step_diff); } } } if (upliftmodel == "rebound") { // ISOSTATIC REBOUND WITH (LINEAR) DECREASING RATE for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { uplift[xpos] [ypos] = Umax * ((TST0P - t) / TSTOP); } } } // FLUVIAL PROCESSES void dnet(Matrix landscape, Matrix slope, Matrix riverchange, Matrix flow, double k_fluvial, float t) { * II Compute the drainage network from elevations using D8 algorithm // each cell in Q is the annual average discharge at that point, // based on the area, integrated upslope, above the point. // Each cell in rivernet is the amount removed by the river (i.e., // negative production of sediment) 139 Matrix acc, riverflux; flowroute(landscape, slope, flow); bool tampered = false; double before_tamper = volume(landscape); while (is_a_zero(flow)) { tampered = true; fillsinks(landscape, flow); // fills small sinks only flowroute(landscape, slope, flow); } if (tampered) { // print a warning message if any closed pits have been filled // as a results of enforcing the flow-routing rules // never happens with usual erosion rules because no // process is capable of producing a pit. double after_tamper = volume(landscape); if ((after.tamper - before_tamper) > (after_tamper * 0.01)) { cout « "lscape (dnet): WARNING! elevation field modified by > 17. to route flow" « endl; } } accumulate(flow, acc); Bagnold(acc, slope, riverflux, k_fluvial); fluxsep(riverchange, riverflux, flow, t); if (t != -1) matrix_add(channelch, riverchange); if (t = TSTOP) { Matrix network; filKnetwork, "nodata"); threshold(acc, network, k_fluvial/(NSIDE * NSIDE)); export(network, diagnose, 96); system("mv lscape96.out lscape.net"); } } void fillsinks(Matrix landscape, Matrix flow) { ./* for every cell, determine its lowest neighbour if the lowest neighbour is higher than the original cell 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 = false; for (int ypos = 2; ypos < NSIDE-2; ypos++) { // skip the edges // where there is always a sink for (int xpos = 2; xpos < NSIDE-2; xpos++) { double low_height = 1000000; // a very large number bool is_sink = true; for (int neighbour = 0; neighbour < 8; neighbour++) { int xj = xpos + xlookup[neighbour]; int yj = ypos + ylookup[neighbour]; if (landscape[xpos] [ypos] < landscape[xj][yj]) { if (landscape[xj][yj] < low_height) low_height = landscape[xj][yj]; } else is_sink = false; } if (is_sink) { landscape[xpos] [ypos] = low_height; sinks = true; > } } } // For all cells with undefined flow direction // find the global pour point double global_min = 10000; // arbitrarily large value for (int ypos = 2; ypos < NSIDE - 2; ypos++) { for (int xpos = 2; xpos < NSIDE - 2; xpos++) { if (flow[xpos][ypos] == 0) { bool has_higher = false; for (int neighbour = 0; neighbour < 8; neighbour++) { int xj = xpos + xlookup[neighbour]; int yj = ypos + ylookup[neighbour]; if ((landscape[xpos][ypos] != landscape[xj][yj]) kk (landscape[xj][yj] < global_min)) 141 global_min = landscape[xj][yj]; > } } } // assign global_min to all cells with undefined flow direction for (int ypos = 2; ypos < NSIDE - 2; ypos++) { for (int xpos = 2; xpos < NSIDE - 2; xpos++) { if (flow[xpos][ypos] == 0) { landscape[xpos] [ypos] = global_min; fillnumber[xpos][ypos]++; } } } } void flowroute(Matrix inaltgrid, Matrix slope, Matrix outdirgrid) { // kernel of D8 flow routing algorithm. // Input must be a depressionless DEM; if it is not, depressions // are returned with zero drainage direction // Make edges flow outward; assumes that we are well inside area of interest // a two-cell strip is 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 interior node compare its elevation with those of its neighbours double diag_div = pow(2 * pow(SIDELENGTH, 2), 0.5); // horizontal // length for points on diagonal double diag_corr = SIDELENGTH / diag.div; // correction for diagonals for (int ypos = 2; ypos < NSIDE-2; ypos++) { for (int xpos = 2; xpos < NSIDE-2; xpos++) { double cellheight = inaltgrid[xpos][ypos]; double local_max = 0; int local_flow = 0; for (int neighbour = 0; neighbour < 8; neighbour++) { int xj = xpos + xlookup[neighbour]; int yj = ypos + ylookup[neighbour]; double n_slope = (cellheight - inaltgrid[xj][yj]) / SIDELENGTH; if (is_even(neighbour)) { // if on a diagonal, divide by sqrt(2) n_slope = n_slope * diag_corr ; // correct for diagonals } if (n_slope > local_max) { local_flow = (int) pow(2, neighbour); local_max = n_slope; } else if (n_slope == local_max) { local_flow += (int) pow(2, neighbour); } } if (local_max == 0) local_flow = 0; slope[xpos][ypos] = local_max; outdirgrid[xpos] [ypos] = (double) local_flow; } > // END OF FIRST PASS // Choose direction randomly for cells that could have several for (int ypos = 2; ypos < NSIDE-2; ypos++) { for (int xpos = 2; xpos < NSIDE-2; xpos++) { 143 int local_flow = (int) outdirgrid[xpos][ypos]; if (local_flow == 0) ; // do nothing: wait for next sweep // assign central cell if a whole side is equally lower else if (local_flow == 7) local_flow = 2; else if (local_flow == 28) local_flow = 8; else if (local_flow == 112) local_flow = 32; else if (local_flow == 193) local_flow = 128; else { // unpack flow directions and make a random choice if // several are equally valid dvec choices; int pips = local_flow; if (pips >= 128) { choices.push_back(128); pips -= 128; } if (pips >= 64) { choices.push_back(64); pips -= 64; > if (pips >= 32) { choices.push_back(32); pips -= 32; } if (pips >= 16) { choices.push_back(16); pips -= 16; } if (pips >= 8) { choices.push_back(8); pips -= 8; } if (pips >= 4) { choices.push_back(4); pips -= 4; } if (pips >= 2) { choices.push_back(2); pips -= 2; } if (pips >= 1) { choices.push_back(l); 144 pips -= 1; } if (choices.size() == 1) local_flow = (int) choices[0]; else { int index = billy.integer(choices.size()-l, 0); local_flow = (int) choices[index]; } > outdirgrid[xpos][ypos] = local_flow; } } // END OF SECOND PASS // Resolve flow for flat nodes by iteration // flow may be to any node that has an assigned direction provided // that node does not flow back to the original. // the while loop may get stuck if there are extensive areas of // perfect flatness if it does, tail_chaser bails it out and dnet // will call fill_sinks. Few landscapes are perfectly flat // so this is not usually a problem. int tail_chaser = 0; while ((is_a_zero(outdirgrid)) kk (tail_chaser < NSIDE*NSIDE)) { tail_chaser++; for (int ypos = 2; ypos < NSIDE-2; ypos++) { for (int xpos = 2; xpos < NSIDE-2; xpos++) { if (outdirgrid[xpos][ypos] == 0) { int flat.flow = 0; dvec flat_choices; for (int neighbour = 0; neighbour .< 8; neighbour++) { int xj = xpos + xlookup[neighbour]; int yj = ypos + ylookup[neighbour]; if ((inaltgrid[xpos][ypos] - inaltgrid[xj][yj] == 0) kk (outdirgrid[xj][yj] != 0) kk (outdirgrid[xj][yj] != flowback[neighbour])) { double dim = pow(2, neighbour); flat_choices.push_back(dirn); } } if (flat_choices.size() == 0) flat_flow =0; //no change else if (flat_choices. sizeO == 1) flat_flow = (int) flat_choices[0]; 145 else { int index = billy.integer(flat_choices.size()-l, 0); flat_flow = (int) flat_choices[index]; > outdirgrid[xpos][ypos] = flat_flow; } } } // END OF THIRD PASS } } int drain(Matrix ingrid, int xi, int yi) { // find upslope area of point xi, yi on the // flowdirection grid ingrid // (recursive) int cellcount = 0; for (int neighbour = 0; neighbour < 8; neighbour++) { int xj = xi + xlookup[neighbour]; int yj = yi + ylookup[neighbour] ; if ((xj > 0) && (xj < NSIDE) && (yj > 0) && (yj < NSIDE)) { if (ingrid[xj][yj] == flowback[neighbour]) { int upstream = drain(ingrid, xj, yj); cellcount += 1 + upstream; } > } return cellcount; } void accumulate(Matrix indirgrid, Matrix outareagrid) { // given a grid of flow-directions, return upslope area // of each cell for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { outareagrid[xpos][ypos] = drain(indirgrid, xpos, ypos); } } > void Bagnold(Matrix acc, Matrix slope, Matrix riverflux, double k_fluvial) { // convert upslope area into (bedload) sediment flux 146 // define parameters that give regional context double k_fl = k_fluvial; 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_fl * pow(upslparea / 1000000, 0.68); if (upslparea >= c_chan_maint) { if (slope[xpos][ypos] < 0.01) { // compute sediment transport rate using // normal Bagnold relation // (Martin and Church, 1999; Knighton, 1998) double pointdiff = 6.44 * pow(q, 0.35) * pow(slope[xpos] [ypos], 1.5); riverflux[xpos][ypos] = pointdiff; } else if (slope[xpos] [ypos] < 0.15) { // compute sediment transport rate using // headward Bagnold relation (Martin, 1998; Day, 1969) double pointdiff = 40.32 * pow(q, 0.34) * pow(slope [xpos] [ypos], 2.22); riverflux[xpos][ypos] = pointdiff; } else { riverflux[xpos][ypos] = 0; } } else riverflux [xpos] [ypos] = 0; // no fluvial processes above channel head } } > void fluxsep(Matrix riverchange, Matrix riverflux, Matrix flow, float t) { 147 // use predicted total flux at each point and input flux from upstream // neighbours to compute amount of erosion or deposition 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]; if ((xj > 0) kk (xj < NSIDE) kk (yj > 0) kk (yj < NSIDE)) { if (flow[xj] [yj] == flowback[neighbour]) { inflow += riverflux[xj][yj]; } } } double netchange = inflow - riverflux[xpos][ypos]; riverchange[xpos] [ypos] = netchange / (SIDELENGTH * SIDELENGTH); } } // write total bedload sediment flux to sedtransp, as a row of // NSIDE*NSIDE+1 columns // (first column is time) // ... records one row every time dnet is run from main.. one row // each 100 years if (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 initial, Matrix final, Matrix production, Matrix uplift, float start, float stop, float dt, string kmodel) { // Explicit 2-D solver for non-liear landscape diffusion equation // with initial surface, from start to stop with interval dt // using kmodel specified and including terms for // production of material (tectonics, rivers etc.) (spatially dependent) // with results written to final (which must exist when fn is called) Matrix U, V; matrix_copy(initial,U); fill(V, "zeros"); int toggle = 1; double Kl = K / 2; for (float time = start; time <= stop; time+=dt) { if (toggle == 1) { // set boundary using method of images for (int ypos = 1; ypos < NSIDE-1; ypos++) { // left and right 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 derivatives 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 flux using non-linear function double slpA = XA + YA double slpB = XB + YB double slpC = XC + YC double slpD = XD + YD double QA = Kl * (1 + tanh( (abs( slpA ) - hxO) / hscale)) * (slpA) double QB = Kl * (1 + tanh( (abs( slpB ) - hxO) / hscale)) * (slpB) double QC = Kl * (1 + tanh( (abs( slpC ) - hxO) / hscale)) * (slpC) double QD = Kl * (1 + tanh( (abs( slpD ) - hxO) / hscale)) * (slpD) // compute elevation change double diffuse = (dt/SIDELENGTH)*(QA-QB+QC-QD); if ((slpA > 1.6) || (slpB > 1.6) II 150 (slpC > 1.6) II (slpD > 1.6)) { // turn off diffusion on slopes > " 60 degrees diffuse = 0; } double change = diffuse + (production[xpos][ypos] * dt) + (uplift[xpos] [ypos] * dt); V[xpos] [ypos] = U[xpos] [ypos] + change; } } } if (toggle == -1) { // set boundary using method of images for (int ypos = 1; ypos < NSIDE-1; ypos++) { // left and right 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 derivatives 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 flux using non-linear function double slpA double slpB double slpC double slpD = XA + YA = XB + YB = XC + YC = XD + YD double QA = Kl * (1 + tanh( (abs( slpA ) - hxO) / hscale)) * (slpA) double QB = Kl * (1 + tanh( (abs( slpB ) - hxO) / hscale)) * (slpB) double QC = Kl * (1 + tanh( (abs( slpC ) - hxO) / hscale)) * (slpC) double QD = Kl * (1 + tanh( (abs( slpD ) - hxO) / hscale)) * (slpD) // compute elevation change double diffuse = (dt/SIDELENGTH)*(QA-QB+QC-QD); if ((slpA > 1.6) II (slpB > 1.6) || (slpC > 1.6) II (slpD > 1.6)) { // turn off diffusion 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; } if (toggle == 1) matrix_copy(U, final); else if (toggle == -1) matrix_copy(V, final); else { cout « "lscape (expsolv): internal error" << endl; exit(l); } > // DEEP-SEATED BEDROCK LANDSLIDES void dsbl(Matrix landscape, Matrix slidescape, Matrix sliderecord, Matrix slope, Matrix flow, float t) { // modify landscape to create slidescape based on the occurrence of // deep-seated bedrock landslides; slides occur at random. // most likely sites are identified using model nodes where gradient // elevation and length of slope are highest // put time in slidefile slidefile « t « endl; // Generate sliding probabilities Matrix slideprob; matrix_copy(landscape, slideprob); matrix_copy(landscape, slidescape); for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { slideprob[xpos] [ypos] *= slope [xpos] [ypos]; } } // Make a list of probabilities 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); if ((xpos <= 5) || (xpos >= NSIDE-5) II (ypos <= 5) || (ypos >= NSIDE-5)) { // zero probability of sliding near edge nodelist[xpos*NSIDE+ypos].elev(O); } else { nodelist[xpos*NSIDE+ypos].elev(slideprob[xpos][ypos]); } } } // Sort the list of probabilities sort(nodelist, (NSIDE*NSIDE)); // Choose those probabilities greater than a critical threshold int limit = -1; bool steep_end = false; for (int probcount = 0; probcount < (NSIDE*NSIDE); probcount++) if (nodelist[probcount] <= crit_prob) limit++; } // cout « limit « "," « (NSIDE*NSIDE)-1 « endl; if (limit == (NSIDE*NSIDE)-1) { // bail out here if we arent having a slide slidefile << "no slide occurred" « endl; } else { //otherwise choose a slide location at random int xpos, ypos, cougar; bool selected = false; cougar = billy.integer((NSIDE*NSIDE)-1, limit); xpos = (int) nodelist[cougar].xcoordO; ypos = (int) nodelist[cougar].ycoordO; slidefile << nodelist[cougar] « endl; slidedata « nodelist[cougar] « " "; sliderecord[xpos] [ypos] = 1; // Calculate the maximum failure size at the chosen site double current_height = slidescape[xpos][ypos]; 154 slidefile << "current height = " « current_height << endl; int nextx = dsxpos(xpos, flow[xpos][ypos]); int nexty = dsypos(ypos, flow[xpos][ypos]); slidefile « "Next cell: " « nextx « "," « nexty « endl; double next_cell_height = slidescape[nextx][nexty]; slidefile « "Next cell height = " « next_cell_height << endl; double biggest_slide = current.height - next_cell_height - (SIDELENGTH * tan(alpha)); if (biggest_slide < 0) biggest_slide =0; //if slope not steep enough slidefile << "biggest slide = " « biggest_slide << endl; double ndmag =1; // all available relief fails slidefile « "weibull proportion = " « ndmag « endl; double slide_size = ndmag*biggest_slide; slidefile << "actual slide size = " << slide_size « endl; // Begin volume check during slide double pre_vol = volume(landscape); // reduce elevation at failure site slidescape[xpos][ypos] -= slide_size; // decide how much talus we have generated double stash = slide_size*(rho_rock / rho_sed); // allow for density change slidefile « "available regolith = " « stash << endl; double ttl_dep_vol = stash; // track ttl amount moved double ds_difference = SIDELENGTH*tan(alpha); // difference to maintain // angle of repose // move the material downslope until it runs out int dscount =0; // how many times have we moved downslope with debris int current_xpos = nextx; // start at node downslope of originating node int current_ypos = nexty; while (stash > 0) { slidefile << "Moved downstream" << endl; slidefile << "at location: " « 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 slidefile << "done backflow" « endl; if ((new_xpos <= 0) II (new_xpos >= NSIDE-1) I I (new_ypos <= 0) I I (new_ypos >= NSIDE-1)) { // in case we go off the edge slidefile « "lscape (dsbl): a freak slide has occurred: " << stash << " units of sediment have left 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; slidefile << "downslope height = " << ds_height « " ;ideal_height = << ideal_height « " ; reqdheight = " « required_height « endl; if (required_height >= 0) { // if we are in a depositional node slidefile << "Depositing: "; if (stash >= required_height) { // if there is enough material deposit to ideal_height slidescape[current_xpos][current_ypos] = ideal_height; stash -= required_height; slidefile << required_height << " units" << endl; } else { // otherwise just leave the entire amount that we have slidescape[current_xpos][current_ypos] += stash; slidefile « stash « " units" « endl; stash = 0; } } else { // otherwise the slide scours the downslope area slidefile « "Scouring "; slidescape[current_xpos][current.ypos] = ideal_height; double surplus = slidescape[current_xpos][current_ypos] -ideal_height; 156 slidefile « surplus << " units" << endl; surplus *= (rho_rock / rho_sed); stash += surplus; ttl_dep_vol += surplus; // add to total slide volume } current_xpos = new_xpos; current_ypos = new_ypos; > } double post_vol = volume(slidescape); // to check conservation of volume slidefile « "Total Volume: " « ttl_dep_vol « endl; slidefile « "Runout length: " << dscount « endl; slidefile « ttl_dep_vol « " " « dscount « endl; slidefile « "Pre-vol: " « pre_vol « " . Post-vol: " « post_vol « endl slidedata « ttl_dep_vol « " " « dscount « endl; cout « "DSBL @ (" « xpos « "," « ypos « "); Vol = " « ttl_dep_vol * (SIDELENGTH * SIDELENGTH) « " m~3 ; Runout = " « dscount * SIDELENGTH « " m." « endl; } } // INPUT-OUTPUT FUNCTIONS void export(Matrix output, bool display, int id) { // write a file containing output, named lscape<id>.out. // if the flag display is true, the output is also written to the screen. string prefix = "lscape"; int tens = (int)(id/10); int units = id - (10*tens); char tnum = tens + 48; char unum = units + 48; string filename = prefix + tnum + unum + ".out"; ofstream outfile(filename.c_str()); for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { outfile « output[xpos][ypos] « " "; 157 } outfile « endl; } if (display == true) { for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { cout « output[xpos][ypos] << " "; } cout « endl; } } } void fill(Matrix input, string file) { // Fill a matrix, input, with values taken from the file // if file is "random" or "zeros" then no file is read // and the matrix is filled with uniform random values or // zeros accordingly if (file == "random") { for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { input[xpos] [ypos] = billy.integer(HMAX,0); > } } else if (file == "zeros") { for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { input[xpos][ypos] = 0; } } } else if (file == "nodata") { for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { input[xpos] [ypos] = -9999; } } } else if (file == "sine") { 158 for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { if ((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 if (file == "comb") { for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { if ((xpos ==0) II (ypos ==0) || (xpos == NSIDE - 1) || (ypos == NSIDE - 1)) input[xpos][ypos] = 0; else i if (is_even(xpos)) input[xpos] [ypos] = (xpos*NSIDE + ypos); else input[xpos] [NSIDE - ypos - 1] = (xpos*NSIDE + ypos); } } } } else if (file == "ramp") { double step_diff = HMAX / (NSIDE / 2); for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { if (xpos < (NSIDE / 2)) input[xpos][ypos] = xpos * step_diff; else input[xpos] [ypos] = (2 * HMAX) - (xpos * step_diff); > } } else if (file == "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 it easier to enter like that above! input[xpos][ypos] = testinput[ypos][xpos]; // / 832; > } } else { ifstream infile(file.c_str()); if (Unfile) { cerr << "lscape (fill): cannot open input file: " << file << endl; exit(1); } cout << "lscape (fill): reading initial values from: " « file « endl; string dataline; for (int ypos = 0; ypos < NSIDE; ypos++) { { get1ine(inf ile, dat aline); double temp; const char* buffer = dataline.c_str(); //dataline -> c-string for streaming istrstream iss(buffer, 10000); for (int xpos = 0; xpos < NSIDE; xpos++) { iss » 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 ttl_elev = 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 titleO { cout « "erode2d: non-linear diffusion with threshold landsliding" « endl; cout « "Diffusive-advective landform evolution model" << endl; cout « "Written by Simon Dadson, April 2000" « endl; } void writelogO { // write log file ofstream logfileClscape.log"); logfile « "ER0DE2D Log File lscape.log: " « endl; logfile « "Initial surface = " << inputfile « endl; logfile « "SIDELENGTH = " « SIDELENGTH « endl; logfile « "TST0P = " « TST0P « endl; logfile « "NSIDE = " « NSIDE « endl; logfile « "NT = " « NT « endl; logfile « "HMAX = " « HMAX « endl; logfile « "D = " « D « endl; logfile « "K, hxO, hscale" « endl; logfile « K « " " « hxO « " " « hscale « endl; logfile « "stability parameter = " « r « endl; logfile « "Pcrit = " « crit_prob « endl; logfile « "alpha = " « alpha « endl; logfile « "K_fluvial = " « k_fluvial « endl; logfile « "hftrans = " « hftrans « endl; logfile « "uplift model = " « upliftmodel « endl; logfile « "max uplift = " « Umax << endl; logfile << "slide test interval = " « slide_test_int « endl; logfile << "rho_rock = " « rho_rock « endl; logfile << "rho_sed = " « rho_sed « endl; logfile << "diagnose = " << diagnose << endl; 161 int dsxpos(int oldx, double flowdir) { // find x coordinate of point downstream of oldx, given flowdir int newx = -1; for (int nc = 0; nc < 9; nc++) { if (flowdir == pow(2, nc)) { newx = oldx + xlookup[nc]; // back to front ! > } if (newx == -1) { cout « "lscape (dsxpos): error routing sediment" « endl; exit(1); } return (newx); > int dsypos(int oldy, double flowdir) { // find y coordinate of point downstream of oldy, given flowdir int newy = -1; for (int nc = 0; nc < 9; nc++) { if (flowdir == pow(2, nc)) { newy = oldy + ylookup[nc]; > } if (newy == -1) {. cout « "lscape (dsypos): error routing sediment" « endl; exit(l); } return (newy); > void matrix_copy(Matrix old, Matrix newmatrix) { // copy old into new // new must exist prior to the call 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 first, Matrix second) { // add second to first, result is in first, // second is unchanged // both must exist prior to the call for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { first[xpos][ypos]+=second[xpos][ypos] ; } } } bool is_even(int testint) { // return true if testint is even; false if odd float teststat = 0; teststat = testint°/,2; if (teststat != 0) return false; else if (teststat == 0) return true; else exit(1); bool is_mult(double a, double b) { // return true if a is a multiple of b double teststat = 0; teststat = int(a)'/,int(b) ; if (teststat != 0) return false; else if (teststat == 0) return true; else exit(1); void getparamsO { cout << "Enter name of input file: " ; cin » inputfile; int location; cout « "Use nonlinear diffusion rule for GVRD (1) or QCI (2) cin >> location; if (location == 1) { K = 0.12; hxO = 0.81; 163 hscale = 0.32; } else if (location == 2) { K = 0.43; hxO = 1.06; hscale = 0.19; } else { cout « "Invalid location" « endl; exit(1); } cout << "Enter uplift style (none, const, thrust, pluton, rebound): "; cin >> upliftmodel; cout « "Enter maximum uplift rate (metres per year): "; cin » Umax; cout « "Enter coefficient for fluvial transport (k): "; cin » k_fluvial; cout « "Enter hillslope-fluvial transition (in upslope area units (m"2)) "; cin » hftrans; cout << "Enter minimum dsbl recurrence interval (years): "; cin » slide_test_int; cout << "Enter critical slide probability (angle * elevation): "; cin » crit_prob; } bool is_a_zero(Matrix grid) { // return true if grid contains any zero, false otherwise, bool found_one = false; for (int ypos = 0; ypos < NSIDE; ypos++) { for (int xpos = 0; xpos < NSIDE; xpos++) { if (grid[xpos][ypos] == 0) found_one = true; } } return found_one; } void threshold(Matrix in, Matrix out, double cutoff) { // produce out from in, 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++) { if (in[xpos] [ypos] > cutoff) out[xpos][yp< else out[xpos][ypos] = -9999; 165 


Citation Scheme:


Usage Statistics

Country Views Downloads
United States 60 2
China 13 4
Russia 9 0
Japan 5 0
France 4 0
Netherlands 3 2
Romania 3 0
Italy 1 0
India 1 0
Philippines 1 0
City Views Downloads
Kansas City 46 0
Unknown 14 3
Hangzhou 8 0
Penza 7 0
Tokyo 5 0
Washington 4 0
Plano 3 0
Mountain View 3 0
Baotou 2 0
Ashburn 2 0
Beijing 2 0
Yekaterinburg 1 0
Redmond 1 0

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



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items