TERRAIN STABILITY ASSESSMENT USING LOGISTIC REGRESSION ANALYSIS FOR THE JAMIESON-ORCHID-ELBOW CREEKS SUBDRAINAGE, SEYMOUR RIVER BASIN, BRITISH COLUMBIA by Gyula Gulyas Diploma in Forest Engineering The University of Forestry and Wood Sciences, Sopron, Hungary, 1992 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (The Department of Forest Resources Management) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA JUNE 1995 © Gyula Gulyas, 1995 In presenting this degree at the thesis in University of partial fulfilment of of department this thesis for or by his or requirements British Columbia, I agree that the freely available for reference and study. I further copying the representatives. an advanced Library shall make it agree that permission for extensive scholarly purposes may be granted her for It is by the understood that head of copying my or publication of this thesis for financial gain shall not be allowed without my written permission. Department of The University of British Columbia Vancouver, Canada T7 DE-6 (2/88) ABSTRACT The purpose of this research project was to develop a procedure for terrain stability assessment by applying case-control sampling and multiple logistic regression analysis, widely used statistical techniques in biomedical research and in epidemiology. The idea of applying statistical methods used in epidemiology to terrain stability assessment was based on the observation that landslides, like some diseases, are rare phenomena. The implementation of a terrain stability assessment based on these statistical techniques was expected to help understand the causeeffect relationships between landsliding and various terrain attributes. In contrast to the currently used approaches, the study procedure provided a quantitative tool to assess the risk of landsliding and to define the most important terrain attributes that contribute to soil mass movements. A case-control study of 20x20 m grid cells with average slope greater than 10 degrees was conducted on the Jamieson-Orchid-Elbow Creeks subdrainage of the Seymour River Basin, British Columbia. All of the 101 landslide cases were compared with 264 control grid cells. Multi-way cross classification tables were constructed to study the relationship between landsliding and several terrain attributes. A possible interaction between slope angle and the drainage condition of the soil was detected. A logistic regression analysis was then performed within a Geographic Information System (GIS) environment to develop a landslide risk model for the Jamieson-Orchid-Elbow Creeks study area. A landslide risk matrix was then constructed based on the landslide risk model. It was found that sites located in the transient snow zone, with slope angle greater than 55 degrees, on bedrock outcrop surficial material type and on shallow soil have the greatest risk of experiencing rapid, shallow soil mass movements. It was also found that holding all the other variables constant, slope angle had the greatest effect on the magnitude of landslide risk. Based on the data, sites with very steep slopes (over 55 degrees) have, on the average, five times the chance of experiencing a landslide event relative to sites with gentle slopes (10-25 degrees). The landslide risk matrix was used to create landslide risk categories. The spatial distribution of landslide risk, categorized as very low, low, moderate, high and very high, is portrayed within 20-m square grid cells on the landslide risk map. The major advantage of using the landslide risk assessment of this study is that it provides the terrain mapper with quantitative information about the relative risk of landsliding. This information can be used as a tool in planning watershed management activities and in an overall risk assessment for a given geographic area. - in - T A B L E OF CONTENTS Pag ABSTRACT ii TABLE OF CONTENTS iv LIST OF TABLES viii LIST OF FIGURES xiii ACKNOWLEDGEMENTS xiv INTRODUCTION 1 1. LITERATURE REVIEW ON TERRAIN STABILITY ASSESSMENT . . . 4 1.1. General problems 5 1.1.1. Complexity 5 1.1.2. Subjectivity 6 1.1.3. Reliability 6 1.1.4. Transferability 7 1.2. Classification of methods of terrain stability assessment 7 1.2.1. Classification by size of the study area 7 1.2.2. Classification by survey intensity level 8 1.2.3. Classification by analytical methods 9 1.2.3.1. Landslide inventories 9 1.2.3.2. Factor overlay methods 11 1.2.3.3. Subjective geotechnical models 13 1.2.3.4. Factor of safety models 16 1.2.3.5. Statistical models 19 iv Page 1.2.3.5.1. Univariate statistical analysis 19 1.2.3.5.2. Multivariate statistical analysis 21 2. METHODS 25 2.1. Introduction 25 2.2. Study area 26 2.2.1. Location 26 2.2.2. Climate 26 2.2.3. Topography 29 2.2.4. History of soil mass movements 31 2.2.5. Vegetation 32 2.2.6. Forest harvesting and hydrology 34 2.3. Theoretical background 35 2.3.1. Case-control studies 35 2.3.1.1. Assessment of risk in case-control studies 37 2.3.1.2. Relative risk and the odds ratio 38 2.3.1.3. The case-control sampling scheme 40 2.3.1.4. Sample size calculation in a case-control study 43 2.3.1.5. Calculation of the power of the study 45 2.3.1.6. Confounding 45 2.3.1.7. Interaction 50 2.3.2. Logistic regression analysis 2.3.2.1. The logistic regression model 2.3.2.1.1. Testing for significance of the coefficients 2.3.3. Logistic regression and case-control studies 3. DATA ANALYSIS AND DISCUSSION 3.1. A case-control study on the study area -v- 54 55 61 62 67 67 Page 3.1.1. Stating the research questions 67 3.1.2. Variable selection 69 3.1.3. Data collection and generation 78 3.1.4. Study design 80 3.1.5. Case definition and selection 81 3.1.5.1. Sources of cases 82 3.1.6. Control definition and selection 83 3.1.6.1. Definition... 83 3.1.6.2. Eligibility criteria and sources 84 3.1.7. Calculation of the sample size 84 3.1.8. The power of the study 85 3.2. Analysis of the case-control study 86 3.2.1. Statistical analyses using 2x2 contingency tables 86 3.2.2. Multi-way cross classification tables 90 3.2.3. Stepwise multiple logistic regression 97 3.2.3.1. Model building strategy 97 3.2.3.1.1. Variable selection 98 3.2.3.1.2. Multiple logistic regression 102 3.2.3.1.3. Assessing the fit of the model 109 3.3. Landslide risk assessment for the study area 110 3.3.1. The logistic regression model 110 3.3.2. Landslideriskassessment 112 3.3.2.1. Choosing the base level of variables 112 3.3.2.2. The landslideriskmatrix 113 3.3.2.3. Landslideriskrating 115 4. SUMMARY AND CONCLUSIONS 123 - vi - Page GLOSSARY OF TERMS 128 REFERENCES 131 Appendices A FORM.FOR FORTRAN PROGRAM 142 B CATCH.FOR FORTRAN PROGRAM 146 CRAW DATA.. 150 - vii -r LIST OF TABLES Table 1.1. Survey intensity levels for terrain stability mapping 8 Table 1.2. Classification of terrain stability assessment by analytical method .9 Table 1.3. Advantages and disadvantages of landslide inventories 11 Table 1.4. Advantages and disadvantages of the factor overlay techniques 13 Table 1.5. Advantages and disadvantages of the subjective geotechnical models 16 Table 1.6. Advantages and disadvantages of the factor of safety models Table 1.7. Advantages and disadvantages of the univariate statistical method. 22 Table 1.8. Advantages and disadvantages of the multivariate statistical method 24 Table 2.1. Annual precipitation and intensity data for the study area 28 Table 2.2. The distribution of slope angle and soil depth in the study area 29 19 31 Table 2.3. Table 2.4. Overview of landslide studies in the study area 32 Table 2.5. Advantages and disadvantages of a case-control study 36 Table 2.6.A Frequency of landsliding among exposed and unexposed sites in the target population 38 - viii - Table 2.6.B Relationship of landsliding and slope angle (SL35) in the JamiesonOrchid-Elbow Creeks study area using SL35 as a dichotomized variable 39 Table 2.7.A Case-control study based on the target population of Table 2.6.B Expected sample outcome for a case-control study Table 2.7.B Numerical example of a case-control study 41 41 Table 2.8. The relationship between bedrock outcrop(SOIL2) and landsliding (SLIDE) in the study area 46 Table 2.9. Multi-way cross classification table for examining the relationship between SOIL2 and SLIDE stratified by TRANS 47 Table 2.10. Distribution of cases and controls by TRANS from Table 2.9 47 Table 2.11. The relationship between transient snow zone (TRANS) and landsliding (SLIDE) in the study area Table 2.12. Calculation of the Woolf s test statistic from Table 2.9 48 49 Table 2.13. Relationship between elevation zones (ELEVZ) and slope angle (SL35) with landsliding (SLIDE) 51 Table 2.14. Odds (in parentheses) and odds ratios for landsliding associated with various combinations of elevation zones and slope angle classes ... 52 Table 2.15. Relationship between drainage conditions (DR01) and slope angle (SL35) with landsliding (SLIDE) 53 Table 2.16. Odds (in parentheses) and odds ratios for landsliding associated with various combinations of drainage conditions and slope angle classes 53 Table 2.17. Frequency table of slope class by SLIDE 58 Table 2.18. Results of fitting a logistic regression model to the data in Table 2.17. 60 Table 2.19. Likelihood ratio test for 3 different models using LOGISTIC Version 3.11Ef of STATOOLS™ free software package 61 - ix - Table 2.20. Assessing the homogeneity of the odds ratios by using the likelihood ratio test 65 Table 3.1.A Variables "measured" at 20x20 m grid cells 74 Table 3.1.B Variables created by transformation 75 Table 3.2. Data types and sources used in the study 79 Table 3.3. Source of cases and method of data collection 83 Table 3.4. Relationship between slope angle and landsliding based on a dichotomized slope variable 85 Table 3.5. Relationship between slope angle (SL35) and landsliding (SLIDE) 87 Table 3.6. Relationship between transient snowzone (TRANS) and landsliding (SLIDE) 87 Table 3.7. Relationship between drainage conditions (DR01) and landsliding (SLIDE) 88 Table 3.8. Relationship between surficial material (SOIL2) and landsliding (SLIDE) 88 Relationship between soil depth (SD) and landsliding (SLIDE) 88 Table 3.9. Table 3.10. Relationship between planform curvature (PL01) and landsliding (SLIDE) 88 Table 3.11. Relationship between profile curvature (PR01) and landsliding (SLIDE) 89 Table 3.12. Relationship between aspect (ASP2) and landsliding (SLIDE) 89 Table 3.13. Relationship between biogeoclimatic zones (BGCZ) and landsliding (SLIDE) 89 Table 3.14. Relationship between local catchment area (CA30) and landsliding (SLIDE) 89 Table 3.15.A. Multi-way cross classification table for variables SLIDE, SL35, SD, SOIL2 and TRANS 92 Table 3.15.B. Stratified breakdown of the data in Table 3.15.A. with SL35 as the main risk factor 92 Table 3.16. Odds ratios for interpreting interaction and confounding 93 Table 3.17. Relation of landsliding (SLIDE) to slope angle (SL35) according to transient snowzone (TRANS), ignoring soil depth (SD) and surficial material (SOIL2) 94 Table 3.18. Relation of landsliding (SLIDE) to slope angle (SL35) according to soil depth (SD), ignoring transient snowzone (TRANS) and surficial material (SOIL2) 95 Table 3.19. Relation of landsliding (SLIDE) to slope angle (SL35) according to surficial material (SOIL2), ignoring transient snowzone (TRANS) and soil depth (SD) 96 Table 3.20. Contingency table for variables SOILS and SLIDE 98 Table 3.21. Contingency table for variables PROFC and SLIDE 99 Table 3.22. Contingency table for variables PLANC and SLIDE 99 Table 3.23. Contingency table for variables DISTUR and SLIDE 99 Table 3.24. Contingency table for variables ELEVZ and SLIDE 100 Table 3.25. Contingency table for variables SLCL and SLIDE 100 Table 3.26. Contingency table for variables ASPECT4 and SLIDE 100 Table 3.27. Contingency table for variables SDEPTH and SLIDE 101 Table 3.28. Univariate logistic regression models for variables considered in this study.. 103 Table 3.29A Results of applying stepwise variable selection using the maximum likelihood method to the study data in the development of model A105 Table 3.29B Results of applying stepwise variable selection using the maximum likelihood method to the study data in the development of model B 106 Table 3.30A Log-likelihood for Model A at each step and likelihood ratio test statistics (G), degrees of freedom (d.f.) and the p-values corresponding to the addition of the variable 106 - xi - Table 3.30B Log-likelihood for Model B at each step and likelihood ratio test statistics (G), degrees of freedom (d.f.) and the p-values corresponding to the addition of the variable 107 Table 3.31. Results of applying stepwise variable selection to interaction terms from the main effects model (B), using the maximum likelihood method 108 Table 3.32. Observed and expected frequencies within each decile of probability for each outcome (SLIDE = 0, SLIDE = 1) using the fitted logistic regression model 110 Table 3.33. Coefficients, standard errors (SE) and odds ratios of the final model Ill Table 3.34. Landslide risk matrix for the Jamieson-Orchid-Elbow Creeks study area 114 Table 3.35. Classification scheme for the study 116 Table 3.36. Cross classification table 117 Table 3.37. Calculation of areal risk based on the two methods 121 - xii - LIST OF FIGURES Figure 2.1. Location of the Jamieson-Orchid-Elbow Creeks study area 27 Figure 2.2. Slope categories in the Jamieson-Orchid-Elbow Creeks subdrainage, Seymour River Basin, coastal British Columbia 30 Figure 2.3. Landslide initiation points in the Jamieson-Orchid-Elbow Creeks subdrainage, Seymour River Basin 33 Figure 2.4. General shape of the logistic regression model 56 Figure 2.5. Scatterplot of SLIDE by SLOPE angle 59 Figure 2.6. Plot of the mean of SLIDE in each SLOPE class 59 Figure 3.1. Flow chart of the study procedure 68 Figure 3.2. Catchment area-slope relationship for small mountainous drainage basins along the Pacific Coast 71 Figure 3.3. Landslide risk assessment for the Jamieson-Orchid-Elbow Creeks subdrainage, Seymour River Basin, British Columbia 118 Figure 3.4. Terrain stability assessment by the Ecological Inventory Pilot Study (Acres International Ltd. 1993) 119 - xiii - ACKNOWLEDGEMENTS I wish to express my grateful thanks to Dr. Douglas L. Golding for his encouragement and guidance during the course of this work. It has been an interesting and rewarding learning experience. I would like to extend my sincere thanks to Dr. Valerie M. LeMay for her help, advice and suggestions concerning the statistical matters of this work. I would also like to express my appreciation to Mr. Lome Gilmour GIS specialist and Mr. Derek Bonin at the Greater Vancouver Water District for all their cooperation with this project. My special thanks go to Dr. Antal Kozak and the Sopron Alumni Fund for making my studies possible at the University of British Columbia. Most importantly, I want to thank my wife, Marianna for her unfailing support and faith in me. Her belief in my ability to overcome any problem has given me the courage to accomplish this work. It is her and our daughter Monika to whom I dedicate this thesis. - xiv - "A new epicycle of erosion has been initiated in this land of ours since the settlement, the like of which has not previously occurred since the glacial epoch of the Pleistocene." Bailey, R. W. 1937. A new epicycle of erosion. J. of Forestry. 35(11): p. 997. INTRODUCTION Resource management activities in many forest lands of the western United States, British Columbia and Alaska are severely restricted by unstable terrain that is susceptible to soil mass movements varying from surface creep to catastrophic landslide events. Removal of forest cover through harvesting and associated engineering operations such as road construction can accelerate mass erosion and cause significant degradation of water quality, fisheries and aesthetics. Forestry operations may severely alter slope geometry and the natural drainage system by improper road construction and may change the hydrologic conditions of the site by the removal of forest cover. In coastal British Columbia and the Pacific Northwest, soil mass movements typically occur on steep slopes where relatively shallow and cohesionless soils are underlain by impermeable bedrock or glacial till (O'Loughlin 1972). These conditions together with high intensity, long-duration rainfall predispose the area to a moderate to high natural level of rapid, shallow soil mass movements. In this study, the terms landslide, soil mass movement, mass wasting, and slope failure will be used interchangeably. These terms describe the displacement and downslope movement of soil, rock and organic material due to gravity and/or excess of water. There are several classification schemes in the literature. For terms used in this study the reader should refer to the "Glossary of terms" section. Detailed discussion on the various types of soil mass movement can be found in Swanston and Swanson (1976). The impacts of soil mass movement on the forest environment can be classified as follows: • damage to fish habitat (spawning sites); • sedimentation resulting in lower water quality and aquatic productivity; • damage to capital investment (property loss, bridges, roads); • site productivity loss and • visual impact (aesthetics). Landslides have their greatest impact on forest management by inflicting damage on anadromous fish habitat. Large landslides from clearcuts or roads may greatly alter spawning sites. Soil mass movement may increase turbidity and sedimentation in streams and rivers which, in turn, reduce water quality for aquatic organisms and human consumption. Large landslides may scour channel banks and destroy streamside vegetation, thereby exposing stream channels to direct sunlight that can increase stream water temperature. The effects of landslides on site productivity have been studied very little because landslides usually occupy less than one percent of the landscape, even in highly unstable clearcut areas. According to Miles (1983), trees grow slower on landslide tracks than on adjacent, undisturbed terrain. Landslides frequently damage roads, bridges, railways etc. and large landslides may even reach valleys where people live, causing injury and death along with property damage. Since it is almost impossible to control soil mass movement once it has started, it is important that forest managers understand causes, recognize areas of instability -2- and adjust management plans accordingly. For successful management of forest lands an accurate assessment of landslide hazard must support management decisions. The main objective of this thesis was to develop a procedure for terrain stability assessment based on case-control sampling and multiple logistic regression, widely used statistical techniques in biomedical research and in epidemiology. The idea of using statistical techniques from biomedical research originates from the observation that landslides, like some diseases, are rare phenomena, thus statistical procedures applied in epidemiology might be useful in terrain stability assessment. The procedure would provide an objective method to define the cause-effect relationship between slope stability and various terrain attributes. To fulfill this objective, the thesis was divided into four main parts: 1. Literature review of existing methods of terrain stability assessment. 2. Review of statistical techniques used in biomedical research to examine how these techniques can be implemented in a terrain stability assessment. 3. Development of a landslide risk model and the production of a landslide risk map for the Jamieson-Orchid-Elbow Creeks study area based on the procedure developed in this study. 4. Summary of the results and suggestions for further research. 1. LITERATURE REVIEW ON TERRAIN STABILITY ASSESSMENT In the past three decades, the problem of landsliding became a very important issue in forest resource management. There are three main reasons. First, valuable timber resources in flat valley bottoms have already been depleted. Second, advanced harvesting technologies (cable yarding, powerful machines) have been introduced providing access to timber on steeper slopes at higher elevation. Third, water quality and environmental concerns have become an integral part of forest resource management activities. A careful evaluation of landslide hazard within any areas proposed for forest development should be conducted to ensure that forestry operations are performed in a manner that will not increase the frequency and magnitude of soil mass movements (Forest Practices Code 1995). Understanding the relationship between landsliding and specific environmental conditions or terrain types can help forest managers predict how forestry operations affect slope stability. A number of different strategies for predicting slope instability, or modeling its likelihood, have been developed. They range from simple inventories of existing soil mass movements, multifactor mapping of slope, geology, soil, hydrology and other variables to complex multivariate statistical analyses. The major problem with any method of terrain stability assessment is to minimize the subjective elements. There is an urgent need for the development and application of more quantitative approaches that can provide sufficient information about the risk of landsliding in -4- the presence or absence of various landscape attributes. 1.1. General problems Of the many issues that must be addressed in developing a useful method of terrain stability assessment, four appear to be particularly important: 1. Complexity 2. Subjectivity 3. Reliability 4. Transferability 1.1.1. Complexity Landsliding is a complex process. The major challenge in the assessment of unstable terrain is to develop reliable techniques for evaluating slope stability hazards under any set of specific conditions. This requires the identification of cause-effect relationships among the numerous variables which are interrelated and multi-dimensioned. Natural variations in climate, slope, soil, geology, hydrology, vegetation and human activities are tremendous and collectively generate literally an infinite number of unique conditions. The large variation in the size, shape, topography and subsurface site conditions among different landslides suggests that each failure results from a unique combination of movement-promoting and movement-resisting forces (O'Loughlin 1972). Probably most soil mass movements result from a combination of many causes with one factor being finally dominant. This is why attempts to assess the relative stability of slopes and predict the location of future slope failures are often unsuccessful. In an ideal case variables included in the slope stability model should: 1. clearly separate stable and unstable sites on a quantitative basis 2. be easily extracted from aerial photographs or topographic maps to reduce costs 3. be geomorphologically meaningful Unfortunately in most cases the slope stability specialists are forced to oversimplify their models to stay within budget. One of the challenges to landslide specialists is to develop and refine slope stability analyses that take advantage of sophisticated techniques (e.g., multivariate modeling) while minimizing data requirements. 1.1.2. Subjectivity Subjectivity cannot be ruled out from the assessment of unstable terrain. Approaches using statistical tools try to reduce the subjective elements of the process. However, local observations and experience will always play a significant role in the assessment procedure. An important factor is the ability of the mapper to use his or her experience along with research findings from outside the area and adapt them to local conditions. 1.1.3. Reliability Many of the disputes concerning reliability revolve around questions of model appropriateness. The primary challenge is to choose models and supporting data that provide reasonable representation of the landsliding process. In examining the datamodel fit and the underlying reliability issue, it is important to understand the limitations and uncertainties of the available data and of the function of the models. Since data collection is the most expensive phase of terrain stability assessment, -6- efforts should be made to plan what kind of information will be needed and how the information will be collected and analyzed. It is very painful to realize that some valuable data are missing while unimportant and costly information has been collected that will not be used in the analysis. 1.1.4. Transferability The soil mechanics part of the landsliding process is quite well understood, but there is much room for improvement to understand this very complex process. A model that has been developed on one area should be used on other areas with extreme caution. The extrapolation of inventory results to other areas is risky because the causal factors may change from one location to the other. 1.2. Classification of methods of terrain stability assessment The different approaches to terrain stability assessment can be classified in many ways, depending on the purpose of the classification. 1.2.1. Classification by the size of the study area Slope stability assessments can be classified by the size of the study area. Local studies are usually conducted on relatively small areas (1-10 km ) They 2 include site specific landslide studies, large-scale landslide inventories (Dyrness 1967; Chatwin and Rollerson 1983; Carrara 1988). These studies usually incorporate very detailed information about the site characteristics. Regional studies are carried out on larger areas ranging from several hundred to several thousand square kilometers (Burroughs 1984; Rood 1984; Gimbarzevsky 1988). Regional -7- studies are used to delineate problem-areas where actions should be taken. The information collected is less detailed. 1.2.2. Classification by survey intensity level Terrain stability assessment is usually carried out at two levels. Reconnaissance mapping is used to map larger areas at a smaller scale (1:50, 000 to 1:250,000) with survey intensity level D (Table 1.1). The assessment is usually made based on limited field checking, usually helicopter overflights. The major objective is to pinpoint areas that require detailed terrain stability assessment. Detailed terrain stability mapping is usually carried out at the scale of 1:20, 000 with survey intensity level B. The objective is to identify those areas that require site-specific terrain stability assessment prior to approval of road construction, layout of cutblock boundaries, timber harvesting methods and silvicultural systems. Table 1.1. Survey intensity levels for terrain stability mapping % of terrain polygons fieldchecked 75-100 < 1:20, 000 1:10, 000 to 30-75 1: 50, 000 1: 20, 000 to 25-50 1: 100, 000 1 : 50, 000 to 0-25 1: 250, 000 any scale 0 Terrain Survey Map Scale Intensity Level A B C D E ModifiedfromGerath et al (1994). Method of Fieldchecking foot traverses foot and vehicle traverses vehicle and helicopter overflights vehicle and helicopter overflights no field work, air photo interpretation only 1.2.3. Classification by analytical methods Terrain stability assessments can be classified based on the analytical method. The following section provides an overview of the different methods placing emphasis on the recognition of their advantages and disadvantages. Terrain stability assessments can be divided into five major groups (Table 1.2). Table 1.2. Classification of terrain stability assessment by analytical method 1. Landslide inventories 2. Factor overlay methods 3. Subjective geotechnical models 4. Factor of Safety models 5. Statistical models Damage inventories Landslide activity analyses Landslide distribution analyses Distribution models Probability models Univariate models Multivariate models 1.2.3.1. Landslide inventories The aim of landslide inventories is to document soil mass movement processes on a given area. Landslide inventory maps are sometimes assembled for the purpose of documenting the distribution and the degree of damage incurred in a region (Kienholz 1978; Ellen and Wieczorek 1988; Mark 1992). Damage inventories usually document the impact of major storms (Schwab 1983; Evans and Lister 1984) or other catastrophic events such as earthquake. A landslide inventory following the 1946 Vancouver Island earthquake was documented by Mathews in 1979. The data collection can be done systematically to provide a good base for further analyses. The Gifford Pinchot National Forest in southern Washington developed a Geologic Resource Database in 1982. This database was used to construct Geological -9- Resources and Conditions Maps. All existing landslides were classified on the basis of the process of failure and mapped along with engineering and geological information related to existing soil and rock conditions. Using this information, the geotechnical personnel could establish the priority of the various slope stability problems on a Forest-wide basis and schedule field investigations in a more systematic manner (Reilly and Powell 1985). Landslide distribution analysis is used for more detailed hazard analysis. It can be used to develop an understanding of the relationship of different land use patterns (e.g., logging) to the occurrence of soil mass movements (Gerath et al. 1994). Many studies have dealt with the impacts of logging and road building on the frequency and yield of landslides from steep cleared slopes and road areas (Swanson and Dyrness 1975; Morrison 1975; Swanson and Swanston 1977; Fiksdahl 1974; O'Loughlin 1972). These landslide distribution analyses have shown that logging accelerates the frequency of landsliding by up to fifteen times (Rood 1984), while road construction may accelerate slope failures by up to 346 times (Morrison 1975). Very detailed landslide inventories and distribution analyses were prepared on the Queen Charlotte Islands (Rood 1984; Gimbarzevsky 1988) as part of the Fish/Forestry interaction program. This program was initiated in 1981 by the British Columbia Ministry of Forests, British Columbia Ministry of Environment, and Canada Department of Fisheries and Oceans as a positive action toward the resolution of conflicts concerning steep slope logging and integrated management of fish and forest resources on the Queen Charlotte Islands. The major problem with landslide inventories and distribution analyses is that they record activity only within a specified area, since they provide no information on the potential instability of areas other than those which experienced slope failures during the inventory time period. Landslide inventories in the past did not include information on landslide rate for forested areas, partly because of the difficult access - 10- to such areas and partly because landslides are very hard to recognize on aerial photographs taken of densely forested terrain. Such inventories provide insufficient information for comparison, since it is almost impossible to evaluate any acceleration of landsliding associated with forest management activities. Another problem of inventories is that shallow failures of colluvial veneer are cyclical in nature. When a failure occurs at a given site, a dormant period follows while the forest cover regenerates and new surficial layer develops. The landslide inventory does not take into account sites that are mature and primed for a new series of soil mass movements (Gerath et al. 1994). In landslide activity analyses, data are obtained from serial aerial photographs. The information is included from several time periods to detect how the landscape has changed over time. This information can also be used to compare landslide frequencies before and after logging activities (Swanson et al. 1982). Table 1.3 summarizes the advantages and disadvantages of landslide inventories. Table 1.3. Advantages and disadvantages of landslide inventories Advantages Disadvantages 1 • provide good database for future use • record within specific time interval • represent simple hazard assessment • some slides are cyclical in nature • do not follow changed conditions • do not distinguish stable/unstable sites | 1.2.3.2. Factor overlay methods One of the most common landslide delineation methods is factor overlay of a combination of landslide producing elements. Once key landslide producing factors are identified they must be combined in a manner that will yield realistic landslide potential designations. Most investigators rely on the basic static factors of slope gradient, geology and soils plus a dynamic factor that represents water (Ward 1976). In the overlay approach, certain factors related to landslide occurrence are individually delineated. These factors are then combined via overlay techniques to yield sets of intersecting factors that can be assigned landslide potentials (Nielsen and Brabb 1977). This approach is subjective in the sense that it assigns an equal weight or value to each factor. The factors can be assigned numerical values in order to develop threshold levels between potential classifications. A procedure in the WRENSS (Water Resources Evaluation of Non-point Silvicultural Sources) was presented as a guide for assessing the stability of natural slopes, the potential impacts of silvicultural activities on slope stability, and forecast sediment contribution to water courses from soil mass movements (Swanston et al. 1980). The procedure involves adding a series of "weighting factors" to estimate landslide hazard. Different factors are considered for shallow, rapid failures (debris avalanches-debris flows) and for deep-seated, slowly moving failures (slumpsearthflows) in assessing natural hazard. In many cases the factor overlay technique is applied to regional level analysis. The reason for this is that the method usually provides acceptable results for broad delineation of potentially unstable areas. More detailed analysis needs to be done on the landslide-prone terrain. Early efforts to rate landslide severity of various physiographic regions of the United States using factor overlay were reported by Baker and Chieruzzi (1959). Their ratings of landslide severity were based on the information gathered for the various physiographic regions with regard to frequency of occurrence, size of moving mass, and dollars expended per year. A more recent study described the magnitude, intensity, and distribution of slope stability problems within the northwestern United States (Burroughs 1985). He overlaid areas with high hazard levels for natural landslides, important fishery-water resources and - 12- areas of major levels of management activity to pinpoint zones with a large potential for environmental damage as a result of harvesting operations. Rood (1990) used a previous landslide inventory (Rood 1984) to identify the site and regional characteristics of a large number of failures on the Queen Charlotte Islands. He described the areal distribution of these characteristics within clearcut areas. The data did not distinguish between stable and unstable sites, but rather described the frequency and magnitude of soil mass movements associated with the different characteristics. There are two major problems with the regional approach. First, it does not provide a clearly delineated landslide potential classification. Second, the analysis is based on fairly large areas, so many of the local factors causing slope failures are integrated or averaged out. This approach is good for delineating potentially hazardous areas that are subject to subsequent detailed stability analysis. Table 1.4 summarizes the advantages and disadvantages of the factor overlay techniques. Table 1.4. Advantages and disadvantages of the factor overlay techniques Advantages • good to pinpoint problem areas • useful in regional development plans Disadvantages • provide broad delineation only • assign equal weight to each factor • oversimplify • danger of wrong factor selection • local factors are averaged out | 1.2.3.3. Subjective geotechnical models Subjective geotechnical models can also be called empirical models since they require a high degree of specific local knowledge and experience. The accuracy of this method depends upon the mapper's experience and ability to make subjective - 13 - observations and associations. The method involves the selection of a classification scheme to divide the heterogeneous landscape into units which respond similarly to various management activities. These homogeneous units (called terrain polygons) are delineated based on certain classification criteria that may change from one study to the other. Ryder and Howes (1984) delineated terrain units by type of surficial material, texture, slope gradient (surface expression) and geological processes. Others used landform, geologic parent materials, soils and vegetative habitat types as mapping criteria (Wilson 1985). Due to regional variations in climate, geology, soils and other factors, there are no specific criteria that can be applied within all physiographic and climatic regions of the province. After delineating the polygons, the mapper subjectively assigns a qualitative hazard rating to the various mapping units on the basis of experience, field observations and informed guessing. It is important to note that the rules of assignment of the hazard classes are not specified and can vary from polygon to polygon. Several studies applied this method using the British Columbia Terrain Mapping System for polygon-delineation (Ryder and MacLean 1980; Howes 1981). The key to any subjective mapping system designed to interpret slope stability is that the mapping system must be capable of extrapolating limited research data to a variety of landscapes on the basis of observable characteristics of those landscapes. The criteria for terrain stability hazard classes are typically qualitative and are highly dependent on the knowledge and the local experience of the terrain mapper. Recognition of unstable terrain is often guided by observations of existing failures. This approach was used to evaluate slope stability in the Clearwater National Forest in northern Idaho. Their mapping system, based on dominant geomorphic erosion processes, allows mappers to stratify the base erosion rate over the landscape (Wilson 1985). Another example of empirical modeling efforts was - 14 - made at the Siuslaw National Forest, Oregon. The method was based on the logic that landslide processes accelerated by specific management activities on land with carefully defined physical characteristics will accelerate in the future at a similar rate providing that the management impacts and climatic parameters are similar (Bush 1985). The main advantage of empirical modeling is that no attempt is made to understand the interaction of the myriad of variables involved in the landsliding process as theoretical models do. The empirical modeler can effectively predict future landslide rates while avoiding some of the pitfalls of having a less than perfect understanding of all the processes involved. The subjective geotechnical approach is quite flexible, but it also has some drawbacks. Due to the subjective elements involved, it is not repeatable and it totally relies on the skills and local experience of the mapper. Lack of repeatability can be overcome by assigning an algorithm for the entire study area based on subjective weighting of relevant factors. This algorithm is then applied to the entire project (Chang 1992; Sidle 1985) The MacMillan Bloedel Limited Subjective Rating System, originally developed by Bourgeois (1974, 1978), defines five landslide susceptibility classes on the basis of set combinations of landform, material texture, slope characteristics, soil type, moisture regime, landscape position and underlying bedrock type. The system is applied to terrain polygons defined with the B.C. Terrain Classification System (Howes and Kenk 1988). The main advantage of this method is its repeatability. The class assignments can be independently checked and a record of the procedure (rating algorithm) exists. Table 1.5 summarizes the advantages and disadvantages of the subjective geotechnical models. - 15 - Table 1.5. Advantages and disadvantages of the subjective geotechnical models Advantages Disadvantages • simple, easy to implement • flexible • • • • subjective local experience is essential not repeatable rating algorithm non-transferable 1.2.3.4. Factor of safety models The assessment of slope failures may involve the limit equilibrium mechanics of slope stability. The analyses compute a safety factor (FS) that defines the ratio of shear strength (S) to the shear stress (T) required to bring the slope into a state of limit equilibrium along a failure plane. This limit equilibrium analysis can be applied to planar slip surfaces, circular slip surfaces and non-circular slip surfaces (Sidle 1991). The most widely applied analytical method for steep forested slopes has been the planar infinite slope analysis. It is based on the assumption that the thickness of the soil mantle is small compared to the slope length and the failure plane is parallel to the slope. Several landslide models have been developed on the basis of the infinite slope equation (Ward 1976; Wu and Swanston 1980; Ward et al. 1981; Sidle 1985, 1991, 1992; Hammond et al. 1992). Two-dimensional, shallow slope failures (debris slides, debris avalanches) on infinite slopes are expected to occur when the available shear strength (S) is less than the applied shear stress (T) FS = S/T < 1 Shear strength can be calculated on the basis of the modified Mohr-Coulomb equation: S = C + A C + (cr-//)tan^; - 16- where S: soil shear strength; C: soil cohesion; AC: cohesion derived from root strength; a: total normal stress; p.: pore water pressure; <j>: effective angle of internal friction. In forested areas root systems contribute to soil strength by providing an artificial cohesion (AC) that can be added to the effective soil cohesion. The applied shear stress (T) is a result of the weight of trees, soil, and moisture on the slope and slope angle. Theoretically, when FS < 1 the slope is unstable and failure is imminent. Soil properties are highly variable over space and time on natural slopes. The important factors affecting shallow soil mass movements are: pore water pressure, which varies with storm patterns and local slope hydrology; soil characteristics; slope steepness; the depth to the potential failure plane; and in clearcut areas the removal of vegetation and the subsequent root system deterioration (Sidle et al. 1985). Other than slope steepness, all these parameters are difficult to measure in the field and vary greatly along and across natural slopes. The use of Factor of Safety models in regional approaches is limited, since it is impractical to directly measure any of the above mentioned parameters other than slope gradient. There are two basic approaches in Factor of Safety models. The deterministic approach deals with the distribution of the Factor of Safety over the study area, while in the probabilistic approach the distribution of probability of a safety factor less than one is examined. The deterministic techniques must be applied at grid points, while the probabilistic approach can be applied to terrain polygons (Hammond et al. 1992). Ward (1976) applied a simplified factor of safety model. In his segmentation-landslide potential delineation model, factor of safety was calculated at grid points. Parameter maps were overlaid to get information for the - 17- grid points. He calculated an average safety factor for each grid cell on the basis of the average input data of each node point representing the cell. In nature, soil conditions are never uniform. The probabilistic approach tries to estimate the probability of failure of a site based on the knowledge of the temporal distribution of the deterministic factors. With the probabilistic approach, the engineer can judge the acceptability of a finite probability of failure instead of the adequacy of a calculated factor of safety (Schroeder 1985). Sidle (1992) presented an infinite slope stability model that incorporated changes in root cohesion and vegetation surcharge through several timber management cycles along with the stochastic influence of rainfall on pore water pressure. These features allow the user to simulate a variety of complex vegetation management scenarios in relation to slope stability. The probability of slope failure can be related to change in root strength and other site variables in any given year, as well as to the probability of occurrence of a landslide triggering storm. The model applies only to shallow, translational failures. Actual values of probability of failures calculated by the model should be viewed with caution due to the assumptions inherent in the infinite slope model and the spatial variability of site parameters. The model provides a viable relative comparison of different silvicultural systems. Soil strength at lateral boundaries resists movement in a landslide (vertical variation). The factor of safety is usually underestimated, because these forces are not taken into account in the usual two-dimensional analysis. Cost and time constraints limit our knowledge of input data. Therefore, on natural slopes where geologic conditions are nearly always complex, where groundwater conditions are hardly ever known reliably, and where the topography must be approximated for the analysis, engineering techniques will never have great precision. Table 1.6 summarizes the advantages and disadvantages of the factor of safety models. - 18- Table 1.6. Advantages and disadvantages of the factor of safety models Advantages • seem to be objective • easy to calculate • repeatable Disadvantages • neglects vertical variation in soil • depend highly on input data • some parameters difficult to obtain 1.2.3.5. Statistical models Statistical analysis is another approach to landslide hazard assessment. Numerical data and statistical techniques are used to produce equations that yield a value related to slope stability. In the assessment of slope stability, a direct statistical correlation is sought between the probability of landslide occurrence and a single variable (univariate models) or a group of variables (multivariate models). Regardless of the number of variables used, there are two basic approaches. The great majority of statistical methods reported in the literature attempt to set up relative correlation based on the assumption that mapping units (usually terrain polygons) that are similar in certain factors to areas which failed in the past are most likely to fail in the future (Brabb et al 1972; Rice and Pillsbury 1982). Only a few studies deal with absolute correlation assuming that future landslide frequency in similar units can be predicted, knowing landslide frequency in failed units, over a given time period (Rollerson and Sondheim 1985; Rollerson 1992). The problem with the latter approach is that landslides are rare phenomena for which it is very difficult to establish a temporal frequency. 1.2.3.5.1 Univariate statistical analysis In univariate methods, the relationship between landslide frequency and slope - 19- attributes is examined separately for each factor. The relationship represents a set of weighting factors which are combined to produce a relative hazard rating. In most cases the analysis is based on terrain polygons that are delineated by overlaying parameter maps (Brabb et al. 1972). The major problem with the polygon-based methods is that these polygons are delineated on the basis of surficial geology (surficial materials, slope gradient, geological processes, soil characteristics etc.) without considering the subsurface conditions. This can bias the resultant landslide frequencies which provide the base of the entire analysis. Most of the studies are restricted to logged areas. They compare terrain characteristics of logged sites where failures have occurred to those characteristics on non-failed sites. In this way, factors that appear to most influence landsliding can be isolated (Rice and Pillsbury 1982; Chatwin and Rollerson 1983; Rollerson and Sondheim 1985; Howes 1987; Rollerson 1992). Rollerson (1992) applied the univariate method to 6-to 15-year old logged areas on the Queen Charlotte Islands. The lower age limit was imposed to allow time for root strength to deteriorate and to provide reasonable opportunity for large storms to act on logged terrain. The upper age limit was set because well-advanced conifer regeneration often masks the presence of smaller slope failures. The study addressed post-logging landslide frequencies as they related to individual landscape attributes and combinations of landscape attributes. He found that different attributes were significant for clearcut landslides and road fill failures. The choice of relevant attributes and their use in formulating the multi-factor classification is guided by a parallel relative analysis of each separate landscape attribute. A highly systematic approach to ranking the different types of terrain for likelihood or frequency of failure has been developed on the basis of this work is known as "terrain attribute study" (Rollerson, T.P. 1994. personal comm.). The objective of the terrain attribute study is to characterize terrain types that are subject to landslides following harvesting operations or road -20- building. The first step involves the terrain mapping and data collection in logged areas using the Terrain Classification System for British Columbia. In the second step data are collected on terrain attributes such as surficial material, soil type, drainage, slope gradient, aspect, slope morphology, bedrock and pre- and post logging landslide activities. In the third step a statistical analysis is carried out on the database to identify relationships between post-logging landslide frequency and terrain attributes. Although the univariate statistical technique gives quite reliable estimates of landslide hazard, there are some problems. It does not account for the interaction (effect modification) and confounding effects that may exist among the different factors. It relies on high quality data that are often difficult to obtain. Most of the analyses are concentrated on logged areas, so the inventory of landslides and terrain attributes may be biased with ease of road access. If there is no data collected on terrain without landslides the study will not be suited to developing predictions of probability of landslides occurring in any particular place of the landscape. Another problem is that the frequency distribution of the various parameters is not known. The advantage of the univariate model is that it allows the mapper to observe the influence of individual slope attributes. It is simple to implement and test. Selection of relevant attributes and definition of terrain classes require careful and thorough work. Table 1.7 summarizes the advantages and disadvantages of the univariate statistical methods. 1.2.3.5.2. Multivariate statistical analysis Multivariate analyses use multivariate regression techniques to establish a correlation between landslide frequency and a group of attributes. This method also depends highly on the quality of the input data. A combination of the factor overlay -21 - method and multivariate statistical techniques was used by De Graff and Romesburg (1984). In their matrix approach, they defined a separate class for each combination of the independent variables. This leads to a large number of combination classes even with a very limited selection of variables. The method requires a detailed database of landslide frequencies to achieve statistically significant correlation. More sophisticated techniques were used to define stable and unstable sites by Furbish and Rice (1983), Rice et al. (1985), and Rice and Lewis (1991). These techniques were based on discriminant function analysis (Fisher 1936). The discriminant function is a statistical procedure used to classify observations into populations based on multiple parameters. The discriminant function provides the maximum separation of populations (stable and unstable sites) on the basis of site characteristics. Discriminant analysis produces an equation that can be used to estimate the probability that a site will become unstable if logged or roaded. Table 1.7. Advantages and disadvantages of the univariate statistical method Advantages • more objective • easy to implement and test • repeatable • relevant attributes can be selected • influence of individual attributes can be observed Disadvantages • rely highly on data quality • do not account for interactions • not transferable to other areas • studies concentrated on logged sites • biased with ease of road access • freq. distribution of parameters is not known • terrain polygons can be biased Discriminant function analysis has been used as an objective method for estimating landslide risk in grass and bush environment (Rice and Foggin 1971). Pillsbury (1976) used discriminant function analysis to identify potential landslide sites after logging. His equation correctly identified 80 percent of the sites in the -22- developmental data. A subset of Pillsbury's data was used by Rice and Pillsbury (1982) to develop a new discriminant function having a logarithmic form. The new equation identified 90 percent of the sites accurately based on the developmental data. In a more recent study, Lewis and Rice (1990) used discriminant functions to separate stable and unstable sites using different discriminant scores for forest roads and logged areas. The resulting equations had a classification accuracy of 78 percent for road plots and 69 percent for logged area plots. In central Italy, Carrara (1983, 1988) applied discriminant function analysis to classify stable and unstable units on the basis of their morphological, geological and vegetational characteristics. He randomly split the dataset into two parts. One part was used to estimate the model (learning or developmental dataset) and the remaining part was used to test the goodness of fit (test dataset). Although discriminant functions seem to be effective in the identification of landslide-prone terrain, their use does not lack problems. Whether they are used for prediction or as an interpretive tool, their reliability depends on how well they describe the operable processes and conditions affecting slide occurrence (Rice et al. 1985). Discriminant function analysis is based on the assumption of multivariate normality. Since geographical variables are often not normally distributed, this assumption will not be met. Certain statistical techniques can be applied to the data to achieve normality, but in most cases the violation of the assumption of normality cannot be avoided. The major difficulty in the use of the predictive equations is establishing levels of "acceptable" risk marking the boundary between stable and unstable sites. Threshold values are extremely difficult to deal with since they may increase the percentage of misclassified sites. The costs of such errors include erosional damages at undetected critical sites and the treatment cost or loss of timber revenue. The choice of threshold value should be based on economic, social and political considerations. Studies in the past tended to use too many variables in the -23- analysis. Recently, more sophisticated statistical techniques have been used to screen variables for inclusion. The main disadvantage of the multivariate approach is that it eliminates human experience and judgement from the analysis, thus the results are highly dependent on the quality of input data. Table 1.8 summarizes the advantages and disadvantages of the multivariate statistical methods. Table 1.8. Advantages and disadvantages of the multivariate statistical method Advantages • objective • repeatable • fairly high classification accuracy Disadvantages • difficult to establish threshold values • depend highly on input data • non-transferable to other areas • ignore mapper's experience • complicated procedures -24- 2. METHODS 2.1. Introduction This section provides an overview of certain statistical techniques used in epidemiology and biomedical research. The case-control sampling procedure and the multiple logistic regression analysis is introduced and solutions are sought for implementing them in terrain stability assessment. The idea of using logistic regression models in landslide hazard assessment originates from two observations: 1. Landslides like some diseases are rare phenomena, thus statistical techniques used in epidemiology may be applied to the problem of terrain stability assessment. 2. Geographic variables are hardly, if ever, normally distributed over the landscape which is an assumption of multivariate statistical procedures such as discriminant function analysis. Examples used in this section are part of this study that was carried out on the Jamieson-Orchid-Elbow Creeks study area. -25 - 2.2. Study area 2.2.1. Location The study has been carried out on the 1201-ha Jamieson-Orchid-Elbow Creeks subdrainage. This mountainous drainage is located within the rugged Coast Mountain physiographic area, in the Seymour River Basin. The Seymour Basin is one of the three major Greater Vancouver Water District (GVWD) community watersheds (Seymour, Capilano and Coquitlam) that provide water to the city of Vancouver (Figure 2.1). The study area has been chosen for two reasons. First, it is located directly north of Vancouver, so there is no problem with accessibility. Second, in 1992 an ecological pilot study was undertaken by an interdisciplinary team on the study area to develop an inventory method based on biogeoclimatic classification principles (Acres International Ltd. 1993). This team developed a digital database that will be used in the study. The database provides information on vegetation, terrain, hydrology and wildlife. This study gives an opportunity to examine the usefulness of such digital databases for terrain stability assessment. 2.2.2. Climate The outstanding features of the climate along the southern coast of British Columbia are: - mild and wet winters; -rain-on-snow events in the winter; - very heavy precipitation and cloudiness; - large accumulation of snow at high elevations, which generally lasts into the -26- -27- summer; - striking differences in temperature and precipitation as a result of changes in elevation, distance from the coast, and complex mountain topography (Beaudry 1984). Several precipitation gauges were installed at Jamieson Creek and Elbow Creek as a part of a research project and maintained jointly by the GVWD and the University of British Columbia, Faculty of Forestry. Average annual precipitation and intensity-duration data are shown for three precipitation gauges in the study area in Table 2.1. Table 2.1. Annual precipitation and intensity data for the study area Station Name Elevation (m) Avg. annual ppt. (mm) Intensity (mm/24hr) for 10 year recurrence interval 21-A 640 3294 (71-87) 167.6 25-B 762 3141 (81-87) 203.2 ECr 300 3427 (71-87) 154.9 Adapted from Hall, R.G. (1989) Historical data show that about 60-70 percent of the annual precipitation falls in the rainy season from October to March. The area is subjected to frequent, high intensity storms with intensities over 152 mm/24 hr at 10-year recurrence interval or less (Table 2.1). Because the greatest amounts of precipitation coincide with the presence of snow in the winter, flood and debris torrent producing events and headwater erosion may occur from the combination of large amounts of rainfall coupled with rapid snowmelt (Beaudry 1984). Precipitation is highly affected by elevation and distance inland, Schaefer and Nikleva (1973) found that the proportion of annual precipitation falling as snow increased by about 20-25 percent per km -28- increase in elevation in the Coast Mountains near Vancouver. 2.2.3. Topography The study area has a well-defined topographic boundary with a watershed area of 1201 ha. The elevation ranges from 255 m to slightly over 1400 m at the height of land dividing the Orchid Creek watershed from the Capilano watershed. The slopes are generally steep with a considerable area of shallow soils. More than 50 percent of the area occurs on steep (over 30 degrees) slopes (Table 2.2 and Figure 2.2). About two-third (65.7 %) of the shallow soils are located on steep slopes while about the same proportion of deep soils (69.4 %) can be found on relatively gentle slopes (less than 30 degrees). Most sites with poor drainage are located on slopes greater than 30 degrees (Acres International Ltd. 1993). Table 2.2. The distribution of slope angle and soil depth in the study area Slope class in degrees 0-10 10-20 20-30 30-40 40-50 50< Total: Percentage % 6.3 12.6 25.7 37.1 12.9 5.4 100.0 Shallow 2.8 9.0 22.5 42.3 18.1 5.3 100.0 Soil depth Medium 6.8 4.7 11.8 41.4 13.7 21.6 100.0 Deep 12.3 21.4 35.7 26.6 3.6 0.4 100.0 Topographic features of the area include: - relatively flat valley bottoms along the Seymour River; - moderately and steeply sloping valley sides dissected by V-shaped gullies; - and mountainous terrain along the western border (Acres International Ltd. 1993). -29- -30- These characteristics along with intense and long-duration rainfalls predispose the area to shallow, rapid soil mass movements such as debris avalanches, debris flows and torrents. 2.2.4. History of soil mass movements There are about one hundred landslides on the study area. These landslides are typically small in volume, representing the debris avalanche-debris flow type soil mass movements (Swanston and Swanson 1976). These landslides are rapid, shallow soil mass movements from hillslope areas, occurring on shallow, non-cohesive soils, where subsurface water may be concentrated by subtle topography on bedrock or glacial till surfaces. Table 2.3. Historical landslides in the study area Year Location Jamieson, Elbow 1930's Cks. 1970's Orchid Creek Orchid Lake area 1980 Jamieson Creek 1983 Jamieson Creek 1990 Triggering event Remarks max. daily rain-fall: Debris torrent 250 mm stable debris torrent deposit Dec. 26, max. daily Debris torrent rainfall: 150 mm Nov. 15 moderately channel bank erosion large flood with 5 year return period Nov. 23, max. daily originated in a cutblock, 1500 m long rainfall: 300 mm debris torrent Most of the information on landsliding was collected from the Ecological Inventory Pilot Study digital database (Acres International Ltd. 1993) and the review of landslide studies that had been carried out in the Seymour River Basin. -31 - Figure 2.3. shows the location of landslide initiation points on the Jamieson-Orchid -Elbow Creeks study area. Table 2.3 provides an overview of historical debris torrents within the study area. There have been several terrain stability studies conducted on the GVWD community watersheds. Most of these studies have been completed at an overview level. Table 2.4 provides a short summary of these studies. Table 2.4. Overview of landslide studies in the study area Author Schultz, C.W. Year 1956 O'Loughlin, C.R. 1972 Briere, D. 1979 Luttmerding, H.A. 1981 Thompson, G. and 1990 Golding, D.L. 1991 Thurber Engineering Ltd. Acres International 1993 Ltd. Remarks to compare relative erodibility of various materials examination of root slope stability analysis strength deterioration and slope stability Aqua-Terra Classification adopted by GVWD in 1984 System comments on potential soil mapping for surface erosion after timber harvesting and on relationship between soil types and biogeoclimatic zones Evaluation of the AquaLandslide inventory Terra Class. System geotechnical assessment after 1990 November rainstorm events of 35 landslides Ecological Inventory Terrain stability assessment for the study area Pilot Study using B.C. Terrain Classification System Nature of study surficial geology 2.2.5. Vegetation The Jamieson Creek watershed study area consists of mature and over mature -32- -33- forest stands. It includes two biogeoclimatic zones, the Coastal Western Hemlock Zone (CWHZ) from the valley floor to about 900 meters above sea level (masl) consisting of Douglas fir (Pseudotsuga menziezii (Mirb) Franco), western hemlock (Tsuga heterophylla (Rafn.) Sarg), western red cedar (Thuja plicata Donn), sitka spruce (Picea sitchensis (Bong.) Carr.), and the Mountain Hemlock Zone above 900 masl dominated by mountain hemlock (Tsuga mertensiana (Bong.) Carr.), yellow cedar (Chaemacyparis nootkatensis (D.Don) Spach) and amabilis fir (Abies amabilis (Dougl.) Forbes). The Mountain Hemlock Zone receives more precipitation, especially snow, than the CWHZ and is subject to cooler and shorter summers (Acres International Ltd. 1993). 2.2.6. Forest harvesting and hydrology Several graduate theses (Ph.D. 8, Masters 3 and Bachelors 8) and scientific papers have dealt with the hydrologic regime of the study area, often focusing on the effects of harvesting operations (Golding 1988). An important characteristic of the Jamieson-Orchid-Elbow watersheds is the almost total absence of surface runoff. This is due to the high permeability of the mineral soil and the negligible presence of frozen soils. Thus, all rainfall and snow melt water percolates through the soil to stream channels (Beaudry 1984). Jamieson and Elbow Creeks were instrumented in 1969 as treatment and control watersheds in a paired watershed experiment to examine the effect of harvesting on stream flow. They were established in response to (1) lack of knowledge of the hydrologic system in South Coastal British Columbia and (2) concern for the effects of timber harvest on stream flow. The Greater Vancouver Water District practices very conservative forest management in order to provide good quality water instead of maximizing the return from the timber resource. The -34- pre-treatment period extended from 1972 to 1978. About nineteen per cent (19.2%) of the Jamieson Creek watershed was harvested during the 1978-84 period. The post-treatment period began in the fall of 1984 (Golding 1988). 2.3. Theoretical background The following section provides an introduction to the case-control sampling method and the logistic regression model. It is shown how a case-control sampling can be analyzed by multiple logistic regression. The discussion is restricted to the techniques considered to be useful in the development of a terrain stability assessment procedure. For detailed discussion of the various topics mentioned in this study the reader should consult any of the biostatistics textbooks. All the examples used in this section are part of this study that was carried out on the Jamieson-Orchid-Elbow Creeks study area. Detailed analysis and discussion will be given in the "Data analysis and discussion" section. 2.3.1. Case-control studies Generally speaking, there are two fundamental research strategies. Observational studies do not involve human intervention, while experimental studies involve planned intervention on factors suspected of altering the phenomenon under study. The objective of both strategies is to examine the cause-effect relationships between a phenomenon and various factors. Clearly, the experimental approach cannot be applied to landslide studies since: 1. human intervention on some factors is limited or impossible 2. ethical considerations restrict experimentation, since causing a disastrous event should never be the objective of any study -35 - 3. required sample size may be exceptionally large to detect a possible increased risk 4. duration of experimental study may be very long Observational studies are designed to identify factors that cause or prevent the occurrence of an event. There are two major approaches in observational studies: 1. Cohort (prospective or follow-up) studies 2. Case-control (retrospective) studies. The major difference between these two approaches is in the selection of study subjects. In a cohort study, subjects free of the hazardous event under study are selected for observation and followed over time to determine the rate at which the event develops. The term "event" in this study will be used as a substitute for hazardous event or landslide event. In a case-control study, subjects are selected on the basis of the presence (case) or absence (control) of the event. Cases and controls are then compared with respect to attributes or "exposures" thought to be relevant to the development of the hazardous event. The latter approach proceeds from effect to cause while the former proceeds from cause to effect. Table 2.5 provides a summary of the advantages and drawbacks of the case-control approach. Table 2.5. Advantages and disadvantages of a case-control study Advantages: Disadvantages 1. Well suited to study rare events 2. Relatively quick to conduct 3. Relatively inexpensive 4. Requires relatively small sample size 5. Allows study of multiple potential causes of an event 6. Analysis can be carried out on a microcomputer 1. Control of extraneous variables may be incomplete 2. Method is unfamiliar to terrain stability experts ModifiedfromSchlesselman (1982). -36- Landslides are rare phenomena. In this situation a cohort study is inefficient, since practically all of the effort would be devoted to follow up of subjects that remain free of the study event. In the case-control method, we study sites that have already developed the event, so there is no need to wait for the occurrence of the event. A case-control study is relatively inexpensive since it requires a small sample size and therefore it is relatively quick to conduct. 2.3.1.1. Assessment of risk in case-control studies The probabilistic approach provides the basis of discussion for evaluating the effect of various factors on the risk of an event. Let y denote the dichotomous response variable indicating the presence (1) or absence (0) of a hazardous event, such as landsliding. Let x denote an "exposure" variable, such as slope angle in degrees. X is said to be the cause of y if the probability (relative frequency) that y occurs is increased as a consequence of the presence of certain x values. The relationship between x and y may be influenced by different levels of other extraneous variables; it may be strengthened, weakened or may entirely disappear. The problem of interpreting the apparent effect of the x "exposure" variable, in view of potential bias from extraneous factors, will be addressed in the discussion of confounding and interaction. The assessment of risk can be based on the number of new cases that occur during a specific period of time or on the number of cases at a given point in time. In epidemiological studies, the former is called incidence and the latter is called prevalence. In landslide hazard assessment both approaches may be used, though in most cases prevalence should be used due to insufficient information on the time of occurrence of landslide events in the past. In this study prevalence will be used, i.e., all landslide events that occurred till the end of 1992 that could be identified -37- from aerial photographs will be used. The ratio of the number of cases occurring over a time period to the number initially at risk is called incidence rate and is expressed by: p=c/n; (Eqn. 2.1), where p: incidence rate c: number of cases n: number of subjects initially at risk This ratio can be interpreted as the average conditional probability that a landslide develops in the given time period. Since each subject may be considered to have a distinct risk of landsliding, p may be interpreted as an estimate of the average risk of landsliding for sites in the given geographical area. 2.3.1.2. Relative risk and the odds ratio Consider two groups of sites that are similar to each other in all attributes relevant to a landslide event, apart from the presence or absence of exposure to slope angle greater than 35 degrees (Tables 2.6.A and 2.6.B). Table 2.6.A. Frequency of landsliding among exposed and unexposed sites in the target population Landslide Present Absent Exposure A B Yes D C No N2 Nl Total B/D Odds of exposure A/C Total Ml M2 N Odds of event A/B C/D Adapted from Schlesselman (1982). The incidence rate among exposed sites can be calculated by : pi = A/Ml (Eqn. 2.2) -38- The incidence rate among the unexposed is: P2=C/M2 (Eqn. 2.3) Table 2.6.B. Relationship of landsliding and slope angle (SL35) in the JamiesonOrchid-Elbow Creeks study area using SL35 as a dichotomized variable SL35 Slope angle > 35° Slope angle < 35° Total Odds of exposure Landslide Absent Present 44* 5474 57 20972 101 26446 0.77 0.26 Total Odds of event 5518 0.0080 21029 0.0027 26547 * the study subjects are 20x20 m grid cells on the Jamieson-Orchid-Elbow Creeks study area The proportion of the two incidence rates is called relative risk (R): R=P1/P2 (Eqn. 2.4) Relative risk represents how much more likely or unlikely a hazardous event occurs in the exposed group as compared to the unexposed group. If R differs from unity, then the study factor is said to be associated with the risk of an event. The odds ratio is also a measure of association that is closely related to the relative risk. If an event occurs with a probability of p, then the ratio p/q is called the odds, where q=l-p. If pi denotes the incidence rate among exposed sites, the odds of event are pi/qi while P2/q2 represents the odds of event among unexposed sites. For rare events, such as landsliding, the risk of an event p and the odds of event p/q are virtually identical since p is a very small number and q approximates unity. The ratio of the odds of event in exposed sites relative to the unexposed sites is called the odds ratio: y=(Pl/qi)/(P2/q2) (Eqn. 2.5) According to Table 2.6. A, the odds ratio can be expressed by the formula: -39- Y=(A/B)/(C/D)= AD/BC (Eqn. 2.6) The odds ratio is particularly important in the risk assessment of rare events since the odds ratio provides a good approximation to the relative risk. Consider the example given in Table 2.6.B: pi = A/Ml =44/5518=0.0080 P2=C/M2 =57/21029 =0.0027 R=Pl/P2=2.96 SV=AD/BC=(44x20972)/(57x5474) =2.96 The risk of landsliding among unexposed sites (p2) is 0.0027 which represents a rate of 27 per 10,000 sites. The risk of landsliding for sites with an average slope of 35.01 degrees or greater is 80 per 10,000 sites. The relative risk of a landslide event for sites with steep slopes as compared with sites with gentle slopes is 2.96. The calculated odds ratio Q¥) is also 2.96. Using either measure of association, one can state that the risk of landsliding for steep slopes is about 3 times higher than for gentle slopes. Among sites with a landslide event, the odds of exposure can be defined as A/C (Table 2.6.A). The odds of exposure among sites without landslides are expressed as B/D. The ratio of the odds of exposure in sites that have experienced landslide events as compared with those without landslides is given by: F=(A/C)/(B/D)= AD/BC (Eqn. 2.7) X Thus the "exposure odds ratio" is equivalent to the event odds ratio defined by Equation 2.6. This is a very important and useful relationship in the consideration of design and analysis of case-control studies. 2.3.1.3. The case-control sampling scheme In a case-control study we take a random sample of cases of the hazardous -40- event under study and a random sample of controls (subjects that are free of the study event). The sampling fractions of cases and controls are denoted as fl and f2, respectively (Table 2.7. A). Table 2.7.A. Case-control study based on the target population of Table 2.6.B Expected sample outcome for a case-control study* Landslide Exposure Present Absent b (f2B) a (flA) Yes d(f2D) c (flC) No nl n2 Total Odds of exposure a/c (A/C) b/d (B/D) Total ml m2 n Odds of event a/b (flA/f2B) c/d (flC/f2D) Adapted from Schlesselman (1982). * Random samples of f1*100 percent of Nl affected and f2*100 percent of N2 unaffected sites are taken from the target population presented in Table 2.6.B. Table 2.7.B. Numerical Example of a case control study Landslide SL35 Present Absent 44 slope angle > 35° 44 57 168** slope angle <35° 212** 101 Total 0.26 Odds of exposure 0.77 Total 88 225 313 1.00 0.34 f1 = 1.0 and f2=0.008 applied to Table 2.6.B. ** rounded to nearest integer The proportion of cases among exposed sites in the sample is: pi =a/(a+b)=flA/(flA+f2B) (Eqn. 2.8) , which is generally not equal to the population incidence rate (pi). p2 '=c/(c+d)=flC/(flC+f2D) (Eqn. 2.9) represents the proportion of cases among unexposed sites which is generally not -41 - equal to p2, thus the ratio Pl7p2' will not estimate the population relative risk (pi/p2) correctly. Unless the ratio of the sampling fractions (fl/f2) is known or unity, the population incidence rates pi and p2 cannot be determined from a casecontrol study. In most instances there are only a few cases compared with the number of potential controls in studies of rare events. Such studies usually choose a value of fl close to unity, whereas the value of f2 is chosen close to zero. Thus the ratio of fl/f2 is usually a very large number (Schlesselman 1982). Random sampling in a case-control study implies that, among cases and controls, an exposed site has the same chance of being selected as an unexposed site. Thus sampling on the basis of a hazardous event does not affect the odds of exposure (Table 2.7.B). The sample estimates of the odds of exposure for the cases and controls agree with the population values. The odds of event are altered by the ratio of fl/f2 in a case-control study compared with the population values. The sampling ratio cancels on division in the calculation of the odds ratio, so the sample estimate of *F will agree with the population odds ratio apart from sampling variability (Schlesselman 1982). ¥ = ( f l A / f 2 B ) / ( f l C / f 2 D ) =AD/BC = 1.00/0.34=2.95 ccs (Eqn. 2.10) In the practice of biomedical research all eligible cases within a geographic area are selected. Since the number of eligible controls greatly exceeds the number of cases available for the study with rare events, sampling considerations are more concentrated on the selection of the control group. Sampling is carried out from a frame, which is a list of all potentially eligible cases and controls in the target population. There are many sampling procedures that may be used. The simplest one is simple random sampling that refers to the method of selecting subjects from a frame such that each possible sample of size of n has an equal probability of being selected. Stratified sampling involves the selection of subject at random from subgroups (strata) of the target population. A -42- sample of different size from each stratum may be selected, so that the total sample size is made up of subgroups of possibly differing size. In epidemiology, matched sampling (matching) is used extensively. This method involves the pairing of one or more controls to each case on the basis of specified variables, the effects of which one wants to eliminate from the case-control comparison. More complex sampling techniques, such as multistage sampling may be considered. Interested readers may consult Cochran (1977) for details on various sampling methods. The purpose of any sampling procedure is to avoid biased selection. In terms of case-control sampling, each eligible case and control in the target population, irrespective of exposure, should have an equal chance of being selected. With a properly designed and implemented sampling procedure, an unbiased selection can be assured. 2.3.1.4. Sample size calculation in a case-control study The number of subjects to be selected for the study is crucial in planning a case-control investigation. A sufficiently large sample helps to avoid two sources of errors: 1. Claiming that a variable is associated with landsliding when in fact it is not (Type I or alpha error) 2. Claiming that the independent variable is not associated with landsliding, when in fact it is (Type II or beta error) Basically, there are four values to specify in order to decide how many subjects should be selected for the case-control study: 1. Relative frequency of exposure among controls (prj) 2. Hypothesized relative risk that is worth detection (R') -43- 3. Level of significance (a) is the probability of accepting the alternative hypothesis (Hg) that the exposure variable(s) is associated with landsliding (H :R*l) when the null hypothesis (HQ:R=1) is true A 4. Power of the study (1-0) (Schlesselman 1982) This section provides the formulae for calculating the required number of cases for a case-control study. For the underlying assumptions, hypotheses and detailed discussion the reader should refer to Schlesselman (1982). The number of cases is calculated as follows: n = p*q*(\ + \lc)*[Z + Z _ fl(p -p f a {x P) x 0 (Eqn. 2.11) with = *R/\l Pl Po + p *(R'-l)] (Eqn. 2.12) 0 p=( +p*c)/(l + c) (Eqn. 2.13) q = \-p (Eqn. 2.14) Pl where R : hypothesized relative risk associated with exposure that is worth detection 1 c: case-control ratio Za and Z(l-(3): values from standard normal distribution p(j: proportion of exposure among controls For a one-sided test (H :R>l or H :R<l) Za is taken to be the value of the A A standard normal distribution that is exceeded by the probability of a. For a twosided test (H :R*l) Z is taken to be the value that is exceeded with probability of A a a/2. -44- 2.3.1.5. Calculation of the power of the study Power calculation in a case-control study is given by the following formulae: Z<h> = >*(P,-^o) /[/>*(l-Jp)*(l + l/c)] - Z 2 a (Eqn. 2.15) where the notation agrees with the previous section. The power of the study is determined from the standard normal distribution by finding the probability with which the calculated value of Z(i_p) is not exceeded: Power = P(Z<Z(i_p)) (Schlesselman 1982). 2.3.1.6. Confounding The term confounder will be used in this study as it is used by epidemiologists, that is, to describe an extraneous variable associated with both the outcome variable of interest (y) and a primary independent variable (x) usually called risk factor. An apparent association between an exposure (independent) variable and the outcome variable may actually be due to an extraneous variable. Suppose we want to investigate the relationship between a specific soil type and landsliding by a casecontrol study. Comparing case and control groups, we might find that the case group contains a greater proportion of sites in the transient snow zone, which is an elevation band where the snowpack is not permanent throughout the winter and rainon-snow events are common, than does the control group. Transient snow zone may be associated with landsliding due to its hydrologic characteristics, and it might also be correlated with certain surficial materials. Thus the apparent increased risk of landsliding found to be associated with soil type in fact may be due to being in the transient snow zone. The determination of the confounder status involves two criteria: -45- 1. the extraneous variable should be associated with the outcome variable and 2. it should also be associated with the exposure variable of interest In the preceding example, transient snow zone is a confounder with respect to an association between soil type and landsliding. In practice, the confounder status of an extraneous variable is ascertained by comparing the estimated effect of an exposure variable from models containing or not containing the extraneous variable. Any logically and geologically important change in the estimated effect of the exposure variable would dictate that the extraneous variable is a confounder and should be included in the model. Table 2.8 presents the result of a case-control study that reports the presence or absence of a certain surficial material (SOIL2bedrock outcrop) among 101 cases of landsliding and 264 controls. Table 2.8. The relationship between bedrock outcrop(SOIL2) and landsliding (SLIDE) in the study area SLIDE SOIL2 Cases Bedrock outcrop (1) 28 Other (0) _73 101 Controls 43 221 264 Estimated Risk I by the odds ratio ¥=1.97 I The odds ratio calculated by Eqn. 2.10. is ¥ = (28x221)/(43x73) = 1.97 which means that the relative risk is estimated to be about two times higher at sites with bedrock outcrop as compared with sites with other surficial materials present. Suppose that now we stratify the cases and controls based on whether the location of the study site is in the transient snow zone (TRANS) or not. The result can be seen in Table 2.9. -46- Table 2.9. Multi-way cross classification table for examining the relationship between SOIL2 and SLIDE stratified by TRANS TRANS 0 1 SOIL2 Bedrock outcrop Other Subtotal Bedrock outcrop Other Subtotal Total SLIDE Cases 19 6 25 9 67 76 101 Estimated Risk by odds ratio (19x85)/(6x39) = = 6.90 Controls 39 85 124 4 136 140 264 (9xl36)/(67x4) = = 4.57 The transient snow zone specific relative risks exceed 1.97 in both instances which indicates that the overall estimate ¥ = 1.97 is spuriously low. Table 2.10 indicates that while controls are approximately uniformly distributed over the two transient snowzone groups (1 and 0), cases are more frequent in the transient zone (75 %). Furthermore, bedrock outcrops are far less frequent in the transient zone, for both cases and controls. Table 2.9 and 2.10 indicates that ignoring the variable TRANS would bias the analysis since there are proportionally more cases than controls in the transient snowzone. TRANS itself appears to be a risk factor for landsliding (Table 2.11). Table 2.10. Distribution of cases and controls by TRANS from Table 2.9 Percent at each instance TRANS 0 1 Cases Controls % % 24.8 (25/101) 75.2 (76/101) 47.0(124/24) 53.0 (140/264) 100.0 100.0 TRANS specifi c percentage of SO [L2 Controls Cases % -47- 76 (19/25) 12 (9/76) % 31 (39/124) 3 (4/140) Table 2.11. The relationship between transient snowzone (TRANS) and landsliding (SLIDE) in the study area SLIDE TRANS 1 0 Cases 76 _25 101 Controls 140 124 264 Estimated Risk by the odds ratio ¥=2.69 The estimated relative risk, adjusted for TRANS can be calculated by the Mantel-Haenszel method (Mantel and Haenszel 1959.). The Mantel-Haenszel estimator (M-H estimator) can be obtained as a weighted average of the stratum specific odds ratios SV[ = (ai/bi)/(q/di), where aj, bi, q and dj are the observed cell frequencies for stratum i. For example in Table 2.9 arj = 19, bo = 39, co = 6 and do = 85. Denote the total number of subjects within stratum i as nj. In the example nQ = 149 (25 + 124). The M-H estimator can then be calculated as follows: V = v (Eqn. 2.16) Using the data in Table 2.9 yields the value of 5.88 for the Mantel-Haenszel odds ratio suggesting a sixfold increased risk among sites with bedrock outcrop as compared with sites where bedrock outcrops are absent. It is important to note that the M-H estimator provides a correct estimate only if the odds ratios are constant across the strata.. The simplest and most easily computed test for the homogeneity of odds ratios is the Woolf s test (Woolf 1955.). This test is based on the weighted sum of squares deviations of the stratum specific log-odds ratios (ln^j) from their weighted mean. Table 2.12. presents the estimated odds ratios for each stratum (^j), the log-odds ratios and their variance (var(ln¥j)) and a weight (w{) which is the inverse of the variance for this example. -48- Table 2.12. Calculation of the Woolf s test statistic from Table 2.9 TRANS = 1 4.57 1.5195 0.3834 2.6083 TRANS = 0 6.90 1.9315 0.2567 3.8955 InCFj) varDn(Ti)]* ** Zwj 6.5038 * var[ln( I )] = 1/3J+ l/bj+ l/q+ 1/dj ** w =l/var[lnCi' )] v / i i i The test statistic by Woolf is calculated as X =S{wiDn(¥i)-ln(¥ )]2} (Eqn. 2.17) ¥ =exp[Iwiln( Pi)/Swi] (Eqn. 2.18) 2 w where x w The statistic x has a chi-square distribution with k-1 degrees of freedom, 2 where k is the number of strata. One would reject the null hypothesis HQ: = ..... = Tfc if x > X^(k-l)- Using the example above: ^ 2 test for homogeneity yields a chi-square value of x w =¥2 = 5.85. The Woolf s = 0.2652 with 1 degree of 2 freedom (k = 2) which is far from significant (p = 0.60). One can state that the odds ratios are homogeneous across strata. It should be noted that the p value calculated from the chi-square distribution will be accurate only when the sample sizes are not too small within each stratum (Hosmer and Lemeshow 1989). This condition holds in the above example (nrj = 149, n\ = 216). The adjustment for the confounder variable may be carried out by stratification, matching or multivariate analysis discussed later in this section. The study design and analysis should be planned to either assess or eliminate the effects of the confounder variables. The degree of confounding can be assessed by the ratio of the "crude" estimate of the relative risk (odds ratio) based on the analysis ignoring the confounding variable (¥c) and the adjusted estimate, using, for example, the Mantel-Haenszel method (^M-H)- F ° -49- rm e previous example ^c = 1.97 and ^M-H = 5 8 8 - The r a t i o 1.97/5.88 = 0.34 indicates that the crude odds ratio underestimates the adjusted value by 66 percent. 2.3.1.7. Interaction Sometimes the association between an exposure variable and the outcome variable differs, or depends in some way on the level of another variable. Epidemiologists use the term effect modifier to describe an extraneous variable that interacts with an exposure variable. Generally, if the joint effect of two factors is significantly higher or lower than the effect of either one of them, then interaction is present. If the joint effect is less than the individual effects (negative interaction) the term antagonism is used. In case of positive interaction, the joint effect exceeds the individual effects, and one may say that there is synergism. (Schlesselman 1982) The assessment of joint effects of two or more variables is usually done by creating subgroups and comparing the various subgroups against a reference group. The reference group is usually the one with minimum levels of exposure on all of the variables under investigation as risk factors for the hazardous event. As an example, consider the investigation of the joint effect of elevation zones and slope angle. For the sake of simplicity, slope angle was dichotomized as 0-1 using the cutpoint of 35 degrees (SL35). This value is generally considered to be a threshold value in the development of debris flow-debris avalanche soil mass movements in the literature (Rollerson 1992; Howes 1987). Both slope angle and elevation zones may be considered as risk factors in the development of a landslide event. Table 2.13 represents the data that provide the basis of further discussion. -50- Table 2.13. Relationship between elevation zones (ELEVZ) and slope angle (SL35) with landsliding (SLIDE) ELEVZ (0-450 m) SL35 1 1 0 (450-900 m) 2 1 0 (900m -) 3 1 0 SLIDE Cases 2 2 29 47 13 8 Odds ratios Controls 2 35 24 115 18 71 ¥1=17.5 ¥ =2.96 2 ¥3=6.41 In Table 2.13, the odds ratio estimates the relative risk of a landslide event associated with slope angle for each of the three elevation groups. The odds of "event" among sites with slope angle greater than 35 degrees is expressed relative to the odds of "event" among sites with slope angle less than or equal to 35 degrees. Since the reference group of slope <35° changes for each of the three elevation categories, the odds ratios are "subgroup-specific". For example ¥ = 6.41 means that in elevation zone 3, the relative risk of landsliding is 6.4 times higher on steep slopes (>35°) than on gentle slopes. Taking as a reference group grid cells which have slope angle<35° and elevation group 1 (0-450 meters), the odds of event for each of the other combinations may be expressed by a ratio relative to the reference group (Table 2.14). This representation allows one to compare the different odds ratios since all of them are based on the same reference group. Sites in the transient snowzone (ELEVZ = 2) with slope angle < 35° are estimated to have 7.2 times greater risk of experiencing a landslide than those sites that are in the rain dominated zone (ELEVZ = 1) and have slope angle < 35°. Sites in the transient snowzone with slope angle >35° are estimated to have 21.2 times increased risk as compared with -51 - the reference group. The analysis of the individual and joint effects of 3 or more risk factors proceed exactly as above using a single subgroup as the reference, expressing the odds of event for all other combinations relative to it. The assessment of the joint effect of individual risk factors can be carried out comparing the subgroup odds ratios. Woolf s test cannot be applied since the odds ratios are not independent, i.e. they all are based on the same reference group. Various statistical tests for interactions are available in the literature. These tests are provided or can easily be implemented with the major statistical packages such as SAS (SAS Institute Inc. 1984) or SYSTAT (Wilkinson 1990). In this section only one test will be presented. For detailed discussion of the various tests the reader may consult with the relevant biostatistical textbooks. Table 2.14. Odds (in parentheses) and odds ratios for landsliding associated with various combinations of elevation zones and slope angle classes ELEVZ 0-450 m 450-900 m 900m- SL35 <35° (2/35) 1.0* (47/115) 7.2** (8/71) 2.0 >35° (2/2) 17.5 (29/24) 21.2 (13/18) 12.6 * Reference group ** Odds ratio calculated as (47/115)/(2/35)=7.2 For testing the significance of the joint effect of two risk factors, consider the following example. Tables 2.15 and 2.16 represent the analysis of slope angle (SL35) and drainage conditions (DR01). Both variables are dichotomized, DR01 = 1 represents poor drainage conditions and DR01 = 0 represents well drained sites. SL35 represents the variable used in previous examples. -52 - Table 2.15. Relationship between drainage conditions (DR01) and slope angle (SL35) with landsliding (SLIDE) SLIDE DR01 1 Cases 43 52 1 5 SL35 1 0 1 0 0 Controls 37 180 7 40 Odds ratio ¥=4.02 ¥ = 1.14 An approximate chi-square test of significance of interaction (synergy) can be carried out as follows. First, the odds of event should be defined for the reference group (a()/bo), for the first variable only (aj/bi), for the second variable (a /b2) 2 and for the joint effect (a^/b^). According to the example: ao/bo=5/40 =0.1250 ai/bi = 1/7=0.1429 a2/b2=52/181 =0.2873 a3/b3 =43/37 = 1.1622 Table 2.16. Odds (in parentheses) and odds ratios for landsliding associated with various combinations of drainage conditions and slope angle classes SL35 DR01 1 0 <35° (5/40) 1.0* (52/180) 2.31 >35° (1/7) 1.14 (43/37) 9.30 * Reference group -53- Let S represent synergy (positive interaction): S=(a3/b3+ao/bo)-(ai/bi +H2^2) (Eqn 2.19) Let vj denote the estimated variance of aj/bj: vi=var(ai/bi)«(ai/bi) (ai+bi)/(aibi) 2 (Eqn. 2.20) and let V denote the sum of the estimated variances: V=vo+vi+v2+v3 (Eqn. 2.21) A chi-square test with 1 degree of freedom can be used to test whether S is significantly greater than zero or not: X (1)=S2/V (Eqn. 2.22) 2 The previous example yields x 7.59 with 1 d.f. (p=0.006) suggesting that there 2= may be a significant joint effect between slope angle and drainage conditions. 2.3.2. Logistic regression analysis In many cases, there are numerous factors to control or adjust for. In such situations a case-control study using cross-classification (sometimes called contingency) tables may not be applicable. If we consider a landslide study that examines 10 variables, a total of 1024 subgroups (210) would result from multiple cross classification tables based on simple dichotomized factors! Thus, if 5 cases and 5 controls per subgroup are desired, a minimum of 5120 cases and 5120 controls would be needed in order to avoid empty cells, and to investigate the individual and joint effects of the ten variables. Creating multiple tables is also a very tedious task especially if there are more than 5 variables involved. It will be shown that multivariate models such as the logistic regression model can be used to analyze case-control studies. -54- 2.3.2.1. The logistic regression model This section describes the logistic regression model and its use in the analysis of case-control studies. In order to follow the discussion of logistic regression, familiarity with linear regression is assumed. Logistic regression has long been used in modern biomedical research and epidemiological studies. The application of logistic regression was first used in biomedical research in the early 1960's for the analysis of the individual and joint effects of a set of variables on the risk of disease (Cornfield et al. 1961). The logistic regression model is being applied in several other areas such as the estimation of probabilities from meteorological data, marketing problems involving the investigation of the probability of purchase, given the values of independent variables characterizing the customers. It is also applied to socio-political problems, involving the estimation of the probability of certain actions taken by an individual, given his/her social, political and economic status (Walker and Duncan 1967). The logistic model specifies that the probability of an event depends on a set of variables xj, x , ... ,Xp in the following way: 2 fii+fi\X +...+fi x l P(y) pfE-W m r p = -gpzix=^ (Eqn. 2.23) The variable E denotes either the presence (E = 1) or absence (E = 0) of an event, and x denotes a set of p variables x = (xj, x , ... ,Xp). The betas (P's) are 2 parameters that represent the effects of the independent variables on the probability of an event. Figure 2.4. shows the general shape of the logistic function which is known as a sigmoid (S-shaped) curve. -55 - The transformation of P(y) is called the logit transformation: g(x) = \n[P(y) I {1 - P(y)}\ = /W , * , +• • -+fi,x f (Eqn. 2.24) The major differences between the logistic model and the linear model can be summarized in three points (Hosmer and Lemeshow 1989): 1. Outcome variable is dichotomous 2. P(y) must be bounded between 0 and 1 3. The distribution of errors is binomial and not normal These properties of the logistic regression model require methods that are different from the linear model. In linear regression the method used most often to estimate the unknown parameters (betas) is called the method of least squares. In the method one chooses those values of betas which minimize the sum of squared deviations of the observed values from the predicted values based on the model. Unfortunately due to the properties of the logistic model the method of least squares cannot be applied. The general method of estimation of parameters in the logistic regression is called maximum likelihood. The method of maximum likelihood yields values of the parameters which maximize the probability of obtaining the observed set of data. This probability is expressed as a function of the unknown parameters and is called the likelihood function. The discussion of the maximum likelihood method can be found in statistical textbooks thus it will not be discussed. However, it is worthwhile to note that the computations require iterative calculations that are best done by computer. Several statistical packages (SAS, SYSTAT etc.) provide routines for maximum likelihood estimation. As an example, consider the data given in Table 2.17. This table lists slope angle in degrees (SLOPE) and the presence or absence of a landslide event (SLIDE) for 365 subjects. The outcome variable (SLIDE) is coded with a value of zero to indicate that SLIDE is absent, or 1 to indicate that it is present in the study subject. -57- The subjects are 20x20 m grid cells from the Jamieson-Orchid-Elbow Creeks study area. The scatterplot of the data is given in Figure 2.5. In this scatterplot, all points fall on one of two parallel lines representing the absence or presence of a landslide event. There is some tendency for the individual grid cells with evidence of a landslide event to be on steeper slopes than those with no evidence of a landslide event. While this scatterplot sufficiently describes the binary (dichotomous) nature of the outcome variable, it does not provide a clear picture of the nature of the relationship between SLIDE and SLOPE angle. One method to describe the functional relationship between the outcome and the independent variable is to create groups for the predictor variable, and compute the mean of the outcome variable within each group. Variable SLOPE has been grouped into slope classes. Table 2.17 contains, for each slope class, the frequency of occurrence of each outcome as well as the mean for each class. Table 2.17. Frequency table of slope class by SLIDE Slope Class n 10°-19.99° 20°-29.99° 30°-39.99° 40°-49.99° 50°-59.99° 60°Total 61 138 121 27 14 4 365 SLIDE Absent 54 113 75 15 5 2 264 Present Mean (Proportion) 7 25 46 12 9 2 101 0.11 0.18 0.37 0.44 0.64 0.50 It appears that as the slope angle increases the proportion of grid cells with evidence of landslide occurrence increases (Figure 2.6). -58 - Figure 2.5. Scatterplot of SLIDE by SLOPE angle W^^r ^^^^^^ LLI Q Z3 (O _0> « (0 0 > < D 0 C m 0 0 0 0.00 20.00 40.00 60.00 80.00 Slope angle in degrees Figure. 2.6. Plot of the mean of SLIDE in each SLOPE class o -i 0 1 1 1 1 1 1 2 3 4 5 SLOPE classes -59- — i 6 While the scatterplot provides a valuable insight into the relationship between SLIDE and SLOPE angle in this example, a functional form for this relationship is needed. The shape of the curve indicates that the relationship between the logit g(x) and the continuous variable (SLOPE) is not linear since a linear logit would produce a curve as presented in Figure 2.4. Thus, a function that takes on values between 0 and 1 and that is capable of assuming such a shape is required to effectively model this kind of relationship. The logistic function quadratic in its argument satisfies these requirements. Under this model the relationship between slope angle (x) and the probability P of occurrence of a landslide event is given by: P{y) = I (q- - > £ gflo+Pix+Pn* 2 + n 2 25 where PQ, Pi and Pi\ are the unknown parameters. Table 2.18 shows the output of LOGISTIC Version 3.11Ef of STATOOLS™ free software package (Dallal 1988) using the continuous independent variable SLOPE and the dichotomous outcome variable SLIDE. Table 2.18. Results of fitting a logistic regression model to the data in Table 2.17. VARIABLE SLOPE SLOPE CONSTANT 2 Estimated coefficient Standard error 0.1858 0.0607 -0.0016 0.0008 -5.0103 1.0889 The maximum likelihood estimates of PQ, Pi and P n are PQ' = -5.0103, 0.1858 and Pn' =-0.0016. The fitted values are given by the equation: ^ ~ i , s.om+o.isss* SLOPE-o.ooi6* SLOPE 2 1+e -60- (Eqn. 2.26) Pi' and the estimated logit g'(x) is given by the equation: g'(x)=-5.0103 +0.1858*SLOPE-0.0016*SLOPE 2 (Eqn. 2.27) 2.3.2.1.1. Testing for significance of the coefficients After estimating the coefficients, assessment of the significance of the variables in the model is needed. This usually involves the assessment of whether the inclusion of a variable significantly contributes to the model or not. The guiding principle is to compare the observed values of the response variable to predicted values from models with and without the variable in question. In logistic regression this is done by the likelihood ratio test: G = -21n[likelihood without the variable(s)/likelihood with the variable(s)] = = -2[loglikelihood without the variable(s)-loglikelihood with the variable(s)] Most statistical packages provide the value of statistic G . Under the hypothesis that HQ: Pi = 0, G follows a chi-square distribution with k degrees of freedom, where k is the number of variables removed from the model. Table 2.19 shows three models. Model 1 does not contain any variables, only the constant term. Model 2 contains SLOPE and Model 3 contains SLOPE and SLOPE . 2 Table 2.19. Likelihood ratio test for 3 different models using LOGISTIC Version 3.11Ef of STATOOLS™ free software package Model Constant SLOPE _ 1 -0.9608 -3.1190 0.0700 2 -5.0103 0.1858 3 SLOPE - 2 -0.0016 Loglikelihood -215.2848 -197.4385 -195.3888 -61 - d.f. - G - p-value - 1 1 35.69 4.10 0.0000 0.0429 Evaluating Model 2, G is calculated by: G= -2[-215.2848-(-197.4385)]=35.6926 Using the symbol x^( ) to denote a chi-square random variable with v degrees of v freedom, the p-value associated with this test is p[x (l) > 35.69] < 0.0000. Thus we 2 can reject the null hypothesis and state that there is evidence that SLOPE significantly contributes to the model (PsLOPE * 0)> an£ ^ * therefore a significant s variable in predicting landslide events (SLIDE). Since adding SLOPE to Model 2 2 is still significant at 95 percent significance level (p = 0.04), we can state, that based on the data, Model 3 provides the "best" fit. However, there are other important factors to consider before concluding the significance of a variable. These involve examining the logical and geological importance of the variable as well as the inclusion of other potentially important variables. 2.3.3. Logistic regression and case-control studies When sampling is performed conditional on the outcome variable, as in a casecontrol study, logistic regression models can easily be applied to obtain the adjusted odds ratios from the estimated coefficients of the variables in the model. This is a tremendous advantage, since instead of creating multiple cross-classification tables, one can use the logistic model to obtain odds ratios and thus estimate relative risk. Farewell (1979) and Prentice and Pyke (1979) showed that only the coefficient Po differs in a logistic model based on the case-control sampling while the remaining parameters are unaffected. In a case control study, the binary response variable is fixed by stratification (1-0). Samples of fixed size are taken from the two strata (cases and controls) defined by the outcome variable with sampling fractions x\ (cases) and %2 (controls. -62 - There are four possible outcomes for a subject with variable x = (xi, X2,..., Xp): 1. Develop a landslide and be in the sample with probability tlP(y); 2. Develop a landslide and not be in the sample with probability (1-Ti)p( ); y 3. Remain free of landsliding and be in the sample with probability t2 l(y)» ( 4. Remain free of landsliding and not be in the sample with probability (1-T2)q( ). y Applying Bayes theorem (Neutra and Drolette 1978), the probability that the 1 response variable (E) = l given the value of x in the sample of cases and controls: P"(y) lP(y) t lP(y)+ 2q(y)] =T / T T w h e r e q(y) = "P(y) ( q - - ) A £ n 2 28 The odds of event associated with variable x=(xi, x2,..., Xp) in the sample can be written as: p"(y)/q"(y)='ciP(y)/t2q(y) where q"( ) = l-p"( ) y y (Eqn. 2.29) Applying the logit transformation gives ln[p"( )/q"( )] =ln(xi/x2)+ln[p( )/q( )] = Po* + P l x i +P2 2+--- + Pp p X y y y x y (Eqn. 2.30) where Po* ( l/ 2)" "PO ==m i; x (Eqn. 2.31) , This shows that a case-control study can easily be analyzed by any statistical package using logistic regression. Parameters P j , P2, Pp can be used to obtain odds ratios. Since parameter PQ depends on the sampling fraction ratio x\/t2> P(y)> which depends on the parameter PQ, cannot be estimated unless we know the value of x\/x2- The odds of event in the sample of cases and controls is p"(y)/q"(y) that is closely related to the odds of event of the target population (P(y)/q(y)) by the factor -63- of ti/t2. (Please recall that the odds of event is p i / q i and the odds ratio is fr oin (Pl^lVC^^ section 2.3.1.2.) In general, the relative odds of event for a subject with x*=(xj*, Xp*) as compared with a subject with x=(xi, Xp) is given by the odds ratio: * ( X *' X ) = f{E P = lx)i% = l\x) ( E c m - 2 3 2 ) p *¥(x*:x)= e x p [ ^ ^ ( ^ ; - *,)] ( E q n 2 3 3 ) /=i Examining Equation 2.33, one can see that Po canceled out in the calculation of odds ratios. ¥ depends only on those factors for which the two subjects differ (Schlesselman 1982). Consider the example given in Table 2.8 where the relationship between surficial material (SOIL2) and landsliding was examined. Using LOGISTIC Version 3.11Ef of STATOOLS™ software (Dallal 1988) the following equation was computed from the data: -1.1077+0.6787*SO7L2 ^ ~ j _i_ g-i.io77+o.6787»so/L2 (Eqn. 2.34) SOIL2 represents the presence (xj = 1) or absence (xj =0) of bedrock outcrop on a site. Assuming that two sites are "equal" on any other factors, then ¥'=exp[Pi(xi*-xi)]=exp[Pi(l-0)]=exp(Pi) (Eqn. 2.35) In our example the crude (unadjusted) odds ratio ¥ ' = exp(Pj) = exp(0.6787) = 1.97 which agrees with the value calculated from the 2x2 cross-classification table in Table 2.8. Using TRANS (transient snowzone) as a stratifier yields the equation: -64- Adjusting for the variable TRANS gives the stratified estimate of the odds ratio 4" = exp(1.7737) = 5.89. This value is the maximum likelihood estimate of the "adjusted" odds ratio, and it agrees with the Mantel-Haenszel estimate (^M-H = 5.88) in Table 2.9. The change in the estimate of the odds ratio from the crude to the adjusted is 1.97 to 5.89 indicating considerable confounding effect due to the variable TRANS. Assessment of the homogeneity of the odds ratios across strata is based on the likelihood ratio test involving the model without and with the SOIL2xTRANS interaction term. Table 2.20 provides the basis of the analysis. Table 2.20. Assessing the homogeneity of the odds ratios by using the likelihood ratio test Model SOIL2 +TRANS +SOIL2xTRANS SOIL2 coefficient * 0.6787 1.7737 1.9318 Loglikelihood -212.3919 -195.6960 -195.5642 G 5.7858** 33.3918 0.2636 d.f. 1 1 1 p value 0.0162 0.0000 0.6101 * Only the coefficient of SOIL2 is presented ** The log-likelihood of the base model containing the constant term only is -215.2848 Note that minus twice the change in the log-likelihood may be used to test for significance of the variables added to the model. The interaction term (SOIL2xTRANS) was created as a product of SOIL2 and TRANS. The G statistic yields the value of 0.2636. This statistic is compared to the chi-square distribution with 1 degree of freedom. The p value (0.61) suggests that the odds ratios are homogeneous across strata which also proves that there is no significant interaction -65- between S0IL2 and TRANS. The confounding effect can be assessed by looking at the model with and without the variable TRANS. The coefficient of SOIL2 changes significantly from 0.68 to 1.77 which suggests a significant confounding between SOIL2 and TRANS. -66- 3. DATA ANALYSIS AND DISCUSSION In this section a case-control study is conducted in the Jamieson-Orchid-Elbow Creeks study area. The data provided by the case-control study are then analyzed by multi-way cross classification tables and multiple logistic regression to set up a landslide risk model. Finally, a landslide risk map is produced that can be used in planning watershed management activities and in the overall risk assessment of the study area. Figure 3.1 shows the flow chart of the study procedure. 3.1. A case-control study on the study area 3.1.1. Stating the research questions The final product of the terrain stability assessment is a landslide risk map that can be used in planning watershed management activities. To assess landslide risk, the following general questions must be answered: (1) What are the causes of landslide events in the Jamieson-OrchidElbow Creeks pilot study area? (2) Is there interaction or confounding between the various factors investigated? (3) If so, what are the possible explanations? -67- Figure 3.1. Flow chart of the study procedure 68 The initial formulation of the general questions serves little purpose beyond pointing to the general problems of interest. These questions will be narrowed and made more precise as the study proceeds to evaluate aspects of one or more specific hypotheses. 3.1.2. Variable selection To answer the research questions, a decision has to be made about the variables (risk factors) to be included and also about the form that these variables should take. Information should also be collected about potential confounding variables to rule out alternative explanations for the apparent presence or absence of an association. Although logistic regression assumes a categorical dependent variable (usually dichotomous), the independent variables may be either categorical or continuous. In this study, these variables are analyzed as either categorical or ordinal in nature. Because the multi-way table of frequencies in a categorical analysis has as many dimensions as variables, one has to take care that the sample size is large enough to cover the tables defined by these variables since using categorical variables can cause empty cells in the multi-way cross classification tables. An attempt was made to keep the total number of cells to the minimum by creating dichotomized variables for the multi-way table analyses. The selection of variables was based on two criteria: 1. All the pertinent risk factors should be included or accounted for in the analysis 2. Variables that are easy to obtain from aerial photographs, digital elevation models and forest cover maps should be considered in the analysis Based on these principles and extensive literature review, the following factors were considered in the study: -69- Slope gradient: this is the most widely used topographic measurement that influences flow rates of water and sediment by controlling the rate of energy expenditure or stream power available to drive the flow. Aspect: defines the slope direction and therefore the direction of the flow. Profile curvature: is the rate of change of slope; it affects flow acceleration and deceleration and therefore influences aggradation and degradation. Irregular shapes with "bumps" on the terrain help to break the energy of downslope movements, therefore reducing the runout zone. Planfonn curvature: defines the curvature of the land surface transverse to the slope direction; it influences flow convergence and divergence. Local catchment area: is of fundamental importance to the operation of fluvial and hydrological processes in the landscape. Erosion by overland flow and shallow, rapid landsliding are common channel initiation mechanisms in humid, soil mantled landscapes (Montgomery 1994). In coastal British Columbia, overland flow rarely occurs on steep slopes due to the high hydraulic conductivity of the low density, colluvial soils and thus shallow landsliding is the primary channel initiation mechanism. The location of the channel head is associated with the relationship between slope and local catchment area (Figure 3.2). The inflection in this relationship reflects a transition from steep debris flow-dominated channels to lower gradient alluvial channels. This inflection occurs in the 50-60 percent ( 25-30 degrees) slope gradient range in the Pacific Northwest and is associated with 1-1.5 ha drainage area based on field mapping of channels in small mountainous drainage basins in Oregon and Washington (Montgomery and Foufoula-Georgiou 1993). Elevation: streamflow in the coastal watersheds is produced by a combination of rainfall, rain-on-snow and radiation snowmelt. Watersheds can be stratified into three major elevation bands: below 300-500 meters, where runoff is most frequently produced as a result of rainfall. The elevation band between 300-500 and 800-900 -70- meters characteristically has a transient snowpack in the winter, that is, the snowpack is not permanent throughout the winter, but tends to come and go depending on the meteorological factors, and rain-on-snow events are common. The most frequent peak flow events are generally produced within this elevation band. At elevations above 800-900 meters, the snowpack tends to accumulate throughout the winter. Runoff from above 800-900 meters is produced during the spring and early summer by the combination of radiation snowmelt and rain-on-snow. In the study area soil mass movements appear to be more frequent in the transient snow zone (Chatwin 1994). Figure 3.2. Catchment area-channel slope relationship for small mountainous drainage basins along the Pacific Coast slope (m/m) Adapted from Montgomery and Dietrich (1994) -71 - Surficial materials (soil texture): are classified according to their mode of origin because there is a close relationship between the process whereby they were formed and their most important physical properties. Texture refers to the size and roundness of particles that constitute a surficial material and to the sorting (uniformity of particle sizes) within a mass of sediment. Texture largely determines the physical properties of soils such as permeability, compressibility, drainage characteristics, stability and erodibility on steep slopes (Ryder and Howes 1984). Bedrock geology: the type and characteristics of bedrock are important factors that determine where bedrock failures, such as rockfalls, rockslides might occur. Bedrock also influences the type and character of terrain, soils and landform which affect the overall stability and erodibility of an area. Soil depth: the depth of impermeable layers provides a quantitative indicator of the depth and type of the failure that may occur at a site. It also indicates principal paths of subsurface water movement or zones of temporary water table development, probable surfaces of failure on the slope, and, in some instances, the depth of root penetration which can be important for identifying areas made unstable by windthrow (Chatwin et al. 1991). Drainage conditions: groundwater plays an important role in the development of slope failures. The relative amount of moisture available within the soil, and the length of time it remains, can have a significant influence on hillslope instability in all types of materials. Soil drainage classes as described in the Canadian System of Soil Classification (Expert Committee on Soil Survey 1987) may be useful as indicators of hydraulic conductivity. Biogeoclimatic zone: at the regional level, biogeoclimatic zones describe changes in climate over distances of up to a hundred kilometers. Regional climatic gradients within an area are delineated and mapped based on observations of changes in vegetation composition and structure (Acres International Ltd. 1993). -72 - Mean annual precipitation: large storms in the winter play a significant role in triggering soil mass movements in Coastal British Columbia. Topography exerts the major influence on rainfall production and distribution in the area. Hetherington (1976) developed a model that was useful for estimating precipitation on windward mountain slopes during stable storms with major precipitation, typically storms involving southwest flows of maritime tropical air masses. These are often the type of events that cause the greatest streamflow peaks when such warm, wet air masses lead to rain-on-snow in the transient snow zone. Although precipitation input plays a very significant role in triggering landslide events, it seems that annual precipitation does not provide enough information for assessing landslide hazard in smaller areas. Neither average annual precipitation nor intensity change appreciably from one location to the other within the study area (Table 2.1). The uncertainties involved make it hard to predict rainfall amount and intensity at a certain location in small mountainous watersheds. Regional terrain stability assessment may be able to utilize isohyetal maps that show zones of equal precipitation resulted from moderately frequent, long-duration storms. These storms are often accompanied by mass movements including sometimes devastating debris torrents. Disturbance history: the distribution of forest stands is a function of timing and nature of stand-replacing disturbance that has occurred. The time since disturbance determines the age and successional stage of the stand. Sites with soil mass movements usually have a history of instability. This instability may be triggered by fire or forest cover removal. The main agents of disturbance responsible for stand replacement in the study area are fire, landslides and harvesting and insect damage (Acres International Ltd. 1993). Forest cover: tree cover influences the amount and intensity of rainfall reaching the surface, the amount of water stored in the overburden, and the strength developed along a potential failure surface. Root system of trees and other -73- vegetation may increase the shear strength of unstable slope by anchoring through the mass into fractures in bedrock, tying the slope together across zones of weakness or instability. Distance from water courses: this characteristic indicates the risk that a slope failure ends up in a stream causing water quality deterioration and affecting anadromous fish habitat. Tables 3.1.A and 3.1.B show landscape attributes that have been considered to be important in the development of soil mass movements. Table 3.I.A. Variables "measured" at 20x20 m grid cells Acronym SLIDE ELEV SOILS GEOL SDEPTH DRAIN Definition Landslide present absent Elevation above sea level Surficial material Bedrock Colluvium Till Fluvial sediment Bedrock geology Plutonic rock Gabbro, quartz diorite Cretaceous rocks (Gambier Group) Soil depth Shallow Moderate Deep Drainage conditions Poorly drained (drainage class 1-2) Moderate (drainage class 3-4) Well drained (drainage class 5-6) -74- Code/Unit 1 0 m 1 2 3 4 1 2 3 1 2 3 1 2 3 Table 3.1.A. cont. Acronym Definition SUCST Successional stages Initial (0-2 yrs) Shrub/Herb regeneration (3-15) Pole-sapling forest (16-40) Young forest (41-100) Mature forest (101-250) Oldgrowth forest (>250) Not vegetated BGCZ Biogeoclimatic zones Mountain Hemlock Coastal Western Hemlock Disturbance history DISTUR Harvesting Landslide Fire Flood Gap replacement Code/Unit 1 2 3 4 5 6 7 1 0 1 2 3 4 5 Table 3.I.B. Variables created by transformation Acronym SLOPE ASPECT PROFC PLANC Definition Slope of ground surf, (from ELEV) Aspect (from ELEV) Profile curvature (frm ELEV) Concave Straight Convex Horizontal curv. (from ELEV) Concave Straight Convex Please see Appendix A -75- Code/Unit degrees degrees Software IDRISI IDRISI form, for 1 1 2 3 form, for 1 2 3 Table 3.I.B. cont... Acronym Definition DISWAT Distance from water bodies ELEVZ Transient snow zone (from ELEV) < 450 metres 450-900 metres > 900 metres Dichotomized variable from ELEV TRANS 0-450 m and > 900m 450-900 m SLCL Slope classes from SLOPE 10-25 degrees 25-40 degrees 40-55 degrees over 55 degrees SLOCL Slope classes from SLOPE 10-25 degrees 25-40 degrees 40-55 degrees over 55 degrees Dichotomized variable from SLOPE SL35 < 35 degrees > 35 degrees ASPECT4 Categorical variable from ASPECT North East South West Dichotomized variable frm ASPECT ASP2 Western Eastern Dichotomized variable frm SD SDEPTH Moderate and deep Shallow -76- Code/Unit m Software IDRISI IDRISI 1 2 3 IDRISI 0 1 IDRISI 1 2 3 4 IDRISI 17.5 32.5 47.5 65.0 IDRISI 0 1 IDRISI 1 2 3 4 IDRISI 0 1 IDRISI 0 1 Table 3.I.B. cont.. Acronym Definition CATCH Catchment area Dichotomized variable frm CATCH CA30 over 30 pixels (1.2 ha) 30 pixels or less Dichotomized variable from DRAIN DR01 Well drained site (Class 2-3) Poorly drained site (Class 1) SOIL2 Dichotomized variable from SOILS Bedrock outcrop Other type Design variable for SOILS SOI SOILS = 2 otherwise Design variable for SOILS S02 SOILS = 3 or 4 otherwise PL01 Dichotomized variable from PLANC PLANC= 1 otherwise Dichotomized variable from PROFC PR01 PROFC= 1 otherwise Design variable for PLANC PL1 PLANC= 2 otherwise Design variable for PLANC PL2 PLANC = 3 otherwise Design variable for PROFC PR1 PROFC = 2 otherwise Please see Appendix B 2 -77- Code/Unit Software no. of cells catch, for IDRISI 0 1 IDRISI 0 1 IDRISI 1 0 IDRISI 1 0 IDRISI 1 0 IDRISI 2 1 0 IDRISI 1 0 IDRISI 1 0 IDRISI 1 0 IDRISI 1 0 Table 3.I.B. cont... Acronym Definition PR2 Design variable for PROFC PROFC= 3 otherwise Dl Design variable for DISTUR DISTUR= 2 otherwise D2 Design variable for DISTUR DISTUR = 3 otherwise Design variable for DISTUR D3 DISTUR = 4 otherwise D4 Design variable for DISTUR DISTUR = 5 otherwise Code/Unit Software IDRISI 1 0 IDRISI 1 0 IDRISI 1 0 IDRISI 1 0 IDRISI 1 0 3.1.3. Data collection and generation Most of the data are available in digital form thanks to the Ecological Inventory Pilot Study (EIPS) that was carried out on the study area in 1992. The purpose of the EIPS was to develop an inventory method based on biogeoclimatic classification principles (Acres International Ltd. 1993). Table 3.2 shows the data types along with their sources. The data were available in TERRASOFT 10.3c Geographical Information Systems (GIS) software. 20x20 m raster layers were created then the information was imported into IDRISI v.4.1. GIS software. IDRISI running on a 486-DX microcomputer was used to integrate spatial information for the project. -78- Table 3.2. Data types and sources used in the study Factor landslide occurrence slope gradient aspect profile curvature horizontal curvature catchment area elevation surficial material bedrock geology soil depth drainage conditions biogeoclimatic zones mean annual precipitation disturbance forest cover Source EIPS, field reconnaissance, literature Digital Terrain Model (DTM) DTM DTM DTM DTM TRIM* 1:20, 000 base map EIPS EIPS, literature EIPS EIPS, field reconnaissance EIPS GVWD** prec. records EIPS GVWD 1:20, 000 Arc/Info file *TRIM: Terrain Resource Inventory Mapping digital topographical database ** GVWD: Greater Vancouver Water District IDRISI is a raster-based geographic analysis software package developed, distributed and supported by The IDRISI Project, a non-profit organization within the Graduate School of Geography at Clark University in Worcester, Massachusetts. IDRISI is designed to be affordable and easy to use, yet provide professionallevel analytical capability on DOS-PC platforms. IDRISI provides an extensive suite of tools for image processing, geographic and statistical analysis, spatial decision support, data display, and import/export and conversion. IDRISI had been chosen for a number of reasons: 1. affordable 2. the analysis is to be carried out at grid cell level -79- 3. easy to use 4. provides sufficient analytical tools on DOS-PC platform 5. modules written in any programming language can easily be integrated into the IDRISI system. These characteristics made IDRISI an excellent choice for the study. In IDRISI, a digital terrain model had been created that was used to calculate topographical indices such as slope, aspect and local catchment area. Programs were written in FORTRAN 77 programming language to get profile and horizontal curvatures as well as catchment area belonging to any grid cell using 3x3 altitude matrices (See Appendices A and B). The information on selected factors was extracted from IDRISI image files using simple routines written in FORTRAN to create a database. 3.1.4. Study design The logistic regression model using case-control sampling was chosen for the analysis. The purpose of the analysis was twofold: 1. to examine the effect of various landscape attributes on the development of shallow, rapid soil mass movements and 2. to develop a procedure that can be used to estimate the relative risk of landsliding as a function of a relatively large number of independent variables The study subjects are 20x20 m grid cells that were created using an equidistant raster placed over the study area. Since most of the data sources were digitized from 1:20,000 maps, 20 m represent a sufficient resolution for the study. The observed values of the dependent variable are 1 or 0 corresponding to the presence or absence of landslide initiation point within the grid cell. Due to the large number of cells (over 30,000) and the relatively small number of cells with landslide initiation points (101), a case-control sampling was carried out. The case-80- control method is suited to the study for a number of reasons. This method has long been used to study rare diseases in epidemiology in biomedical research. Since landslides, like some diseases, are rare phenomena, the case-control method seemed to be applicable. Most of the variables were chosen in such a way that they might easily be obtained from aerial photographs, topographical maps, soil maps, forest cover maps and digital elevation models. Most of these sources are available in digital form in British Columbia, so there was no need for digitizing which would increase the cost of the analysis. With the smaller sample obtained from a casecontrol study, the analysis could be carried out on a microcomputer. 3.1.5. Case definition and selection Defining a case is crucial to the case-control study. It involves two distinct specifications: 1. Establishment of objective criteria for the definition of the hazardous event under study 2. Selection of study subjects based on eligibility criteria The Jamieson-Orchid-Elbow Creeks subdrainage is characterized by shallow, rapid soil mass movements. Debris avalanches and debris flows dominate mass erosion processes on terrain characterized by steep slopes, cohesionless soils and relatively competent bedrock. The term landslide used throughout this study describes the displacement and downslope movement of soil, rock and organic material due to gravity and/or excess of water. A study subject (20x20 m grid cell) was considered to be a potential "case" if there was a sign of landslide initiation recognized from aerial photographs or field reconnaissance. All of the 101 cases within the study -81- area were selected. Complete sampling of cases is a common practice in casecontrol studies dealing with rare events. Eligibility criteria were established to restrict the study to subjects that were potentially at risk of "exposure". These subjects define the target population, that is a set of the general population that is both at risk of the study factors and the development of a landslide event. To restrict the study to sites that were potentially at risk, all grid cells that have slope angle less than or equal to 10 degrees were excluded from the target population. This decision was made because slope angle is considered the most important factor in the development of the debris avalanche-type soil mass movement. The data also supported this decision since none of the 101 cases could be found on slopes less than 15 degrees. 3.1.5.1. Sources of cases Potential cases can be identified from different sources: 1. Existing landslide inventories 2. Aerial photographs (satellite imagery) 3. Existing slope stability maps 4. Field reconnaissance 5. Maps, reports, and surveys Table 3.3 presents the sources of data and the method of data collection used in this study to select eligible cases. -82- Table 3.3. Source of cases and method of data collection Source EIPS digital database Method Aerial photograph interpretation + ground truthing, AP-190 analytical plotter was used to transfer points and lines into Terrasoft GIS (Acres International Ltd. 1993.) Thurber Engineering Ltd. 1991. Literature review Geotechnical Assessment of 1990-91 landslide events in GVWD watersheds Briere, D. 1979. The stratifica- Literature review tion of forested landscapes for intensive management: development and application Thompson, G. and Golding, Literature review D.L. 1990. Mass wasting in the GVWD watersheds: an inventory of events and an evaluation of the Aqua-Terra Classification System. Field reconnaissance ground truthing 3.1.6. Control definition and selection 3.1.6.1. Definition A control group was used to compare the history of exposure to the risk factors in the cases with that in subjects that were free of landsliding. - 83 - Controls are intended to provide an estimate of the exposure rate to the risk factors that would be expected to occur in the cases if there were no association between the event under study and the risk factors (Schlesselman 1982). 3.1.6.2. Eligibility criteria and sources Controls should be similar to the cases in regard to potential exposure of risk factors. All grid cells that have slope angle less than or equal to 10 degrees were excluded from the analysis. This criterion restricted the number of potential controls to 26,430 grid cells. A decision was made to use one percent of the eligible controls in the study. I selected 264 controls randomly from the study area using IDRISI v.4.1. sample module. This provided an approximate 1:2.6 case-control ratio for this study. Statistical tests were used to check whether the sample size and the power of the study is sufficient to carry out the analysis. 3.1.7. Calculation of sample size Suppose that approximately 35 percent of the grid cells free of landsliding have an exposure to the main risk factor (slope angle over 35 degrees). Take this as an estimate of the expected rate of exposure among controls (po = 0.35), assuming 2.6 controls per case. Let R = 2 be the minimum relative risk that one would like 1 the study to detect with a = 0.05 (1-sided test) and |3 = 0.10. Using equations 2.11 through 2.14: pi =0.35*2/[l+0.35*(2-l)] =0.5185 p=(0.5185+2.6*0.35)/(l+2.6)=0.3968 q=l_p=0.6032 n=0.3968*0.6032*(l + l/2.6)*(1.64+1.28) /(0.5185-0.35) =99.5 2 -84- 2 Thus, n=100 cases and 2.6*n=260 controls are required to carry out the study with a power of 90 percent. Thus, the 101. cases and a random sample of 264 controls selected for the study seem to be efficient for carrying out the analysis. 3.1.8. The power of the study Table 3.4 shows a dichotomized variable on slope angle which takes the value of one in grid cells with slope greater than 35 degrees, and zero otherwise. The variable SLIDE represents the presence (1) or absence (0) of landslide initiation points in a grid cell. Table 3.4. Relationship between slope angle and landsliding based on a dichotomized slope variable SLOPE Exposure > 35 degrees < 35 degrees Total SLIDE Cases 69 32 101 Controls 97 167 264 a=0.05 n=101 Za=1.64 (1-sided) The power of the study for detecting a relative risk of R' = 2 is calculated as follows: c= 264/101= 2.6139 P0= 97/264= 0.3674 Pl= 0.3674*2.0/[l+0.3674*(2.0-l)]= 0.5374 p= (0.5374+2.6139*0.3674)/(l+2.6139)= 0.4144 Z _ = Vl01*(0.5374 - 0.3674) /0.4144*(l - 0.4144)*(1 +1/2.6139) -1.64 = 1.31 2 n A -85- P(Z<1.31) = 0.90, thus the study has an estimated 90 percent chance of detecting a 2 times increase in the risk of landsliding. Since the purpose of the study was to detect an increase in the relative risk of landsliding, the alternative hypothesis Ha:R> 1 (one-sided test) was used. The two-sided test (Ha: R * 1) yields p = 0.84 that indicates a 84 percent chance of detecting a two times increase or decrease in the risk of landsliding. 3.2. Analysis of the case-control study The analysis of the case-control study was carried out in three steps: 1. Statistical display using 2x2 contingency tables 2. Multi-way cross classification tables 3. Stepwise multiple logistic regression First, 2x2 contingency tables were used to provide an overview of the relative risk of the individual variables. Next, relationships between variables were explored by simple cross tabulations. Measures of association, such as crude and adjusted odds ratios as estimates of relative risk were computed. Finally, a stepwise multiple logistic regression was used to control for many potential confounding factors, that could not be achieved by using cross-classification tables. 3.2.1. Statistical analyses using 2x2 contingency tables In the first step of the analysis of the data shown in Appendix C, twodimensional tables were created. In each table a dichotomized factor was paired - 86- with the dependent variable SLIDE that represents the presence (1) or absence (0) of landsliding in the particular grid cell. The odds ratios were calculated to get a crude estimate of the relativeriskof landsliding for the individual factors (Tables 3.5-3.14). The description of the variables is in Table 3.1 in section 3.1.2 "Variable selection". Table 3.5. Relationship between slope angle (SL35) and landsliding (SLIDE) SL35 1 0 Total SLIDE Cases 44 57 101 Controls 44 220 264 Chi-square* p-value 28.89 0.0000 odds ratio: 3.86 lower limit:** 2.32 upper limit: * * * 6.42 Total 88 277 365 The Pearson Chi-square value was computed by Epi Info v.6.0 free software that was developed by the Centers for Disease Control and Prevention, Delaware Division of Public Health (Columbier et al. 1994). It can be calculated by = [(ad-bc) n]l(p n m m^) using the notation given in Table 2.7.A. The value 28.89 indicates that the sample odds ratio CP=3.86) is significantly greater than unity. 2 x The lower confidence limit was calculated by ¥ confidence level. 2 x = ¥ e x p [ - Z / v a r ( l n F ) ] using 95 percent v L aX *** The upper confidence limit was calculated by % = Y exp[ + Z -Jvar(ln ¥ ) ] using 95 percent confidence level. a Table 3.6; Relationship between transient snowzone (TRANS) and landsliding (SLIDE) TRANS 1 0 Total SLIDE Cases 76 25 101 Controls 140 124 264 Total 216 149 365 - 87- Chi-square 14.93 odds ratio: lower limit: upper limit: p-value 0.0001 2.69 1.61 4.50 | Table 3.7. Relationship between drainage conditions (DR01) and landsliding (SLIDE) DR01 1 0 Total SLIDE Cases 95 6 101 Controls 217 47 264 Total 312 53 365 Chi-square 8.28 odds ratio: lower limit: upper limit: p-value 0.0040 3.43 1.42 8.29 Table 3.8. Relationship between surficial material (SOIL2) and landsliding (SLIDE) SOIL2 1 0 Total SLIDE Cases 28 73 101 Controls 43 221 264 Total 71 294 365 Chi-square 6.08 odds ratio: lower limit: upper limit: p-value 0.0137 1.97 1.10 3.54 Table 3.9. Relationship between soil depth (SD) and landsliding (SLIDE) SD 1 0 Total SLIDE Cases 78 23 101 Controls 160 104 264 Total 238 127 365 Chi-square 8.90 odds ratio: lower limit: upper limit: p-value 0.0029 2.20 1.30 3.73 Table 3.10. Relationship between planform curvature (PL01) and landsliding (SLIDE) PL01 1 0 Total SLIDE Cases 47 54 101 Controls 87 177 264 Total 134 231 365 -88- Chi-square 5.80 odds ratio: lower limit: upper limit: p-value 0.0160 1.77 1.11 2.83 Table 3.11. Relationship between profile curvature (PR01) and landsliding (SLIDE) SLIDE Cases 39 62 101 PR01 1 0 Total Controls 115 149 264' Total 154 211 365 Chi-square 0.73 odds ratio: lower limit: upper limit: p-value 0.3919 0.82 0.51 1.30 Table 3.12. Relationship between aspect (ASP2) and landsliding (SLIDE) ASP2 1 0 Total SLIDE Cases 83 18 101 Controls 200 64 264 Total 283 82 365 Chi-square 1.73 odds ratio: lower limit: upper limit: p-value 0.1885 1.48 0.82 2.64 Table 3.13. Relationship between biogeoclimatic zones (BGCZ) and landslidin (SLIDE) BGCZ 1 0 Total SLIDE Cases 73 28 101 Controls 165 99 264 Total 238 127 365 Chi-square 3.08 odds ratio: lower limit: upper limit: p-value 0.0794 1 1.56 0.95 2.58 Table 3.14. Relationship between catchment area (CA30) and landsliding (SLIDE) CA30 1 0 Total SLIDE Cases 17 84 101 Controls 19 245 264 Total 36 329 365 Chi-square 7.63 odds ratio: lower limit: upper limit: p-value 0.0057 2.61 1.30 5.25 From the Tables 3.5 to 3.14, it seems that all but three variables are significant at a 95 percent confidence interval. It means that they are positively associated with the generation of soil mass movements. PR01 (profile curvature, -89- dichotomized variable), BGCZ (biogeoclimatic zones) and ASP2 (aspect, dichotomized variable) are not significant at 95 % confidence level, but they will be included for further analysis due to their possible confounding effects on other variables. Looking at the calculated odds ratios for the variables, one can see that variables SL35, TRANS, DR01, SOIL2, SD, PL01 and CA30 all have a positive association with the hazardous event under study since the lower odds ratios are all higher than one. (Note that 1.00 represents no association between the independent variable and the dependent variable in a study.) Significant associations were found (p values are less than 0.05) for factors: slope angle, transient snow zone (elevation), drainage conditions, surficial material, soil depth and planform (across slope) curvature. Although these results are informative, they do not allow one to study the relationship between various factors and the hazardous event, while controlling for possible confounding effects. This may be studied by multi-way cross classification tables and logistic regression. 3.2.2. Multi-way cross classification tables With 365 subjects, an attempt was made to keep the total number of cells to a minimum. This was needed since cells with zero frequency (empty cells) in a multiway table would create numerical problems in calculating the odds ratios. The estimation of the effect of various landscape attributes on landsliding was the major objective of the multi-way analysis. The main risk factor was considered to be slope angle since the study area can be characterized by rapid, shallow soil mass movements that are primarily controlled by slope angle. The continuous variable SLOPE was treated as a dichotomous variable (SL35) in order to carry out the multi-way table analysis. SL35 has the value of one when slope angle is less than 35 degrees and zero otherwise. The response variable SLIDE takes the value of -90- one in the presence of landsliding in the grid cell and zero otherwise. The above configuration of the dependent variable SLIDE by the dichotomized risk factor SL35 produced a table with four cells. The introduction of additional variables would multiply these original four cells by the product of the number of categories in each of the variables. In choosing between the number of additional variables to be studied and the number of categories each variable should contain, it was decided that the diversification of variables should take precedence over the preciseness of categories. For this reason, each of the variables was dichotomized. This permitted the study of as many as three variables since the dichotomous scheme yielded a factor of 8 (2^) new cells per case. The total number of cells was then (8)(4)=32, which allowed for an average of 11.4 (365/32) frequencies per cell. The three additional variables selected for the analysis were soil depth (SD ) surficial material (SOIL2) and transient snowzone (TRANS). Due to the empty cell restrictions, drainage characteristics, slope shape (profile and across slope curvature), biogeoclimatic zone and aspect could not be included in the multi-way table analysis. This obviously affects the outcome, since important variables may be missed. The soil depth variable (SD) was dichotomized with Moderate-Deep = 127 and Shallow = 238 grid cells, while the transient snow zone variable (TRANS) had 150 in the TRANS = 0 group and 215 in the TRANS = 1 group. Surficial material (SOIL2) was dichotomized in a manner that bedrock outcrop (SOILS = 1) was assigned the value of 1 and all the other surficial materials were assigned the value of 0. The group SOIL2 = 1 had 71 subjects while the group SOIL2 = 0 had 294 subjects. The variables are presented in the multi-way framework of Table 3.15. A. -91 - Table 3.15.A. Multi-way cross classification table for variables SLIDE, SL35, SD, SOIL2 and TRANS SD Shallow Deep SOIL2 slope angle > 35 degrees slope angle < 35 degrees SL35 Shallow Deep Other Bedrk Other Bedrk Other Bedrk Other Bedrk SLIDE TRANS Absent No 23 14 56 12 2 8 4 5 Yes 53 0 62 0 2 2 19 2 No 1 0 2 7 0 6 3 6 Yes 11 1 35 0 4 0 17 8 Present Interaction can be assessed by testing the homogeneity of odds ratios at different levels of the third factor. Mantel and Haenszel (1959) suggested stratification as a means of reducing confounding effects. Table 3.15.B contains a breakdown of the event-exposure (SLIDE/SL35) tables obtained by stratifying on three factors, namely soil depth, surficial material, and transient snow zone. Given this scheme, one can compute the overall odds ratios for landsliding. This is accomplished by weighting the crossproducts in each stratum by the number of individuals in that stratum. Any biases due to these factors are controlled for in the computation of the overall (Mantel-Haenszel) odds ratio. Table 3.15.B. Stratified breakdown of the data in Table 3.15.A with SL35 as the main risk factor SD SOIL2 Deep Other Deep BR* SLIDE 1 0 1 0 SH ** Other SH BR 0 0 TRANS SL35 No 1 0 2 6 8 3 4 6 5 0 1 23 0 14 2 56 7 12 1 4 2 0 2 17 19 8 2 0 11 53 1 0 35 62 0 0 Yes Bedrock outcrop ** Shallow soil -92- Using the Mantel-Haenszel procedure on the strata for SOIL2=l (bedrock outcrop), the overall odds ratio O^m-H) o v e r m four tables would include control e for biases in the factors of soil depth (SD) and transient snow zone (TRANS). The computation would be: (6)(14) m ^ (0X0) 28 3 (6X12) , (8X0) 30 10 5 (8w om (ZM (°M + 28 + 3 + 30 A 0 = = i.84 2 t yi • 10 Using this method, we can obtain the weighted odds ratios for all the factors. Moreover, using the criterion for interaction, one can check to see which factors are potential effect modifiers (interactors). These values are presented in Table 3.16, where the odds ratios are also computed on tables where landsliding (SLIDE) is first paired with the three factors (SD, SOIL2, TRANS) and then where the main risk factor (SL35) is paired with these variables. These latter odds ratios can be used to check out the degree to which the variables are confounding. Table 3.16. Odds ratios for interpreting interaction and confounding Soilclepth Surficial material SO][L2 S]D Shallow Other Bedrock Deep 2.05 Interaction 2.46 2.93 5.69 Overall odds ratio = 3.86 1.97 2.20 SLIDE 1.58 5.19 SL35 Transient snowzone 1 TRANS No Yes 5.82 1.83 2.69 1.13 Based on these results, one can see that ¥ does not differ appreciably at the two levels of surficial material (SOIL2), but it does differ for the two levels of variables soil depth and transient snowzone indicating possible interaction. Surficial material and soil depth seem to be related to both landsliding and slope angle -93 - suggesting that they are potential confounders. Transient snow zone does not seem to be a possible confounder if we consider slope angle as the main risk factor since its odds ratio is close to unity = 1.13). Its lower and upper limits are 0.67 and 1.91, respectively. Suppose that one concentrates on the relationship between slope angle (SL35) and landsliding (SLIDE) and ignores information on the other variables. A simple 2X2 table can be created showing slope angle among cases and controls (Table 3.5). The crude odds ratio (i.e., not adjusted for soil depth, transient snowzone and surficial material) of SLIDE-SL35 is given by *F = 3.86. An approximate 95 percent confidence interval (2.32, 6.42) is computed by Woolf s method from Table 3.5. This suggests that slope angle (SL35) is associated with a 3.9 times increased risk of landsliding. Table 3.17. Relation of landsliding (SLIDE) to slope angle (SL35) according to transient snowzone (TRANS), ignoring soil depth (SD) and surficial material (SOIL2) SL35 <35° > 35° ^i n; TRANS=1 SLIDE Control 25 29 47 115 2.84 216 TRANS=0 SLIDE Control 19 15 105 10 8.29 149 15*105 29*115 M-H 'o'19 47 4 5 =3.87(2.30;6.52)« 4 + 2 6 1 149 216 *95 % confidence interval was calculated using the Mantel-Haenszel method (1959). If one stratifies on transient snowzone (TRANS) as shown in Table 3.17, the adjusted odds ratio ^ M H = 3.87 does not differ from the crude odds ratio -94- = 3.86. The two odds ratios being virtually identical suggests that there is no confounding effect between SL35 and TRANS. However, looking at the individual odds ratios in the two strata (TRANS = 1, TRANS = 0) one would find that these odds ratios might differ for the two levels of variable TRANS. The odds ratio for the TRANS = 0 group, ¥ = 8.29, means that on sites outside the transient snowzone the relative risk of landsliding is about 8.3 times higher on steep slopes (slope angle >35 degrees) than on gentle slopes. The Woolf s test for the homogeneity of odds ratios (Section 2.3.1.6.) yields the chisquare value of x 2 = 3.44 with 1 degree of freedom, which is not significant at 95 % confidence level (p = 0.0636), providing no evidence of interaction between slope angle and transient snow zone. Based on Table 3.9, soil depth (SD) itself is a risk factor for landsliding. The estimated odds ratio, ¥ = 2.20 (1.30, 3.73), suggests that shallow soils are associated with a 2.2 times increase in the risk of a landslide event. Table 3.18 shows that the SL35/SLIDE odds ratio adjusted for soil depth is ^M-H = 3.68. Table 3.18. Relation of landsliding (SLIDE) to slope angle (SL35) according to soil depth (SD), ignoring transient snowzone (TRANS) and surficial material (SOIL2) SL35 <35° > 35° ¥; n; SD=1 SLIDE Control 34 30 44 130 3.35 238 SD =0 SLIDE Control 14 10 90 13 4.95 127 M-H 10*90 34*130 127 13*14 238 = 3.68(2.20;6.15)* 44*30 127 •+ 238 *95 % confidence interval was calculated using the Mantel-Haenszel method (1959). -95- The crude odds ratio W = 3.86 for SL35 and the soil depth adjusted odds ratio ^SD = 3.68 does not indicate confounding effect between SL35 and SD. Using Woolf s test for the homogeneity of the odds ratios yields the chi-square value of 0.43 (p = 0.5120) with 1 degree of freedom. This value suggests that the soil depth-specific odds ratios does not differ significantly, providing no evidence of interaction between slope angle (SL35) and soil depth (SD). Table 3.19. Relation of landsliding (SLIDE) to slope angle (SL35) according to surficial material (SOIL2), ignoring transient snowzone (TRANS) and soil depth (SD) SL35 <35° > 35° n i SOIL2=l SLIDE Control 17 20 8 26 3.82 71 SOIL2=0 SLIDE Control 24 27 194 49 3.52 294 24*194 20*26 | = 11*21 294 8*17 7 = 3 6 1 ( 2 ' U ; 6 - 1 9 ) * '95 % confidence interval was calculated using the Mantel-Haenszel method (1959). A similar analysis was carried out with surficial material (SOIL2) and slope angle yielding ¥ § o i L 2 = 3.61 (Table 3.19). This value does not provide evidence of confounding effect between slope angle and surficial material. The Woolf s test yields x 2 = 0.02 with 1 degree of freedom which is far from significant (p = 0.8875), so there is no interaction between SOIL2 and SL35 based on the data. As indicated earlier (Table 3.8), surficial material itself is a risk factor for landsliding ( ¥ = 1.97). Taking variable SOIL2 as the main risk factor of interest, -96- that there is a significant confounding effect between SOIL2 and TRANS, since the calculated Mantel-Haenszel odds ratio OFM-H = 5.87) is about three times higher than the crude odds ratio (Table 2.9). This example shows that for a set of variables each variable can be chosen as the main risk factor of interest. By forming the appropriate subgroups, any of the adjusted odds ratios using the Mantel-Haenszel method can be estimated. The choice of which variable should be used for adjustment and in what order they should be introduced are not statistical issues. Knowledge on the landsliding process and its known associations with various landscape attributes should guide the analysis. 3.2.3. Stepwise multiple logistic regression Multi-way cross classification table analysis is an excellent tool for becoming familiar with the data. It provides valuable insight into the relationship between the various terrain attributes. Unfortunately, if one needs to control or adjust for numerous variables, the number of levels in each variable or the number of variables should be reduced to avoid empty cells. The logistic regression procedure not only accomplishes what the multi-way table analysis does but it goes one step further. Not only can one test for effect modification (interaction) but one is able "to quantify" the degree to which each factor is an effect modifier or a confounder, thus gaining important insight into the landsliding process. 3.2.3.1. Model building strategy The purpose of model building is to select those variables that result in the "best" model for assessing landslide risk. To achieve this goal, three major steps were carried out: -97- 1. Variable selection 2. Building the logistic regression model 3. Assessment of the adequacy and overall fit of the model 3.2.3.1.1. Variable selection Each variable was examined carefully in univariate analysis by constructing contingency tables for variables with categorical nature. The Pearson chi-square test with k-1 degrees of freedom was used to test the significance of a categorical variable with k levels. Creating k-1 design (dummy) variables and carrying out univariate analysis provided exactly the same value for the likelihood ratio test. Tables 3.20-3.27 provide the univariate analyses for categorical and ordinal variables using the output of Epi Info v.6. In addition to the overall test, for those variables exhibiting at least a moderate level of association with the dependent variable (SLIDE), the individual odds ratios were calculated for one of the levels as a reference group. Table 3.20. Contingency table for variables SOILS and SLIDE SLIDE Surficial material Present Absent SOILS 1- Bedrock* 28 43 2-Colluvium 26 79 47 3- Till 139 3 4-Fluvial sediment** 0 X =7.04 with 3 degrees of freedom ^Colluvium *Till= 2 I p=0.0707 * Reference group ** Fluvial sediment was lumped with category 3 for the analysis -98- 0.51 Since the odds ratios (¥rj u vium 0 U a n £ * ^Till) do not differ when using bedrock outcrop surficial material type as the reference level, it was decided to use the dichotomous variable SOIL2 (Table 3.1.B) in the further analysis. Table 3.21. Contingency table for variables PROFC and SLIDE SLIDE Profile curvature PROFC Present Absent 104 1- Concave 49 2- Straight 26 79 3- Convex* 47 139 X =2.67 with 2 degrees of freedom p=0.2626 ^Concave 1-39 = ^Straight 0.85 = 2 Reference group Table 3.22. Contingency table for variables PLANC and SLIDE SLIDE Planform curvature Present Absent PLANC 36 118 1- Concave 18 59 2- Straight 3- Convex* 47 87 X =5.80 with 2 degrees of freedom p=0.0551 ^Concave 0.56 - ^Straight 0.56 = 2 * Reference group Table 3.23. Contingency table for variables DISTUR and SLIDE SLIDE Disturbance history Present Absent DISTUR 1- Harvesting* 26 5 37 2- Landslides 15 50 3- Fire 20 142 53 4- Flood 9 5- Gap replacement 8 X =5.34 with 4 degrees of freedom 2 * Reference group -99- ^Landslides 2.11 ¥ i r e = 2.08 = F ^Flood= 1-94 ^ G a = 4.62 P p=0.2538 Table 3.24. Contingency table for variables ELEVZ and SLIDE Elevation zones SLIDE ELEVZ Present Absent 4 1- 0-450 m* 36 140 2- 450-900 m 76 21 3- 900m88 X = 16.18 with 2 degrees of freedom p=0.0003 4.89 ¥ 3 = 2.75 2 * Reference group Table 3.25. Contingency table for variables SLCL and SLIDE SLIDE slope classes SLCL Present Absent 17 1- 10-25 degrees* 109 2- 25-40 degrees 61 133 21 18 3- 40-55 degrees 2 4 4- 55 degrees X =27.48 with 3 degrees of freedom p=0.0000 ¥ = 2.94 2 ¥ 3 = 7.48 ¥ 4 = 3.21 2 Reference group Table 3.26. Contingency table for variables ASPECT4 and SLIDE SLIDE Aspect ASPECT4 Absent Present 67 1- North* 26 114 41 2- East 34 79 3- South 4 0 4- West** X =1.98 with 3 degrees of freedom 2 I p=0.5767 Reference group Ignored due to zero cell - 100- 0.93 ¥ 3 = 1.11 Table 3.27. Contingency table for variables SDEPTH and SLIDE Soil depth SLIDE SDEPTH Present Absent 1- Shallow 78 160 7 24 2- Medium 3- Deep* 16 80 X =9.31 with 2 degrees of freedom p=0.0095 *! = 2.44 1.46 2 Reference group Attention was paid to empty cells in any contingency tables during the univariate analysis. This could cause problems in the calculation of odds ratios yielding zero or an undefined value. A solution to avoiding empty cells is to collapse categories of the independent variable, or eliminate the category completely. For continuous variables, the easiest way of analysis involves fitting an univariate logistic regression model to obtain the estimated coefficient, standard error and the likelihood ratio test for the significance of the coefficient. It is also a good practice to provide a scatterplot of the continuous variable as was shown earlier for the variable SLOPE. This plot when done on the logit scale can provide valuable information on the scale of the variable. Such a scatterplot can help to decide whether the variable is linear in the logit [g(x)] or whether higher order terms should be included. Assuming linearity in the logit in the variable selection stage is a common practice and is consistent with the goal of ascertaining if a variable should be in the model (Hosmer and Lemeshow 1989.). The variables identified in Table 3.1 have been shown to be associated with landsliding in the literature. The goal of the current study was to ascertain if these variables are important in the Jamieson-Orchid-Elbow Creeks study area. - 101 - The result of fitting univariate logistic regression models to the data is given in Table 3.28. In this table, for purposes of preliminary analysis only, all of the variables shown in Tables 3.1.A and 3.1.B were included. The analyses were done using Logistic V.3.11E freeware (Dallal 1988). The continuous variable SLOPE that was shown to be quadratic in the logit, was analyzed assuming linearity in the logit for the sake of simplicity. In Table 3.28 the following information is presented: (1) The estimated coefficients for the univariate logistic regression model (2) The estimated standard error for the estimated coefficients (3) The estimated odds ratio which was obtained by exponentiating the estimated coefficients. Please note that for variables ELEV and SLOPE the odds ratios were calculated for 20 metres and 10 degrees, respectively. This was done since the change of one meter and 1 degree would not be as meaningful. (4) The 95 % confidence interval for the odds ratios (p ±Z,_ t a/2 *SE(J? )) i (5) The value of the log-likelihood (LL) for the model and (6) the likelihood ratio test statistic G, for the hypothesis that the coefficient is zero (7) The degrees of freedom related to the test statistic G (8) The p-value indicating the probability that the estimated coefficient is zero 3.2.3.1.2. Multiple logistic regression After completing the univariate analyses, variables need to be selected for the multivariate analysis. Based on Hosmer and Lemeshow (1989), any variable whose univariate likelihood ratio test yields a p value less than 0.25 should be considered a candidate for the multivariate model. The problem with the univariate approach is - 102- Table 3.28. Univariate logistic regression models for variables considered in this study Variable Constant ELEV SLOPE SLCL* SLl SL2 SL35 SDPTH* SDl SD2 SD DR01 501 502 SOIL2 BGCZ Dl D2 D3 D4 FOREST ASP2 CA30 ELEVZ* TRANS PR1 PR2 PL1 PL2 3 -0.9608 0.0155 0.0700 0.8352 1.2129 2.1042 1.3506 -0.4502 -0.5137 -0.8910 0.7904 1.2324 -0.6824 -0.6767 -0.6787 -0.4474 0.7458 0.7324 0.6631 1.5309 0.0429 0.3890 0.9592 -0.0769 0.9905 -0.4891 -0.3288 0.6748 0.5714 SE(P) 0.1170 0.0447 0.0126 0.1792 0.2561 0.5206 0.2599 0.1506 0.4512 0.3067 0.2686 0.4506 0.3318 0.2954 0.2778 0.2561 0.5763 0.5554 0.5142 0.6889 0.2412 0.2970 0.3570 0.1916 0.2615 0.3594 0.2357 0.3298 0.2627 ¥ 1.36 2.01 2.31 3.36 8.20 3.86 0.64 0.60 0.41 2.20 3.43 0.51 0.51 1.97 0.64 2.11 2.08 1.94 4.62 1.04 1.48 2.61 0.93 2.69 0.61 0.72 1.96 1.77 95 % CI 0.24 1.57 1.62 2.04 2.96 2.32 0.47 0.25 0.22 1.30 1.42 0.26 0.28 1.14 0.39 0.68 0.70 0.71 1.20 0.65 0.82 1.30 0.64 1.61 0.30 0.44 1.03 1.06 7.86 2.58 3.28 5.56 22.8 6.42 0.86 1.45 0.75 3.73 8.29 0.97 0.91 3.40 1.06 6.52 6.18 5.32 17.8 1.67 2.64 5.25 1.35 4.50 1.24 1.18 3.75 2.96 LL -215.3 -215.2 -197.4 -203.5 -198.9 G df p value 0.12 1 35.69 1 23.58 1 32.80 2 0.7279 0.0000 0.0000 0.0000 -201.8 -210.4 -210.4 26.97 9.82 9.84 1 1 2 0.0000 0.0017 0.0073 -210.6 -210.5 -212.4 9.31 9.61 5.79 1 1 2 0.0023 0.0019 0.0554 -212.4 -213.7 -212.7 5.79 3.15 5.27 1 1 4 0.0162 0.0759 0.2609 -215.3 -214.4 -211.8 -215.2 -207.5 -213.9 0.03 1.79 6.96 0.16 15.58 2.68 1 1 1 1 1 2 0.8586 0.1810 0.0083 0.6883 0.0001 0.2626 -212.5 5.70 2 0.0578 * SLCL, SDPTH and ELEVZ were treated as continuous variables since they all have an ordinal scale. Please refer to Table 3.1.A and B. - 103 - that it ignores the possibility that a collection of variables that may be weakly associated with the outcome, can become an important predictor when taken together. Due to the complexity of landsliding, it was decided that the stepwise logistic regression would be used to select variables for the final model. The technique used in the stepwise logistic regression is forward variable selection with a test for backward elimination. The stepwise approach is useful and intuitively appealing in that it builds the models in a sequential fashion based solely on statistical criteria. The logistic regression model was fitted to the case-control study using the stepwise variable selection procedure. Those variables that showed at least a moderate significance (p<0.25) were selected for the stepwise procedure. Based on this decision, the following variables were considered as possible predictors: 1. slope angle 2. soil depth 3. drainage characteristics 4. surficial material 5. biogeoclimatic zone 6. aspect 7. local catchment area 8. elevation (transient snowzone) 9. slope shape (horizontal and across slope curvature) The stepwise procedure for selection or deletion of variables from a model is based on a statistical algorithm which checks for the significance of the variables, and either includes or excludes them on the basis of a fixed decision rule (Hosmer and Lemeshow 1989). This decision rule is based on the p value of the likelihood ratio chi-square test of the models with and without the variable in question. The forward variable selection continues until all the p values of the remaining variables - 104- exceed a predefined p value (PE). In this study the value of P E was set to 0.05. A backward checking for elimination was carried out at each step of the stepwise procedure. If a variable became insignificant after including another variable into the model, it was taken out of the model. The significance of the variables were assessed by the p value of the likelihood ratio chi-square test. This p value is usually called pR and was set to 0.10 in the present study. Two models were developed on the basis of two different forms of the variable SLOPE. In the development of model "A", SLOPE was used as a continuous variable. Model "B" was developed using the variable SLOCL which is a categorical variable that was created from variable SLCL using the slope class mid-points. Tables 3.29.A and 3.29.B present the stepwise variable selection for the two models on the basis of the p values calculated from the likelihood ratio test described previously. The order of the variables given columnwise in the tables is the order in which they were entered in the models. In each row the p values to the left of the vertical line are used in the backward checking and the values to the right of the vertical line are used to enter the next variable into the model (forward selection). Table 3.29.A Results of applying stepwise variable selection using the maximum likelihood method to the study data in the development of model A. # SLOPE TRANS SOIL2 SD CA30 SLOPE2 PR1 PL1 PR2 PL2 BGCZ ASP2 DR01 0 0.000 0.000 0.016 0.002 0.008 0.263 0.058 0.076 0.181 0.002 1 0.000 0.001 0.502 0.021 0.006 0.144 0.081 0.009 0.292 0.010 2 0.000 0.001* 0.005 0.020 0.007 0.184 0.030 0.584 0.435 0.036 3 0.000 0.000 0.005* 0.007 0.016 0.248 0.069 0.197 0.647 0.177 4 0.000 0.000 0.002 0.007* 0.023 0.173 0.078 0.010 0.475 0.071 5 0.000 0.000 0.004 0.010 0.023* 0.273 0.266 0.092 0.354 0.090 * the maximum p value to remove at each step. - 105 - Table 3.29.B Results of applying stepwise variable selection using the maximum likelihood method to the study data in the development of model B. ft SLOCL TRANS SOIL2 CA30 SD SLOCL2 PR1 PL1 PR2 PL2 BGCZ ASP2 DR01 0 0.000 0.000 0.016 0.002 0.008 0.263 0.058 0.076 0.181 0.002 1 0.000 0.000 0.160 0.008 0.009 0.195 0.091 0.021 0.230 0.003 2 0.000 0.000* 0.000 0.007 0.010 0.247 0.081 0.999 0.359 0.023 3 0.001* 0.000 0.000 0.003 0.021 0.312 0.086 0.279 0.603 0.199 4 0.004* 0.000 0.000 0.003 0.030 0.224 0.102 0.124 0.377 0.066 5 0.003 0.000 0.000 0.004 0.030* 0.338 0.317 0.086 0.306 0.087 | * the maximum p value to remove at each step. Table 3.30.A and 3.30.B present the steps needed to develop the 2 models. Logistic 3.11Ef of STATOOLS™ was used in the analysis. The columns of the tables respectively contain the label, the effect being tested, the log-likelihood (LL), the likelihood ratio test statistic (G), the degrees of freedom of the model and the pvalue. Table 3.30.A Log-likelihood for Model A at each step and likelihood ratio test statistics (G), degrees of freedom (d.f.) and the p-values corresponding to the addition of the variable Step# 0 1 2 3 4 5 Effect tested Constant + SLOPE + SLOPE + TRANS + SOIL2 +SD +CA30 2 LL -215.2848 -195.3887 -189.7926 -185.7605 -182.0933 -179.4884 G 39.7921 11.1923 8.0642 7.3344 5.2098 d.f. 2 3 4 5 6 p-value 1 0.0000 0.0008 0.0045 0.0068 0.0225 Final model: g(x)= -5.337+0.103*SZ,OP£ -0.001 *SLOPE + 1AS3*TRANS + l.229*SOIL2+0J63*SD +0.923*G430 7 Log-likelihood = - 179.4884 - 106- G=71.5928 with six degrees of freedom (p=0.0000) Table 3.30.B Log-likelihood for Model B at each step and likelihood ratio test statistics (G), degrees of freedom (d.f.) and the p-values corresponding to the addition of the variable Step# 0 1 2 3 4 5 Effect tested Constant + SLOCL* + SLOCL + TRANS -1- SOIL2 +SD +CA30 LL 2 -215.2848 -202.0228 -195.6395 -189.0507 -184.4677 -182.1237 G 26.5241 12.7665 13.1776 9.1660 4.6880 d.f. 2 3 4 5 6 p-value 0.0000 0.0004 0.0003 0.0025 0.0304 * SLOCL was created using the categorical variable SLCL (Table 3.1.B). The slope class midpoints were used to create SLOCL as a categorical variable with ordinal scale. This variable may easily be used in terrain stability mapping due to its categorical nature while using SLOPE as a continuous variable may be inefficient for such purpose. Final model: g(x}= -4.974 + 0.W3*SLOCL -0.00l*SWCL 2 + l.629*TRANS +1.499*5(9/12 +0.843*SD+0.865*C430 Log-likelihood: - 182.1237 G =66.3221 with six degrees of freedom (p=0.0000) The same set of variables were identified for models A and B as important predictors namely, slope angle, transient snowzone, soil type, soil depth and catchment area. The scale of the continuous variable SLOPE was identified as quadratic in the logit, so the quadratic term was included in the models. Since both models identified the same variables and the coefficients are similar, I decided to use model B for further analysis. This decision was made since the variable SLOCL may easily be used for mapping due to its categorical nature, while using SLOPE as a continuous variable is inefficient for such purpose. Once the main effect model has been decided, we may consider stepwise selection to identify possible interactions. The candidate interaction terms are those that seem logical and geologically reasonable given the main effect variables in the - 107- model. Since landsliding is a very complex process and there are only five variables in the main effect model, it was decided to create and test all of the ten possible interaction terms. Table 3.31 shows the stepwise variable selection process for the possible interaction terms. Table 3.31. Results of applying stepwise variable selection to interaction terms from the main effects model (B), using the maximum likelihood method Interaction term SLOCLxTRANS SLOCLxSOIL2 SLOCLxSD SLOCLxCA30 TRANSxSOIL2 TRANSxSD TRANSxCA30 SOIL2xSD SOIL2xCA30 SDxCA30 p-value 0.1897 0.5220 0.7290 0.6468 0.2943 0.6315 0.3994 0.1726 0.1398 0.5777 None of the 10 interactions specified in Table 3.31 were chosen since none of the calculated p values were higher than the value of p£ = 0.05. This means that the data do not provide sufficient evidence for effect modification among the variables that are included in the main effect model. It is very important to note at this point that the fundamental reason for developing a model was to provide as clear description as possible of the associations between the outcome variable and the independent variables based on the available data. Entering any of the ten possible interaction terms into the model does not seem to improve the estimates of the relevant associations, so none of the interaction terms were used in the final model. - 108- 3.2.3.1.3. Assessing the fit of the model In the assessment of the goodness-of-fit of the logistic regression model one would like to know how effective the developed model is in describing the response variable. Generally speaking a model is considered to be well fitted if 1. the summary measures of the distance between the observed and expected values are small and 2. the contribution of each pair of the observed and expected values to the summary measures is unsystematic and is small relative to the error structure of the model (Hosmer and Lemeshow 1989) In this study, the Hosmer-Lemeshow goodness-of-fit test was used to assess the fit of the model (Hosmer and Lemeshow 1989). This test involves the grouping of subjects based on the estimated probabilities, usually into g = 10 groups. These groups represent the deciles of probabilities from zero to one. The HosmerLemeshow goodness-of-fit statistic, xC> * obtained by calculating the Pearson chi2 s square statistic from a 2xg table of observed and expected frequencies. Most statistical packages that can be used for developing a logistic regression model provide the x C 2 t e s t statistic that is approximated by a chi-square distribution with (g-2) degrees of freedom under the hypothesis (Hg) that there is a significant lack of fit of the model at 95 percent confidence level. The results of applying the deciles of probabilities grouping strategy to the estimated probabilities computed from the model developed for the study are given in Table 3.32. The value of the Hosmer-Lemeshow goodness-of-fit statistic computed from the frequencies in Table 3.32 is xC> 2 = 7.12 and the corresponding p value computed from the chi-square distribution with 8 degrees of freedom was 0.52. This - 109- indicates that we can accept the hypothesis that there is no evidence for lack of fit of the model at 95 percent confidence level. Table 3.32. Observed and expected frequencies within each decile of probability for each outcome (SLIDE = 0, SLIDE = 1) using the fitted logistic regression model Deciles of probabilities 0.0-0.1 0.1-0.2 0.2-0.3 0.3-0.4 0.4-0.5 0.5-0.6 0.6-0.7 0.7-0.8 0.8-0.9 0.9-1.0 Total SLIDE = 1 SLIDE=0 Observed 62 58 58 68 4 4 6 3 1 0 264 Expected 60.8 57.6 62.2 64.6 5.6 6.1 4.1 2.1 0.8 0.1 264 Observed 1 8 23 37 6 11 5 5 4 1 101 Expected 2.2 8.4 18.8 40.4 4.4 8.9 6.9 5.9 4.2 0.9 101 Total 63 66 81 105 10 15 11 8 5 1 365 3.3. Landslide risk assessment for the study area 3.3.1. The logistic regression model The logistic regression model developed in the previous section indicated that the most important factors in landslide development in the Jamieson-Orchid-Elbow Creeks study area were: 1. slope angle 2. transient snowzone (elevation) - 110- 3. surficial material (bedrock outcrops) 4. soil depth 5. local catchment area The variables selected for the model are logically and geologically meaningful as they are known to be associated with shallow, rapid soil mass movements. The variables can easily be obtained from aerial photographs, topographical maps, soil maps and forest cover maps. Table 3.33 presents the coefficients, standard errors and odds ratios for the final model. The coefficient of a variable measures potential confounding since it represents the degree of association between landsliding and the variable in question, once the other factors are controlled for. Based on the information in Table 3.33, we can state, for example, that there is a 5.1 times higher risk of landsliding on sites located in the transient snowzone compared with sites outside this zone when the other important variables are controlled for. One can see that all the variables in the model indicate a positive association with landsliding since the lower odds ratios are still higher than one. Table 3.33. Coefficients, standard errors (SE) and odds ratios of the final model Variable Constant SLOCL SLOCL2 TRANS SOIL2 SD CA30 Coefficient - 4.9738 0.0834 - 0.0006 1.6286 1.4987 0.8425 0.8649 SE 0.9749 0.0558 0.0008 0.3729 0.4230 0.2985 0.3958 0.77 2.17 6.83 2.45 1.95 1.29 1.09 5.10 4.48 2.32 2.37 10.59 10.26 4.16 5.16 Since one degree increment in slope angle does not provide sufficient information about the increasedriskof landsliding, a ten-degree increment was used - Ill - in calculating the odds ratio ¥ = 2.17 (0.77, 6.83). This indicates that for every increase of 10 degrees in slope angle, the risk of landsliding increases 2.17 times. The validity of such a statement is questionable without considering the appropriate scale of the continuous variable SLOPE, since the additional risk of landsliding for a site with an average slope angle of 30 degrees compared to a site with slope angle of 20 degrees may be quite different from the additional risk of landsliding for sites with 50 degrees compared to 40 degrees. Since it was shown that slope angle is not linear in the logit, the higher order term (quadratic) was included in the model. The confidence interval is shifted to the right providing evidence that there is a positive association between slope angle and landsliding. The next step in the terrain stability assessment procedure is to estimate the relative risk (odds ratios) of landsliding for the 20x20 m grid cells in the study area. 3.3.2. Landslide risk assessment A way of illustrating the effect of the variables in the model is to estimate the "relative" risk of landsliding as compared to a base level. In this way, one could compare all variable combinations to a common base and in the process compare the relative risk of landsliding for different variable combinations. 3.3.2.1. Choosing the base level of variables For the study area the logistic regression model predicted the expected frequencies to be minimum for subjects (20x20 m grid cells) that 1. have gentle slopes (SLOCL = 17.5°) 2. are outside the transient snowzone (TRANS = 0) 3. have surficial material type other than bedrock outcrop (SOIL2 = 0) - 112- 4. have deep soils (SD = 0) and 5. have local catchment area greater than 1.2 ha (30 grid cells) (CA30 = 0). This selection also ensures that the calculated relative risk estimates for the other variable combinations will result in values greater than one. 3.3.2.2. The landslide risk matrix Given the base level of variables, one can compute the "relative" risk in terms of odds ratios for any variable combinations using Equation 2.33: V(x:x)= x &/3 (x;-x )] e V i l The odds ratio depends only on those factors for which the two subjects differ. Suppose that one wished to estimate the relative risk for grid cells with slope category 4 (SLOCL = 65.0°), located in the transient snowzone (TRANS = 1), with bedrock outcrop surficial material type (SOIL2 = 1) and shallow soils (SD=1). The relative risk would then be the odds ratio calculated for the above variable combination as compared with the base level. Using Equation 2.33: F=exp[(0.0834)(65.0-17.5)-(0.0006)(65.0 -17.5 )+(1.6286)(l-0)+(1.4987)(l- x 2 2 0)+ +(0.8425)(l-0)+(0.8649)(l-0)]= exp(6.44495) = 629.5 This means that there is 630 times increased risk of landsliding on sites with the maximum levels of the variables relative to the base level. One should note that because the variables are "standardized" to the base level, the relative risk for the base variable combination is unity. Table 3.34 shows the relative risks of - 113- landsliding for the various combinations of factors. This type of representation of the relative risk will be called the Landslide Risk Matrix. Table 3.34. Landslide risk matrix for the Jamieson-Orchid-Elbow Creeks Study Area CA30 0 - catchment area > 1.2 ha SD Deep soils SOIL2 Ot ler Shallow soils Bedrock OC* Ot ler Bedrock OC TRANS 0 1 0 1 0 1 0 1 10 -25 1.0 2.2 3.8 5.5 5.1 11.4 19.3 25.5 4.5 10.0 17.0 22.4 22.8 50.8 86.4 114.2 2.3 5.2 8.8 11.6 11.8 26.4 44.8 59.2 10.4 23.2 39.4 52.0 53.0 118.0 200.7 265.1 c 25°-40° 40--55 • 55°- 1 - catchment area < 1.2 ha CA30 SD Shallow soils Deep soils SOIL2 Ot ler Bedrock OC Ot ler Bedrock OC TRANS 0 1 0 1 0 1 0 1 I0 -25 25°-40° 40*'-55 55°- 2.4 5.3 9.0 11.9 12.1 27.0 45.0 60.6 10.6 23.7 40.3 53.2 54.2 120.7 205.2 271.1 5.5 12.3 20.0 27.6 28.1 62.6 106.5 140.6 24.7 55.0 93.5 123.5 125.8 280.3 476.6 629.5 4 C c Bedrock outcrop If one were to rank order the "relative" risks in the Risk Matrix, one would find that, holding all the other variables constant, slope angle has the greatest effect on the magnitude of the relative risk. However, transient snowzone, surficial material, soil depth and catchment area are so related to landsliding that, in most cases, they override the effect of slope angle. Hence, sites with gentle slopes (10°-25°) in the transient snowzone that have shallow soils with bedrock outcrops and with a catchment area of 1.2 ha or less, have greater relative risk of landsliding than sites with steep slopes located outside the transient snowzone with relatively deep - 114- soils. Given the results in Table 3.34, one can see that study subjects with very steep slopes (over 55°) have, on the average, five times the chance of experiencing a landslide event relative to sites with gentle slope (10°- 25°) Moreover, this is constant over the sixteen variable groups because of the lack of effect modification (interaction). With interaction, the ratio of relative risk would vary over the variable groups. 3.3.2.3. Landslide risk rating After the analysis is performed to obtain the Landslide Risk Matrix for a particular area, it can be used as a planning tool. The overall risk value for a particular site may be increased by calculating the effects of additional land use risk values for planned activities. There have been some attempts to carry out overall risk assessment of a given watershed using GIS to evaluate the watershed sensitivity to land use changes (Lull et al. 1995). The landslide risk matrix can provide valuable information in the overall risk assessment of a particular watershed. Since such overall risk assessment for watersheds is still under development, this study can only suggest the use of the Landslide Risk Matrix in watershed risk analysis in the future. Given the information in the landslide risk matrix, a classification scheme for terrain stability mapping can be established. For the sake of comparison with the terrain stability classes currently used in British Columbia, 5 risk categories were set up in the following manner: - 115 - Table 3.35. Classification scheme for the study Risk category 1 2 3 4 5 Description Very low risk Low risk Moderate risk High risk Very high risk Relative Risk Class boundaries 1.0-10.0 10.1-20.0 20.1-30.0 30.1-60.0 60.1- The 5 categories represent increasing level of "relative" risk. Risk category 4, for example, represents an approximate 6 times greater risk than category 1 (60.0/10.0) based on the upper class boundaries. The class boundaries should be set on the basis of local knowledge, experience and individual judgement of the terrain mapper. Using IDRISI v.4.1. GIS software, a landslide risk map was created for the Jamieson-Orchid-Elbow Creeks study area. This was done by creating and manipulating the base maps for the individual factors. Subsequently, these maps were combined to produce a map with the odds ratios for the individual pixels. Finally, this map was classified using the classification scheme given in Table 3.35. In addition to this classification, all areas within a 30-m buffer zone around the stream channels were considered as high risk areas. This decision was guided by the principle that the protection and enhancement of a quality water supply are priorities in the management of the Greater Vancouver Water District community watersheds. According to Thurber Engineering Ltd. (1991.) the most obvious source of turbidity-generating sediment is the erosion of stream banks that carry silt-clay and organic fines into the reservoirs. Assigning the highest risk category to the 30-m buffer zone around the stream channels helps to avoid the acceleration of channel bank erosion due to forest management activities. Figure 3.3. shows the final landslide risk map that was created for the study area. Generalization of the map -116- was achieved by IDRISI using its filter module. This module generalizes the risk map based on mathematical algorithms. The "median" option was used for filtering the image, which calculates the median for a grid cell from the surrounding 8 cells. Cross tabulation of the original image and the filtered image showed that filtering provided a "smoother" image without significantly modifying the original map. This helped to create larger, generalized areas that could be used in forest management planning. A terrain stability assessment according to the British Columbia Terrain Classification System was carried out on the study area in 1992 as a part of the Ecological Inventory Pilot Study (EIPS) (Acres International Ltd. 1993). This procedure involves the creation of terrain units (usually referred as terrain polygons) on the basis of surficial material, texture, slope gradient (surface expression) and geological processes. After delineating the terrain polygons, the mapper subjectively assigns a qualitative hazard rating to these polygons. This procedure will be referred to as "terrain polygon method" or TPM later in this section. The terrain stability map based on the TPM can be seen in Figure 3.4. Table 3.36 presents a cross classification table of the map created in the EIPS (Acres International Ltd. 1993) compared with the map created in this study. Table 3.36. Cross classification table Classes s T U D Y 1 2 3 4 5 Total 1 1467* 0 0 0 454 1921 Terrain polygon method 4 2 3 2527 162 1959 652 3096 86 219 2276 1028 2 312 18 394 958 2975 1427 11186 4051 * The values represent the number of 20x20 m grid cells - 117- 5 Total 250 1038 1550 978 7463 11279 6365 4872 5073 1310 12244 29864 -118- -119- At first glance, there seem to be a great difference between the two maps. Only 39.7 percent of the 29,864 grid cells have been identified identically. The two maps differ for a number of reasons. First of all, though both maps are based on 5 classes, they are not based on the same classification scheme. Five risk categories were established in this study, but these risk categories may not coincide well with the currently used terrain stability classification system. There is a great deal of subjectivity involved in the assignment of a hazard rating to the individual terrain polygons in the TPM. This can be problematic since the quality of the hazard map produced by the TPM depends greatly on the mapper's experience. This subjectivity is partly eliminated by the quantitative risk assessment used in this study, although setting the risk categories does entail an arbitrariness as to the specific cutpoints used to define the categories. The 454 grid cells classified as risk category 5 by this study were identified as class 1 by the TPM (Table 3.36). This difference is because in this study a 30-m buffer zone was created around the stream system and was classified as risk category 5 to avoid possible channel bank erosion. Areas around the stream channels on the lower part of the drainage were rated as class 1 or 2 by the TPM (Figure 3.4). Of 11,279 grid cells classified as terrain stability class 5 in the TPM, about 38 percent (3816 cells) were classified other than category 5 by this study, because the 20-m resolution used by the current study provided a "finer" areal classification. This is mainly because the most important variable, slope angle, was calculated for each individual grid cell in this study. Some of the terrain polygons in the TPM are 15-20 hectares in size. Since slope angle is clearly one of the most important variables in the generation of shallow, rapid soil mass movements, the "extra" information on the average slope angle for the 20x20 m cells provided a "finer" classification. The incorporation of this information was found to be the major cause of differences in the two maps. - 120- There seem to be no agreement in the mid-classes of the two classification systems. This is expected, since subjectivity cannot be ruled out from the TPM, and it is most apparent in class 2, 3 and 4 where class boundaries are hard to define. Both maps identified sites with active or inactive landsliding as class 4 or 5, though the average class for "landslide sites" was a little bit lower in this study (4.4) than in the TPM (4.8). This is mainly because the TPM tends to assign class 4 or 5 to polygons with ongoing or past landslide activities. Table 3.37 provides information on the areal risk (weighted average) produced by the two methods. The areal risk is similar for the two methods (3.3 and 3.5) (Table 3.37). This can be explained by the fact that the fractional areas of the risk categories are very similar with the exception of category 1 and 3. It appears that the study procedure identified a greater proportion (22 %) of the watershed as risk category 1 and a smaller proportion (17 %) as risk category 3, whereas the terrain polygon method identified a very small portion of the watershed as category 1 (7 %) and a much higher percentage (37 %) as category 3. Table 3.37. Calculation of areal risk based on the two methods Study procedure Risk category Percentage FWS* x Risk 0.2174 21.74 1 16.22 2 0.3244 16.90 0.5070 3 0.1744 4 4.36 40.78 2.039 5 3.2622 100.00 Total Terrain polygon method FWS x Risk Percentage 0.0695 6.95 0.2700 13.50 1.1175 37.25 0.1900 4.75 1.8875 37.55 3.5245 100.00 fractional area of watershed The major advantage of using the landslide risk matrix for terrain stability classification purposes is that it gives quantified information about the relative risk - 121 - of landsliding while letting the mapper decide on the class boundaries based on local knowledge, experience and the planned forest management activities. Altogether, we can state that the maps produced by two methods are different because of: 1. subjective elements involved in the classification schemes 2. incorporation of "extra" information on slope angle in this study 3. arbitrariness of the cutpoints of the risk categories - 122- 4. SUMMARY AND CONCLUSIONS This study was designed to investigate the usefulness of logistic regression models and the case-control method in terrain stability assessment. A case-control study was carried out on the Jamieson-Orchid-Elbow Creeks study area and it was analyzed by multi-way cross classification tables and multiple logistic regression. A landslide risk model was created that was used to produce a landslide risk map for the study area. The statistical procedures were shown to be a useful tool in determining the cause-effect relationships between landsliding and various landscape attributes. It was found that slope angle, transient snowzone, surficial material type, soil depth and local catchment area play an important role in the development of rapid, shallow soil mass movements in the Jamieson-Orchid-Elbow Creeks study area. These attributes are logically and geologically meaningful and they are known to be associated with rapid, shallow soil mass movements in the literature. Based on the landslide risk model developed in this study, slope angle has the greatest effect on the magnitude of landslide risk, holding all the other variables constant. Slope angle influences flow rates of water and sediment by controlling the rate of energy expenditure or stream power available to drive the flow. Sites in the transient snow zone have approximately five times increased risk of landsliding compared to sites outside the transient snow zone in the study area. This can be explained by the hydrologic characteristics of this zone, where rain-on-snow events are common in - 123 - the winter. The most frequent peak flow events are generally produced within this elevation band (Chatwin 1994). Soil depth provides a quantitative indicator of the depth and type of the slope failure that may occur at a site. Debris flow-debris avalanche landslides are generated on soils where impermeable bedrock or glacial till is close to the surface. The texture of surficial materials largely determines the physical properties (permeability, drainage conditions, erodibility etc.) of the soils. It was found that sites with bedrock outcrop surficial material have a 4.5 times increased risk of landsliding as compared to sites with other surficial material types. This can be attributed to soil depth and also to the physical properties of the soil determined by the soil texture. Local catchment area was found to be an important factor in recognizing landslide-prone terrain. Sites with a local catchment area of 1.2 ha or less have a 2.5 times increased risk of landsliding as compared to sites with catchment area greater than 1.2 ha based on this study. In the study area, channel heads seem to coincide with small-scale debris flow scars in topographic hollows. Local catchment area is a good indicator of the location of channel heads (Montgomery and Dietrich 1994), thus it is also associated with the occurrence of shallow landslides in the Jamieson-Orchid-Elbow Creeks study area. During the preliminary analysis, a possible interaction between drainage conditions and slope angle was detected. Since drainage condition was not found significant at a 95 percent confidence level in the logistic regression analysis, the drainage condition-slope angle interaction term was not included in the final landslide risk model. There was no evidence of interaction between the model variables based on the stepwise multiple logistic regression analysis. The minimum risk of landsliding was found to be on sites with gentle slopes, that are located outside the transient snow zone, with deep soils and with a local catchment area of 1.2 ha or greater. - 124- A simple, raster-based GIS was shown to be capable of managing and manipulating the spatial, quantitative information and creating color maps and images. Geographical Information System, in combination with terrain stability assessment and overall risk assessment, is a valuable tool in planning forest management activities. The study provided a method for quantitative landslide risk assessment that has long been awaited by landslide hazard specialists. There are two major reasons why risk assessment procedures, rather than subjective evaluations, will be used in terrain stability assessment in the near future: 1. geotechnical consulting firms with legal liability will discourage qualitative landslide hazard assessments due to their inherent subjectivity 2. due to the current legal framework, terrain stability assessment will require quantified definitions of hazard and risk that allow approval authorities to evaluate risk logically (Gerath 1995) The results of the landslide risk assessment for the study area are encouraging, but additional research is required to refine the mechanics of the model. There are two basic types of uncertainties in the risk assessment and decision making procedure that should be addressed. There is an inherent uncertainty in the underlying database on which the analysis is based and there is uncertainty in the decision rule. Basically all the topographical indices (slope, aspect, curvatures) were collected from the Terrain Resource Inventory Mapping (TRIM) 1:20,000 digital topographical database. The TRIM digital data sets conform to the British Columbia mapping specifications and guidelines (B.C. Ministry of Environment and Parks 1988). The assessment of the accuracy of the TRIM map was beyond the scope of this research but it would be desirable to assess the accuracy of such maps in the field. Establishing cutpoints for the categories based on the Landslide Risk Matrix involves decision rule uncertainty. There is a need for studies in the field of GIS - 125 - that incorporate the uncertainties of the database and the decision rule applying recently developed decision and probability theories. In the assessment of the risk of shallow, rapid soil mass movements, slope angle and pore water pressure are known as key variables. Information on pore water pressure is difficult to collect for large areas. Since channel initiation by landsliding is controlled by critical pore water pressure in colluvial soils, it can be approximated by local catchment area, drainage conditions of the soil and slope shape (planform curvature) derived from a Digital Terrain Model (DTM) (Montgomery and Foufoula-Georgiou 1993). Local catchment area was computed from a grid based DTM by using a simple flow-routing algorithm, directing flow from a grid cell to one of its eight neigbours. This method does not allow for the representation of divergent flow that may occur in some areas of the drainage basin. A number of newer algorithms for representing divergent flow are being developed (Montgomery 1994). A better assessment may be achieved by incorporating information from rainfall frequency-intensity-duration maps developed for small mountain drainage basins. Unfortunately such maps are hard to produce due to the need for long-term precipitation records. It is felt that a 40-m resolution may be sufficient to carry out the study procedure. This would reduce the target population from which the case-control sample is taken. Another area that may be worth exploration is the use of a combination of terrain polygons created by the B.C. Terrain Classification System (Howes and Kenk 1988) and a slope map to generate terrain units for landslide risk assessment in areas where shallow, rapid soil mass movements dominate. The logistic regression model could then be applied to these units to create the landslide risk matrix for the area under investigation. This study did not differentiate between the two major types of soil mass movements, namely the debrisflow-debrisavalanche landslides and the slump-126- earthflow landslides. Steep slopes with coarse, shallow soils tend to generate debris flows and debris avalanches. Slump-earthflow landslides are deep seated, rotational failures that do not depend heavily on slope angle. Since the study area is characterized by shallow soils underlain by impermeable bedrock, debris flowdebris avalanche landslides prevail. In areas where both landslide types are present, the analysis should be done by differentiating between the two types. The statistical analysis can be improved by the introduction of more "sophisticated" sampling techniques, such as stratified random sampling, multistage sampling, or frequency sampling. For areas where the number of landslides (cases) may not provide sufficient study power, sampling with replacement could be a solution. Matched case-control studies may provide a better understanding of the cause-effect relationships between landsliding and various factors, but it is still to be explored. The theories for different sampling methods in other sciences, such as modern biomedical research are under development. The monitoring of epidemiological literature would be very helpful, but difficult due to the "distance" between forestry and biomedical research. - 127- GLOSSARY OF TERMS Biogeoclimatic ecosystem classification: A system of ecosystem classification for developing ecologically based forest management strategies. Case-control study: subjects with a particular condition (i.e., a landslide event) (the cases) are selected for comparison with series of subjects without the particular condition (the controls). Cases and controls are compared with respect to certain attributes or exposures thought to be relevant to the development of the condition. Debris avalanche-debris flow landslides: Rapid, shallow soil mass movements from hillslope areas, occurring on shallow, non-cohesive soils, where subsurface water may be concentrated by subtle topography on bedrock or glacial till surfaces. Debris torrent: the rapid movement of water-charged soil, rock and organic debris in a channel. Movement is very rapid because of the high water content. These slope failures are commonly initiated by debris slides on steep gully headwalls and side walls. Digital Terrain Model (DTM): A quantitative model of the landscape in digital form. It is also called a Digital Elevation Model (DEM). Geographical Information System (GIS): An integrated system of hardware, software, and procedures designed to support the capture, management, analysis, modeling and display of spatially represented data for solving complex management and planning problems. Hazard: The probability of a potentially damaging event occurring within a given time and a given area. - 128- Isohyetal method: lines of equal precipitation, or isohyets, drawn on the map from the observer's knowledge of the basin topography, storm patterns, and from the amount of precipitation measured at each gauge. Matched case-control study: the pairing of one or more controls to each case on the basis of their "similarity" with respect to the selected variables. Pore water pressure: a measure of the pressure produced by the head of water in a saturated soil and transferred to the base of the soil through the pore water. Porewater pressure is a key factor in failure of a steep slope soil, and operates primarily by reducing the weight component of soil shear strength. Risk: represents the outcome of a hazard and is expressed as the product of the probability of a hazard, the probability of exposure to the hazard, and the probability of a particular outcome given the exposure. Slump-earthflow landslides: Slowly moving, deep seated rotational movements of a block of soil material over a broadly concave slip surface that usually occur on deep, poorly drained cohesive soils. Surface expression: Refers to the form and pattern of forms expressed by a surficial material at the land surface. The surface expression of surficial materials is classified according to slope, geometric shape and spatial pattern. Terrain stability assessment: the assessment of the susceptibility of the landscape to soil mass movement. Stability is expressed on a relative scale ranging from stable to unstable. Terrain polygons: Terrain is mapped and presented in homogeneous units based on texture, surficial materials, surface expression, geological processes and qualifying descriptors. Transient snowzone: An elevation band where the snowpack is not permanent throughout the winter, but tends to come and go depending on the meteorological conditions, and rain-on-snow events are common. The most frequent peak flow events are produced in this elevation band that spans from about 300-500 m up to 800-1000 m in coastal British Columbia. - 129- TRIM: Terrain Resource Inventory Mapping refers to B.C. Ministry of Crown Lands digital mapping of B.C. at 1:20,000 scale from aerial photography. Turbidity: The measure of the "cloudiness" of water produced by the scattering of light by suspended material. - 130- REFERENCES Acres International Ltd. 1993. GVWD Watershed Ecological Inventory Pilot Study (Jamieson-Orchid-Elbow Drainage) Final Report. Prepared by Acres International Ltd., Blackwell and Associates, Northwest Hydraulic Consultants, Oikos Ecological Consulting, Phero Tech Inc. and Remtech Inc. Vancouver, B.C. 211 pp. Bailey, R.W. 1937. A new epicycle of erosion. J. of Forestry. 35(11): 997-999. Baker, R.F. and Chieruzzi, R. 1959. Regional concept of landslide occurrence. Highway Research Board Bull. 216, pp. 1-16. B.C. Ministry of Environment and Parks. 1988. Specifications and guidelines 1:20,000 digital mapping. Release 3.0, Surveys and Resource Mapping Branch. 281 pp. Beaudry, P.G. 1984. Effect of forest harvesting on snowmelt during rainfall in Coastal British Columbia. M.F. Thesis. Faculty of Forestry, University of British Columbia, Vancouver, B.C. 183 pp. Bourgeois, W.W. 1974. Alliford Bay area geotechnic inventory. Land Use Planning Advisory Team, Woodlands Services, MacMillan Bloedel Limited, Nanaimo, British Columbia. Bourgeois, W.W. 1978. Timber harvesting activities on steep Vancouver Island terrain. In Procs. 5th Nort American Forest Soils Conference. Colorado State University. Fort Collins, Colorado, pp. 393-409. r - 131 - Brabb, E.E. , E.H. Pampeyan and M.G. Bonilla 1972. Landslide susceptibilityin San Mateo County, California. U.S. Geological Survey Miscellaneous Field StudiesMap MF-360. Briere, D. 1979. The stratification of forested landscapes for intensive management: development and application. Unpublished Ph.D. Thesis. University of British Columbia. 254 pp. Burroughs, E.R. 1984. Landslide hazard rating for portions of the Oregon Coast Range. Procs., Symposium of Effects of Forest Land Use on Erosion and Slope Stability (711 May 1984), University of Hawaii, Honolulu. Burroughs, E.R. 1985. Survey of slope stability problems of forest lands in the West. In Procs., Workshop of Slope Stability: Problems and Solutions in Forest Management General Technical Report PNW-180, US Forest Service, Portland, Oregon, pp. 516. Bush, G. 1985. Harvest planning and layout on steep terrain - The Siuslaw model. In Procs., Workshop of Slope Stability: Problems and Solutions in Forest Management General Technical Report PNW-180, US Forest Service, Portland, Oregon, pp. 6467. Carrara, A. 1983. Multivariate models for landslide hazard evaluation. Mathematical Geology, Vol. 15, pp. 403-426. Carrara, A. 1988. Drainage and divide networks derived from high fidelity digital terrain models. Quantitative Analysis of Mineral and Energy Resources, Chung, C.F.Ed., D. Riedel Publishing Co., pp. 581-597. Chang, S. 1992. The Simprecise mapping and evaluation system for engineering geological and landslide hazard zonation. Procs., 6th International Symposium on Landslides, Christchurch, New Zealand, pp. 905-911. Chatwin, S. and Rollerson. T. 1983. Landslide study: TFL 39, Block 6. Woodland Services Division. MacMillan Bloedel Ltd., Nanaimo, British Columbia 18 pp. - 132- Chatwin, S.C., Howes, D.E., Schwab, J.W. and Swanston, D.N. 1991. A guide for management of landslide-prone terrain in the Pacific Northwest. B.C. Ministry of Forests, Land Management Handbook No. 18. Victoria, B.C. 220 pp. Chatwin, S.C. 1994. The Coastal Watershed Assessment Procedure (Draft). Research Branch, Ministry of Forests. Victoria, B.C. 62 pp. Cochran, W.G. 1977. Sampling Techniques. Wiley, New York. 428 pp. Columbier, D., Hathcock, L . , Smith, C. and Fagan, R. 1994. Epi Info Executive Health Information Shell. Centers for Disease Control and Prevention, Delaware Division of Public Health. Cornfield, J., Gordon, T. and Smith, W.W. 1961. Quantal response curves for experimentally uncontrolled variables. Bulletin of the International Statistical Institute. 38:97-115. Dallal, G.E. 1988. LOGISTIC: A Logistic Regression Program for the IBM PC. The American Statistician. 42:272-273. DeGraff, J.V., H.C. Romesburg 1984. Regional landslide-susceptibility assessment for wildland management: a matrix approach. Thresholds in Geomorphology, Coates, D.R. & Vitak, J. Ed., Allen & Unwin, Boston, pp. 401-414. Dyrness, C.T. 1967. Mass soil movements in the H.J. Andrews Experimental Forest. U.S. Dep. Agric. For. Serv., Pac. Northwest For. Range Exp. Stn., Res. Paper PNW-42, 12 pp. Ellen, S.D. and Wieczorek, G.F. 1988. Landslides, floods and marine effects of the storm of January 3-5, 1982 in the San Francisco Bay Region, California. USGS Professional Paper 1434, 310 pp. Evans, S.G. and Lister, D.R. 1984. The geomorphic effects of the July, 1983 rainstorms in the southern Cordirella. Geol. Survey of Canada Paper 84-1B, pp. 223-235. - 133- Expert Committee on Soil Survey (Canada) 1987. The Canadian System of Soil Classification. Dept. of Agriculture, Research Branch. Ottawa. 164 pp. Farewell, V.T. 1979. Some results on the estimation of logistic models based on retrospective data. Biometrika. 66:27-32. Fiksdahl, A.J. 1974. A landslide survey of the Stequalekop Creek watershed. Supplement to Final Report FRI-UW-7404. Fisheries Research Institute, Univ. of Washington, Seattle, wash. 8 pp. Fisher, R.A., 1936: The use of multiple measurements in taxonomic problems. Annals of Eugenics, 7:179-188. Forest Practices Code of British Columbia 1995. Mapping and assessing terrain stability guidebook. Ministry of Forests, Victoria, B.C. 34. pp. Furbish, D.J. and Rice, R.M. 1983. Predicting landslides related to clearcut logging, Nortwestern California, USA. Mountain Research and Development, 3(3):253-259. Gerath, R.F., Hungr, O. and VanDine D.F. 1994. Landslide hazard mapping guidelines for British Columbia (Draft). Report to the Earth Sciences Task Force British Columbia Government Resources Inventory Committee, Vancouver, B.C. 56 pp. Gerath, R.F. 1995. Landslide hazards in BC's forests: The Case for Formal Risk Assessment. The BC Professional Engineer. 3:10-13. Gimbarzevsky, P. 1988. Mass wasting in the Queen Charlotte Islands; regional inventory. British Columbia Ministry of For. Land Management Report 29., 96 pp. Golding, D.L. 1988. Jamieson Creek experimental watershed in the Greater Vancouver municipal catchments. In Proceedings of the Canadian Hydrology Symposium No. 17. pp. 229-236. Hall, R.G. 1989. Precipitation intensity duration frequency curves for Jamieson Creek and Elbow Creek watersheds in the Seymour Basin. University of British Columbia B.S.F. Graduate Essay. Forest Research Management. 67 pp. - 134- Hammond, C , Hall, D., Miller, S. and Swetik, P. 1992. Level I Stability Analysis (LISA), Documentation for version 2.0. Intermountain Research Station General Technical Report INT-285, USDA Forest Service, Moscow, Idaho. 190 pp. Hetherington, E.D. 1976. Investigation of orographic rainfall in the South Coastal Mountains of British Columbia. Unpublished PhD Thesis, University of British Columbia. 130 pp. Hosmer, D.W. and Lemeshow, S. 1989. Applied logistic regression. John Wiley & Sons. New York. 307 pp. Howes, D.E. 1981. Terrain inventory and geological hazards: Northern Vancouver Island. Assessment and planning Division Bulletin 5, Terrestrial Studies Branch, British Columbia Ministry of Environment, Victoria. 105 pp. Howes, D.E. 1987. A terrain evaluation method for predicting terrain susceptible to post-logging landslide activity. MOEP Technical Report 28, British Columbia Ministry of Environment and Parks, Victoria, British Columbia, 38 pp. Howes, D.E. and Kenk, E. 1988. Terrain classification system for British Columbia (revised edition). Manual 10, British Columbia Ministry of Environment, Victoria, B.C. 90 pp. Kienholz, H. 1978. Maps of geomorphology and natural hazards of Grindelwald, Switzerland. Arctic and Alpine Research, 10(2): 169-184. Lewis, J. and Rice, R.M. 1990. Estimating erosion risk on forest lands using improved methods of disriminant analysis. Water Resources Research, 26(8): 1721-1733. Lull, J.K., Tindall, J.A. and Potts, D.F. 1995. Assessing nonpoint-source pollution risk. Journal of Forestry. 1:35-40. Luttmerding, H.A. 1981. Soil Survey of the Langley-Vancouver area; Resources and Analysis Branch, B.C. Ministry of Environment, RAB Bulletin 18. Volume 2: Soil maps, southern Sunshine Coast and southern Coast Mountains, scale 1:50,000. Volume 3.: Description of soils. - 135 - Mantel, N. and Haenszel, W. 1959. Statistical aspects of the analysis of data from retrospective studies of disease. Journal of the National Cancer Institute 22. pp. 719-748. Mark, R.K. 1992. Map of debris flow probability San Mateo County, California. Map I-1257-M, US Geological Survey. Martz, L.W. and Jong, E. 1988. CATCH: A FORTRAN program for measuring catchment area from digital elevation models. Computers and Geosciences. 14(5):627-640. Mathews, W.H. 1979. Landslides of central Vancouver Island and the 1946 earthquake. Bull., Seismological Society of America, 69:445-450. Miles, D.W.R. 1983. Impact of landslide erosion on the growth and vegetation in the Western Oregon Cascades. Oregon State University, M.S. Thesis, Corvallis, Oregon, pp. 46. Montgomery, D.R. and W.E. Dietrich 1988. Where do channels begin? Nature. 336: 232-234. Montgomery, D.R. and E. Foufoula-Georgiou 1993. Channel network source representation using Digital Elevation Models. Water Resources Research. 29(12): 3925-3934. Montgomery, D.R. 1994. Road surface drainage, channel initiation, and slope instability. Water Resources Research. 30(6): 1925-1932. Montgomery, D.R. and W.E. Dietrich 1994. Landscape dissection and drainage areaslope thresholds. In Process models and theoretical geomorphology. M.J. Kirkby (ed.) John wiley & Sons Ltd., Baffins Lane, Chichester,England. pp. 221-246. Morrison, P.H. 1975. Ecological and geomorphological consequences of mass movements in the Alder Creek watershed and implications of forest land management. Univ. Oreg., B.Sc. Thesis, Eugene, Oreg. - 136- Neilsen, T.H. and E E . Brabb. 1977. Slope-stability studies in the San Francisco Bay region, California. Reviews in Engineering Geology, vol. 3, Landslides, Geological Society of America, pp. 235-243. Neutra, R.R. and Drolette, M.E. 1978. Estimating exposure-specific disease rates from case-control studies using Bayes' Theorem. American Journal of Epidemiology. 108(3):214-222. Niemann, K.O. and Howes, D.E. 1992. Slope stability evalutions using digital terrain models. Land Management Report No. 74., British Columbia, Ministry of Forests, Victoria, B.C. 28 pp. O'Loughlin, C.L. 1972. An investigation of the stability of the steepland forest soils in the Coast Mountains, southwest British Columbia. Ph.D. thesis, Univ. of British Columbia, Vancouver, British Columbia. 147 pp. Pillsbury, N.H. 1976. A system for landslide evaluation on igneous terrain. Colorado State Univ., Ph. D. Thesis, Fort Collins, Colo., 109 pp. Prentice, R.L. and Pyke, R. 1979. Logistic disease incidence models and case-control studies. Biometrika. 66:403-411. Reilly, T. and Powell, B. 1985. Applications of geotechnical data to forest management. In Procs., Workshop of Slope Stability: Problems and Solutions in Forest Management General Technical Report PNW-180, US Forest Service, Portland, Oregon, pp. 87-93. Rice, R.M. and Foggin, G.T. 1971. Effect of high intensity storms on soil slippage on mountainous watersheds in Southern California. Water Resources Research, 7(6): 1485-1495. Rice, R.M. and Lewis, J. 1991. Estimating erosion risks associated with logging and forest roads in Nortwestern California, Water Resources Bull., 27(5):809-818. Rice, R.M. and Pillsbury, N.H. 1982. Predicting landslides in clearcut patches. Recent Developments in the Explanation and Prediction of Erosion and Sediment Yield. In: - 137- Proceedings Exeter Symposium 1982. Int. Assoc. Hydrol. Sciences. Publication No. 137. IASH-AISH. pp. 303-311. Rice, R.M., Pillsbury, N.H. and Schmidt, K.W. 1985. Risk analysis approach for using discriminant functions to manage logging related landslides on granitic terrain. For. Sci., 31(3):772-784. Rollerson, T.P. 1992. Relationships between landscape attributes and landslide frequencies after logging-Skidegate Plateau, Queen Charlotte Islands. Land Management Report No.76, Forest Science Research Branch, British Columbia Ministry of Forests, Victoria, B.C. 11 pp. Rollerson, T.P. and Sondheim, M. 1985. Predicting post-logging terrain stability - a statistical-geographical approach. In Procs. Joint Symposium of the IUFRO Mountain Logging Section and the 6th Pacific Northwest Skyline Symposium, Vancouver, British Columbia. Rood, K.M. 1984. An Aerial Photograph Inventory of the Frequency and Yield of Mass Wasting on the Queen Charlotte Islands, British Columbia. Land Management. Report No. 34. Ministry of Forest. British Columbia. 48 pp. 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, B.C. 46 pp. Ryder, J.M. and D. E. Howes. 1984. Terrain information, a user's guide to terrain maps in British Columbia. British Columbia Ministry of Environment, Victoria, B.C. 16 pp. Ryder, J.M. and MacLean, B. 1980 Guide to the preparation of a geological hazards map. Resources Analysis Branch Report 1980-04-17, British Columbia Ministry of Environment, Victoria, B.C. SAS Institute Inc., 1984. SAS/STAT User's Guide, Version 6, Fourth Edition. Vol. 1, Cary, NC: SAS Institute Inc., 943 pp. - 138- Schaefer, D.G. and Nikleva, S.N. 1973. Mean precipitation and snowfall maps for a mountainous area of potential urban development. In Forty-first Annual Western Snow Conference. Denver. Colorado, pp. 58-67. Schroeder, W.L. 1985. The engineering approach to landslide risk analysis. In Procs., Workshop of Slope Stability: Problems and Solutions in Forest Management General Technical Report PNW-180, US Forest Service, Portland, Oregon, pp. 4350. Schlesselman, J.J. 1982 Case-Control Studies: Design, Conduct, Analysis. Monographs in epidemiology and biostatistics. Oxford University Press. New York. 354 pp. Schwab, J.W. 1983. Mass wasting: October-November 1978 storm, Rennell Sound, Queen Charlotte Islands , British Columbia. Min. For., Res. Note 91. 23 pp. Sidle, R.C. 1985. Factors influencing the stability of slopes. In Procs.. Workshop of Slope Stability: Problems and Solutions in Forest Management General Technical Report PNW-180. US Forest Service, Portland, Oregon, pp. 17-25. Sidle, R.C. 1991. A conceptual model of change in root cohesion in response to vegetation management. J. Environ. Qual., Vol. 20. 52 pp. Sidle, R.C. 1992. A theoretical model of the effects of timber harvesting on slope stability. Water Resour. Res., 28(7): 1897-1910. Sidle, R . C , Pearce, A.J. and O'Loughlin, C.L. 1985. Hillslope stability and land use. Water Resources Monograph 11., American Geophysical Union, 140 pp. Swanson, F.J. and Dyrness C.T. 1975. Impact of clearcutting and road construction on soil erosion by landslides in the Western Cascade Range, Oregon. U.S. Dep. Agric. For. Serv. Reprinted from Geology, 3:393-396. Swanson, F.J. and Swanston, D.N. 1977. Complex mass movement terrains in the western Cascade Range, Oregon. In Rev. Eng. Geol. vol. 3, Landslides. Geol. Soc. Am., Boulder, Colorado, US, pp. 113-124. - 139- Swanson, F.J. Janda, R.J., Dunne, T. and Swanston, D.N. 1982. Sediment budgets and routing in forested drainage basins. General Technical Report PNW-141, Pacific Nortwest Forest and Range Exp, Stn., USD A Forest Service. Swanston, D.N. and F.J. Swanson. 1976. Timber harvesting, mass erosion, and steepland forest geomorphology in the Pacific Northwest. In Geomorphology and engineering. D.R. Coates (editor). Dowden, Hutchinson and Ross, Inc., Stroudsburg, Pa., pp. 199-221. Swanston, D.N., Swanson, F.J. and Rosgen, D. 1980. Soil mass movement. In An approach to Water Resources Evaluation of Non-point Silvicultural Sources. Chapter V. Mulkey, L.A. (ed). USDA Forest Service, Washington, D.C., V.I.V.49 pp. Thompson, G. and Golding, D.L. 1990. Mass wasting in GVWD watersheds: an inventory of events and an evaluation of the Aqua-Terra Classification System. Faculty of Forestry, University of British Columbia. 85 pp. Thurber Engineering Ltd. 1991. Geotechnical assessment of 1990-1991 landslide events in Greater Vancouver Water District watersheds. Report to Greater Vancouver Water district. Vancouver. B.C. 28 pp. Walker, H.S. and Duncan, D.B. 1967. Estimation of the probability of an event as a function of several independent variables. Biometrika. 54(1-2): 167-178. Ward, T.J. 1976. Factor of Safety approach to landslide potential delineation. Colorado State Univ., Ph.D Thesis, Fort Collins, Colorado, 119 pp. Ward, T.J. 1985. Computer-based landslide delineation andriskassessment procedures for management planning. In Procs., Workshop of Slope Stability: Problems and Solutions in Forest Management General Technical Report PNW-180, US Forest Service, Portland, Oregon, pp. 51-57. Ward, T.J., Li, R.M. and Simons, D.B. 1981. Use of a mathematical model for estimating potential landslide sites in steep forested drainage basins. IAHS Publ. 132., pp. 21-41. - 140- Wilkinson, L. 1990. SYSTAT: The System for Statistics, SYSTAT, Inc., Evanston, IL. Wilson, D. 1985. Subjective techniques for identification and hazard assessment of unstable terrain. In Procs.. Workshop of Slope Stability: Problems and Solutions in Forest Management General Technical Report PNW-180. US Forest Service, Portland, Oregon, pp. 36-42. Woolf, B. 1955. On estimating the relation between blood group and disease. Annals and Human Genetics. 19:251-253. Wu, T.H. and Swanston, D.N. 1980. Risk of landslides in shallow soils and its relation to clearcutting in southeastern Alaska. Forest Science 26(3):495-510. -141 - APPENDIX A - FORM.FOR FORTRAN PROGRAM Appendix A: Flow chart of the FORM.FOK program INPUT DATA * READ IDRISI ELEVATION MATRIX • CALCULATE LAGRANGE POLYNOMIALS 4 CALCULATE SPATIAL DERIVATIVES • WRITE IDRISI IMAGE FILES - 142 - APPENDIX A - FORM.FOR FORTRAN PROGRAM Source code of the FORM.FOR PROGRAM c c C C C C C C C C C FORM.FOR — • A FORTRAN-77 PROGRAM TO CALCULATE THE SPATIAL DERIVATIVES OF AN ELEVATION MATRIX USING LaGRANGE POLYNOMIALS IT TAKES IDRISI ELEVATION IMAGE IN ASCII FORMAT AND CALCULATES TERRAIN ATTRIBUTES SUCH AS SLOPE, ASPECT, PROFILE AND PLANFORM CURVATURE ORIGINAL PROGRAM BY D.E. HOWES K.O. WRITTEN NIEMANN AND C C MODIFIED BY C GYULA GULYAS C DEPT. OF FOREST RESOURCES MANAGEMENT C UNIVERSITY OF BRITISH COLUMBIA C C REFERENCE: C NIEMANN, K.O. AND HOWES, D.E. 1992. SLOPE STABILITY C EVALUATIONS USING DIGITAL TERRAIN MODELS. B.C. C MINISTRY OF FORESTS. VICTORIA, B.C. 28 P. C C C C DIMENSION ARRAYS PARAMETER (M=310) INTEGER ELEV(M,M),NROW,NCOL,I,J CHARACTER*15 INPUT,OUTP1,OUTP2,OUTP3,OUTP4 REAL D,E,F,G,H,GS,SQG,SQH,TMP,SLP(M,M),ASP(M,M),DWNC(M,M),ACRC(M,M) $,TMP1 C C C OPEN FILES AND GET INFORMATION C PRINT*,'NAME OF INPUT ELEVATION IMAGE (M):' READ(*,1000) INPUT PRINT*,'WHAT IS THE RESOLUTION OF THE IMAGE?(M): READ(*,*) GS 1000 FORMAT(Al5) 7 PRINT*,'ENTER THE NUMBER OF COLUMNS:' READ(*,*) NCOL IF(NCOL.GT.M) THEN PRINT*,'SORRY, THE MAX. NUMBER IS:',M GO TO 7 END IF 8 PRINT*,'ENTER THE NUMBER OF ROWS:' READ(*,*) NROW IF(NROW.GT.M) THEN PRINT*,'SORRY, THE MAX. NUMBER IS:',M GO TO 8 END IF C OUTPl='SLOPE.IMG' OUTP2='ASPECT.IMG' OUTP3='DOWNSL.IMG' OUTP4='ACRSL.IMG' OPEN(3,FILE=INPUT) OPEN(5,FILE=OUTPl) OPEN(7,FILE=OUTP2) OPEN(9,FILE=OUTP3) 1 - 143- APPENDIX A - FORM.FOR FORTRAN PROGRAM OPEN(11,FILE=OUTP4) C C READ IDRISI ELEVATION MATRIX IN ROW AND C COLUMN ORDER C DO 10 I=1,NR0W READ(3,*)(ELEV(I,J),J=l,NCOL) 10 CONTINUE C C C START MAIN LOOP C C THE EDGE ELEMENTS WILL B E LEFT OUT OF THE CALCULATION C DO 20 I=2,NROW-l DO 30 J=2,NCOL-l D=(((ELEV(I,J-1)+ELEV(I,J+l))/2)-ELEV(I,J))/GS*GS E=(((ELEV(1-1,J)+ELEV(I+1,J))/2)-ELEV(I,J))/GS*GS F=(-1*ELEV(1-1,J-l)+ELEV(1-1,J+l)+ELEV(1+1,J-l)-ELEV(1=1,J+l))/ $4*GS*GS G=(-1*ELEV(I,J-1)+ELEV(I,J+l))/2*GS H=(ELEV(I-1,J)-ELEV(I+1,J))/2*GS C C SOG=G*G SQH=H*H C CALCULATE THE SPATIAL DERIVATIVES C C S L P ( I , J ) = (SQG+SQH)* *0.5 C IF(H.EQ.O.O.AND.G.EQ.O.O) THEN ASP(I,J)=361.0 ELSE TMP1=ABS(H/G) TMP=(ATAN(TMPl))*57.958 IF(TMP.LT.O.O) THEN ASP(I,J)=361.0 ELSE IF(H.GT.O.O.AND.G.GT.O.O) ASP(I,J)=270.0-ABS(TMP) IF(H.GT.O.O.AND.G.LT.O.O) ASP(I,J)=90.0+ABS(TMP) IF(H.LT.O.O.AND.G.LT.O.O) ASP(I,J)=90.0-ABS(TMP) IF(H.LT.O.O.AND.G.GT.O.O) ASP(I,J)=270.0+ABS(TMP) END IF END IF C C IF(G.EQ.O.O.AND.H.EQ.O.O) THEN DWNC(I,J)=0.00001 ACRC(I,J)=0.00001 ELSE DWNC(I,J)=-2.0*(D*SQG+E*SQH+F*G*H)/(SQG+SQH) ACRC(I,J)=2.0*(D*SQH+E*SQG-F*G*H)/(SQG+SQH) END IF 30 CONTINUE 20 CONTINUE C C WRITE IDRISI IMAGE FILES C DO 40 I=l,NROW DO 50 J=l,NCOL WRITE(5,1111) S L P ( I , J ) - 144- APPENDIX A - FORM.FOR FORTRAN PROGRAM WRITE(7,1111) ASP(I,J) WRITE{9,1111) DWNC(I,J) WRITE(11,1111) ACRC(I,J) FORMAT(F10.2) CONTINUE CONTINUE 1111 50 40 C C CLOSE FILES CLOSE(3) CLOSE(5) CLOSE(7) CLOSE(11) STOP END - 145 - APPENDIX B - CATCH.FOR FORTRAN PROGRAM Appendix B: Flow chart of the CATCH.FOR program START INPUT DATA READ IDRISI ELEVATION MATRIX START NEW FLOWLINE IS CURRENT POSITION THE LOWER RIGHT CORNER? Yes WRITE IDRISI IMAGE FILE No IS CURRENT POSITION AN EDGE ELEMENT? No CAN FLOWLINE ADVANCE? Yes INCREASE CATCHMENT AREA BY ONE STOP - 146- Yes No INCREASE CATCHMENT AREA BY ONE DECREASE CATCHMENT AREA BY ONE APPENDIX B - CATCH.FOR FORTRAN PROGRAM Source code of the CATCH.FOR PROGRAM C c C C C C C C C C C C C C C C C C C C C C C C C CATCH.FOR A FORTRAN-77 PROGRAM TO DETERMINE THE LOCAL CATCHMENT AREA AT EVERY ELEMENT OF AN ELEVATION MATRIX. IT TAKES IDRISI ELEVATION IMAGE IN ASCII FORMAT AND CALCULATES CATCHMENT AREA BELONGING TO ANY PIXEL IN THE IMAGE ORIGINAL PROGRAM WRITTEN BY DR. LAWRENCE W. MARTZ DEPT. OF GEOGRAPHY, UNIV. OF SASKATCHEWAN SASKATOON, CANADA S7N 0W0 MODIFIED BY GYULA GULYAS DEPT. OF FOREST RESOURCES MANAGEMENT UNIVERSITY OF BRITISH COLUMBIA REFERENCE: MARTZ, L.W. AND JONG, E . 1988. CATCH: A FORTRAN PROGRAM FOR MEASURING CATCHMENT AREA FROM DIGITAL ELEVATION MODELS. COMPUTERS & GEOSCIENCES 5(14): PP. 627-640. DIMENSION ARRAYS PARAMETER (M=310) INTEGER ELEV(M,M),IROW,ICOL,NROW,NCOL,IR1,IC1,DEPR(M,M) CHARACTER*15 INPUT,OUTPl,OUTP2 REAL CA(M,M) C C OPEN FILES AND GET INFORMATION C PRINT*,'NAME OF INPUT ELEVATION IMAGE (M):' READ(*,1000) INPUT 1000 FORMAT(A15) 7 PRINT*,'ENTER THE NUMBER OF COLUMNS:• READ(*,*) NCOL IF(NCOL.GT.M) THEN PRINT*,'SORRY, THE MAX. NUMBER IS:',M GO TO 7 END IF 8 PRINT*,'ENTER THE NUMBER OF ROWS:' READ(*,*) NROW IF(NROW.GT.M) THEN PRINT*,'SORRY, THE MAX. NUMBER IS:',M GO TO 8 END IF OUTPl='catch.img * OUTP2='depr.img' OPEN(3,FILE=INPUT) OPEN(5,FILE=OUTPl) OPEN(7,FILE=OUTP2) C C READ IDRISI ELEVATION MATRIX IN ROW AND C COLUMN ORDER C DO 105 IROW=l,NROW READ(3,*)(ELEV(IROW,ICOL),ICOL=l,NCOL) C 105 CONTINUE - 147- APPENDIX B - CATCH.FOR FORTRAN PROGRAM C C C SCAN ELEVATION MATRIX... C DO 5 IROW=l,NROW DO 4 ICOL=l,NCOL C SET THE CURRENT POSITION OF THE FLOWLINE AT THIS ELEMENT IRl=IROW ICl=ICOL C C CHECK IF CURRENT FLOWLINE POSITION IS EDGE ELEMENT 2 CALL EDGE(*10,IR1,ICl,NROW,NCOL) C C I F NOT EDGE ELEMENT, FIND NEIGHBOUR TO WHICH SLOPE IS C GREATEST C CALL ADVANC(*3,IR1,ICl,ELEV,M,IR2,IC2) C C IF FLOWLINE ADVANCES, INCREMENT CATCHMENT AREA BY ONE C AND CHANGE CURRENT POSITION OF FLOWLINE C IF(IR1.NE.IROW.OR.ICl.NE.ICOL) CA(IR1,ICl)=CA(IR1,ICl)+1.0 IR1=IR2 IC1=IC2 GO TO 2 C C I F FLOWLINE CANNOT ADVANCE ( I . E . IN A DEPRESSION), DECREASE C CATCHMENT AREA BY ONE AND BEGIN NEW FLOWLINE C 3 CA(IR1,IC1)=CA(IR1,IC1)-1.0 DEPR(IR1,IC1)=1 GO TO 4 C C I F CURRENT FLOWLINE POSITION WAS AN EDGE ELEMENT, C INCREMENT CATCHMENT AREA BY ONE AND BEGIN NEW C FLOWLINE 10 IF(IRl.NE.IROW.OR.ICl.NE.ICOL) CA(IR1,ICl)=CA(IR1,ICl)+1.0 4 CONTINUE 5 CONTINUE C C WRITE LOCAL CATCHMENT AREA MATRIX TO OUTPUT FILE DO 13 IROW=l,NROW DO 12 ICOL=l,NCOL WRITE(5,103) ABS(CA(IROW,ICOL)) WRITE(7,104) DEPR(IROW,ICOL) 12 CONTINUE 13 CONTINUE 103 FORMAT(F10.4) 104 FORMAT(12) CLOSE(3) CLOSE(5) CLOSE (7) STOP END C C************************* SUBROUTINES * ************** *********** SUBROUTINE EDGE(*,IROW,ICOL,NROW,NCOL) C CHECKS I F AN ELEMENT IS AN EDGE ELEMENT INTEGER IROW,ICOL,IR,IC,NROW,NCOL DO 2 IR=-1,1 DO 1 IC=-1,1 IF(IROW+IR.LT.1.0R.IROW+IR.GT.NROW.OR.ICOL+IC.LT.l.OR.ICOL+IC $.GT.NCOL) RETURN 1 1 CONTINUE - 148 - APPENDIX B - CATCH.FOR FORTRAN PROGRAM 2 CONTINUE RETURN END C*************************************************************** SUBROUTINE ADVANC(*,IR1,IC1,ELEV,M,IR2,IC2) C FINDS GREATEST SLOPE BETWEEN ELEMENT AND ITS NEIGHBOURS C AND CHECKS I F SLOPE IS POSITIVE OR NEGATIVE INTEGER IR,IC,IR1,IC1,IR2,IC2,ELEV(M,M) REAL SMAX,SLP DO 2 IR=-1,1 DO 1 IC=-1,1 IF(IR.EQ.O.AND.IC.EQ.O) GO TO 1 CALL SLOPE(IR1,IC1,IR,IC,ELEV,M,SLP) IF(IR.EQ.-l.AND.IC.EQ.-l) THEN SMAX=SLP IR2=IR1-1 IC2=IC1-1 GO TO 1 END IF IF(SLP.GT.SMAX) THEN SMAX=SLP IR2=IR1+IR IC2=IC1+IC END IF 1 CONTINUE 2 CONTINUE IF(SMAX.LE.O.O) RETURN 1 RETURN END C*************************************************************** SUBROUTINE SLOPE(IROW,ICOL,IR,IC,ELEV,M,SLP) C C CALCULATES RELATIVE SLOPE MAGNITUDE BETWEEN ELEMENTS C INTEGER IROW,ICOL,IR,IC,ELEV(M,M) REAL SLP IF(ABS(IR).EQ.ABS(IC)) THEN SLP=(ELEV(IROW,ICOL)-ELEV(IROW+IR,ICOL+IC))/(2 * *0.5) ELSE SLP=(ELEV(IROW,ICOL)-ELEV(IROW+IR,ICOL+IC)) END IF RETURN END C* ****** * * ****************************************************** - 149- OIO(OOOT-(<)Oi-0(00>T-i-CM(MT-T-0)OT-100(Min<0«-T-tMOlOi-CON t- u < u C M T - T- C M T - •<- — I (0 (O^-nnnnnnnnnnnnT-nnnT-ni-csinnnncoT-T-nMnnn o OT © • " - O O T - O O T - T - O O O T - T - T - T - T - O T - O T - O O O O O O T - T - O O O I - O 8 o CD Q. UJ a C M C M C O C O I O C M C M I O I O C M C O C O I O C O C M C O I O C O C M C M C M C N J C M C M C O C O C M C M C M C ^ OT 2 O NM(OT-ronnT-(Ni*)T-i-T-i-T-tNMMT-MT-i-i-T-r«j(M'-'-T-T-rOT-(on O Q. E-t < Q r-n«N»-ntMT-(OT-T-Nn(Ni'-tOT-MnT-T-nT-cv)'-(N(M(oni-T-T-T- 1 Q. CN^co^T-^^^^^<o^^^^TrcocNiT-^cocococ\icO'r-CNi'«j-'*'*CNi-*'i3 X Q (0 5 w C M 0. T- O O T - T ~ O O •<- O O O O O T - O - r - T — O Y— ^— V - T — O O T — T— O O T — (0 UJ o u. co m co co > m CO oo CM co •* UJ co CM CO o> _l o T— ^~ UJ o <o UJ 0. CO 00 O CM to o o o> o> oo CD o CM m CM 1c—o o> co m m Oo oo o> O) O "o* - o> 00 * " m CO o> o> o CO m CM co c00o 00 CO oo CM CO CM o co m CM TlCO lO m m CO CM o c- co m l»CD O CD •>*• o CO m o a> "* ~ T - ) CM m CO o O) co m CM o O co O CO oo m co m O) lO CO O CM 00 (O coo C ^— r-CO co o o o o o C M o O C O C O C O C O m o > C M C M oo o oo o> in T— o o •* ' a ' r-«° C OC O in CM CO oo' •<* 1^ *— ^— T— '* '" co CM CM CM CO CM CM CM co CM CM CM CM CM CO CO CM CM m CM 00 00 00° I D — _l OT o o o o o o o o o o O o o o o o o — o ^_ o o o O O o o o o UJ Q -j W UJ >< 0. •• U — o m o o> 1C—O o o CM CO CO ^— CM CO CO CO oo CM o T— CM co O) CM o r-O) CO oo T— CM T - O) o CM CO oo T— o M co co CO uo m CO CO CD CO CO 00 o> oo CM co m CO CO oo CN C CM CM CM CM CM CM CM CM CM CM CM CM CM CM CO N N CO •<- O CM CD m o O N i n O W O l N * ^ ' o o i - i - M o m i O f l o r a00o o T - 150- 00(Oa>Sinr^N!OOi-0)(D(0'*rtr-0)rOT-<00(DN'<l (0'<t-mMT-T-rt(DOO , Mi- O N C O CO iCM T- T~ O T- < o rtneonnT-nnrtnnnnNT-i-MT-NnnnnnnNNnrtnncvinn CO O CO T— T— O T O — T— T— OOi— O O O T T— O — T — T - O O O O O T O o — o o o o o o 8 o m CO CO CO T— I O Q. UJ Q CO inin^c\icNic\iiocococococococ\i^csiincN4csicococococococ\ic^ C T - C O C M C O C O T - C O C M C O T T C M O T - C O C M C O C N C O C O C N J C O - T - C C O O u. O D£ Q. EH < C O C \ I C \ I T - T - C O C O ^ T - C O C O C O C O C O C O C O C N C O T - ^ T - C O C N I T - C N I C N T - C O T - T - T - C O T - Q U z Q. ^^co^^^*T-T-^^r\iincNi^^^^^^Tr^T-coTj-TT'<i-'*'t-T-r>i'«*-T- u X H Q CO Q o CM T- T-OT-T-T-t-T-Ol-T-OOT-O o O O T- W ft T- O T - cm < OL 2 O O ' ' — T - T — O O ' - O T — O T — T - T - O O T — O T — - ^ - T — T — T — T — T — O T — T — T — T — T - T - T — O C0 UJ K O (DiDNOOiDooneoi-i-oionN^i-iooiDinnino stMiosooi-nfflO) inooiaoMsnoiN 4O 4Ou)io<M(Dno)ioo)too)in(oo)(NnNcoo>ONOi<) > ,,iO Tr^-'»-oo)r^cDO>oocD'*T-T-ooo>oioo)o»i ~ > > > UJ CO T-T-nO)NNO)CO(00(DW'*T-om( )NO)NCO«OT-*(<)fCMT-tOO)(N<0(0 , S° x-nc»owNin«NC»(Oi^<qMC»c\ici><DCNNro*NNONcq co'Tr^^coior^odioo'oo^odoco^cOTf O J C s | T - T - C 0 C N C 0 C N | C \ | c 0 C O C N C N C N C 0 C N c 0 C N C 0 C 0 T f ^ C N I C \ I M CO o o o o o o o o o o OOOOOOOOOT-t-'r-OOOT-T-OT-OOi-T- UJ Q (0 CO O CO CO N m o o f f l i n o o N n m N O hOntMO(0NONtT-0)*nN O ooo5i--<-coior-~o)OT-co O) O — I CO o .) K)(OUlNCM(OONnN(MlOMn<0 0000C»C»O>O>O>C»^-'T-T-T-T-CNCS|C 0^t O oot-T-cMC\icococo-<j-ir>iomcocD -_ _ . „-__._._._._._. ...... — cococococococococococococococococococococococococo^TfTfTfTfTr^TtTf UJ , 151 X o nNNOW^coofflinrin-oinNNNnoiooTtT-T-nmon^oincoT-iD (M T - O O T - T - C O CM T - T - CM CM CM < CO -J o CO O O O O T O - I - O O O O O O O O O O O T O - 0 0 0 0 0 0 0 0 0 0 0 0 N o o m T-T-nnT-T-rtnT-T-nnT-i-T-T-T-wnrtrtnonrt OL 111 Q CO Z o ( O T - r t n o T - n n i - T - T - n T - T - n T - T - n T - r t ^ - n M n T - N N T - e o i - T - T - M O OL O z nr-NT-i-T-nT-conNNNnT-noimpiCNT-NT-i-^-N^-NT-'-nNroN 1 Q. tf)^^^T-CO^CNCN^*-COCN^^COCN*^T-COCNTrCN^^-*-*''+''*'*'*CM'* •• U X H CO Q H Q W ft ft <« CM Q. 2 T - T - T - T - O I - T - O O T - O T - T - T - T - T - T - T - X - O T - O O O CO UJ OC oo CO O O) h- o CO o o m CM CO C MV O) O C o CM rc>o CM O 00 » co o> 00 m co 00 m o> co _l UJ 1231 u. CM r*- oo CO o o> o o CM CO o co oo CD o o> o> oo m 1132 1049 O Si Q o M o co CM CM CM co CM CM CO o CM 1— C CM ^— CD CM 001 O oo oo 00 oo 00 OO oo 00 oo CO o> oo oo M Oi O) CO co Ti- o U) CO C o o CO o> CO m o> CO oo m m CM CO O co CM CO r- CO o> CM m ^— e n CO 00 cq CM o co o> co m CO C cq ur>CO o m o> co c o C M o o > Tt T— c o iri iri d oi c o C O C O d « r i d C O CM co d CO O CM co CO oi CO CO CO •<* CO CM d co co co co CM co d co co CO CO co co CM co CM co co C CM ^— O C O c o c o c o C M c o c o _i CO ,_ o o o o ^— ^ _ o "»— o o O o o o o o o o o o o o o o o UJ Q. UJ o -i CO Jo*(MONO>T-(D***^'-0010iniOO)*N( )OnO)(ONO)*( )fM(OnW S ^ - I O I O C D C O C O C D I ^ ^ C O ^CMCMlOU^COCDCOr^COCOOOCOOOOOOO^COCO^lO , O O ftT-T-T-T-^T-X-T-T-T-T-T-r-T-T-T--I-T-'T-T-T--r-^V-T-T--I-T--I--r-T-T-' -TI - 152- , NN<0*NM'OncNC0O'i-OT-(0Oir)inT-N(0OV0T-OtNOU)l0T-(0T-(0N _ CM * ~ I T ~ CM CM o < o CO Tj- T~ T - T~ * ~ T - niNn^-nnnnNrOrtnnT-NT-T-^nriironn'-NnnnconNNtNi- O (0 O O I O - O O O O O O O O O O O T O - O O O O O O O O O O O O O O O T O - N u 0 m O T - T — C O C O ^ ^ ^ C O C O C O T — T— C O T - T — coco^ T - C O C O C O C O T — C O T - T - O OL HI Q CO C M C M C N C M ( N C M C M C M C M C M C > J C 0 < \ I C M ( N C O C M * C M C M C M C N C M W u u. 0 COCOT— CO CM COCOT— C M C O C O ^ C O T - T - C O C O T - C M C O C O C O C M ^ — T - C O C O - I - C M ' I — co OC 01 CO C O C O C O C O C M C O C O - r - C M T - T - C O T - C M T - T - C M C O C O CO CO CM Q u z OL u OS X ICO Q H Q 55 w 0^-T-T-T-OOT-r-^-T-T-^T-T-^T-O^T-OOt-OT-OT-^>r-^0 CM OL 3 O T - O O T O - O O T OOl-0'<--<-'<--<-T— - OH T - O O l - O T - T - T - T - T - T - T — co HI Of o CM CM I*- oinMntMnincoot-ooiMOffl^BrtCMW O (ONCO'*T-NCOOO'<f(MMOO'*NST-N'*00 co r~- CM r~- T J -c o ^ T - n i o c o i o i o i o t M c o o o s N i o o i w n ' UJ moOOOOCOCOmcOOOh-COCOCOO — I UJ T - to m M CO T - c O T r i n c M c o c o r ~ o o T - c D c o r - - i n c M C D T j - h ~ o t ^ c o c o o > T - r - ~ c D h~ co co co m Ul C ^'r^oicbuiuSTtcocoinuicor^cOT-^ i- C D 00 1^ 0. o _1 CO o o O O O O O O O O O O O O O T - O O T - T - T - T - O O O O O O T - T - O * - UJ a -I CO .1 —I 1-- in cOT-coo>T-cNO)r-CMcor-oomin cor>~o>cOT-coo)tDO>T-i^o)'*ooioo UJooo5CNo>o)o>T-mcoooo)T-r-cMcomcoco iflifllOMOOaoOOIOOOOOrrnn i^i^r^r^r^f^r^oooooooooooooooooo XW(DOOOOIO)OOrT-T-NNnOf " - 153 m^rOONMDO^OIO^OOOIN^OT-OOi-onOO^O^tDOT-ifl CO O CN T - C O ti• > * • C M < U co N N « r t n N e o e o < i N n n n T - n < o < O M i O T - * n e o N ^ - i - ( \ i N n < < ) i - T - N i O CO O O O O O O O O O O O O O I - O O O O O I - O O O T O - T - O O O O O I - O O N o o ta co^— co^ cococo^— - T— C O T — T— T— C M C O T — T— T— C O C M C O C O C O T — CM CM CO Q. UJ 0 CO z 1 C O T - ^ C O C N C N C N C M T - C M C O C O T - C M C O ^ ^ C O T - C O C M C O C N T - C O C O ' I - ' I - ' I - T - C O T - C N o u. o EH Q. o z < ( O T - T - T - T - ( M C N T - c 1 ( N r g < - » - r t T - i - N C M T - T - T - r > | T - r - » - T - T - M T - N C O T - r t Q 1 Q. •• U 01 3 X H Q H CO Q W ft ft CM 0. < 3 O O O T - 1 - O O T- co UJ O > UJ _l co CO CD cn OJ CM 00 cn co CO •<* oo o r- CN lO CN CO m ^— oo co CO co O J 1— r*- oo co OJ co oo in CO •«*• OJ CN CD CO co CO 0CM CM OJ OJ r•<r- CM o O o oo CO or » CD CM CM ^— oo co m 0 l~- lO m m •* C CO co lO CM h- CO CO hco in oo T— l>CO CO oo co " T ^ CM o> CO CO T™ co CO CM CO CO co t>co CN ^ _ co 00 CM CO CO CO OJ CM CN mm 00 o o m m 00 m o o q C O co C O co co O I co C O C O C M C N M C O C M co C D co UJ d ^_* T— d d d CO oi CO d d cri CM d iri CM CM OJ CO co ^— OJ CO m iri C O co 0. CM co CM co co co lO CM CM CM CN CN CD co CM CM CM T— CM co co •* CO o co co CM co co _l CO ^_ ^_ t— o T— o o o o o o o o ,_ o o o o o o o ^— UJ O O O O O UJ g -i CO J ^_ o OJ OJ o o o OJ m co m CN CO CO m OJ OJ m h^— o ^— OJ co OJ CM oin o^— co O co i n CO OJ oo 00 ^— OJ oo 00 CO CM CO m r>» •* 00 OJ CM C D co C O C O o o O J O J C M co co co co C O C O O J C N C M co co o o o o C oO — co co oo co oo oo oo oo co OJ OJ OJ OJ OJ OJ OJ OJ OJ OJ OJ OJ OJ OJ o O O o o o co CN CN CM CM CM CM CM CM CM CM CM UJ 0. T" - 154- S*COn(\|*(')<t(OOcOO<ONO>0(<)*Oi-i-COCOT-*N*0)COOC\l'-lOTC M O O T - T - C M C o < D CM T - CM i- CO C M r- T - T - in N i- 3 o CO O O O O O O O O O T - O O O O O O T O - O O T O - O O O I O - T - T - O O O N o O m nT-nncocoT-nt-nnNnT-nconnT-i-nT-nNn^-'i-'i-'-N'-T-nn 0. HI O CO (OCM^rtCMi-tMNCOCM CMCOCOCMCOCMCOCMCMcOCO CMCMCOCOCOCMCMCOCOCMCO 2 Q T-rMr-ni-i-T-rti-T-(<)T-T-n?-nf\|T-ni-n(<)n u u. o on a. En n n r t T - ( N M n c \ i n c \ i n c o ( N n i \ ' - ( N i r t ( N ( N c o ^ i ( N T - T - r t ' r - T - r o n ( 0 ' i - T - T - 0 z •• 01 u OC 3 I- X H co Q 55 Q W cu a, CM OL 3 O O t T " T - - T— T— O O O T - T - O - T - O T - O O T • < - O O T- T- O - CO HI Of o O CM CO T - CO COIOtOO>(DO)MOrtfM'*CM«>'*CO'*0)OC»<DN*NlO^-0)C»eO* C m T J - r~- c o o o HI o o c o c D t n m o o o o c o c o t ^ c o O I S N O I T - O N C O W _l r 5t o _l CO HI oo co o T I - r~ oo m co C M «o C O C O O) o in eo ed * T— CO CM T- O O O i- CM co CD •*' T— CM CM CM _ o oo O *— o co d m oo c o CM CM oo CM CO CM 00 co' d co t— CM T™ CM CM o o o o O o o o o o oq r c o co' o co o o CD iri oo m ' m CM m O O o CM rco CM CM CO o o oo oo c o oo o> CO o CM CM CM c o CM CM CM CM CM o oo CO c o o> CM in c o o> o o> o> CO T— c o c o CO iri CM 00 oo d — CO co co CM c o CM o x — o o o o oo m ^— CO CM T— rCO c o oo I - oo Y— a> Y — Y— •f m CD CO c o c o c o c o c o CO c o CM CM CM CM CM CM CM g _l CO —I HI a I 5 HI < CO CM CD 00 o> o roCM CM >< g co m o o m ^ - CM Y— oo o> o CM c o m o o o CM CM CM CM CM CM CM co o 00 CO o> o> oo CM m CM co CM oo CM CD c o r - r- 00 o> c o o> CM o 1— T— ^— CM CM CM CM CM CM CM CM CM CM CN CM CN CM - 155 - o o co CM CM s CO O ^ i - I O N < 0 ( 0 O 0 1 0 0 O i - N * M D O O ( 0 m ' t n ( D T - N ' - N M 0 O ( 0 * * C N C M t - t - C M C N ^ - T - CO •«-•<- •<- < O C O C M C M C M C N C M C O C M C N C N C N C M C M C O T - ^ C M C O ^ C N C N C N C O C O C M C O C O C N C O C N N CO -J O o O O O O O O O O O T - O O O T - T - O O O T - O O O O O CO 8 o m C O C M O J T - T — C O T - T — T— C O C O T — T - T — T — 1— T — T— CN X \— 0. UJ a CO ^ C O C O C N ^ * ^ C O ^ C O C N C O C N * T - * - C N C M ^ C N C O C O ^ ^ C M ^ ^ C O - * C N C O C O C N 2 o u u. 0 01 C M ^ ^ C N T - ^ ^ C N C O C O T - T - C N T - C N T - C N C O C O T - C O C O C O T - C O T - C N C O T - C O C O C O ^ C O C N ^ C O ^ C O C O T - r - T - e O C O C N C N C N C O C O C N ^ T - C O ^ T - T - C N T - C O ' T - T - C N C N ' r - T - C O T - o z Q. ^ ^ C O C N ^ ^ ^ C O ^ C N C N - * ' * C N - * ' * , * ' * C N C N C O C O C N C N c O ^ - ' * ' * 1 * - * C N ' * l O - * 0£ 3 H CO Q 55 o w ft ft o < o UJ -I UJ c o 00 CO o co ^; o fco CO o o CO CO 0 0 hCO CM CM 00 co c o h - co o o co CO OJ o CO oo oo OJ CO co m c o CM o CN O CO m CN o OJ 0 0 OJ c o O J CO o CO CO in *~ o> Q. co oo o oo m i— o co — t oi oi co d co CN CN CO CO c o O) oo co •* m oq o OJ iri ^— CM CO od d c o c o c o CO CM CM c o oo oo m CN oo CM co T _ iri CO o o -I CO o o ^— o o o ^— o o o — i OJ m co d d CM CM *— ^— O co co m co c o CO CO •<r O J c o CO oo T— CO CM CM CO CO q co OJ T— co co iri co CM -<* c o CO CM o o o t«00 CO m CM T— rrOJ m m m CM CN CM T— o o CM •<r m CM rco m m CM m 00 m co m m CM CN d O oo lCO O— t OJ OJ T— o ^— o o o T— CN q co CO c o c o CM CN c o oi o o r~ CO 0 0 o ^— T— CO o CM ^— o CO o co IO OO CO d co •* o o o 00 ^— in CO CM OJ 00 r-oo c o CO CM CM o o m iri iri CO in o UJ o -I CO Ui X •• u H CO UJ u. En X o CM Q. < < a ^> OJ roo co co co CM CM CM OJ OJ 00 co CM CN O OJ CO CN m CM c o c o CO O J OJ CN co CM CM CM i CO — m o CM m m CM CM CM CM CM CM T— m OJ oo CO oo CO co o OJ CO o CN m CM CN 00 CM m CN - 156- T— oo — i CD CN 00 O CO co co CO CO CN CM CO CO c o OJ OJ o CO c o CM CM CM CO lO CO CM O lO CO N NCM T - T - CO CO CM O h- CM K) r f- CM T}- o) m o co coK) N i - IO CM CM T co CO - < O COCMCMCMCNCMCOCM CM CM CM CM •«- C N C M ' « - T - T - C M C M C M C M ' « - C M C O C M C M C M T - C M r- 3 o CO N O O m CM CM 0. HI Q CO z i Q O T-n(\inr-t-T-NWi-ncinnnT-nT->-c\i^rtn»-nNNT-n>-Nni-N U. 0 OC EH 01 C O T - C N C N C O C O C O T - ^ T - ^ ^ C O C O T - C O C O C O C O T - C O C O T - T - T - C M C O T - T - C O C N T - ' - Q 0 Z 5 01 OC 3 ICO Q ^^^^^CM^^CM^^^kO-*'*-*-*CM *'*-*'*-* *'*'*'*'* *'*'*CO-*'* , , , X H Q 55 O O T - - « - O T - * - O T - T - O T - 0 T - T - T — 0 O O T - T - T - T - O O T - T — T - T - T — O CM ^ ^ - ^ O O O T - T - O O T - 0 0 ' < - T — O O O O O O T - O O T - O O I - 0 0 0 - < — o o Ico UJ OC o 0>«Ot-lO'<4-COCMO«OCOCOinCM'*OOCMlO o > r r o u ) T - c M O ' * r - i - T - - * i n c o o o NO)T-COCOtNlOO)<000)(MO)f»10* > •uOJtwoi-oh-CROcoooi-roOTOT-roco , s s 0 » 0 0 0 > T - T - C O O T - O I - T - ' » - C O C M ' « - Ul UJ _ i f l M T - N T f o o o i i q q N ^ T - T f N N r t q c o co u ) N M- o o c o c o u ) n e d o ) r g u> e n o ) N . t S i i c O C O C O C O C O C O C M C M M ' C M N t i f C O C O C O C O C O o _l w CH CL CO < u T - O O T - I - O Y - O O O O O T - O O O T - 00VO'«tT-CMlOCMT-COT-<OCOC000O>O> aNcocooo(Ooqqr;oqncoiosoq iri CN ^ O) (O (D * Ol d O N N N CO •*COCO^T-CMCMCMCM'*COCMCM'*COCO T- T - O O O O O O O O O O O O O T - O O CO UJ o _l CO —j' (O CO O) O) CO O UJ CM CO CC OO C T CO CO «_» «U if o w (M m I OO CNW COO CM t- m o) mo roo - C N CM CO O) r T 1- | ^ U I " -«*- I— 1_» VI i«_» C 1M - >C 1O ! ) i-r- lO ^ N CO CO i - N s T J - VJ CO in W oo »/ in c o o w co yTi-T-nnTtTtNMocBOiooooi-T-MTfTnoffliotooT-T-T-T-Mn^ifl O.CNCMCMCMCMCMCMCMCMCMCMCMCMCMCMCMCMCMCMCMCMCMCMCMCMCM - 157- ~< N ( o n n n o i n ' - i o o ^ - N ^ o o n n i o n o ' t f f l o ^ O T - O T - f f l r t T - c O ' t i - i n _ • X « - • « - CO CM o T- co I< O C O C O C N I C N C M C O T - C N J C O C O C N C O C O C M C M C O ^ C O C M C O C M T - T - C M C O 3 o CO O O O O ^ - O O O O O T - O O O O T - O O O T- O O 8 o CD C M C M N C M N M N C N I M T - n r - n r t T - T - i - T - T - ' - ' - O ' - n i - T - r - T - T - n T - T - T - T - £L UJ Q CO C M C O C M C O C O C O C M C M 2 a CM CM C N C O C O C N C O C N C N C M C N C N C O C N C N T - C N N T - n n e o n n n n n i - N n n T - n r ' n N i - n n n M T - r t T - r t r t ' - n n N M u U. 0 0u1 z En C M T - C N C O ^ ^ C O C N I T - C N C O T - T - T - T - C O C O C O T - C O T - T - C M ^ C O C O X - T - T - T - C O C O T - C O Q 1 Q. O X CO Q Q 55 CM < w ft ft 0- CO < o o o o o o o o O o O O H CO UJ 0£ O UJ - J CD in CO CO -* O o o CD 00 m CO m m o m ^— — i C M o o C M o o o O C M C M 0 0 o> C co m co CO CD CO CO CM r- CD CO rCM C C~M N co CMm oo C O CO m o -xj- co CO oo oCO CO w CD CO o co CM D CO o rCM C CMCM co CMC M CM CM CMC CM O C CO O •* ^— T _ UJ CM CO CM CO m oco CO co co r»O co CM oo r- „_ O o CO CO o Oi O) rCO CM •* lO l O cq co o o co C D c n oo •<Jco co" CO co ai iri iri d 1^ CO CM — o M CO CM CM d d t o' C C M C M C M C M C M C M C O C O C M m m C M C M m C M CN m o UJ Q. ^—• •»— _l CO UJ Q o o O o o o o o o o o o o o T- CO CN ^— — t o CM m O - T— o O T CO O i— m CN o CO co o> •* co O CD co CO CO d iri 1^ CM co CN " CM m CO CM CO CN CM co •*— O o o o O o ,- to m CO 00 00 CO T— CO 00 r» CO CO CO CO CO CM CM CO O CO CO hCO CO CO oo CO CO CO Ol •f T— CO co a> Oi CO Oi CO O) CO Oi O CM CN CM CM CN CM CM CN CM CO co CO - 158- co o Oi oo r-i- oo oi CM co Oi r-- co _o co o oi m co co co Tf < ^ co i*- t o t- m o o o in m m CONN^CONOOOCOT-MOONfOrgiOT-N^OlONlO^lOlOtOCBT-ON TT T N rTjT— C M (Nl CM t - CM ^ T i- UJ - i - T - T - < O co co T- t - C O C O ( O N N ( O C O » - N N N N C O ( O N C O C O ( O C O T - N N N ( O N C O CO o CO O O O T - T - t - T - O O O O O O O t - O O T - T - O O T - T - O O O T - O O i - O T - O T - N o o a COCOCOCMCOT-*-T-f-COCM co OL UJ o CO CMCOCMCMCMCMCMCMCOCNCMCMCOCOCNCMCNCN Q CO CO CM C O T - C O C O C M C O T - T - ^ T - C O T - C O C O C M C M C M C O C O T - C O T - ^ C O C M C M ^ C M C O C O 0 01 OC o z E-t T - N ( 0 ( O T - N C O C O T - T - C O N C O N C O N N C O C O T - ( 0 ( O N C O T - C O C O C O C O T - C O N ^ - C O Q 5 CL CJ DC 3 •- X H Q 55 CO Q o< T - O O O O T - 0 0 ' « - T - T - T - Y - T - T - O O T - O T - O T - T - O T - T - T - 0 ' < - O I - T - CM OL 3 T - T - T - O O O O T - T - T - T - T - T - T - O T - T - O O T - T - O O T - T - T - O T - T - O T - O T - O CO UJ OC O cococococococoo>o>T-T- co co ^ at co o) )M i \C iOuT)r ^l Ool O) iCM iD C O I - C \ | 0 > { N l - C M(' o * -o *C M U>)n Ca Mn CM C Or -o* iC fOtOaOo T - M O ( 0 N C 0 1 0 S * » I > l . | C O ( O C O T - ( \ | T - N C O ( D S * * N N C > l * T f T - T - C O ( O T - N N ( 0 ( D ( \ I C O U O s s , • UJ - " " O CO O N CM _l UJ ( O N O I M W C O T - I O C O CO O) Ol CD CO I-CO VO CD CM CD COoCO 00 C D (M C O O) C O C O 0 0 CD I - 00 C Mo>CMOCOCOCNI^-lOCM"* O CM O CO CM O 0> o co i n ^ - ' N CO O) CO •T -TtcNi^T-O'^o^'r^iriTtcD-'i-r^r^-' o o o to c o c o c o c o m •lOCNCOCOCMCOCNCOCNCMCOCOCN^COCNCNCNCOCOCNC^ ^ r ^ ^ o o c o c o o c M o o o T - co o> s s UJ OL O - J CO I - O O ^ O O O O O O O O O T - T - O O T - O O O O O O O I — O O O T - O O - " — UJ o _l CO J 111 NNT-000)0)C0(0N«O l O CM lO CMCOCMC0'<crC\ICOI-O>C0 (M CO O O ) N o wcococM^*OT-cor -o>ooo>cNcocoT-T-cococMcoocoo>i>-ooiOT-r-(\ir -ococo S s s o inoui io«i o i uf fi lwi o( D o oi o ) at oo t)tTi u -T l2uonoiN n cwouni w u iuuiiuuinNi N flO i n)iaoci folci n i ooi D D-iN o si fslsisnNi o S sScNooOoT co - 159- X n o i \ i n ^ o < > | O i - N N a ) o i o N c o T - r - a ) i n i n ( o o n o T - T - CO C M C M • > CN co o < u 3 C O C O C O C O C O C O C N C O C O C O C N C O C O C O C O C O C O C O C O T - C O * - C N C O C O o to T - T — O O O T — O T - T - O O O O * — T— O O T — O * — T— T— O 8 o m OL til Q (0 I O I O C D C N C O C N C N C O C D C D C N C O C O C O I O C O C O C O C O C N C O T - C M C D C N 2 a <OCOCOT-COi-'<r-T-eOCOC0COCOCNT-CNCNC0COCO'^T-i- o a c O T - e O T - e O ' « - T - < o c o c O f - c N T - c N c o c N c O T - T - c o c N c o c o c N ' » U z Q. ' * ' * - * ' * C O T - C N - * ' * * C O I J - ' * C O C N ' * ' * - * ' * ' * C O - * - * ' « » C O 1 H (0 a CN Q. o o 3 O O O T - T - T - O T - O O T - T - T - T - O O T - T - T - O T - O O O T - (0 UJ o o UJ _J UJ CN CO CN CN CN c o CO r — CO O ^— CO m CO o co' o _l (0 o o *~ CM CO c o o r - CO m co CO r~- o> CM CM CO CM c o ^— CO co' CO i r i CN i r i CN i— CN CO CN CM c o CO CD CN c n m CO CO c o 111 oi a. c o o CM t» co o CD m co o o o o o m co C M oo CD hCN CO m CO m C M CO C M CO co oo oo o o o CO CO m oo m co m oo CO CN 00 00 m m co CM CN o oo m CO o CD o m CD r - CO ^— CN c o O O) m CO m CO CO CM i r i c o CM d CO CO CM CM x— O o o o CM CO CO 0O CM r - CN Oi m cq iri co o co o o •«* T— o o> '' 00 m — CN CN oo oo oo oo c o CM c o o T— m ^— o cq CN CO cd i— CM *~ " o o O UJ 2 _1 to M CO C CO C M CO CO •* co Oi 1*- oo o o CO CO o co oo co co C O 1— o C M lO co T— T— CO co co o CO oo CO m •«* •<* in a> co CO oo oo 0 0 OJ CO co o> C O o> CO OJ - 160-
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Terain stability assessment using logistic regression...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Terain stability assessment using logistic regression analysis for the Jamieson-Orchid-Elbow Creeks subdrainage,.. Gulyás, Gyula 1995
pdf
Page Metadata
Item Metadata
Title | Terain stability assessment using logistic regression analysis for the Jamieson-Orchid-Elbow Creeks subdrainage, Seymour River Basin, British Columbia |
Creator |
Gulyás, Gyula |
Date | 1995 |
Date Issued | 2009-01-20T17:24:59Z |
Description | The purpose of this research project was to develop a procedure for terrain stability assessment by applying case-control sampling and multiple logistic regression analysis, widely used statistical techniques in biomedical research and in epidemiology. The idea of applying statistical methods used in epidemiology to terrain stability assessment was based on the observation that landslides, like some diseases, are rare phenomena. The implementation of a terrain stability assessment based on these statistical techniques was expected to help understand the causeeffect relationships between landsliding and various terrain attributes. In contrast to the currently used approaches, the study procedure provided a quantitative tool to assess the risk of landsliding and to define the most important terrain attributes that contribute to soil mass movements. A case-control study of 20x20 m grid cells with average slope greater than 10 degrees was conducted on the Jamieson-Orchid-Elbow Creeks subdrainage of the Seymour River Basin, British Columbia. All of the 101 landslide cases were compared with 264 control grid cells. Multi-way cross classification tables were constructed to study the relationship between landsliding and several terrain attributes. A possible interaction between slope angle and the drainage condition of the soil was detected. A logistic regression analysis was then performed within a Geographic Information System (GIS) environment to develop a landslide risk model for the Jamieson-Orchid-Elbow Creeks study area. A landslide risk matrix was then constructed based on the landslide risk model. It was found that sites located in the transient snow zone, with slope angle greater than 55 degrees, on bedrock outcrop surficial material type and on shallow soil have the greatest risk of experiencing rapid, shallow soil mass movements. It was also found that holding all the other variables constant, slope angle had the greatest effect on the magnitude of landslide risk. Based on the data, sites with very steep slopes (over 55 degrees) have, on the average, five times the chance of experiencing a landslide event relative to sites with gentle slopes (10-25 degrees). The landslide risk matrix was used to create landslide risk categories. The spatial distribution of landslide risk, categorized as very low, low, moderate, high and very high, is portrayed within 20-m square grid cells on the landslide risk map. The major advantage of using the landslide risk assessment of this study is that it provides the terrain mapper with quantitative information about the relative risk of landsliding. This information can be used as a tool in planning watershed management activities and in an overall risk assessment for a given geographic area. |
Extent | 9526325 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2009-01-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0103800 |
URI | http://hdl.handle.net/2429/3807 |
Degree |
Master of Science - MSc |
Program |
Forestry |
Affiliation |
Forestry, Faculty of |
Degree Grantor | University of British Columbia |
Graduation Date | 1995-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- ubc_1995-0352.pdf [ 9.09MB ]
- Metadata
- JSON: 1.0103800.json
- JSON-LD: 1.0103800+ld.json
- RDF/XML (Pretty): 1.0103800.xml
- RDF/JSON: 1.0103800+rdf.json
- Turtle: 1.0103800+rdf-turtle.txt
- N-Triples: 1.0103800+rdf-ntriples.txt
- Original Record: 1.0103800 +original-record.json
- Full Text
- 1.0103800.txt
- Citation
- 1.0103800.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
United States | 18 | 2 |
United Kingdom | 14 | 0 |
Japan | 5 | 0 |
Bulgaria | 4 | 0 |
Czech Republic | 4 | 0 |
India | 3 | 0 |
France | 2 | 0 |
Russia | 2 | 0 |
China | 2 | 6 |
South Africa | 2 | 0 |
Germany | 1 | 3 |
Mauritius | 1 | 0 |
Ukraine | 1 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 28 | 4 |
Ashburn | 7 | 0 |
Tokyo | 5 | 0 |
Sofia | 4 | 0 |
Boardman | 3 | 0 |
Seattle | 3 | 0 |
Saint Petersburg | 2 | 0 |
Pietermaritzburg | 2 | 0 |
Rose Hill | 1 | 0 |
London | 1 | 0 |
Los Angeles | 1 | 0 |
Shenzhen | 1 | 5 |
Beijing | 1 | 1 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0103800/manifest